Untouchable numbers: Difference between revisions
Content added Content deleted
m (→{{header|Pascal}}: reduced to only calc sum of divisors nearly halves time on TIO.RUN.Longer list of count of untouchable numbers) |
(→{{header|Pascal}}: changed search range to n^(4/3) an unmark all odd numbers and correct 5.Now the count untouchable numbers fits up to 1E8) |
||
Line 893: | Line 893: | ||
Nigel is still right, that this procedure is not the way to get proven results. |
Nigel is still right, that this procedure is not the way to get proven results. |
||
<lang pascal>program UntouchableNumbers; |
<lang pascal>program UntouchableNumbers; |
||
program UntouchableNumbers; |
|||
{$IFDEF FPC} |
{$IFDEF FPC} |
||
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON} |
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON} |
||
{$CODEALIGN proc=16,loop=4} |
|||
{$ELSE} |
{$ELSE} |
||
{$APPTYPE CONSOLE} |
{$APPTYPE CONSOLE} |
||
Line 903: | Line 905: | ||
; |
; |
||
const |
const |
||
MAXPRIME = 1742537; |
|||
//100*1000*1000;//10*1000*1000;//5*1000*1000;//1000*1000; |
|||
//sqr(MaxPrime) = 3e12 |
|||
LIMIT = 5*1000*1000; |
LIMIT = 5*1000*1000; |
||
LIMIT_mul = trunc(exp(ln(LIMIT)/3))+1; |
|||
//384;//152;//104;//64; |
|||
LIMIT_mul = 104; |
|||
//###################################################################### |
|||
//gets sum of divisors of consecutive integers fast |
|||
const |
const |
||
SizePrDeFe = 16*8192;//*size of(tprimeFac) =16 byte 2 Mb ~ level 3 cache |
SizePrDeFe = 16*8192;//*size of(tprimeFac) =16 byte 2 Mb ~ level 3 cache |
||
Line 921: | Line 921: | ||
tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac; |
tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac; |
||
tPrimes = array[0..65535] of Uint32; |
|||
tPrimes = array[0..1 shl 17-1] of Uint32; |
|||
var |
var |
||
Line 929: | Line 930: | ||
SmallPrimes: tPrimes; |
SmallPrimes: tPrimes; |
||
pdfIDX,pdfOfs: NativeInt; |
pdfIDX,pdfOfs: NativeInt; |
||
TD : Int64; |
|||
procedure OutCounts(pUntouch:pByte); |
|||
var |
|||
n,cnt,lim,deltaLim : NativeInt; |
|||
Begin |
|||
n := 0; |
|||
cnt := 0; |
|||
deltaLim := 100; |
|||
lim := deltaLim; |
|||
repeat |
|||
cnt += 1-pUntouch[n]; |
|||
if n = lim then |
|||
Begin |
|||
writeln(Numb2USA(IntToStr(lim)):13,' ',Numb2USA(IntToStr(cnt)):12); |
|||
lim += deltaLim; |
|||
if lim = 10*deltaLim then |
|||
begin |
|||
deltaLim *=10; |
|||
lim := deltaLim; |
|||
writeln; |
|||
end; |
|||
end; |
|||
inc(n); |
|||
until n > LIMIT; |
|||
end; |
|||
function OutN(n:UInt64):UInt64; |
|||
begin |
|||
write(Numb2USA(IntToStr(n)):15,' dt ',(GettickCount64-TD)/1000:5:3,' s'#13); |
|||
TD := GettickCount64; |
|||
result := n+LIMIT; |
|||
end; |
|||
//###################################################################### |
|||
//gets sum of divisors of consecutive integers fast |
|||
procedure InitSmallPrimes; |
procedure InitSmallPrimes; |
||
//get primes. |
//get primes. Sieving only odd numbers |
||
const |
|||
MAXLIMIT = (821641-1) shr 1; |
|||
var |
var |
||
pr : array[0.. |
pr : array[0..MAXPRIME] of byte; |
||
p,j,d,flipflop :NativeUInt; |
p,j,d,flipflop :NativeUInt; |
||
Begin |
Begin |
||
Line 946: | Line 981: | ||
until pr[p]= 0; |
until pr[p]= 0; |
||
j := (p+1)*p*2; |
j := (p+1)*p*2; |
||
if j> |
if j>MAXPRIME then |
||
BREAK; |
BREAK; |
||
d := 2*p+1; |
d := 2*p+1; |
||
Line 952: | Line 987: | ||
pr[j] := 1; |
pr[j] := 1; |
||
j += d; |
j += d; |
||
until j> |
until j>MAXPRIME; |
||
until false; |
until false; |
||
Line 958: | Line 993: | ||
SmallPrimes[2] := 5; |
SmallPrimes[2] := 5; |
||
j := 3; |
j := 3; |
||
d := 7; |
|||
flipflop := (2+1)-1;//7+2*2->11+2*1->13 ,17 ,19 , 23 |
flipflop := (2+1)-1;//7+2*2->11+2*1->13 ,17 ,19 , 23 |
||
p := 3; |
p := 3; |
||
Line 964: | Line 998: | ||
if pr[p] = 0 then |
if pr[p] = 0 then |
||
begin |
begin |
||
SmallPrimes[j] := |
SmallPrimes[j] := 2*p+1; |
||
inc(j); |
inc(j); |
||
end; |
end; |
||
d += 2*flipflop; |
|||
p+=flipflop; |
p+=flipflop; |
||
flipflop := 3-flipflop; |
flipflop := 3-flipflop; |
||
until (p > |
until (p > MAXPRIME) OR (j>High(SmallPrimes)); |
||
end; |
end; |
||
Line 1,129: | Line 1,162: | ||
if pdfIDX >= SizePrDeFe then |
if pdfIDX >= SizePrDeFe then |
||
if Not(NextSieve) then |
if Not(NextSieve) then |
||
Begin |
|||
writeln('of limits '); |
|||
EXIT(NIL); |
EXIT(NIL); |
||
end; |
|||
result := @PrimeDecompField[pdfIDX]; |
result := @PrimeDecompField[pdfIDX]; |
||
inc(pdfIDX); |
inc(pdfIDX); |
||
Line 1,141: | Line 1,177: | ||
result := SieveOneSieve(PrimeDecompField); |
result := SieveOneSieve(PrimeDecompField); |
||
end; |
end; |
||
//gets sum of divisors of consecutive integers fast |
|||
//###################################################################### |
|||
procedure CheckRest(n: Uint64;pUntouch:pByte); |
procedure CheckRest(n: Uint64;pUntouch:pByte); |
||
var |
var |
||
k : Uint64; |
k,lim : Uint64; |
||
begin |
begin |
||
lim := 2*LIMIT; |
|||
repeat |
repeat |
||
k := GetNextPrimeDecomp^.pfSumOfDivs-n; |
k := GetNextPrimeDecomp^.pfSumOfDivs-n; |
||
inc(n); |
inc(n); |
||
if (k <= |
if Not(ODD(k)) AND (k<=LIMIT) then |
||
pUntouch[k |
pUntouch[k] := 1; |
||
// showing still alive not for TIO.RUN |
|||
// if n >= lim then lim := OutN(n); |
|||
until n >LIMIT_mul*LIMIT; |
until n >LIMIT_mul*LIMIT; |
||
end; |
|||
function CheckPrime(n:Uint64;prmEndIdx:NativeInt;pUntouch:pByte):NativeInt; |
|||
var |
|||
i,k: NativeInt; |
|||
Begin |
|||
//n= prime,n+1 would be marked by n*n with proper factors 1,n |
|||
//here n is aready n+1 |
|||
pUntouch[n] := 1; |
|||
//marked by prime*n with proper factors 1,(prime),n |
|||
For i := 0 to prmEndIdx do |
|||
begin |
|||
k := smallprimes[i]+n; |
|||
If k > LIMIT then |
|||
Begin |
|||
dec(prmEndIdx); |
|||
BREAK; |
|||
end; |
|||
pUntouch[k] := 1; |
|||
end; |
|||
result := prmEndIdx; |
|||
end; |
end; |
||
Line 1,178: | Line 1,198: | ||
Untouch : array of byte; |
Untouch : array of byte; |
||
pUntouch: pByte; |
pUntouch: pByte; |
||
puQW : pQword; |
|||
T0:Int64; |
T0:Int64; |
||
n,k |
n,k : NativeInt; |
||
Begin |
Begin |
||
if sqrt(LIMIT_mul*LIMIT) >=MAXPRIME then |
|||
setlength(untouch,LIMIT+1); |
|||
Begin |
|||
writeln('Need to extend count of primes > ', |
|||
trunc(sqrt(LIMIT_mul*LIMIT))+1); |
|||
HALT(0); |
|||
end; |
|||
setlength(untouch,LIMIT+8+1); |
|||
pUntouch := @untouch[0]; |
pUntouch := @untouch[0]; |
||
//Mark all odd as touchable |
|||
puQW := @pUntouch[0]; |
|||
For n := 0 to LIMIT DIV 8 do puQW[n] := $0100010001000100; |
|||
InitSmallPrimes; |
InitSmallPrimes; |
||
T0 := GetTickCount64; |
T0 := GetTickCount64; |
||
prmEndIdx := 0; |
|||
repeat |
|||
inc(prmEndIdx); |
|||
until smallprimes[prmEndIdx] > trunc(sqrt(Limit)); |
|||
writeln('LIMIT = ',Numb2USA(IntToStr(LIMIT))); |
writeln('LIMIT = ',Numb2USA(IntToStr(LIMIT))); |
||
writeln(prmEndIdx:10,smallprimes[prmEndIdx]:10); |
|||
writeln('factor beyond LIMIT ',LIMIT_mul); |
writeln('factor beyond LIMIT ',LIMIT_mul); |
||
n := 0; |
n := 0; |
||
Init_Sieve(n); |
Init_Sieve(n); |
||
pUntouch[0] := 1; |
|||
pUntouch[1] := 1;//all primes |
pUntouch[1] := 1;//all primes |
||
repeat |
repeat |
||
Line 1,206: | Line 1,232: | ||
pUntouch[k] := 1 |
pUntouch[k] := 1 |
||
else |
else |
||
begin |
|||
//n-1 is prime p |
|||
prmEndIdx := CheckPrime(n,prmEndIdx,pUntouch); |
|||
//mark p*p |
|||
pUntouch[n] := 1; |
|||
//mark 2*p |
|||
//5 marked by prime 2 but that is p*p, but 4 has factor sum = 3 |
|||
pUntouch[n+2] := 1; |
|||
end; |
|||
end; |
end; |
||
until n >LIMIT; |
until n > LIMIT-2; |
||
//unmark 5 and mark 0 |
|||
puntouch[5] := 0; |
|||
pUntouch[0] := 1; |
|||
//n=limit-1 |
|||
k := GetNextPrimeDecomp^.pfSumOfDivs-n; |
|||
inc(n); |
|||
If (k <> 1) AND (k<=LIMIT) then |
|||
pUntouch[k] := 1 |
|||
else |
|||
pUntouch[n] := 1; |
|||
//n=limit |
|||
k := GetNextPrimeDecomp^.pfSumOfDivs-n; |
|||
If Not(odd(k)) AND (k<=LIMIT) then |
|||
pUntouch[k] := 1; |
|||
n:= limit+1; |
|||
writeln('runtime for n<= LIMIT ',(GetTickCount64-T0)/1000:0:3,' s'); |
writeln('runtime for n<= LIMIT ',(GetTickCount64-T0)/1000:0:3,' s'); |
||
writeln('Check the rest ',Numb2USA(IntToStr((LIMIT_mul-1)*Limit))); |
writeln('Check the rest ',Numb2USA(IntToStr((LIMIT_mul-1)*Limit))); |
||
TD := GettickCount64; |
|||
CheckRest(n,pUntouch); |
CheckRest(n,pUntouch); |
||
writeln('runtime ',(GetTickCount64-T0)/1000:0:3,' s'); |
|||
T0 := GetTickCount64-T0; |
T0 := GetTickCount64-T0; |
||
writeln('runtime ',T0/1000:0:3,' s'); |
|||
OutCounts(pUntouch); |
|||
k := 0; |
|||
lim :=10; |
|||
for n := 0 to LIMIT DIV 10 do |
|||
Begin |
|||
if n = lim then |
|||
Begin |
|||
writeln(lim:10,k:10); |
|||
lim *= 10; |
|||
end; |
|||
k += 1-pUntouch[n]; |
|||
end; |
|||
lim := 2*LIMIT DIV 10; |
|||
for n := LIMIT DIV 10+1 to LIMIT do |
|||
Begin |
|||
if n = lim then |
|||
Begin |
|||
writeln(lim:10,k:10); |
|||
lim += LIMIT DIV 10; |
|||
end; |
|||
k += 1-pUntouch[n]; |
|||
end; |
|||
end.</lang> |
end.</lang> |
||
{{out}} |
{{out}} |
||
Line 1,243: | Line 1,273: | ||
TIO.RUN |
TIO.RUN |
||
LIMIT = 5,000,000 |
LIMIT = 5,000,000 |
||
factor beyond LIMIT 171 |
|||
331 2237 |
|||
runtime for n smaller than LIMIT 0.204 s |
|||
factor beyond LIMIT 104 |
|||
Check the rest 850,000,000 |
|||
runtime for n<= LIMIT 0.272 s |
|||
runtime 32.643 s |
|||
Check the rest 515,000,000 |
|||
100 5 |
|||
runtime 19.052 s |
|||
200 10 |
|||
300 22 |
|||
400 30 |
|||
500 38 |
|||
600 48 |
|||
700 55 |
|||
800 69 |
|||
900 80 |
|||
2500000 383422 |
|||
3000000 462110 |
|||
3500000 540769 |
|||
4000000 619638 |
|||
4500000 698504 |
|||
5000000 777672 |
|||
1,000 89 |
|||
2,000 196 |
|||
3,000 318 |
|||
4,000 443 |
|||
5,000 570 |
|||
6,000 689 |
|||
7,000 801 |
|||
8,000 936 |
|||
9,000 1,082 |
|||
10,000 1,212 |
|||
20,000 2,566 |
|||
30,000 3,923 |
|||
40,000 5,324 |
|||
50,000 6,705 |
|||
60,000 8,153 |
|||
70,000 9,586 |
|||
80,000 10,994 |
|||
90,000 12,429 |
|||
100,000 13,863 |
|||
200,000 28,572 |
|||
300,000 43,515 |
|||
400,000 58,459 |
|||
500,000 73,565 |
|||
600,000 88,828 |
|||
700,000 104,062 |
|||
800,000 119,302 |
|||
900,000 134,758 |
|||
1,000,000 150,232 |
|||
2,000,000 305,290 |
|||
3,000,000 462,110 |
|||
4,000,000 619,638 |
|||
5,000,000 777,672 |
|||
Real time: 32.827 s CPU share: 99.30 % |
|||
//url=https://math.dartmouth.edu/~carlp/uupaper3.pdf |
//url=https://math.dartmouth.edu/~carlp/uupaper3.pdf |
||
100000 13863 |
100000 13863 |
||
Line 1,292: | Line 1,353: | ||
90000000 14606549 |
90000000 14606549 |
||
100000000 16246940 |
100000000 16246940 |
||
... at home up to 1E8 |
|||
50,000,000 8,060,163 |
|||
60,000,000 9,694,467 |
|||
70,000,000 11,330,312 |
|||
80,000,000 12,967,239 |
|||
90,000,000 14,606,549 |
|||
100,000,000 16,246,940 |
|||
real 18m51,234s |
|||
</pre> |
</pre> |
||