|
=={{header|Phix}}==
{{trans|GoAppleScript}}
Certainly puts my previous rubbish attempts ([[Self_numbers\Phix|archived here]]) to shame.<br>
Replacing the problematic sv[a+b+... line with a bit of dirty inline assembly improved performance by 90%<br>
The precise nature of the difference-pattern eludes me, I will admit.
(Of course you lose bounds checking, type checking, negative subscripts, fraction handling, and all that jazz.)<br>
<lang Phix>integer startIndex, endIndex, counter
We use a string of Y/N for the sieve to force one byte per element ('\0' and 1 would be equally valid).
atom currentSelf
<lang Phix>if machine_bits()=32 then crash("requires 64 bit") end if
sequence output
function sievedoneAfterAdding(integer interval)
-- Advance to the next self number in the sequence, append it to the output if required, indicate if finished.
string sv = repeat('N',2*1e9+9*9+1) -- (1.86GB)
integer ncurrentSelf += 0interval
forcounter a+=0 to 1 do
if counter < startIndex forthen b=0return tofalse 9end doif
output &= currentSelf
for c=0 to 9 do
return (counter = endIndex)
for d=0 to 9 do
for e=0 to 9 do
for f=0 to 9 do
for g=0 to 9 do
for h=0 to 9 do
for i=0 to 9 do
for j=0 to 9 do
-- n += 1
-- sv[a+b+c+d+e+f+g+h+i+j+n] = 'Y'
#ilASM{
[32] -- (allows clean compilation on 32 bit, before crash as above)
[64]
mov rax,[a]
mov r12,[b]
mov r13,[c]
mov r14,[d]
mov r15,[e]
add r12,r13
add r14,r15
add rax,r12
mov rdi,[sv]
add rax,r14
mov r12,[f]
mov r13,[g]
mov r14,[h]
mov r15,[i]
add r12,r13
shl rdi,2
mov rcx,[n]
mov r13,[j]
add r14,r15
add rax,r12
add rax,r14
add r13,rcx
add rax,r13
add rcx,1
mov byte[rdi+rax],'Y'
mov [n],rcx }
end for
end for
end for
end for
end for
end for
end for
end for
printf(1,"%d,%d\r",{a,b}) -- (show progress)
end for
end for
return sv
end function
atom t0 = time()
string sv = sieve()
printf(1,"sieve build took %s\n",{elapsed(time()-t0)})
integer count = 0
printf(1,"The first 50 self numbers are:\n")
for i=1 to length(sv) do
if sv[i]='N' then
count += 1
if count <= 50 then
printf(1,"%3d ",i-1)
if remainder(count,10)=0 then
printf(1,"\n")
end if
end if
if count == 1e8 then
printf(1,"\nThe %,dth self number is %,d (%s)\n",
{count,i-1,elapsed(time()-t0)})
exit
end if
end if
end for</lang>
{{out}}
<pre>
sieve build took 6.6s
The first 50 self numbers are:
1 3 5 7 9 20 31 42 53 64
75 86 97 108 110 121 132 143 154 165
176 187 198 209 211 222 233 244 255 266
277 288 299 310 312 323 334 345 356 367
378 389 400 411 413 424 435 446 457 468
The 100,000,000th self number is 1,022,727,208 (11.0s)
</pre>
=== generator dictionary ===
While this is dog-slow (see shocking estimate below), it is interesting to note that even by the time it generates the
10,000,000th number, it is only having to juggle a mere 27 generators.
Just a shame that we had to push over 10,000,000 generators onto that stack, and tried to push quite a few more.
Memory use is pretty low, around ~4MB.<br>
[unlike the above, this is perfectly happy on both 32 and 64 bit]<br>
Long story short: this works much the same as a prime sieve, in which you only need to eliminate multiples
of previous primes. Here, you only need to eliminate digital additions of previous safe numbers.
Only after writing this did I understand how to write a sliding sieve, which it turns out is a much better idea (see below).
Still, this is pretty interesting and quite neat.
<lang Phix>integer gd = new_dict(), gdhead = 2, n = 0
function nggetShorterRunLength(integer n)
-- Work out a shorter run length based on the most significant digit change about to happen.
integer r = n
integer shorterRunLength = 9,
while r do
n + temp = remainderfloor(r,10currentSelf/1000)
while r = floorremainder(r/temp,10)=9 do
shorterRunLength -= 1
temp = floor(temp/10)
end while
return nshorterRunLength
end function
function selfselfNumbers(integersequence /*i*/indexRange)
startIndex = indexRange[1]
-- note: assumes sequentially invoked (arg i unused)
nendIndex += 1indexRange[$]
whilecounter n=gdhead do0
gdheadcurrentSelf = pop_dict(gd)[-1]
output = {}
setd(ng(gdhead),0,gd)
n += (n!=gdhead)
-- Main process. Start with the single-digit odd numbers and first run.
end while
for i=1 to 5 do
setd(ng(n),0,gd)
if doneAfterAdding(2) then return output end if
return n
end functionfor
for i=1 to 9 do
if doneAfterAdding(11) then return output end if
end for
-- If necessary, fast forward to the last self number before the
atom t0 = time()
-- lowest-order block containing the first number required.
printf(1,"The first 50 self numbers are:\n")
if counter!=startIndex then
pp(apply(tagset(50),self),{pp_IntFmt,"%3d",pp_IntCh,false})
-- The highest-order blocks whose ends this is thought to handle correctly contain 97,777,777,778 self numbers.
-- The difference between equivalently positioned numbers in these blocks is 1,000,000,000,001.
constant limit = 10000000
-- The figures for successively lower-order blocks are obtained by successively removing 7s and 0s!
integer chk = 100
atom indexDiff = 97777777778,
printf(1,"\n i n size time\n")
numericDiff = 1000000000001
for i=51 to limit do
while indexDiff >= 98 and counter != startIndex do
n = self(i)
atom test = counter + indexDiff
if i=chk then
if (test > startIndex) then
printf(1,"%,11d %,11d %6d %s\n",{i,n,dict_size(gd),elapsed(time()-t0)})
chk * indexDiff = floor(indexDiff/10) + 1
numericDiff = floor(numericDiff/10) + 1
else
counter = test
currentSelf += numericDiff
end if
end while
end if
end for
-- Sequencing loop, per lowest-order block.
printf(1,"\nEstimated time for %,d :%s\n",{1e8,elapsed((time()-t0)*1e8/limit)})</lang>
{{out}}
<pre>
The first 50 self numbers are:
{ 1, 3, 5, 7, 9, 20, 31, 42, 53, 64, 75, 86, 97,108,110,121,132,143,
154,165,176,187,198,209,211,222,233,244,255,266,277,288,299,310,312,323,
334,345,356,367,378,389,400,411,413,424,435,446,457,468}
i n size time
100 973 18 0.1s
1,000 10,188 13 0.2s
10,000 102,225 10 1.0s
100,000 1,022,675 20 9.3s
1,000,000 10,227,221 17 1 minute and 37s
10,000,000 102,272,662 27 16 minutes and 40s
Estimated time for 100,000,000 :2 hours, 46 minutes and 37s
</pre>
For the record, I would not be at all surprised should a translation of this beat 20 minutes (for 1e8)
=== sliding sieve ===
Mid-speed, perhaps helped by a slightly smarter way of calculating/updating the digit sums<br>
Similar to some other entries, we (only) need a sieve of 9*9, +1 here as I test an entry after the slide.
<lang Phix>--sequence sieve = repeat(0,82), -- (~25% slower)
sequence sieve = repeat(0,8192),
digits = repeat(0,10)
integer offset = 0,
digit_sum = 0,
n = 0
procedure next_self()
while true do
-- Eight ten-number runs, each at a numeric interval of 2 from the end of the previous one.
n += 1
for i=length(digits)1 to 1 by -18 do
integerif ddoneAfterAdding(2) =then digits[i]return output end if
iffor d!j=1 to 9 thendo
digits[i]if =doneAfterAdding(11) d+1then return output end if
end digit_sum += 1for
exit
end if
digits[i] = 0
digit_sum -= 9
end for
-- Two shorter runs, the second at an interval inversely related to their length.
integer k = n+digit_sum-offset
integer shorterRunLength = getShorterRunLength(),
if k>length(sieve) then
integer j interval = 12
for i=n-offset1 to length(sieve)2 do
if doneAfterAdding(interval) then return sieve[j]output =end sieve[i]if
for j +=1 to shorterRunLength-1 do
if doneAfterAdding(11) then return output end if
end for
sieve[j..$]interval += 0(10-shorterRunLength)*13
end offset = n-1for
k = digit_sum+1
end if
sieve[k] = 1
if sieve[n-offset]=0 then exit end if
end while
--return (result in n)output
end procedurefunction
constant limit = 100000000
atom t0 = time()
printf(1,"The first 50 self numbers are:\n")
pp(selfNumbers({1, 50}),{pp_IntFmt,"%3d",pp_IntCh,false})
integer chk = 100
for ip=18 to limit9 do
next_selfinteger n = power(10,p)
printf(1,"The %,dth safe number is %,d\n",{n,selfNumbers({n})[1]})
if i<=50 then
printf(1," %3d%s",{n,iff(mod(i,25)=0?"\n":"")})
elsif i=chk then
if chk=100 then
printf(1,"\n i n time\n")
end if
printf(1,"%,12d %,13d %s\n",{i,n,elapsed(time()-t0)})
chk *= 10
end if
end for
printf(1,"\nEstimatedcompleted timein for %,d :%s\n",{1e9,elapsed((time()-t0)*1e9/limit)})</lang>
{{out}}
<pre>
The first 50 self numbers are:
{ 1 , 3 , 5 , 7 , 9 , 20 , 31 , 42 , 53 , 64 , 75 , 86 , 97 ,108 ,110 ,121 ,132 ,143 154 165 176 187 198 209 211,
154,165,176,187,198,209,211,222 ,233 ,244 ,255 ,266 ,277 ,288 ,299 ,310 ,312 ,323 334 345 356 367 378 389 400 411 413 424 435 446 457 468,
334,345,356,367,378,389,400,411,413,424,435,446,457,468}
The 100,000,000th safe number is 1,022,727,208
i n time
The 1,000,000,000th safe number is 10,227,272,649
100 973 0.1s
completed in 0.1s
1,000 10,188 0.1s
10,000 102,225 0.1s
100,000 1,022,675 0.1s
1,000,000 10,227,221 0.5s
10,000,000 102,272,662 4.5s
100,000,000 1,022,727,208 44.5s
Estimated time for 1,000,000,000 :7 minutes and 25s
</pre>
|