Untouchable numbers: Difference between revisions

Content added Content deleted
(added pascal fast but with problems after 40,000,000 count -> 6,430,223 < wrong by -1)
m (→‎{{header|Pascal}}: corrected SieveOneSieve Uint64(1) shl.. not 1 shl.., the 1 is Uint32 ...#@!... Now only use TIO.RUN)
Line 892: Line 892:
see also math.dartmouth.edu/~carlp/uupaper3.pdf with list count up to 1e9
see also math.dartmouth.edu/~carlp/uupaper3.pdf with list count up to 1e9
<lang pascal>program UntouchableNumbers;
<lang pascal>program UntouchableNumbers;
// gets factors of consecutive integers fast
//limited to 6.75E11 (Max smallprimes^2)
// limited to 1.2e11 ( = sqr(821641)) aka max(smallprimes)
{$IFDEF FPC}
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
Line 899: Line 900:
{$ENDIF}
{$ENDIF}
uses
uses
sysutils
sysutils,strutils {Numb2USA commatize}
{$IFDEF WINDOWS},Windows{$ENDIF}
{$IFDEF WINDOWS},Windows{$ENDIF}
;
;
const
//100*1000*1000;//10*1000*1000;//5*1000*1000;//1000*1000;
LIMIT = 1000*1000;
//384;//164;//110;//64;
LIMIT_mul = 64;

//######################################################################
//######################################################################
//prime decomposition
//prime decomposition
//http://rosettacode.org/wiki/Factors_of_an_integer#using_Prime_decomposition
const
const
SizePrDeFe = 6*8192;//*size of(tprimeFac) =40 byte level I or 2 Mb ~ level 2 cache
SizePrDeFe = 6*8192;//*size of(tprimeFac) =56 byte level I or 2 Mb ~ level 2 cache
type
type
tdigits = array [0..31] of Uint32;
tdigits = array [0..31] of Uint32;
//the first number with 12 different prime factors =
//the first number with 12 different prime factors =
//2*3*5*7*11*13*17*19*23*29*31*37 = 7,42E12
//2*3*5*7*11*13*17*19*23*29*31*37 = 7,42e12
//40 byte
//56 byte
tprimeFac = packed record
tprimeFac = packed record
pfSumOfDivs,
pfSumOfDivs,
Line 1,067: Line 1,073:
var
var
dgt:tDigits;
dgt:tDigits;
i,j,k,pr,fac,n,MaxP : NativeUInt;
i,j,k,pr,fac,n,MaxP : Uint64;
begin
begin
n := pdfOfs;
n := pdfOfs;
Line 1,098: Line 1,104:
pfpotMax[0] := j;
pfpotMax[0] := j;
pfRemain := (n+i) shr j;
pfRemain := (n+i) shr j;
pfSumOfDivs := (1 shl (j+1))-1;
pfSumOfDivs := (Uint64(1) shl (j+1))-1;
// writeln(n+i,' # ',j,' ## ',pfRemain,' ### ',pfSumOfDivs);
pfDivCnt := j+1;
pfDivCnt := j+1;
end;
end;
Line 1,197: Line 1,204:
end;
end;


procedure CheckRest(n: Uint64;pUntouch:pByte);
const
var
LIMIT = 100*1000*1000;
k : Uint64;
LIMIT_mul = 6 * (1 shl (trunc(ln(Limit)/ln(10))-2));
begin
repeat
k := GetNextPrimeDecomp^.pfSumOfDivs-n;
inc(n);
if (k <= LIMIT) then
pUntouch[k ] := 1;
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;

var
var
Untouch : array of byte;
Untouch : array of byte;
pUntouch: pByte;
pUntouch: pByte;
pPrimeDecomp :tpPrimeFac;
T0:Int64;
T0:Int64;
n,k,lim,prmIdx : NativeInt;
n,k,lim,prmEndIdx : NativeInt;
Begin
Begin
setlength(untouch,LIMIT+1);
setlength(untouch,LIMIT+1);
pUntouch := @untouch[0];
pUntouch := @untouch[0];

InitSmallPrimes;
InitSmallPrimes;
T0 := GetTickCount64;
T0 := GetTickCount64;
prmIdx := 0;
prmEndIdx := 0;
REPEAT
repeat
inc(prmIdx);
inc(prmEndIdx);
until smallprimes[prmIdx]> 2*sqrt(LIMIT);
until smallprimes[prmEndIdx] > trunc(exp(ln(Limit)/2.32));
writeln(prmIdx:10,smallprimes[prmIdx]:10);
writeln(prmEndIdx:10,smallprimes[prmEndIdx]:10);

writeln(LIMIT_mul);
writeln(LIMIT_mul);

n := 0;
n := 0;
Init_Sieve(0);
Init_Sieve(n);
pUntouch[0] := 1;
pUntouch[1] := 1;//all primes
repeat
repeat
pPrimeDecomp:= GetNextPrimeDecomp;
k := GetNextPrimeDecomp^.pfSumOfDivs-n;
inc(n);//n-> n+1
k := pPrimeDecomp^.pfSumOfDivs-n;
If k <> 1 then
if k <= LIMIT then
Begin
if k > 1 then
if k <= LIMIT then
pUntouch[k] := 1;
end
else
begin
begin
If k <> 1 then
//n= prime, would be marked by n*n with proper factors 1,n
if n<= LIMIT+1 then
pUntouch[k] := 1
Begin
else
pUntouch[n+1] := 1;
if n>3 then
prmEndIdx := CheckPrime(n,prmEndIdx,pUntouch);
//marked by prime*n with proper factors 1,(prime),n
For lim := 0 to prmIdx do
begin
k := smallprimes[lim]+n+1;
If k > LIMIT then
begin
prmIdx:= lim-1;
BREAK;
end;
pUntouch[k] := 1;
end;
end;
end;
end;
until n >LIMIT;
writeln('runtime for n<= LIMIT ',(GetTickCount64-T0)/1000:0:3,' s');
writeln('Check the rest ',Numb2USA(IntToStr((LIMIT_mul-1)*Limit)));
CheckRest(n,pUntouch);


inc(n);
until n >LIMIT_mul*LIMIT;
T0 := GetTickCount64-T0;
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
writeln('runtime ',T0/1000:0:3,' s');
Line 1,259: Line 1,283:
if n = lim then
if n = lim then
Begin
Begin
writeln(lim:10,k:10);
writeln(Numb2USA(IntToStr(lim)):10,Numb2USA(IntToStr(k)):10);
lim *= 10;
lim *= 10;
end;
end;
Line 1,269: Line 1,293:
if n = lim then
if n = lim then
Begin
Begin
writeln(lim:10,k:10);
writeln(Numb2USA(IntToStr(lim)):10,Numb2USA(IntToStr(k)):10);
lim += LIMIT DIV 10;
lim += LIMIT DIV 10;
end;
end;
Line 1,278: Line 1,302:
{{out}}
{{out}}
<pre>
<pre>
TIO.RUN
runtime 2.486 s
runtime for n<= LIMIT 0.070 s
Check the rest 63,000,000
runtime 4.214 s
10 2
10 2
100 5
100 5
1000 89
1,000 89
10000 1212
10,000 1,212
100000 13863
100,000 13,863
200000 28572
200,000 28,572
300000 43514
300,000 43,514
400000 58459
400,000 58,459
500000 73565
500,000 73,565
600000 88828
600,000 88,828
700000 104061
700,000 104,061
800000 119302
800,000 119,302
900000 134757
900,000 134,757
1,000,000 150,232
1000000 150232
####
823 6329
192 -> test til 192x10,000,000= 1.92e9
runtime 75.065 s
10 2
100 5
1000 89
10000 1212
100000 13863
1000000 150232
2000000 305290
3000000 462110
4000000 619638
5000000 777672
6000000 936243
7000000 1095710
8000000 1255015
9000000 1414783
10000000 1574973
####
2262 20011
384 -> test til 384x100,000,000= 38.4e9

runtime 1750.728 s
10 2
100 5
1000 89
10000 1212
100000 13863
1000000 150232
10000000 1574973
20000000 3184111
30000000 4804331
40000000 6430223 < wrong -1
50000000 8060162 < wrong -1
60000000 9694467
70000000 11330312
80000000 12967238 < wrong -1
90000000 14606549
100000000 16246941 < wrong +1


real 29m10,824s
//url=https://math.dartmouth.edu/~carlp/uupaper3.pdf
//url=https://math.dartmouth.edu/~carlp/uupaper3.pdf
6000000 936244
6000000 936244
Line 1,351: Line 1,337:
100000000 16246940
100000000 16246940
</pre>
</pre>

=={{header|Perl}}==
=={{header|Perl}}==
{{libheader|ntheory}}
{{libheader|ntheory}}
Retrieved from "https://rosettacode.org/wiki/Untouchable_numbers"