Rosetta Code

Untouchable numbers: Difference between revisions

Page Discussion

Untouchable numbers (view source)

Revision as of 14:00, 2 September 2021

136 bytes removed ,  2 years ago
m
→‎{{header|Pascal}}: corrected SieveOneSieve Uint64(1) shl.. not 1 shl.., the 1 is Uint32 ...#@!... Now only use TIO.RUN
(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:
see also math.dartmouth.edu/~carlp/uupaper3.pdf with list count up to 1e9
<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}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
Line 899 ⟶ 900:
{$ENDIF}
uses
sysutils,strutils {Numb2USA commatize}
{$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
//http://rosettacode.org/wiki/Factors_of_an_integer#using_Prime_decomposition
const
SizePrDeFe = 6*8192;//*size of(tprimeFac) =4056 byte level I or 2 Mb ~ level 2 cache
type
tdigits = array [0..31] of Uint32;
//the first number with 12 different prime factors =
//2*3*5*7*11*13*17*19*23*29*31*37 = 7,42E1242e12
//4056 byte
tprimeFac = packed record
pfSumOfDivs,
Line 1,067 ⟶ 1,073:
var
dgt:tDigits;
i,j,k,pr,fac,n,MaxP : NativeUIntUint64;
begin
n := pdfOfs;
Line 1,098 ⟶ 1,104:
pfpotMax[0] := j;
pfRemain := (n+i) shr j;
pfSumOfDivs := (Uint64(1) shl (j+1))-1;
// writeln(n+i,' # ',j,' ## ',pfRemain,' ### ',pfSumOfDivs);
pfDivCnt := j+1;
end;
Line 1,197 ⟶ 1,204:
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
Untouch : array of byte;
pUntouch: pByte;
pPrimeDecomp :tpPrimeFac;
T0:Int64;
n,k,lim,prmIdxprmEndIdx : NativeInt;
Begin
setlength(untouch,LIMIT+1);
pUntouch := @untouch[0];
 
InitSmallPrimes;
T0 := GetTickCount64;
prmIdxprmEndIdx := 0;
REPEATrepeat
inc(prmIdxprmEndIdx);
until smallprimes[prmIdxprmEndIdx] > 2*sqrttrunc(LIMITexp(ln(Limit)/2.32));
writeln(prmIdxprmEndIdx:10,smallprimes[prmIdxprmEndIdx]:10);
 
writeln(LIMIT_mul);
 
n := 0;
Init_Sieve(0n);
pUntouch[0] := 1;
pUntouch[1] := 1;//all primes
repeat
pPrimeDecompk := GetNextPrimeDecomp^.pfSumOfDivs-n;
inc(n);//n-> n+1
k := pPrimeDecomp^.pfSumOfDivs-n;
Ifif k <>= 1LIMIT then
Begin
if k > 1 then
if k <= LIMIT then
pUntouch[k] := 1;
end
else
begin
If k <> 1 then
//n= prime, would be marked by n*n with proper factors 1,n
if n< pUntouch[k] := LIMIT+1 then
Beginelse
pUntouch[n+1]if :=n>3 1;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;
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;
writeln('runtime ',T0/1000:0:3,' s');
Line 1,259 ⟶ 1,283:
if n = lim then
Begin
writeln(Numb2USA(IntToStr(lim)):10,Numb2USA(IntToStr(k)):10);
lim *= 10;
end;
Line 1,269 ⟶ 1,293:
if n = lim then
Begin
writeln(Numb2USA(IntToStr(lim)):10,Numb2USA(IntToStr(k)):10);
lim += LIMIT DIV 10;
end;
Line 1,278 ⟶ 1,302:
{{out}}
<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
100 5
10001,000 89
10,000 10000 12121,212
100,000 100000 1386313,863
200,000 200000 2857228,572
300,000 300000 4351443,514
400,000 400000 5845958,459
500,000 500000 7356573,565
600,000 600000 8882888,828
700,000 700000 104061104,061
800,000 800000 119302119,302
900,000 900000 134757134,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
6000000 936244
Line 1,351 ⟶ 1,337:
100000000 16246940
</pre>
 
=={{header|Perl}}==
{{libheader|ntheory}}
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/307846"