Rosetta Code

Zumkeller numbers: Difference between revisions

Page Discussion

Zumkeller numbers (view source)

Revision as of 17:43, 15 July 2021

1,812 bytes removed ,  2 years ago
→‎{{header|Pascal}}: improved sieve for prime decompostion runtime for 1e9 down from 474 s to 170 s tested til 1e11
m (→‎{{header|Pascal}}: inserted TIO.RUN for counting 1..10,000,000)
(→‎{{header|Pascal}}: improved sieve for prime decompostion runtime for 1e9 down from 474 s to 170 s tested til 1e11)
Line 2,755:
1145529 1162161 1198197 1224531 1270269 1307691 1324323 1378377</pre>
=={{header|Pascal}}==
Using sieve for primedecomposition<BR>
modified [[Abundant,_deficient_and_perfect_number_classifications#Pascal]] so that I only need to check abundant numbers.<BR>
1E9 needs 8*1E9 bytes.Couldn't try 2e9.Size for speed ( factor 4(1e4)..6(1e9)<BR>
[[Factors_of_an_integer#using_Prime_decomposition]] which is the most time consuming.<BR>
prime decomposition included to get it run on TIO.RUN<BR>
Now using the trick, that one partition sum must include n and improved recursive search.<BR>
Limit is ~1.2e11
 
<lang pascal>program zumkeller;
//https://oeis.org/A083206/a083206.txt
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
// {$O+,I+}
{$ELSE}
{$APPTYPE CONSOLE}
Line 2,774 ⟶ 2,772:
//######################################################################
//prime decomposition
const
//HCN(86) > 1.2E11 = 128,501,493,120 count of divs = 4096 7 3 1 1 1 1 1 1 1
HCN_DivCnt = 4096;
//stop never ending recursion
RECCOUNTMAX = 100*1000*1000;
DELTAMAX = 1000*1000;
type
tItem = Uint32Uint64;
tDivisors = array [0..1920HCN_DivCnt-1] of tItem;
tpDivisor = pUint32pUint64;
const
 
SizePrDeFe = 12697;//*72 <= 1 or 2 Mb ~ level 2 cache -32kB for DIVS
tPot = record
type
potSoD : Uint64;
tdigits = packed potPrim,record
potMax dgtDgts : array [0..31] of Uint32;
end;
 
//the first number with 11 different divisors =
tprimeFac = record
// 2*3*5*7*11*13*17*19*23*29*31 = 2E11
pfPrims : array[0..15] of tPot;
tprimeFac = packed record
pfSumOfDivs :Uint64;
pfCntpfSumOfDivs,
pfRemain : Uint64; //n div (p[0]^[pPot[0] *...) can handle primes <=821641^2 = 6.7e11
pfDivCnt,
pfNum,pfdummypfpotPrim : array[0..9] of : Uint32UInt32;//+10*4 = 56 Byte
pfpotMax : array[0..9] of byte; //10 = 66
pfMaxIdx : Uint16; //68
pfDivCnt : Uint32; //72
end;
 
tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac;
tPrimes = array[0..65535] of Uint32;
 
type
tSmallPrimes = array[0..6541] of Word;
var
SmallPrimes: tSmallPrimestPrimes;
//######################################################################
isSmallPrimesInited : boolean = false;
//prime decomposition
 
procedure InitSmallPrimes;
//only odd numbers
const
MAXLIMIT = (821641-1) shr 1;
var
pr : array[0..MAXLIMIT] of byte;
pr,testPr,j,maxprimidx: Uint32;
p,j,d,flipflop :NativeUInt;
isPrime : boolean;
Begin
maxprimidx := 0;
SmallPrimes[0] := 2;
fillchar(pr[0],SizeOf(pr),#0);
pr := 3;
p := 0;
repeat
isprime := true;
j := 0;
repeat
testPrp :+= SmallPrimes[j];1
IF sqr(testPr) >until pr[p]= then0;
j := break(p+1)*p*2;
if If pr mod testPr = 0j>MAXLIMIT then
BeginBREAK;
isprimed := false2*p+1;
break;repeat
endpr[j] := 1;
j += d;
until j>MAXLIMIT;
until false;
 
SmallPrimes[1] := 3;
SmallPrimes[2] := 5;
j := 3;
d := 7;
flipflop := 3-1;
p := 3;
repeat
if pr[p] = 0 then
begin
SmallPrimes[j] := d;
inc(j);
until falseend;
d += 2*flipflop;
p+=flipflop;
flipflop := 3-flipflop;
until (p > MAXLIMIT) OR (j>High(SmallPrimes));
end;
 
function OutPots(const pD:tprimeFac;n:NativeInt):Ansistring;
if isprime then
var
s: String[31];
Begin
str(n,s);
result := s+' :';
with pd do
begin
str(pfDivCnt:3,s);
result += s+' : ';
For n := 0 to pfMaxIdx-1 do
Begin
inc(maxprimidx);if n>0 then
SmallPrimes[maxprimidx]: result += pr'*';
str(pFpotPrim[n],s);
result += s;
if pfpotMax[n] >1 then
Begin
str(pfpotMax[n],s);
result += '^'+s;
end;
end;
If pfRemain >1 then
inc(pr,2);
Begin
until pr > 1 shl 16 -1;
str(pfRemain,s);
isSmallPrimesInited:= true;
result += '*'+s;
end;
str(pfSumOfDivs,s);
result += '_SoD_'+s+'<';
end;
end;
 
function CnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):NativeInt;
function DivCount(const pD:tprimeFac):NativeUInt;inline;
//n must be multiple of base
begin
var
result := pD.pfDivCnt;
q,r: Uint64;
end;
i : NativeInt;
Begin
with dgt do
Begin
fillchar(dgtDgts,SizeOf(dgtDgts),#0);
i := 0;
// dgtNum:= n;
n := n div base;
result := 0;
repeat
r := n;
q := n div base;
r -= q*base;
n := q;
dgtDgts[i] := r;
inc(i);
until (q = 0);
 
result := 0;
function SumOfDiv(const pD:tprimeFac):Uint64;inline;
while (result<i) AND (dgtDgts[result] = 0) do
begin
inc(result);
result := pD.pfSumOfDivs;
inc(result);
end;
end;
 
function IncByBaseInBase(var dgt:tDigits;base:NativeInt):NativeInt;
procedure PrimeDecomposition(n:Uint32;var res:tprimeFac);
Label
Finished;
var
q :NativeInt;
DivSum,fac:Uint64;
i,pr,cnt,DivCnt,quot{to minimize divisions} : Uint32;
Begin
with dgt do
if Not(isSmallPrimesInited) then
InitSmallPrimes;
//already done
if res.pfNum = n then
EXIT;
res.pfNum := n;
cnt := 0;
DivCnt := 1;
DivSum := 1;
i := 0;
if n <= 1 then
Begin
withresult res.pfPrims[0]:= do0;
q := dgtDgts[result]+1;
Begin
// potPrim := ninc(dgtNum,base);
if q potMax= base := 1;then
begin
repeat
dgtDgts[result] := 0;
inc(result);
q := dgtDgts[result]+1;
until q <> base;
end;
cntdgtDgts[result] := 1q;
result +=1;
end
elseend;
end;
Begin
 
pr := 2;
procedure SieveOneSieve(var pdf:tPrimeDecompField;n:nativeUInt);
IF 2*2>n then
var
GOTO Finished;
dgt:tDigits;
quot := n div 2;
i, j, k,pr,fac : NativeUInt;
IF 2*quot = n then
begin
//init
for i := 0 to High(pdf) do
with pdf[i] do
Begin
withpfDivCnt res:= do1;
pfSumOfDivs := 1;
with pfPrims[Cnt] do
pfRemain := Beginn+i;
potPrimpfMaxIdx := 20;
potMax := 0;
fac := 2;
repeat
n := quot;
quot := quot div 2;
inc(potMax);
fac *= 2;
until pr*quot <> n;
DivCnt *= (potMax+1);
DivSum *= (fac-1)DIV (2-1);
potSoD := DivSum;
end;
inc(Cnt);
end;
 
//first 2 prmake :=n+i 3;even
i := IF 3*3>n thenAND 1;
repeat
GOTO Finished;
quotwith :=pdf[i] n div 3;do
IF 3*quot =if n+i > 0 then
begin
with res doj := BsfQWord(n+i);
Begin pfMaxIdx := 1;
with pfPrimspfpotPrim[Cnt0] do:= 2;
BeginpfpotMax[0] := j;
potPrimpfRemain := 3(n+i) shr j;
potMaxpfSumOfDivs := 0(1 shl (j+1))-1;
facpfDivCnt := 3j+1;
repeat
n := quot;
quot := quot div 3;
inc(potMax);
fac *= 3;
until pr*quot <> n;
DivCnt *= (potMax+1);
DivSum *= (fac-1) DIV (3-1);
potSoD := DivSum;
end;
end;
i += inc(Cnt)2;
until i end>High(pdf);
 
// i now index in SmallPrimes
pr := 5;
i := 0;
IF 5*5>n then
repeat
GOTO Finished;
//search next prime that is in bounds of sieve
quot := n div 5;
repeat
IF 5*quot = n then
begin inc(i);
if i >= High(SmallPrimes) then
with res do
Begin BREAK;
pr := with pfPrimsSmallPrimes[Cnti] do;
k := Beginpr-n MOD pr;
if (k = pr) potPrimAND :=(n>0) 5;then
potMax k:= 0;
if k < SizePrDeFe fac := 5;then
repeatbreak;
until n := quotfalse;
if i >= High(SmallPrimes) then
quot := quot div 5;
inc(potMax)BREAK;
//no need to use higher fac *= 5;primes
untilif pr*quotpr <> n;+SizePrDeFe then
DivCnt *= (potMax+1)BREAK;
DivSum *= (fac-1) DIV (5-1);
potSoD := DivSum;
end;
end;
inc(Cnt);
end;
 
pr// :=j 7;is power of prime
j := CnvtoBASE(dgt,n+k,pr);
IF 7*7>n then
repeat
GOTO Finished;
quot := nwith divpdf[k] 7;do
IF 7*quot = n then
begin
with res do
Begin
with pfPrimspfpotPrim[CntpfMaxIdx] do:= pr;
BeginpfpotMax[pfMaxIdx] := j;
pfDivCnt potPrim :*= 7j+1;
potMaxfac := 0pr;
fac := 7;repeat
repeatpfRemain := pfRemain DIV pr;
n := quotdec(j);
fac quot :*= quot div 7pr;
until j<= inc(potMax)0;
facpfSumOfDivs *= 7(fac-1)DIV(pr-1);
until pr*quot <> ninc(pfMaxIdx);
DivCnt *= (potMax+1);
DivSum *= (fac-1) DIV (7-1);
potSoD := DivSum;
end;
end;
inc(Cnt)k += pr;
j := IncByBaseInBase(dgt,pr);
end;
until k >= SizePrDeFe;
until false;
 
//correct sum of & count of divisors
pr := 11;
for i := 0 to High(pdf) do
IF 11*11>n then
Begin
GOTO Finished;
quotwith :=pdf[i] n div 11;do
IF 11*quot = n then
begin
withj res:= dopfRemain;
Beginif j <> 1 then
with pfPrims[Cnt] dobegin
BeginpfSumOFDivs *= (j+1);
pfDivCnt potPrim :*= 112;
potMax := 0;
fac := 11;
repeat
n := quot;
quot := quot div 11;
inc(potMax);
fac *= 11;
until pr*quot <> n;
DivCnt *= (potMax+1);
DivSum *= (fac-1) DIV (11-1);
potSoD := DivSum;
end;
end;
inc(Cnt);
end;
end;
end;
//prime decomposition
//######################################################################
procedure Init_Check_rec(const pD:tprimeFac;var Divs,SumOfDivs:tDivisors);forward;
 
var
pr := 13;
{$ALIGN 32}
IF 13*13>n then
PrimeDecompField:tPrimeDecompField;
GOTO Finished;
{$ALIGN 32}
quot := n div 13;
Divs :tDivisors;
IF 13*quot = n then
SumOfDivs : tDivisors;
begin
DivUsedIdx : tDivisors;
with res do
Begin
with pfPrims[Cnt] do
Begin
potPrim := 13;
potMax := 0;
fac := 13;
repeat
n := quot;
quot := quot div 13;
inc(potMax);
fac *= 13;
until pr*quot <> n;
DivCnt *= (potMax+1);
DivSum *= (fac-1) DIV (13-1);
potSoD := DivSum;
end;
end;
inc(Cnt);
end;
i := 6;
 
pDiv :tpDivisor;
repeat
T0: Int64;
pr := SmallPrimes[i];
count,rec_Cnt: NativeInt;
depth : Int32;
finished :Boolean;
 
procedure Check_rek_depth(SoD : Int64;i: NativeInt);
IF pr*pr>n then
var
Break;
sum : Int64;
begin
if finished then
EXIT;
inc(rec_Cnt);
 
WHILE (i>0) AND (pDiv[i]>SoD) do
quot := n div pr;
dec(i);
IF pr*quot = n then
 
with res do
while i >= 0 Begindo
Begin
with pfPrims[Cnt] do
DivUsedIdx[depth] := pDiv[i];
Begin
potPrimsum := prSoD-pDiv[i];
if sum potMax := 0; then
begin
fac := pr;
finished := repeattrue;
n := quotEXIT;
end;
quot := quot div pr;
incdec(potMaxi);
fac *= princ(depth);
if (i>= 0) AND (sum <= until pr*quot <>SumOfDivs[i]) n;then
DivCnt *= Check_rek_depth(potMax+1sum,i);
if finished then
DivSum *= (fac-1)DIV(pr-1);
potSoD := DivSumEXIT;
// DivUsedIdx[depth] := 0;
end;
incdec(Cntdepth);
end;
end;
inc(i);
 
until false;
procedure Out_One_Sol(const pd:tprimefac;n:NativeUInt;isZK : Boolean);
Finished:
var
//a big prime left over?
IFsum n: > 1 thenNativeInt;
Begin
with res do
if n< 7 then
exit;
with pd do
begin
writeln(OutPots(pD,n));
if isZK then
Begin
Init_Check_rec(pD,Divs,SumOfDivs);
with pfPrims[Cnt] do
Check_rek_depth(pfSumOfDivs shr 1-n,pFDivCnt-1);
write(pfSumOfDivs shr 1:10,' = ');
sum := n;
while depth >= 0 do
Begin
potPrimsum :+= nDivUsedIdx[depth];
potMax := 1write(DivUsedIdx[depth],'+');
inc dec(Cntdepth);
DivCnt *= 2;
DivSum *= n+1;
potSoD := DivSum;
end;
write(n,' = ',sum);
end;
end
else
write(' no zumkeller ');
end;
res.pfCnt:= cnt;
res.pfDivCnt := DivCnt;
res.pfSumOfDivs := DivSum;
end;
 
Line 3,100 ⟶ 3,117:
end;
 
procedure GetDivs(const pD:tprimeFac;var Divs,SumOfDivs:tDivisors);
var
pDivs : tpDivisor;
i,len,j,l,p,pPot,k : NativeIntUInt64;
i,len,j,l,p,k: Int32;
Begin
i := DivCount(pD).pfDivCnt;
pDivs := @Divs[0];
pDivs[0] := 1;
len := 1;
l := 1;
For i := 0 towith pD.pfCnt-1 do
Begin
with pD.pfPrims[i] do
For i := 0 to pfMaxIdx-1 do
Begin
begin
//Multiply every divisor before with the new primefactors
//and append them to the list
k := potMax-1pfpotMax[i];
p := potPrimpfpotPrim[i];
pPot :=1;
repeat
Line 3,126 ⟶ 3,145:
end;
dec(k);
until k<=0;
len := l;
end;
p := pfRemain;
//Sort. Insertsort much faster than QuickSort in this special case
If p >1 then
InsertSort(pDivs,0,len-1);
begin
end;
For j := 0 to len-1 do
 
function OutPots(const pD:tprimeFac):Ansistring;
var
s: String;
i: NativeInt;
Begin
str(pD.pfNum:10,s);
result := s+' :';
str(DivCount(pD):5,s);
result += s+' : ';
For i := 0 to pD.pfCnt-1 do
with pD.pfPrims[i] do
Begin
if i>0 then
result += '*';
str(potPrim,s);
result += s;
if potMax >1 then
Begin
str(potMax,s)pDivs[l]:= p*pDivs[j];
result += '^'+sinc(l);
end;
len := l;
end;
// str(SumOfDiv(pD):12,s);result := ' '+s+' ';
For i := 0 to pD.pfCnt-1 do
with pD.pfPrims[i] do
Begin
str(potSod,s);
result += ';'+s;
end;
end;
//prime decomposition
//######################################################################
type
tValue = QWord;//needed to get beyond 845404560
tpValue = pQWord;
tPower = array[0..31] of Uint32;//2^32
 
type
tSol = array of Uint32;
var
sol : tSol;
DivUsedIdx : array[0..1920] of Uint32;
SumOfDivs : array[0..1920] of Uint64;
power,
PowerFac : tPower;
ds : array of tValue;
pDS : tpValue;
Divs :tDivisors;
 
pDiv :tpDivisor;
T0: Int64;
MAX : tValue;
depth : NativeInt;
count,rek_count: NativeInt;
finished :Boolean;
//######################################################################
//calc sum of divs
Procedure CalcPotfactor(prim:tValue);
//PowerFac[k] = (prim^(k+1)-1)/(prim-1) == Sum (i=0..k) prim^i
var
k: tValue;
Pot, //== prim^k
PFac : Int64;
begin
Pot := prim;
PFac := 1;
For k := 0 to High(PowerFac) do
begin
PFac := PFac+Pot;
IF (POT > MAX) then
BREAK;
PowerFac[k] := PFac;
Pot := Pot*prim;
end;
//Sort. Insertsort much faster than QuickSort in this special case
end;
InsertSort(pDivs,0,len-1);
 
pPot := 0;
procedure InitPW(prim:tValue);
For i := 0 to len-1 do
begin
fillchar(power,SizeOf(power),#0);
CalcPotfactor(prim);
end;
 
function NextPotCnt(p: tValue):tValue;
//return the first power <> 0
//power == n to base prim
var
i : tValue;
begin
result := 0;
repeat
i := power[result];
Inc(i);
IF i < p then
BREAK
else
begin
i := 0;
power[result] := 0;
inc(result);
end;
until false;
power[result] := i;
end;
 
procedure Sieve(pDS :TpValue;max: tValue);
var
prim,
actNumber,idx : tValue;
begin
prim := 2;
pDS[1]:= 1;
idx :=1;
actNumber := prim+1;
while 2*idx <= max do
Begin
move(pDS[1],pDS[idx+1],SizeOF(pDS[1])*idx);
idx *=2;
pDS[idx] := actNumber;
actnumber := 2*actnumber+1;
end;
if 2*idx > max then
move(pDS[1],pDS[idx+1],SizeOF(pDS[1])*(max-idx+1));
 
prim := 3;
//sieve with "small" primes
while prim*prim <= MAX do
begin
InitPW(prim)pPot += pDivs[i];
SumOfDivs[i] := pPot;
Begin
//actNumber = actual number = n*prim
actNumber := prim;
idx := prim;
while actNumber <= MAX do
begin
dec(idx);
IF idx > 0 then
pDS[actNumber] *= PowerFac[0]
else
Begin
pDS[actNumber] *= PowerFac[NextPotCnt(prim)+1];
idx := Prim;
end;
inc(actNumber,prim);
end;
end;
//next prime
repeat
inc(prim);
until pDS[prim]= 1;
end;
 
//sieve with "big" primes, only one factor is possible
while 2*prim <= MAX do
begin
actNumber := prim;
idx := prim+1;
while actNumber <= MAX do
begin
pDS[actNumber] *= idx;
inc(actNumber,prim);
end;
 
repeat
inc(prim);
until pDS[prim]= 1;
end;
end;
//calc sum of divs
//######################################################################
 
procedure OutSolInit_Check_rec(const solpD:tsoltprimeFac;colWidthvar Divs,ColCountSumOfDivs:NativeInttDivisors);
var
i,col: NativeInt;
Begin
col := ColCount;
For i := 0 to High(sol) do
begin
write(' ',Sol[i]:colWidth-1);
dec(col);
if col = 0 then
Begin
writeln;
col := ColCount;
end;
end;
writeln;
end;
 
procedure Out_One_Sol(const pd:tprimefac;isZK : Boolean);
var
i : NativeInt;
Begin
i := DivCount(pD);
writeln(OutPots(pD));
if isZK then
Begin
// write(SumOfDiv(pD)shr 1 - pd.pfNum:10,' = ');
write(SumOfDiv(pD)shr 1:10,' = ');
i:= 0;
while DivUsedIdx[i] > 0 do
Begin
write(DivUsedIdx[i],'+');
inc(i);
end;
writeln;
end;
end;
 
procedure out_(count,n:NativeUint);
begin
GetDivs(pD,Divs,SUmOfDivs);
writeln(count:10,n:13,' time ',(gettickcount64-t0)/1000:8:3,' s');
finished := false;
depth := 0;
pDiv := @Divs[0];
end;
 
Line 3,350 ⟶ 3,182:
sum : Int64;
begin
if finished then
inc(rek_Count);
EXIT;
if rec_Cnt >RECCOUNTMAX then
begin
rec_Cnt := -1;
finished := true;
exit;
end;
inc(rec_Cnt);
 
WHILE (i>0) AND (pDiv[i]>SoD) do
dec(i);
 
inc(depth);
while i >= 0 do
Begin
DivUsedIdx[depth] := pDiv[i];
sum := SoD-pDiv[i];
if sum = 0 then
begin
finished := true;
DivUsedIdx[depth+1] := 0;
EXIT;
end;
Line 3,369 ⟶ 3,208:
if finished then
EXIT;
end;
DivUsedIdx[depth] := 0;
end;
dec(depth);
end;
 
function GetZumKeller(n: Uint32NativeUint;var primeDecomppD:tPrimefac): boolean;
var
SoD,sum : Int64;
Div_cnt,i,pracLmt: NativeInt;
begin
rek_countrec_Cnt := 0;
SoD:= pd.pfSumOfDivs;
PrimeDecomposition(n,primeDecomp);
SoD := SumOfDiv(primeDecomp);
 
//sum must be even and n not deficient
if Odd(SoD) or (SoD<2*n) THEN
Line 3,391 ⟶ 3,226:
If SoD < 2 then //0,1 is always true
Exit(true);
//idx [0..DivCount(primeDecomp)-1];
 
Div_cnt := DivCount(primeDecomp)pD.pfDivCnt;
 
if Not(odd(n)) then
Line 3,400 ⟶ 3,234:
 
//Now one needs to get the divisors
GetDivsInit_check_rec(primeDecomppD,Divs,SumOfDivs);
sum := 0;
i := 0;
repeat
sum += Divs[i];
SumOfDivs[i] := sum;
inc(i);
until i >= Div_Cnt;
 
pracLmt:= 0;
Line 3,418 ⟶ 3,245:
Begin
pracLmt := i;
// writeln(n,'_',i,':',Div_Cnt:10,Divs[i+1]:10,SumOfDivs[i]:10);
BREAK;
end;
Line 3,426 ⟶ 3,252:
Exit(true);
end;
// number is practical followed by one big prime
if pracLmt = (Div_Cnt-1) shr 1 then
begin
i := SoD mod Divs[pracLmt+1];
with primeDecomppD do
begin
EXIT((pfPrims[pfCnt-1].potPrim<=i)OR (i<=sum));
if pfRemain > 1 then
EXIT((pfRemain<=i) OR (i<=sum))
else
EXIT((pfpotPrim[pfMaxIdx-1]<=i)OR (i<=sum));
end;
end;
 
Begin
IF Div_cnt <= 1920HCN_DivCnt then
Begin
finished := false;
DivUsedIdx[0]:= 0;
depth := -1;
pDiv := @Divs[0];
Check_rek(SoD,Div_cnt-1);
IF rec_Cnt = -1 then
exit(true);
exit(finished);
end;
Line 3,449 ⟶ 3,278:
 
var
primeDecompOfs,i,n : tPrimeFacNativeUInt;
Max: NativeUInt;
n,i,limit,tot_rek_count: NativeInt;
 
isZK : boolean;
procedure Init_Sieve(n:NativeUint);
BEGIN
//Init Sieve i,oFs are Global
writeln(#10,'First 220 zumkeller');
begin
T0 := GetTickCount64;
ni := 1n MOD SizePrDeFe;
Ofs := (n DIV SizePrDeFe)*SizePrDeFe;
count :=0;
SieveOneSieve(PrimeDecompField,Ofs);
setlength(sol,220);
end;
MAX := 1000;// must be known in advance ...
 
setlength(ds,0);
procedure GetSmall(MaxIdx:Int32);
setlength(ds,Max+1);
var
pDS := @ds[0];
ZK: Array of Uint32;
Sieve(pDS,max);
idx: UInt32;
Begin
If MaxIdx<1 then
EXIT;
writeln('The first ',MaxIdx,' zumkeller numbers');
Init_Sieve(0);
setlength(ZK,MaxIdx);
idx := Low(ZK);
repeat
if pDSGetZumKeller(n,PrimeDecompField[ni] >2*n) then
Begin
if GetZumKeller(n,primeDecomp) then
beginZK[idx] := n;
sol[count] := ninc(idx);
inc(count)end;
endinc(i);
inc(n);
If i > High(PrimeDecompField) then
until count >= length(sol);
begin
OutSol(sol,4,20);
dec(i,SizePrDeFe);
T0 := GetTickCount64-T0;
inc(ofs,SizePrDeFe);
writeln('runtime ',T0/1000:0:3,' s');
SieveOneSieve(PrimeDecompField,Ofs);
end;
until idx >= MaxIdx;
For idx := 0 to MaxIdx-1 do
begin
if idx MOD 20 = 0 then
writeln;
write(ZK[idx]:4);
end;
setlength(ZK,0);
writeln;
writeln;
end;
 
procedure GetOdd(MaxIdx:Int32);
T0 := GetTickCount64;
var
ZK: Array of Uint32;
idx: UInt32;
Begin
If MaxIdx<1 then
EXIT;
writeln('The first odd 40 zumkeller numbers');
n := 1;
Init_Sieve(n);
count :=0;
setlength(solZK,40MaxIdx);
idx := Low(ZK);
writeln(#10,'First 40 odd zumkeller numbers');
MAX := 20000;// must be known in advance ...
setlength(ds,0);
setlength(ds,Max+1);
pDS := @ds[0];
Sieve(pDS,max);
repeat
if pDSGetZumKeller(n,PrimeDecompField[ni] >2*n) then
Begin
if GetZumKeller(n,primeDecomp) then
beginZK[idx] := n;
sol[count] := ninc(idx);
inc(count)end;
endinc(i,2);
inc(n,2);
If i > High(PrimeDecompField) then
until count>=length(sol);
begin
OutSol(sol,10,8);
dec(i,SizePrDeFe);
T0 := GetTickCount64-T0;
inc(ofs,SizePrDeFe);
writeln('runtime ',T0/1000:0:3,' s');
SieveOneSieve(PrimeDecompField,Ofs);
end;
until idx >= MaxIdx;
For idx := 0 to MaxIdx-1 do
begin
if idx MOD (80 DIV 8) = 0 then
writeln;
write(ZK[idx]:8);
end;
setlength(ZK,0);
writeln;
writeln;
end;
 
procedure GetOddNot5(MaxIdx:Int32);
T0 := GetTickCount64;
var
n := 5;
ZK: Array of Uint32;
count :=0;
idx: UInt32;
setlength(sol,40);
Begin
writeln(#10,'First 40 odd zumkeller numbers not multiple of 5');
If MaxIdx<1 then
MAX := 1400000;// must be known in advance ...
EXIT;
setlength(ds,0);
writeln('The first odd 40 zumkeller numbers not ending in 5');
setlength(ds,Max+1);
pDSn := @ds[0]1;
SieveInit_Sieve(pDS,maxn);
setlength(ZK,MaxIdx);
idx := Low(ZK);
repeat
if pDSGetZumKeller(n,PrimeDecompField[ni] >2*n) then
Begin
if GetZumKeller(n,primeDecomp) then
beginZK[idx] := n;
sol[count] := ninc(idx);
inc(count)end;
endinc(i,2);
inc(n,2);
ifIf n mod 5 = 0 then
begin
inc(i,2);
inc(n,2);
end;
until count>=length(sol);
If i > High(PrimeDecompField) then
OutSol(sol,10,8);
begin
dec(i,SizePrDeFe);
inc(ofs,SizePrDeFe);
SieveOneSieve(PrimeDecompField,Ofs);
end;
until idx >= MaxIdx;
For idx := 0 to MaxIdx-1 do
begin
if idx MOD (80 DIV 8) = 0 then
writeln;
write(ZK[idx]:8);
end;
setlength(ZK,0);
writeln;
writeln;
end;
BEGIN
InitSmallPrimes;
 
T0 := GetTickCount64-T0;
GetSmall(220);
writeln('runtime ',T0/1000:0:3,' s');
setlengthGetOdd(sol,040);
GetOddNot5(40);
 
writeln;
//speed test
MAXn := 101;//8996229720;//1;
Init_Sieve(n);
writeln('Start ',n,' at ',i);
T0 := GetTickCount64;
MAX := (n DIV DELTAMAX+1)*DELTAMAX;
count := 0;
repeat
writeln('Count of zumkeller numbers up to ',MAX:12);
T0 := GetTickCount64;
repeat
//no need to preserve data
if GetZumKeller(n,PrimeDecompField[i]) then
setlength(ds,0);
setlength inc(ds,Max+1count);
pDS := @ds[0]inc(i);
Sieve inc(pDS,maxn);
If i > High(PrimeDecompField) then
writeln('Count of zumkeller numbers up to ',MAX:9);
tot_rek_count :=0; begin
dec(i,SizePrDeFe);
count :=0;
i := 0 inc(ofs,SizePrDeFe);
SieveOneSieve(PrimeDecompField,Ofs);
For n := 1 to MAX do
begin
limit := pDS[n]-2*n;
if (limit>=0) and Not Odd(Limit) then
Begin
inc(i);
isZK :=GetZumKeller(n,primeDecomp);
inc(tot_rek_count,rek_count);
if isZK then
inc(count);
end;
enduntil n > MAX;
writeln(n-1:10,' tested ',i,' found ',count:10,' ratio ',count/in:10:7);
MAX += DELTAMAX;
writeln('runtime ',(GetTickCount64-T0)/1000:8:3,' s');
until MAX >10*= 10DELTAMAX;
writeln('runtime ',(GetTickCount64-T0)/1000:8:3,' s');
until MAX > 1000*1000*1000 ;
writeln;
setlength(ds,0);
writeln('Count of recursion 59,641,327 for 8,996,229,720');
setlength(Divs,0);
n := 8996229720;
END.</lang>
Init_Sieve(n);
T0 := GetTickCount64;
Out_One_Sol(PrimeDecompField[i],n,true);
writeln;
writeln('runtime ',(GetTickCount64-T0)/1000:8:3,' s');
END.
</lang>
{{out}}
<pre>
TIO.RUN only 10,000,000
CountThe offirst 220 zumkeller numbers up to 10,000,000
 
10000000 tested 2474422 found 2287889
6 12 20 24 28 30 40 42 48 54 56 60 66 70 78 80 84 88 90 96
runtime 1.860 s
102 104 108 112 114 120 126 132 138 140 150 156 160 168 174 176 180 186 192 198
204 208 210 216 220 222 224 228 234 240 246 252 258 260 264 270 272 276 280 282
294 300 304 306 308 312 318 320 330 336 340 342 348 350 352 354 360 364 366 368
372 378 380 384 390 396 402 408 414 416 420 426 432 438 440 444 448 456 460 462
464 468 474 476 480 486 490 492 496 498 500 504 510 516 520 522 528 532 534 540
544 546 550 552 558 560 564 570 572 580 582 588 594 600 606 608 612 616 618 620
624 630 636 640 642 644 650 654 660 666 672 678 680 684 690 696 700 702 704 708
714 720 726 728 732 736 740 744 750 756 760 762 768 770 780 786 792 798 804 810
812 816 820 822 828 832 834 836 840 852 858 860 864 868 870 876 880 888 894 896
906 910 912 918 920 924 928 930 936 940 942 945 948 952 960 966 972 978 980 984
 
The first odd 40 zumkeller numbers
 
945 1575 2205 2835 3465 4095 4725 5355 5775 5985
Real time: 2.058 s CPU share: 99.07 %
6435 6615 6825 7245 7425 7875 8085 8415 8505 8925
9135 9555 9765 10395 11655 12285 12705 12915 13545 14175
14805 15015 15435 16065 16695 17325 17955 18585 19215 19305
 
The first odd 40 zumkeller numbers not ending in 5
...at home AMD 2200G linux fpc 3.2.2
First 220 zumkeller
12 20 24 30 40 42 48 54 56 60 66 70 78 80 84 88 90 96 102 104
108 112 114 120 126 132 138 140 150 156 160 168 174 176 180 186 192 198 204 208
210 216 220 222 224 228 234 240 246 252 258 260 264 270 272 276 280 282 294 300
304 306 308 312 318 320 330 336 340 342 348 350 352 354 360 364 366 368 372 378
380 384 390 396 402 408 414 416 420 426 432 438 440 444 448 456 460 462 464 468
474 476 480 486 490 492 498 500 504 510 516 520 522 528 532 534 540 544 546 550
552 558 560 564 570 572 580 582 588 594 600 606 608 612 616 618 620 624 630 636
640 642 644 650 654 660 666 672 678 680 684 690 696 700 702 704 708 714 720 726
728 732 736 740 744 750 756 760 762 768 770 780 786 792 798 804 810 812 816 820
822 828 832 834 836 840 852 858 860 864 868 870 876 880 888 894 896 906 910 912
918 920 924 928 930 936 940 942 945 948 952 960 966 972 978 980 984 990 992 996
 
81081 153153 171171 189189 207207 223839 243243 261261 279279 297297
runtime 0.004 s
351351 459459 513513 567567 621621 671517 729729 742203 783783 793611
812889 837837 891891 908523 960687 999999 1024947 1054053 1072071 1073709
1095633 1108107 1145529 1162161 1198197 1224531 1270269 1307691 1324323 1378377
 
Start 1 at 1
First 40 odd zumkeller numbers
Count of zumkeller numbers up to 1000000
945 1575 2205 2835 3465 4095 4725 5355
1000000 tested found 229026 ratio 0.2290258
5775 5985 6435 6615 6825 7245 7425 7875
Count of zumkeller numbers up to 2000000
8085 8415 8505 8925 9135 9555 9765 10395
2000000 tested found 457658 ratio 0.2288289
11655 12285 12705 12915 13545 14175 14805 15015
Count of zumkeller numbers up to 3000000
15435 16065 16695 17325 17955 18585 19215 19305
3000000 tested found 686048 ratio 0.2286826
Count of zumkeller numbers up to 4000000
4000000 tested found 914806 ratio 0.2287014
Count of zumkeller numbers up to 5000000
5000000 tested found 1143521 ratio 0.2287042
Count of zumkeller numbers up to 6000000
6000000 tested found 1372208 ratio 0.2287013
Count of zumkeller numbers up to 7000000
7000000 tested found 1600977 ratio 0.2287110
Count of zumkeller numbers up to 8000000
8000000 tested found 1829932 ratio 0.2287415
Count of zumkeller numbers up to 9000000
9000000 tested found 2058883 ratio 0.2287648
Count of zumkeller numbers up to 10000000
10000000 tested found 2287889 ratio 0.2287889
runtime 1.268 s
//zumkeller number with highest recursion count til 1e11
Count of recursion 59,641,327 for 8,996,229,720
8996229720 : 96 : 2^3*3^2*5*2237*11171_SoD_29253435120<
14626717560 = 36+45+60+72+90+120+180+360+201330+804312+805320+2010780+4021560+1124528715+4498114860+8996229720 = 14626717560
runtime 7.068 s // at home 9.5s
 
Real time: 8.689 s CPU share: 98.74 %
runtime 0.000 s
 
First//at 40home oddtil zumkeller1e11 with 85 numbers notwith recursion multiplecount of> 51e8
9900000000 tested found 2262797501 ratio 0.2285654 recursion 10.479
81081 153153 171171 189189 207207 223839 243243 261261
runtime 48.805 s
279279 297297 351351 459459 513513 567567 621621 671517
Count of zumkeller numbers up to 10000000000
729729 742203 783783 793611 812889 837837 891891 908523
rek_ -1 @ 9998443080 : 96 : 2^3*3^2*5*3041*9133_SoD_32509184760<
960687 999999 1024947 1054053 1072071 1073709 1095633 1108107
1145529 1162161 1198197 1224531 1270269 1307691 1324323 1378377
 
10000000000 tested found 2285655276 ratio 0.2285655 recursion 10.520
runtime 0.028 s
runtime 28.976 s
 
real 40m7,478s user 40m7,039s sys 0m0,057s
Count of zumkeller numbers up to 10
only 4 til 4,512,612,672
10 tested 1 found 1 ratio 1.0000000
out_1e10.txt:104: rek_ -1 @ 584818848 : 72 : 2^5*3^2*1423*1427_SoD_1665413568<
runtime 0.001 s
out_1e10.txt:105: rek_ -1 @ 589754016 : 72 : 2^5*3^2*1429*1433_SoD_1679457780<
Count of zumkeller numbers up to 100
out_1e10.txt:174: rek_ -1 @ 1956249450 : 72 : 2*3^2*5^2*2083*2087_SoD_5260832928<
100 tested 20 found 20 ratio 1.0000000
out_1e10.txt:291: rek_ -1 @ 4512612672 : 84 : 2^6*3^2*2797*2801_SoD_12943833396<
runtime 0.000 s
Count of zumkeller numbers up to 1000
1000 tested 229 found 224 ratio 0.9781659
runtime 0.000 s
Count of zumkeller numbers up to 10000
10000 tested 2420 found 2294 ratio 0.9479339
runtime 0.000 s
Count of zumkeller numbers up to 100000
100000 tested 24570 found 23051 ratio 0.9381766
runtime 0.006 s
Count of zumkeller numbers up to 1000000
1000000 tested 246816 found 229026 ratio 0.9279220
runtime 0.087 s
Count of zumkeller numbers up to 10000000
10000000 tested 2474422 found 2287889 ratio 0.9246155
runtime 1.312 s
Count of zumkeller numbers up to 100000000
100000000 tested 24753373 found 22879037 ratio 0.9242796
runtime 21.881 s
Count of zumkeller numbers up to 1000000000
1000000000 tested 247588034 found 228620758 ratio 0.9233918
runtime 454.140 s // before 2763.437 s
</pre>
 
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/292838"