Untouchable numbers: Difference between revisions

Content added Content deleted
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: Line 889:
There are 13863 untouchable numbers ≤ 100000.</pre>
There are 13863 untouchable numbers ≤ 100000.</pre>
=={{header|Pascal}}==
=={{header|Pascal}}==
modified [[Factors_of_an_integer#using_Prime_decomposition]]<BR>
modified [[Factors_of_an_integer#using_Prime_decomposition]] to calculate only sum of divisors<BR>
see also math.dartmouth.edu/~carlp/uupaper3.pdf with list count up to 1e9
Appended a list of count of untouchable numbers out of math.dartmouth.edu/~carlp/uupaper3.pdf<BR>
Nigel is still right, that this procedure is not the way to get proven results.
<lang pascal>program UntouchableNumbers;
<lang pascal>program UntouchableNumbers;
// gets factors of consecutive integers fast
// 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 900: Line 899:
{$ENDIF}
{$ENDIF}
uses
uses
sysutils,strutils {Numb2USA commatize}
sysutils,strutils
{$IFDEF WINDOWS},Windows{$ENDIF}
{$IFDEF WINDOWS},Windows{$ENDIF}
;
;
const
const
//100*1000*1000;//10*1000*1000;//5*1000*1000;//1000*1000;
//100*1000*1000;//10*1000*1000;//5*1000*1000;//1000*1000;
LIMIT = 1000*1000;
LIMIT = 5*1000*1000;
//384;//164;//110;//64;
//384;//152;//104;//64;
LIMIT_mul = 64;
LIMIT_mul = 104;


//######################################################################
//######################################################################
//gets sum of divisors of consecutive integers fast
//prime decomposition
const
const
SizePrDeFe = 6*8192;//*size of(tprimeFac) =56 byte level I or 2 Mb ~ level 2 cache
SizePrDeFe = 16*8192;//*size of(tprimeFac) =16 byte 2 Mb ~ level 3 cache
type
type
tdigits = array [0..31] of Uint32;
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
tprimeFac = packed record
pfSumOfDivs,
pfSumOfDivs,
pfRemain : Uint64;
pfRemain : Uint64;
pfpotPrimIdx : array[0..11] of word;
pfDivCnt : Uint16;
pfpotMax : array[0..12] of byte;
pfMaxIdx : byte;
end;
end;
tpPrimeFac = ^tprimeFac;
tpPrimeFac = ^tprimeFac;
Line 932: Line 924:


var
var
{$ALIGN 8}
{$ALIGN 16}
SmallPrimes: tPrimes;
{$ALIGN 32}
PrimeDecompField :tPrimeDecompField;
PrimeDecompField :tPrimeDecompField;
{$ALIGN 16}
SmallPrimes: tPrimes;
pdfIDX,pdfOfs: NativeInt;
pdfIDX,pdfOfs: NativeInt;


Line 979: Line 971:
flipflop := 3-flipflop;
flipflop := 3-flipflop;
until (p > MAXLIMIT) OR (j>High(SmallPrimes));
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;
end;


Line 1,054: Line 998:
end;
end;


function IncByBaseInBase(var dgt:tDigits;base:NativeInt):NativeInt;
function IncByOneInBase(var dgt:tDigits;base:NativeInt):NativeInt;
var
var
q :NativeInt;
q :NativeInt;
Line 1,068: Line 1,012:
dgt[result] := q;
dgt[result] := q;
result +=1;
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;
end;


Line 1,073: Line 1,042:
var
var
dgt:tDigits;
dgt:tDigits;
i,j,k,pr,fac,n,MaxP : Uint64;
i,j,k,pr,n,MaxP : Uint64;
begin
begin
n := pdfOfs;
n := pdfOfs;
Line 1,083: Line 1,052:
with pdf[i] do
with pdf[i] do
Begin
Begin
pfDivCnt := 1;
pfSumOfDivs := 1;
pfSumOfDivs := 1;
pfRemain := n+i;
pfRemain := n+i;
pfMaxIdx := 0;
pfpotPrimIdx[0] := 0;
pfpotMax[0] := 0;
end;
end;
end;
end;
Line 1,100: Line 1,065:
begin
begin
j := BsfQWord(n+i);
j := BsfQWord(n+i);
pfMaxIdx := 1;
pfpotPrimIdx[0] := 0;
pfpotMax[0] := j;
pfRemain := (n+i) shr j;
pfRemain := (n+i) shr j;
pfSumOfDivs := (Uint64(1) shl (j+1))-1;
pfSumOfDivs := (Uint64(1) shl (j+1))-1;
// writeln(n+i,' # ',j,' ## ',pfRemain,' ### ',pfSumOfDivs);
pfDivCnt := j+1;
end;
end;
i += 2;
i += 2;
Line 1,139: Line 1,099:


//no need to use higher primes
//no need to use higher primes
if pr*pr > n+SizePrDeFe then
if pr > maxP then
BREAK;
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;
until false;


Line 1,171: Line 1,112:
j := pfRemain;
j := pfRemain;
if j <> 1 then
if j <> 1 then
begin
pfSumOFDivs *= (j+1);
pfSumOFDivs *= (j+1);
pfDivCnt *=2;
end;
end;
end;
end;
end;
Line 1,251: Line 1,189:
repeat
repeat
inc(prmEndIdx);
inc(prmEndIdx);
until smallprimes[prmEndIdx] > trunc(exp(ln(Limit)/2.32));
until smallprimes[prmEndIdx] > trunc(sqrt(Limit));
writeln('LIMIT = ',Numb2USA(IntToStr(LIMIT)));
writeln(prmEndIdx:10,smallprimes[prmEndIdx]:10);
writeln(prmEndIdx:10,smallprimes[prmEndIdx]:10);
writeln(LIMIT_mul);
writeln('factor beyond LIMIT ',LIMIT_mul);


n := 0;
n := 0;
Line 1,271: Line 1,210:
end;
end;
until n >LIMIT;
until n >LIMIT;

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)));
Line 1,283: Line 1,223:
if n = lim then
if n = lim then
Begin
Begin
writeln(Numb2USA(IntToStr(lim)):10,Numb2USA(IntToStr(k)):10);
writeln(lim:10,k:10);
lim *= 10;
lim *= 10;
end;
end;
Line 1,293: Line 1,233:
if n = lim then
if n = lim then
Begin
Begin
writeln(Numb2USA(IntToStr(lim)):10,Numb2USA(IntToStr(k)):10);
writeln(lim:10,k:10);
lim += LIMIT DIV 10;
lim += LIMIT DIV 10;
end;
end;
k += 1-pUntouch[n];
k += 1-pUntouch[n];
end;
end;
end.
end.</lang>
</lang>
{{out}}
{{out}}
<pre>
<pre>
TIO.RUN
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
10 2
100 5
100 5
1,000 89
1000 89
10,000 1,212
10000 1212
100,000 13,863
100000 13863
200,000 28,572
1000000 150232
300,000 43,514
1500000 227297
400,000 58,459
2000000 305290
500,000 73,565
2500000 383422
600,000 88,828
3000000 462110
700,000 104,061
3500000 540769
800,000 119,302
4000000 619638
900,000 134,757
4500000 698504
5000000 777672
1,000,000 150,232


//url=https://math.dartmouth.edu/~carlp/uupaper3.pdf
//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
6000000 936244
7000000 1095710
7000000 1095710
Retrieved from "https://rosettacode.org/wiki/Untouchable_numbers"