Zumkeller numbers: Difference between revisions

Content added Content deleted
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: Line 2,755:
1145529 1162161 1198197 1224531 1270269 1307691 1324323 1378377</pre>
1145529 1162161 1198197 1224531 1270269 1307691 1324323 1378377</pre>
=={{header|Pascal}}==
=={{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>
Now using the trick, that one partition sum must include n and improved recursive search.<BR>
Limit is ~1.2e11

<lang pascal>program zumkeller;
<lang pascal>program zumkeller;
//https://oeis.org/A083206/a083206.txt
//https://oeis.org/A083206/a083206.txt
{$IFDEF FPC}
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
// {$O+,I+}
{$ELSE}
{$ELSE}
{$APPTYPE CONSOLE}
{$APPTYPE CONSOLE}
Line 2,774: Line 2,772:
//######################################################################
//######################################################################
//prime decomposition
//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
type
tItem = Uint32;
tItem = Uint64;
tDivisors = array [0..1920] of tItem;
tDivisors = array [0..HCN_DivCnt-1] of tItem;
tpDivisor = pUint32;
tpDivisor = pUint64;
const

SizePrDeFe = 12697;//*72 <= 1 or 2 Mb ~ level 2 cache -32kB for DIVS
tPot = record
type
potSoD : Uint64;
potPrim,
tdigits = packed record
potMax :Uint32;
dgtDgts : array [0..31] of Uint32;
end;
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;
pfCnt,
pfSumOfDivs,
pfRemain : Uint64; //n div (p[0]^[pPot[0] *...) can handle primes <=821641^2 = 6.7e11
pfDivCnt,
pfNum,pfdummy : Uint32;
pfpotPrim : array[0..9] of UInt32;//+10*4 = 56 Byte
pfpotMax : array[0..9] of byte; //10 = 66
pfMaxIdx : Uint16; //68
pfDivCnt : Uint32; //72
end;
end;


tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac;
tPrimes = array[0..65535] of Uint32;


type
tSmallPrimes = array[0..6541] of Word;
var
var
SmallPrimes: tSmallPrimes;
SmallPrimes: tPrimes;
//######################################################################
isSmallPrimesInited : boolean = false;
//prime decomposition

procedure InitSmallPrimes;
procedure InitSmallPrimes;
//only odd numbers
const
MAXLIMIT = (821641-1) shr 1;
var
var
pr : array[0..MAXLIMIT] of byte;
pr,testPr,j,maxprimidx: Uint32;
p,j,d,flipflop :NativeUInt;
isPrime : boolean;
Begin
Begin
maxprimidx := 0;
SmallPrimes[0] := 2;
SmallPrimes[0] := 2;
fillchar(pr[0],SizeOf(pr),#0);
pr := 3;
p := 0;
repeat
repeat
isprime := true;
j := 0;
repeat
repeat
testPr := SmallPrimes[j];
p +=1
IF sqr(testPr) > pr then
until pr[p]= 0;
break;
j := (p+1)*p*2;
If pr mod testPr = 0 then
if j>MAXLIMIT then
Begin
BREAK;
isprime := false;
d := 2*p+1;
break;
repeat
end;
pr[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);
inc(j);
until false;
end;
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
Begin
inc(maxprimidx);
if n>0 then
SmallPrimes[maxprimidx]:= pr;
result += '*';
str(pFpotPrim[n],s);
result += s;
if pfpotMax[n] >1 then
Begin
str(pfpotMax[n],s);
result += '^'+s;
end;
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;
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;
end;


function IncByBaseInBase(var dgt:tDigits;base:NativeInt):NativeInt;
procedure PrimeDecomposition(n:Uint32;var res:tprimeFac);
Label
Finished;
var
var
q :NativeInt;
DivSum,fac:Uint64;
i,pr,cnt,DivCnt,quot{to minimize divisions} : Uint32;
Begin
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
Begin
with res.pfPrims[0] do
result := 0;
q := dgtDgts[result]+1;
Begin
potPrim := n;
// inc(dgtNum,base);
potMax := 1;
if q = base then
begin
repeat
dgtDgts[result] := 0;
inc(result);
q := dgtDgts[result]+1;
until q <> base;
end;
end;
cnt := 1;
dgtDgts[result] := q;
result +=1;
end
else
end;
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
Begin
with res do
pfDivCnt := 1;
pfSumOfDivs := 1;
with pfPrims[Cnt] do
Begin
pfRemain := n+i;
potPrim := 2;
pfMaxIdx := 0;
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;
end;


pr := 3;
//first 2 make n+i even
IF 3*3>n then
i := n AND 1;
repeat
GOTO Finished;
quot := n div 3;
with pdf[i] do
IF 3*quot = n then
if n+i > 0 then
begin
begin
with res do
j := BsfQWord(n+i);
Begin
pfMaxIdx := 1;
with pfPrims[Cnt] do
pfpotPrim[0] := 2;
Begin
pfpotMax[0] := j;
potPrim := 3;
pfRemain := (n+i) shr j;
potMax := 0;
pfSumOfDivs := (1 shl (j+1))-1;
fac := 3;
pfDivCnt := j+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;
end;
inc(Cnt);
i += 2;
end;
until i >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;
with pfPrims[Cnt] do
pr := SmallPrimes[i];
Begin
k := pr-n MOD pr;
potPrim := 5;
if (k = pr) AND (n>0) then
potMax := 0;
k:= 0;
fac := 5;
if k < SizePrDeFe then
repeat
break;
n := quot;
until false;
if i >= High(SmallPrimes) then
quot := quot div 5;
inc(potMax);
BREAK;
fac *= 5;
//no need to use higher primes
until pr*quot <> n;
if pr*pr > n+SizePrDeFe then
DivCnt *= (potMax+1);
BREAK;
DivSum *= (fac-1) DIV (5-1);
potSoD := DivSum;
end;
end;
inc(Cnt);
end;


pr := 7;
// j is power of prime
j := CnvtoBASE(dgt,n+k,pr);
IF 7*7>n then
repeat
GOTO Finished;
quot := n div 7;
with pdf[k] do
IF 7*quot = n then
begin
with res do
Begin
Begin
with pfPrims[Cnt] do
pfpotPrim[pfMaxIdx] := pr;
Begin
pfpotMax[pfMaxIdx] := j;
potPrim := 7;
pfDivCnt *= j+1;
potMax := 0;
fac := pr;
fac := 7;
repeat
repeat
pfRemain := pfRemain DIV pr;
n := quot;
dec(j);
quot := quot div 7;
fac *= pr;
inc(potMax);
until j<= 0;
fac *= 7;
pfSumOfDivs *= (fac-1)DIV(pr-1);
until pr*quot <> n;
inc(pfMaxIdx);
DivCnt *= (potMax+1);
DivSum *= (fac-1) DIV (7-1);
potSoD := DivSum;
end;
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;
quot := n div 11;
with pdf[i] do
IF 11*quot = n then
begin
begin
with res do
j := pfRemain;
Begin
if j <> 1 then
with pfPrims[Cnt] do
begin
Begin
pfSumOFDivs *= (j+1);
potPrim := 11;
pfDivCnt *=2;
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;
end;
inc(Cnt);
end;
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
Begin
while i >= 0 do
Begin
with pfPrims[Cnt] do
DivUsedIdx[depth] := pDiv[i];
Begin
potPrim := pr;
sum := SoD-pDiv[i];
potMax := 0;
if sum = 0 then
begin
fac := pr;
repeat
finished := true;
n := quot;
EXIT;
end;
quot := quot div pr;
inc(potMax);
dec(i);
fac *= pr;
inc(depth);
until pr*quot <> n;
if (i>= 0) AND (sum <= SumOfDivs[i]) then
DivCnt *= (potMax+1);
Check_rek_depth(sum,i);
if finished then
DivSum *= (fac-1)DIV(pr-1);
potSoD := DivSum;
EXIT;
// DivUsedIdx[depth] := 0;
end;
inc(Cnt);
dec(depth);
end;
end;
end;
inc(i);

until false;
procedure Out_One_Sol(const pd:tprimefac;n:NativeUInt;isZK : Boolean);
Finished:
var
//a big prime left over?
IF n > 1 then
sum : NativeInt;
Begin
with res do
if n< 7 then
exit;
with pd do
begin
writeln(OutPots(pD,n));
if isZK then
Begin
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
Begin
potPrim := n;
sum += DivUsedIdx[depth];
potMax := 1;
write(DivUsedIdx[depth],'+');
inc(Cnt);
dec(depth);
DivCnt *= 2;
DivSum *= n+1;
potSoD := DivSum;
end;
end;
write(n,' = ',sum);
end;
end
else
write(' no zumkeller ');
end;
end;
res.pfCnt:= cnt;
res.pfDivCnt := DivCnt;
res.pfSumOfDivs := DivSum;
end;
end;


Line 3,100: Line 3,117:
end;
end;


procedure GetDivs(const pD:tprimeFac;var Divs:tDivisors);
procedure GetDivs(const pD:tprimeFac;var Divs,SumOfDivs:tDivisors);
var
var
pDivs : tpDivisor;
pDivs : tpDivisor;
i,len,j,l,p,pPot,k: NativeInt;
pPot : UInt64;
i,len,j,l,p,k: Int32;
Begin
Begin
i := DivCount(pD);
i := pD.pfDivCnt;
pDivs := @Divs[0];
pDivs := @Divs[0];
pDivs[0] := 1;
pDivs[0] := 1;
len := 1;
len := 1;
l := 1;
l := 1;
For i := 0 to pD.pfCnt-1 do
with pD do
Begin
with pD.pfPrims[i] do
For i := 0 to pfMaxIdx-1 do
Begin
begin
//Multiply every divisor before with the new primefactors
//Multiply every divisor before with the new primefactors
//and append them to the list
//and append them to the list
k := potMax-1;
k := pfpotMax[i];
p := potPrim;
p := pfpotPrim[i];
pPot :=1;
pPot :=1;
repeat
repeat
Line 3,126: Line 3,145:
end;
end;
dec(k);
dec(k);
until k<0;
until k<=0;
len := l;
len := l;
end;
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
Begin
str(potMax,s);
pDivs[l]:= p*pDivs[j];
result += '^'+s;
inc(l);
end;
end;
len := l;
end;
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;
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
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;
end;
end;
//calc sum of divs
//######################################################################


procedure OutSol(const sol:tsol;colWidth,ColCount:NativeInt);
procedure Init_Check_rec(const pD:tprimeFac;var Divs,SumOfDivs:tDivisors);
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
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;
end;


Line 3,350: Line 3,182:
sum : Int64;
sum : Int64;
begin
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
WHILE (i>0) AND (pDiv[i]>SoD) do
dec(i);
dec(i);

inc(depth);
while i >= 0 do
while i >= 0 do
Begin
Begin
DivUsedIdx[depth] := pDiv[i];
sum := SoD-pDiv[i];
sum := SoD-pDiv[i];
if sum = 0 then
if sum = 0 then
begin
begin
finished := true;
finished := true;
DivUsedIdx[depth+1] := 0;
EXIT;
EXIT;
end;
end;
Line 3,369: Line 3,208:
if finished then
if finished then
EXIT;
EXIT;
end;
DivUsedIdx[depth] := 0;
end;
dec(depth);
end;
end;


function GetZumKeller(n: Uint32;var primeDecomp:tPrimefac): boolean;
function GetZumKeller(n: NativeUint;var pD:tPrimefac): boolean;
var
var
SoD,sum : Int64;
SoD,sum : Int64;
Div_cnt,i,pracLmt: NativeInt;
Div_cnt,i,pracLmt: NativeInt;
begin
begin
rek_count := 0;
rec_Cnt := 0;
SoD:= pd.pfSumOfDivs;
PrimeDecomposition(n,primeDecomp);
SoD := SumOfDiv(primeDecomp);

//sum must be even and n not deficient
//sum must be even and n not deficient
if Odd(SoD) or (SoD<2*n) THEN
if Odd(SoD) or (SoD<2*n) THEN
Line 3,391: Line 3,226:
If SoD < 2 then //0,1 is always true
If SoD < 2 then //0,1 is always true
Exit(true);
Exit(true);
//idx [0..DivCount(primeDecomp)-1];


Div_cnt := DivCount(primeDecomp);
Div_cnt := pD.pfDivCnt;


if Not(odd(n)) then
if Not(odd(n)) then
Line 3,400: Line 3,234:


//Now one needs to get the divisors
//Now one needs to get the divisors
GetDivs(primeDecomp,Divs);
Init_check_rec(pD,Divs,SumOfDivs);
sum := 0;
i := 0;
repeat
sum += Divs[i];
SumOfDivs[i] := sum;
inc(i);
until i >= Div_Cnt;


pracLmt:= 0;
pracLmt:= 0;
Line 3,418: Line 3,245:
Begin
Begin
pracLmt := i;
pracLmt := i;
// writeln(n,'_',i,':',Div_Cnt:10,Divs[i+1]:10,SumOfDivs[i]:10);
BREAK;
BREAK;
end;
end;
Line 3,426: Line 3,252:
Exit(true);
Exit(true);
end;
end;
// number is practical followed by one big prime
//number is practical followed by one big prime
if pracLmt = (Div_Cnt-1) shr 1 then
if pracLmt = (Div_Cnt-1) shr 1 then
begin
begin
i := SoD mod Divs[pracLmt+1];
i := SoD mod Divs[pracLmt+1];
with primeDecomp do
with pD 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;
end;


Begin
Begin
IF Div_cnt <= 1920 then
IF Div_cnt <= HCN_DivCnt then
Begin
Begin
finished := false;
DivUsedIdx[0]:= 0;
depth := -1;
pDiv := @Divs[0];
Check_rek(SoD,Div_cnt-1);
Check_rek(SoD,Div_cnt-1);
IF rec_Cnt = -1 then
exit(true);
exit(finished);
exit(finished);
end;
end;
Line 3,449: Line 3,278:


var
var
primeDecomp : tPrimeFac;
Ofs,i,n : NativeUInt;
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;
n := 1;
i := n 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
repeat
if pDS[n] >2*n then
if GetZumKeller(n,PrimeDecompField[i]) then
Begin
if GetZumKeller(n,primeDecomp) then
begin
ZK[idx] := n;
sol[count] := n;
inc(idx);
inc(count);
end;
end;
inc(i);
inc(n);
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;
n := 1;
Init_Sieve(n);
count :=0;
setlength(sol,40);
setlength(ZK,MaxIdx);
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
repeat
if pDS[n] >2*n then
if GetZumKeller(n,PrimeDecompField[i]) then
Begin
if GetZumKeller(n,primeDecomp) then
begin
ZK[idx] := n;
sol[count] := n;
inc(idx);
inc(count);
end;
end;
inc(i,2);
inc(n,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);
pDS := @ds[0];
n := 1;
Sieve(pDS,max);
Init_Sieve(n);
setlength(ZK,MaxIdx);
idx := Low(ZK);
repeat
repeat
if pDS[n] >2*n then
if GetZumKeller(n,PrimeDecompField[i]) then
Begin
if GetZumKeller(n,primeDecomp) then
begin
ZK[idx] := n;
sol[count] := n;
inc(idx);
inc(count);
end;
end;
inc(i,2);
inc(n,2);
inc(n,2);
if n mod 5 = 0 then
If n mod 5 = 0 then
begin
inc(i,2);
inc(n,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;
T0 := GetTickCount64;
GetSmall(220);
writeln('runtime ',T0/1000:0:3,' s');
setlength(sol,0);
GetOdd(40);
GetOddNot5(40);


writeln;
//speed test
MAX := 10;
n := 1;//8996229720;//1;
Init_Sieve(n);
writeln('Start ',n,' at ',i);
T0 := GetTickCount64;
MAX := (n DIV DELTAMAX+1)*DELTAMAX;
count := 0;
repeat
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(ds,Max+1);
inc(count);
pDS := @ds[0];
inc(i);
Sieve(pDS,max);
inc(n);
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;
end;
end;
until n > MAX;
writeln(n:10,' tested ',i,' found ',count:10,' ratio ',count/i:10:7);
writeln(n-1:10,' tested found ',count:10,' ratio ',count/n:10:7);
MAX += DELTAMAX;
writeln('runtime ',(GetTickCount64-T0)/1000:8:3,' s');
MAX *= 10;
until MAX>10*DELTAMAX;
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}}
{{out}}
<pre>
<pre>
TIO.RUN only 10,000,000
TIO.RUN
Count of zumkeller numbers up to 10,000,000
The first 220 zumkeller numbers

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 40 odd zumkeller numbers not multiple of 5
//at home til 1e11 with 85 numbers with recursion count > 1e8
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>
</pre>


Retrieved from "https://rosettacode.org/wiki/Zumkeller_numbers"