Jump to content

Untouchable numbers: Difference between revisions

m
→‎{{header|Pascal}}: reduced to only calc sum of divisors nearly halves time on TIO.RUN.Longer list of count of untouchable numbers
m (→‎{{header|Pascal}}: corrected SieveOneSieve Uint64(1) shl.. not 1 shl.., the 1 is Uint32 ...#@!... Now only use TIO.RUN)
m (→‎{{header|Pascal}}: reduced to only calc sum of divisors nearly halves time on TIO.RUN.Longer list of count of untouchable numbers)
Line 889:
There are 13863 untouchable numbers ≤ 100000.</pre>
=={{header|Pascal}}==
modified [[Factors_of_an_integer#using_Prime_decomposition]] to calculate only sum of divisors<BR>
seeAppended a list of count of untouchable numbers out of also math.dartmouth.edu/~carlp/uupaper3.pdf with list count up to 1e9<BR>
Nigel is still right, that this procedure is not the way to get proven results.
<lang pascal>program UntouchableNumbers;
// gets factors of consecutive integers fast
// limited to 1.2e11 ( = sqr(821641)) aka max(smallprimes)
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
Line 900 ⟶ 899:
{$ENDIF}
uses
sysutils,strutils {Numb2USA commatize}
{$IFDEF WINDOWS},Windows{$ENDIF}
;
const
//100*1000*1000;//10*1000*1000;//5*1000*1000;//1000*1000;
LIMIT = 5*1000*1000;
//384;//164152;//110104;//64;
LIMIT_mul = 64104;
 
//######################################################################
//gets sum of divisors of consecutive integers fast
//prime decomposition
const
SizePrDeFe = 616*8192;//*size of(tprimeFac) =5616 byte level I or 2 Mb ~ level 23 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,42e12
//56 byte
tprimeFac = packed record
pfSumOfDivs,
pfRemain : Uint64;
pfpotPrimIdx : array[0..11] of word;
pfDivCnt : Uint16;
pfpotMax : array[0..12] of byte;
pfMaxIdx : byte;
end;
tpPrimeFac = ^tprimeFac;
Line 932 ⟶ 924:
 
var
{$ALIGN 816}
SmallPrimes: tPrimes;
{$ALIGN 32}
PrimeDecompField :tPrimeDecompField;
{$ALIGN 16}
SmallPrimes: tPrimes;
pdfIDX,pdfOfs: NativeInt;
 
Line 979 ⟶ 971:
flipflop := 3-flipflop;
until (p > MAXLIMIT) OR (j>High(SmallPrimes));
end;
 
function smplPrimeDecomp(n:Uint64):tprimeFac;
var
pr,i,pot,fac,q :NativeUInt;
Begin
with result do
Begin
pfDivCnt := 1;
pfSumOfDivs := 1;
pfRemain := n;
pfMaxIdx := 0;
pfpotPrimIdx[0] := 1;
pfpotMax[0] := 0;
 
i := 0;
while i < High(SmallPrimes) do
begin
pr := SmallPrimes[i];
q := n DIV pr;
//if n < pr*pr
if pr > q then
BREAK;
if n = pr*q then
Begin
pfpotPrimIdx[pfMaxIdx] := i;
pot := 0;
fac := pr;
repeat
n := q;
q := n div pr;
pot+=1;
fac *= pr;
until n <> pr*q;
pfpotMax[pfMaxIdx] := pot;
pfDivCnt *= pot+1;
pfSumOfDivs *= (fac-1)DIV(pr-1);
inc(pfMaxIdx);
end;
inc(i);
end;
pfRemain := n;
if n > 1 then
Begin
pfDivCnt *= 2;
pfSumOfDivs *= n+1
end;
end;
end;
 
Line 1,054 ⟶ 998:
end;
 
function IncByBaseInBaseIncByOneInBase(var dgt:tDigits;base:NativeInt):NativeInt;
var
q :NativeInt;
Line 1,068 ⟶ 1,012:
dgt[result] := q;
result +=1;
end;
 
procedure CalcSumOfDivs(var pdf:tPrimeDecompField;var dgt:tDigits;n,k,pr:Uint64);
var
fac,s :Uint64;
j : Int32;
Begin
//j is power of prime
j := CnvtoBASE(dgt,n+k,pr);
repeat
fac := 1;
s := 1;
repeat
fac *= pr;
dec(j);
s += fac;
until j<= 0;
with pdf[k] do
Begin
pfSumOfDivs *= s;
pfRemain := pfRemain DIV fac;
end;
j := IncByOneInBase(dgt,pr);
k += pr;
until k >= SizePrDeFe;
end;
 
Line 1,073 ⟶ 1,042:
var
dgt:tDigits;
i,j,k,pr,fac,n,MaxP : Uint64;
begin
n := pdfOfs;
Line 1,083 ⟶ 1,052:
with pdf[i] do
Begin
pfDivCnt := 1;
pfSumOfDivs := 1;
pfRemain := n+i;
pfMaxIdx := 0;
pfpotPrimIdx[0] := 0;
pfpotMax[0] := 0;
end;
end;
Line 1,100 ⟶ 1,065:
begin
j := BsfQWord(n+i);
pfMaxIdx := 1;
pfpotPrimIdx[0] := 0;
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;
i += 2;
Line 1,139 ⟶ 1,099:
 
//no need to use higher primes
if pr*pr > n+SizePrDeFemaxP then
BREAK;
 
CalcSumOfDivs(pdf,dgt,n,k,pr);
//j is power of prime
j := CnvtoBASE(dgt,n+k,pr);
repeat
with pdf[k] do
Begin
pfpotPrimIdx[pfMaxIdx] := i;
pfpotMax[pfMaxIdx] := j;
pfDivCnt *= j+1;
fac := pr;
repeat
pfRemain := pfRemain DIV pr;
dec(j);
fac *= pr;
until j<= 0;
pfSumOfDivs *= (fac-1)DIV(pr-1);
inc(pfMaxIdx);
k += pr;
j := IncByBaseInBase(dgt,pr);
end;
until k >= SizePrDeFe;
until false;
 
Line 1,171 ⟶ 1,112:
j := pfRemain;
if j <> 1 then
begin
pfSumOFDivs *= (j+1);
pfDivCnt *=2;
end;
end;
end;
Line 1,251 ⟶ 1,189:
repeat
inc(prmEndIdx);
until smallprimes[prmEndIdx] > trunc(exp(lnsqrt(Limit)/2.32));
writeln('LIMIT = ',Numb2USA(IntToStr(LIMIT)));
writeln(prmEndIdx:10,smallprimes[prmEndIdx]:10);
writeln('factor beyond LIMIT ',LIMIT_mul);
 
n := 0;
Line 1,271 ⟶ 1,210:
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)));
Line 1,283 ⟶ 1,223:
if n = lim then
Begin
writeln(Numb2USA(IntToStr(lim)):10,Numb2USA(IntToStr(k)):10);
lim *= 10;
end;
Line 1,293 ⟶ 1,233:
if n = lim then
Begin
writeln(Numb2USA(IntToStr(lim)):10,Numb2USA(IntToStr(k)):10);
lim += LIMIT DIV 10;
end;
k += 1-pUntouch[n];
end;
end.</lang>
</lang>
{{out}}
<pre>
TIO.RUN
LIMIT = 5,000,000
runtime for n<= LIMIT 0.070 s
331 2237
Check the rest 63,000,000
factor beyond LIMIT 104
runtime 4.214 s
runtime for n<= LIMIT 0.272 s
Check the rest 515,000,000
runtime 19.052 s
10 2
100 5
1,000 1000 89
10,000 10000 1,212 1212
100,000 100000 13,863 13863
200,0001000000 28,572150232
300,0001500000 43,514227297
400,0002000000 58,459305290
500,0002500000 73,565383422
600,0003000000 88,828462110
700,0003500000 104,061 540769
800,0004000000 119,302 619638
900,0004500000 134,757 698504
5000000 777672
1,000,000 150,232
 
//url=https://math.dartmouth.edu/~carlp/uupaper3.pdf
100000 13863
200000 28572
300000 43515
400000 58459
500000 73565
600000 88828
700000 104062
800000 119302
900000 134758
1000000 150232
2000000 305290
3000000 462110
4000000 619638
5000000 777672
6000000 936244
7000000 1095710
Anonymous user
Cookies help us deliver our services. By using our services, you agree to our use of cookies.