Rosetta Code

Untouchable numbers: Difference between revisions

Page Discussion

Untouchable numbers (view source)

Revision as of 05:47, 6 September 2021

1,785 bytes added ,  2 years ago
→‎{{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
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:
Nigel is still right, that this procedure is not the way to get proven results.
<lang pascal>program UntouchableNumbers;
program UntouchableNumbers;
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$CODEALIGN proc=16,loop=4}
{$ELSE}
{$APPTYPE CONSOLE}
Line 903 ⟶ 905:
;
const
MAXPRIME = 1742537;
//100*1000*1000;//10*1000*1000;//5*1000*1000;//1000*1000;
//sqr(MaxPrime) = 3e12
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
SizePrDeFe = 16*8192;//*size of(tprimeFac) =16 byte 2 Mb ~ level 3 cache
Line 921:
 
tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac;
 
tPrimes = array[0..65535] of Uint32;
tPrimes = array[0..1 shl 17-1] of Uint32;
 
var
Line 929 ⟶ 930:
SmallPrimes: tPrimes;
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;
//get primes. #0..65535.Sieving only odd numbers
const
MAXLIMIT = (821641-1) shr 1;
var
pr : array[0..MAXLIMITMAXPRIME] of byte;
p,j,d,flipflop :NativeUInt;
Begin
Line 946 ⟶ 981:
until pr[p]= 0;
j := (p+1)*p*2;
if j>MAXLIMITMAXPRIME then
BREAK;
d := 2*p+1;
Line 952 ⟶ 987:
pr[j] := 1;
j += d;
until j>MAXLIMITMAXPRIME;
until false;
 
Line 958 ⟶ 993:
SmallPrimes[2] := 5;
j := 3;
d := 7;
flipflop := (2+1)-1;//7+2*2->11+2*1->13 ,17 ,19 , 23
p := 3;
Line 964 ⟶ 998:
if pr[p] = 0 then
begin
SmallPrimes[j] := d2*p+1;
inc(j);
end;
d += 2*flipflop;
p+=flipflop;
flipflop := 3-flipflop;
until (p > MAXLIMITMAXPRIME) OR (j>High(SmallPrimes));
end;
 
Line 1,129 ⟶ 1,162:
if pdfIDX >= SizePrDeFe then
if Not(NextSieve) then
Begin
writeln('of limits ');
EXIT(NIL);
end;
result := @PrimeDecompField[pdfIDX];
inc(pdfIDX);
Line 1,141 ⟶ 1,177:
result := SieveOneSieve(PrimeDecompField);
end;
//gets sum of divisors of consecutive integers fast
//######################################################################
 
procedure CheckRest(n: Uint64;pUntouch:pByte);
var
k,lim : Uint64;
begin
lim := 2*LIMIT;
repeat
k := GetNextPrimeDecomp^.pfSumOfDivs-n;
inc(n);
if Not(ODD(k)) AND (k<= LIMIT) then
pUntouch[k ] := 1;
// showing still alive not for TIO.RUN
// if n >= lim then lim := OutN(n);
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;
 
Line 1,178 ⟶ 1,198:
Untouch : array of byte;
pUntouch: pByte;
puQW : pQword;
T0:Int64;
n,k,lim,prmEndIdx : NativeInt;
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];
//Mark all odd as touchable
puQW := @pUntouch[0];
For n := 0 to LIMIT DIV 8 do puQW[n] := $0100010001000100;
 
InitSmallPrimes;
T0 := GetTickCount64;
prmEndIdx := 0;
repeat
inc(prmEndIdx);
until smallprimes[prmEndIdx] > trunc(sqrt(Limit));
writeln('LIMIT = ',Numb2USA(IntToStr(LIMIT)));
writeln(prmEndIdx:10,smallprimes[prmEndIdx]:10);
writeln('factor beyond LIMIT ',LIMIT_mul);
 
n := 0;
Init_Sieve(n);
 
pUntouch[0] := 1;
pUntouch[1] := 1;//all primes
repeat
Line 1,206 ⟶ 1,232:
pUntouch[k] := 1
else
if n>3 thenbegin
//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;
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('Check the rest ',Numb2USA(IntToStr((LIMIT_mul-1)*Limit)));
TD := GettickCount64;
CheckRest(n,pUntouch);
writeln('runtime ',(GetTickCount64-T0)/1000:0:3,' s');
 
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>
{{out}}
Line 1,243 ⟶ 1,273:
TIO.RUN
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 200 2 10
100 300 5 22
1000 400 89 30
10000 500 1212 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
100000 13863
Line 1,292 ⟶ 1,353:
90000000 14606549
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>
 
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/307848"