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. #0..65535.Sieving only odd numbers
//get primes. Sieving only odd numbers
const
MAXLIMIT = (821641-1) shr 1;
var
var
pr : array[0..MAXLIMIT] of byte;
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>MAXLIMIT then
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>MAXLIMIT;
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] := d;
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 > MAXLIMIT) OR (j>High(SmallPrimes));
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 <= LIMIT) then
if Not(ODD(k)) AND (k<=LIMIT) then
pUntouch[k ] := 1;
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,lim,prmEndIdx : NativeInt;
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
if n>3 then
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
10 2
200 10
100 5
300 22
1000 89
400 30
10000 1212
500 38
100000 13863
600 48
1000000 150232
700 55
1500000 227297
800 69
2000000 305290
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>


Retrieved from "https://rosettacode.org/wiki/Untouchable_numbers"