Zumkeller numbers: Difference between revisions

Content added Content deleted
(added extra oeis link. Show some partitions see 60 with 17 or 180 with 375 possibilities)
(→‎{{header|Pascal}}: New pre checked gets 95% .Odd til 140,000,000 287073 95.985 % 133.663 s instead of 3933.063 s)
Line 2,724: Line 2,724:
=={{header|Pascal}}==
=={{header|Pascal}}==
modified [[practical numbers]] and [[Factors_of_an_integer#using_Prime_decomposition]]
modified [[practical numbers]] and [[Factors_of_an_integer#using_Prime_decomposition]]
Now using the trick, that one partition sum must include n.
Now using the trick, that one partition sum must include n and that in most cases summing up from top to buttom gets the right sum of divs with no mental acrobatics.
<lang pascal>program zumkeller;
<lang pascal>program zumkeller;
//https://oeis.org/A083206/a083206.txt
{$IFDEF FPC}
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
Line 2,814: Line 2,815:
if Not(isSmallPrimesInited) then
if Not(isSmallPrimesInited) then
InitSmallPrimes;
InitSmallPrimes;
//already done
if res.pfNum = n then
EXIT;
res.pfNum := n;
res.pfNum := n;
cnt := 0;
cnt := 0;
Line 3,091: Line 3,095:
i: NativeInt;
i: NativeInt;
Begin
Begin
str(pD.pfNum,s);
str(pD.pfNum:10,s);
result := s+' : ';
result := s+' : ';
For i := 0 to pD.pfCnt-1 do
For i := 0 to pD.pfCnt-1 do
Line 3,111: Line 3,115:
type
type
tSol = array of Uint32;
tSol = array of Uint32;
tpHasSum =pByte;
var
var
HasSum: array of byte;
HasSum: array of byte;
Line 3,117: Line 3,122:
T0: Int64;
T0: Int64;
count: NativeInt;
count: NativeInt;
checked :Boolean;


function isZumKeller(const Divs:tDivisors;HalfSumOfDivs_N: NativeInt):boolean;
function isZumKeller(HalfSumOfDivs_N: NativeUInt;const Divs:tDivisors):boolean;
//mark sum and than shift by next divisor == add
//mark sum and than shift by next divisor == add
//for practical numbers every sum must be marked
//for practical numbers every sum must be marked
//modified for zumkeller
//modified for zumkeller
var
var
hs0,hs1: pByte;
hs0: tpHasSum;
pU64 :Uint64;
idx, j, i,maxlimit : NativeInt;
pDivs : tpDivisor;
j,maxlimit : NativeInt;
i: NativeUInt;
begin
begin
idx :=0;
pDIvs := @Divs[0];
maxlimit := 0;
maxlimit := 0;
j:= HalfSumOfDivs_N;

repeat
while idx <=High(Divs) do
IF maxLimit>=j then
begin
IF maxLimit>=HalfSumOfDivs_N then
EXIT(true);
EXIT(true);
i := Divs[idx];
i := pDivs[0];
IF i > maxlimit+1 then
IF i > maxlimit+1 then
break;
break;
maxlimit += i;
maxlimit += i;
inc(idx);
inc(pDivs);
end;
until false;
if Length(HasSum) < HalfSumOfDivs_N+10 then
setlength(HasSum,HalfSumOfDivs_N+10);


i := ((HalfSumOfDivs_N-1) DIV 8 +1)*8;
if High(HasSum) < i then
Begin
setlength(HasSum,0);
setlength(HasSum,i+8);
end;
hs0 := @HasSum[0];
hs0 := @HasSum[0];
fillchar(hs0[0],maxLimit+1,#1); //minimal 0+1
j := (maxLimit SHR 3) shl 3;
fillchar(hs0[maxLimit+1],HalfSumOfDivs_N+10-(maxLimit+2),#0); //minimal 0
fillQword(HasSum[(maxLimit SHR 3) shl 3],(i-j) shr 3+1,0);
fillchar(HasSum[0],maxLimit+1,#1); //minimal 0


repeat
while idx <= High(Divs) do
j:= HalfSumOfDivs_N;
begin
if hs0[HalfSumOfDivs_N] <> 0 then
if hs0[j] <> 0 then
EXIT(true);
EXIT(true);
i := Divs[idx];
i := pDivs[0];
if maxLimit+i > HalfSumOfDivs_N then
if maxLimit > j-i then
BREAK;
BREAK;
hs1 := @hs0[i];
inc(pDivs);
//next sum
//next sum
For j := maxlimit downto 0 do
j := MaxLimit;
hs1[j] := hs1[j] OR hs0[j];
hs0 := @hs0[j+1];
repeat
dec(hs0);
dec(j);
hs0[i] := hs0[i] OR hs0[0];
until j< 0;
maxlimit += i;
maxlimit += i;
inc(idx);
until false;
end;


repeat
while idx <= High(Divs) do
j:= HalfSumOfDivs_N;
begin
i := Divs[idx];
if hs0[j] <> 0 then
maxlimit := HalfSumOfDivs_N-i;
hs1 := @hs0[i];
//next sum
For j := maxlimit downto 0 do
hs1[j] := hs1[j] OR hs0[j];
if hs0[HalfSumOfDivs_N] <> 0 then
EXIT(true);
EXIT(true);
inc(idx);
i := pDivs[0];
end;
j := j-i;
if j <0 then
result := false;
Exit(false);
inc(pDivs);
//next sum
hs0 := @hs0[j+1];
repeat
dec(hs0);
dec(j);
hs0[i] := (hs0[0] OR hs0[i]);
until j< 0;
until false;
end;
end;


function GetZumKeller(n: Uint32;var primeDecomp:tPrimefac): boolean;
function GetZumKeller(n: Uint32;var primeDecomp:tPrimefac): boolean;
var
var
pDivs:tpDivisor;
sum,le: NativeUInt;
sum,tmp,SoD : Int64;
h,i,j,min_j: NativeInt;

begin
begin
checked := true;
PrimeDecomposition(n,primeDecomp);
PrimeDecomposition(n,primeDecomp);
sum := SumOfDiv(primeDecomp);
sum := SumOfDiv(primeDecomp);
//sum must be even and n not deficient
//sum must be even and n not deficient
if(sum<2*N) THEN
if Odd(sum) or (sum<2*n) THEN
EXIT(false);
EXIT(false);
if odd(sum) then
//if Odd(n) then Exit(Not(odd(sum)));// to be tested

EXIT(false);
if odd(n) then
if Not(odd(n)) then
EXIT(Not(ODD(sum)));
if ((n mod 18) in [6,12]) then
EXIT(true);
//only even numbers go here

if (n mod 18) in[6,12] then
Exit(true);
//Now one needs to get the divisors
//Now one needs to get the divisors
GetDivs(primeDecomp,Divs);
GetDivs(primeDecomp,Divs);
pDivs := @Divs[0];

//erase divisor = n, one partition includes n
//erase divisor = n, one partition includes n
SoD := SumOfDiv(primeDecomp) shr 1-n;
setlength(Divs,Length(DIVS)-1);
If SoD < 2 then //0,1 is always true
sum := sum shr 1 -n;
Exit(true);
sum := SoD;

//search minmal index to reach sum of divisors
min_j := 0;
if sum > 0 then
Begin
repeat
sum -= pDivs[min_j];
inc(min_j);
until sum <= 0;
IF sum = 0 then
EXIT(true);
dec(min_j,2);
end;

//only 2 loops as start point of search
j := DivCount(primeDecomp)-2;
repeat
i := j;
while (i >= min_j) AND (pDivs[i] > SoD) do
dec(i);
while i >= min_j do
Begin
sum :=SoD-pDivs[i];
if sum = 0 then
EXIT(TRUE);
h:= i-1;
while (h>=0) AND (sum < pDivs[h]) do
dec(h);
while h >= 0 do
Begin
tmp := sum-pDivs[h];
if tmp = 0 then
EXIT(TRUE);
if tmp > 0 then
sum :=tmp;
dec(h);
end;
dec(i);
end;
dec(j);
until j > min_j;


checked := false;
result := isZumKeller(Divs,sum);
sum := SoD;
result := isZumKeller(sum,Divs);
end;
end;


Line 3,207: Line 3,278:
isZK:=GetZumKeller(n,primeDecomp);
isZK:=GetZumKeller(n,primeDecomp);
writeln(n:10,' SumOfDivs ',SumOfDiv(primeDecomp):11,' CountOfDivs ',DivCount(primeDecomp),' isZk ',isZK);
writeln(n:10,' SumOfDivs ',SumOfDiv(primeDecomp):11,' CountOfDivs ',DivCount(primeDecomp),' isZk ',isZK);
OutPots(primeDecomp);
writeln(OutPots(primeDecomp));
end;
end;


Line 3,230: Line 3,301:
var
var
primeDecomp : tPrimeFac;
primeDecomp : tPrimeFac;
n: NativeInt;
n,checkCnt: NativeInt;
BEGIN
BEGIN
setlength(HasSum,1);
setlength(HasSum,1);
T0 := GetTickCount64;
T0 := GetTickCount64;
writeln('Count of zumkeller numbers up to 100,000');
writeln('Count of zumkeller numbers up to 1,000,000');
n := 1;
n := 1;
count :=0;
count :=0;
Line 3,240: Line 3,311:
if GetZumKeller(n,primeDecomp) then
if GetZumKeller(n,primeDecomp) then
begin
begin
IF checked then
inc(checkCnt);
inc(count);
inc(count);
end;
end;
inc(n,1);
inc(n,1);
until n > 100*1000;
until n > 1000*1000;
writeln(count:10,n-1:10);
writeln(count:10,n-1:10,' pre checked ',checkCnt);
T0 := GetTickCount64-T0;
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
writeln('runtime ',T0/1000:0:3,' s');
Line 3,268: Line 3,341:
n := 1;
n := 1;
count :=0;
count :=0;
checkCnt := 0;
setlength(sol,40);
setlength(sol,40);
writeln(#10,'First 40 odd zumkeller numbers');
writeln(#10,'First 40 odd zumkeller numbers');
Line 3,275: Line 3,349:
sol[count] := n;
sol[count] := n;
inc(count);
inc(count);
if checked then
inc(checkCnt);
end;
end;
inc(n,2);
inc(n,2);
until count >= length(sol);
until count>=length(sol);
OutSol(sol,10,8);
OutSol(sol,10,8);
writeln(count:10,n-2:10,' pre checked ',checkCnt);

T0 := GetTickCount64-T0;
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
writeln('runtime ',T0/1000:0:3,' s');
Line 3,318: Line 3,396:
{{out}}
{{out}}
<pre>TIO.RUN
<pre>TIO.RUN
Count of zumkeller numbers up to 100,000
Count of zumkeller numbers up to 1,000,000
23051 100000
229026 1000000 pre checked 228154
runtime 0.064 s
runtime 3.351 s


First 220 zumkeller
First 220 zumkeller
Line 3,353: Line 3,431:
1145529 1162161 1198197 1224531 1270269 1307691 1324323 1378377
1145529 1162161 1198197 1224531 1270269 1307691 1324323 1378377


runtime 0.145 s
runtime 0.147 s
3492 SumOfDivs 8918 CountOfDivs 18 isZk FALSE
3492 SumOfDivs 8918 CountOfDivs 18 isZk FALSE
3492 : 2^2*3^2*97
14184 SumOfDivs 38610 CountOfDivs 24 isZk FALSE
14184 SumOfDivs 38610 CountOfDivs 24 isZk FALSE
14184 : 2^3*3^2*197
58896 SumOfDivs 165230 CountOfDivs 30 isZk FALSE
58896 SumOfDivs 165230 CountOfDivs 30 isZk FALSE
58896 : 2^4*3^2*409
236448 SumOfDivs 673218 CountOfDivs 36 isZk FALSE
236448 SumOfDivs 673218 CountOfDivs 36 isZk FALSE
236448 : 2^5*3^2*821
954432 SumOfDivs 2737358 CountOfDivs 42 isZk FALSE
954432 SumOfDivs 2737358 CountOfDivs 42 isZk FALSE
954432 : 2^6*3^2*1657
2549700 SumOfDivs 7994714 CountOfDivs 54 isZk FALSE
2549700 SumOfDivs 7994714 CountOfDivs 54 isZk FALSE
2549700 : 2^2*3^2*5^2*2833
10884600 SumOfDivs 36560160 CountOfDivs 72 isZk FALSE
10884600 SumOfDivs 36560160 CountOfDivs 72 isZk FALSE
10884600 : 2^3*3^2*5^2*6047
runtime 0.079 s
runtime 0.087 s
Real time: 0.446 s CPU share: 99.21 %
Real time: 3.738 s CPU share: 99.54 %
</pre>
NEW: Much faster: tested count of odd zumkeller numbers ti 1,000,000,000
<pre>count of odd zumkeller numbers up to 1,000,000,000
n count count/ time
prechecked
10000000 20637 96.860 % 2.839 s
20000000 41463 96.698 % 7.519 s
30000000 62002 96.637 % 13.391 s
40000000 82513 96.474 % 20.381 s
50000000 103001 96.351 % 28.329 s
60000000 123482 96.283 % 37.077 s
70000000 143933 96.232 % 46.573 s
80000000 164412 96.202 % 56.877 s
90000000 184816 96.165 % 67.944 s
100000000 205283 96.152 % 79.329 s
150000000 307454 95.956 % 148.857 s
200000000 409569 95.835 % 239.997 s
250000000 511807 95.796 % 344.792 s
300000000 613721 95.792 % 478.399 s
400000000 817456 95.743 % 840.497 s
500000000 1021438 95.720 % 1291.123 s
600000000 1225716 95.708 % 1820.230 s
700000000 1430227 95.683 % 2461.841 s
800000000 1634217 95.625 % 3263.611 s
900000000 1838186 95.569 % 4161.572 s
1000000000 2042196 95.526 % 5134.790 s
2042196 999999999 pre checked 1950822
</pre>
</pre>
Tested odd zumkeller numbers up to 500,000,000 aka 5*10^8 , which takes 14h <BR>
Tested odd zumkeller numbers up to 500,000,000 aka 5*10^8 , which takes 14h <BR>
Retrieved from "https://rosettacode.org/wiki/Zumkeller_numbers"