Jump to content

Zumkeller numbers: Difference between revisions

Zumkeller numbers (view source)

Revision as of 12:27, 17 May 2021

4,209 bytes added ,  3 years ago
→‎{{header|Pascal}}: using shortcut for even n MOD 18 = 6,12 and for odd numbers for speed up
(→‎{{header|Pascal}}: using shortcut for even n MOD 18 = 6,12 and for odd numbers for speed up)
Line 2,705:
<lang pascal>program zumkeller;
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF}
 
uses
sysutils
{$IFDEF WINDOWS},Windows{$ENDIF}
,Windows
{$ENDIF}
;
//######################################################################
//prime decomposition
type
tItem = Uint32;
tDivisors = array of tItem;
tpDivisor = pUint32;
 
tPot = record
potPrim,
Line 2,730 ⟶ 2,733:
end;
 
 
type
tSmallPrimes = array[0..6541] of Word;
tItem = NativeUint;
tDivisors = array of tItem;
tpDivisor = pNativeUint;
//##########################################################
//factors of integer via primedecomposition
var
SmallPrimes: tSmallPrimes;
isSmallPrimesInited : boolean = false;
primeDecomp: tprimeFac;
 
procedure InitSmallPrimes;
Line 2,770:
inc(pr,2);
until pr > 1 shl 16 -1;
isSmallPrimesInited:= true;
end;
 
function DivCount(const primeDecomppD:tprimeFac):NativeUInt;inline;
begin
result := primeDecomppD.pfDivCnt;
end;
 
function SumOfDiv(const primeDecomppD:tprimeFac):Uint64;inline;
begin
result := primeDecomppD.pfSumOfDivs;
end;
 
procedure PrimeDecomposition(n:Uint32;var res:tprimeFac);
Label
Finished;
var
DivSum,fac:Uint64;
i,pr,cnt,DivCnt,quot{to minimize divisions} : NativeUintUint32;
Begin
if Not(isSmallPrimesInited) then
InitSmallPrimes;
res.pfNum := n;
cnt := 0;
Line 2,802 ⟶ 2,807:
end
else
Begin
pr := 2;
IF 2*2>n then
GOTO Finished;
quot := n div 2;
IF 2*quot = n then
Begin
with res do
with pfPrims[Cnt] do
Begin
potPrim := 2;
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);
end;
inc(Cnt);
end;
 
pr := 3;
IF 3*3>n then
GOTO Finished;
quot := n div 3;
IF 3*quot = n then
begin
with res do
Begin
with pfPrims[Cnt] do
Begin
potPrim := 3;
potMax := 0;
fac := 3;
repeat
n := quot;
quot := quot div 3;
inc(potMax);
fac *= 3;
until pr*quot <> n;
DivCnt *= (potMax+1);
DivSum *= (fac-1) DIV (3-1);
end;
end;
inc(Cnt);
end;
 
pr := 5;
IF 5*5>n then
GOTO Finished;
quot := n div 5;
IF 5*quot = n then
begin
with res do
Begin
with pfPrims[Cnt] do
Begin
potPrim := 5;
potMax := 0;
fac := 5;
repeat
n := quot;
quot := quot div 5;
inc(potMax);
fac *= 5;
until pr*quot <> n;
DivCnt *= (potMax+1);
DivSum *= (fac-1) DIV (5-1);
end;
end;
inc(Cnt);
end;
 
pr := 7;
IF 7*7>n then
GOTO Finished;
quot := n div 7;
IF 7*quot = n then
begin
with res do
Begin
with pfPrims[Cnt] do
Begin
potPrim := 7;
potMax := 0;
fac := 7;
repeat
n := quot;
quot := quot div 7;
inc(potMax);
fac *= 7;
until pr*quot <> n;
DivCnt *= (potMax+1);
DivSum *= (fac-1) DIV (7-1);
end;
end;
inc(Cnt);
end;
 
pr := 11;
IF 11*11>n then
GOTO Finished;
quot := n div 11;
IF 11*quot = n then
begin
with res do
Begin
with pfPrims[Cnt] do
Begin
potPrim := 11;
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);
end;
end;
inc(Cnt);
end;
 
pr := 13;
IF 13*13>n then
GOTO Finished;
quot := n div 13;
IF 13*quot = n then
begin
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);
end;
end;
inc(Cnt);
end;
i := 6;
 
repeat
pr := SmallPrimes[i];
 
IF pr*pr>n then
Break;
Line 2,829 ⟶ 2,991:
inc(i);
until false;
Finished:
//a big prime left over?
IF n > 1 then
Line 2,842 ⟶ 3,005:
DivSum *= n+1;
end;
end;
res.pfCnt:= cnt;
res.pfDivCnt := DivCnt;
Line 2,865 ⟶ 3,029:
end;
 
procedure GetDivs(varconst primeDecomppD:tprimeFac;var Divs:tDivisors);
var
pDivs : tpDivisor;
i,len,j,l,p,pPot,k: NativeInt;
Begin
i := DivCount(primeDecomppD);
setlength(Divs,0);
setlength(Divs,i);
Line 2,877 ⟶ 3,041:
len := 1;
l := 1;
For i := 0 to primeDecomppD.pfCnt-1 do
with primeDecomppD.pfPrims[i] do
Begin
//Multiply every divisor before with the new primefactors
Line 2,899 ⟶ 3,063:
InsertSort(pDivs,0,len-1);
end;
//factors of integer via primedecomposition
//##########################################################
 
 
function OutPots(const pD:tprimeFac):Ansistring;
var
s: String;
i: NativeInt;
Begin
str(pD.pfNum,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);
result += '^'+s;
end;
end;
end;
//prime decomposition
//######################################################################
type
tSol = array of Uint32;
var
HasSum: array of byte;
Divs:tDivisors;
sol : tSol;
CheckDivCount: array[0..255] of Uint32;
DIVCOUNTLIMIT T0: Int32Int64;
count: NativeInt;
 
function isZumKeller(varconst Divs:tDivisors;MiddleHalfSumOfDivs_N: NativeInt):boolean;
//mark sum and than shift by next divisor == add
//for practical numbers every sum must be marked
Line 2,917 ⟶ 3,104:
idx, j, i,maxlimit : NativeInt;
begin
hs0idx := @HasSum[0];
hs0[0] := 1; //empty set
maxlimit := 0;
 
for idx := 0 to High(Divs) do
while idx <=High(Divs) do
begin
ifIF hs0[middle] <maxLimit> 0=HalfSumOfDivs_N then
EXIT(true);
i := Divs[idx];
IF i :=> Divs[idx];maxlimit+1 then
break;
if maxLimit+i > middle then
maxlimit :+= middle-i;
hs1 := @hs0[i]inc(idx);
end;
if Length(HasSum) < HalfSumOfDivs_N+10 then
setlength(HasSum,HalfSumOfDivs_N+10);
 
hs0 := @HasSum[0];
fillchar(hs0[0],maxLimit+1,#1); //minimal 0+1
fillchar(hs0[maxLimit+1],HalfSumOfDivs_N+10-(maxLimit+2),#0); //minimal 0
 
while idx <= High(Divs) do
begin
if hs0[HalfSumOfDivs_N] <> 0 then
EXIT(true);
i := Divs[idx];
if maxLimit+i > HalfSumOfDivs_N then
BREAK;
hs1 := @hs0[i];
//next sum
For j := maxlimit downto 0 do
hs1[j] := hs1[j] OR hs0[j];
maxlimit += i;
inc(idx);
end;
 
while idx <= High(Divs) do
begin
i := Divs[idx];
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);
inc(idx);
end;
result := false;
end;
 
function GetZumKeller(n: Uint32;var primeDecomp:tPrimefac): boolean;
var
sum,le: NativeUInt;
Line 2,944 ⟶ 3,161:
sum := SumOfDiv(primeDecomp);
//sum must be even and n not deficient
if Odd(sum) or (sum<2*N) THEN
EXIT(false);
if odd(sum) then
EXIT(false);
if odd(n) then
EXIT(Not(ODD(sum)));
//only even numbers go here
if (n mod 18) in[6,12] then
Exit(true);
//Now one needs to get the divisors
le := DivCount(primeDecomp);
GetDivs(primeDecomp,Divs);
//erase divisor = n, one partition includes n
setlength(Divs,Length(DIVS)-1);
sum := sum shr 1 -n;
 
le := sum+10;
if le > length(HasSum) then
Begin
setlength(HasSum,0);
setlength(HasSum,le);
end
else
fillChar(HasSum[0],le,#0);
result := isZumKeller(Divs,sum);
end;
 
procedure CheckSpecial(n:Uint32;var primeDecomp:tPrimeFac);
var
isZK : boolean;
Begin
isZK:=GetZumKeller(n,primeDecomp);
writeln(n:10,' SumOfDivs ',SumOfDiv(primeDecomp):11,' CountOfDivs ',DivCount(primeDecomp),' isZk ',isZK);
OutPots(primeDecomp);
end;
 
procedure OutSol(const sol:array of Uint32tsol;colWidth,ColCount,limit:NativeInt);
var
i,col: NativeInt;
Begin
col := 0ColCount;
For i := 0 to limit-1High(sol) do
begin
write(' ',Sol[i]:colWidth-1);
incdec(col);
if col = ColCount0 then
beginBegin
writeln;
col := 0ColCount;
end;
end;
writeln;
end;
 
const
ColCnt = 20;
MAX = 220;
 
var
T0primeDecomp : Int64tPrimeFac;
sol n: array of Uint32NativeInt;
BEGIN
n,limit,count: NativeInt;
setlength(HasSum,1);
begin
InitSmallPrimes;
setlength(HasSum,31);
setlength(sol,MAX+1);
DIVCOUNTLIMIT := 3;
T0 := GetTickCount64;
writeln('Count of zumkeller numbers up to 100,000');
count := 0;
Limitn := MAX1;
ncount := 20;//99500
writeln('The first ',MAX,' zumkeller numbers');
repeat
if GetZumKeller(n,primeDecomp) then
Beginbegin
inc(count);
end;
inc(n,1);
until n > 100*1000;
writeln(count:10,n-1:10);
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
 
writeln(#10,'First 220 zumkeller');
T0 := GetTickCount64;
n := 1;
count :=0;
setlength(sol,220);
repeat
if GetZumKeller(n,primeDecomp) then
begin
sol[count] := n;
inc(count);
end;
inc(n,1);
until count >= Limitlength(sol);
OutSol(sol,84,10,Limit20);
For nT0 := 0 to 255 doGetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
if CheckDivCount[n] <> 0 then
writeln(n:4,CheckDivCount[n]:10);
 
T0 := GetTickCount64;
writeln('The first odd 40 zumkeller numbers');
n := 1;
count :=0;
Limit := setlength(sol,40);
writeln(#10,'First 40 odd zumkeller numbers');
repeat
if GetZumKeller(n,primeDecomp) then
begin
sol[count] := n;
inc(count);
end;
inc(n,2);
until count >= Limitlength(sol);
OutSol(sol,6,10,Limit8);
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
 
T0 := GetTickCount64;
n := 1;
count :=0;
Limit := setlength(sol,40);
writeln(#10,'TheFirst first40 odd 40 zumkeller numbers not endingmultiple inof 5');
repeat
if GetZumKeller(n,primeDecomp) then
begin
sol[count] := n;
inc(count);
end;
inc(n,2);
if n MODmod 5 = 0 then
inc(n,2);
until count >= Limitlength(sol);
OutSol(sol,810,8,Limit);
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
{count of divisors | Number
18 3492
24 14184
30 58596
36 236448
42 954432
54 2549700
72 10884600
}
CheckSpecial(3492);
CheckSpecial(14184);
CheckSpecial(58896);
CheckSpecial(236448);
CheckSpecial(954432);
CheckSpecial(2549700);
CheckSpecial(10884600);
 
T0 := GetTickCount64;
CheckSpecial(3492,primeDecomp);
CheckSpecial(14184,primeDecomp);
CheckSpecial(58896,primeDecomp);
CheckSpecial(236448,primeDecomp);
CheckSpecial(954432,primeDecomp);
CheckSpecial(2549700,primeDecomp);
CheckSpecial(10884600,primeDecomp);
T0 := GetTickCount64-T0;
writeln('runtime ',T0/1000:0:3,' s');
 
setlength(sol,0);
setlength(HasSum, 0);
//{$IFNDEF UNIX} readln; {$ENDIF}
endEND.</lang>
{{out}}
<pre>TIO.RUN
Count of zumkeller numbers up to 100,000
23051 100000
runtime 0.064 s
 
First 220 zumkeller
6 12 20 24 28 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 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
 
runtime 0.000 s
The first 220 zumkeller numbers
6 12 20 24 28 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 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
 
TheFirst first40 odd 40 zumkeller numbers
945 1575 2205 2835 3465 4095 4725 5355 5775 59855355
5775 5985 6435 6615 6825 7245 7425 7875 80857425 8415 8505 89257875
8085 8415 8505 8925 9135 9555 9765 10395
9135 9555 9765 10395 11655 12285 12705 12915 13545 14175
11655 12285 12705 12915 13545 14175 14805 15015
14805 15015 15435 16065 16695 17325 17955 18585 19215 19305
15435 16065 16695 17325 17955 18585 19215 19305
 
runtime 0.001 s
The first odd 40 zumkeller numbers not ending in 5
81081 153153 171171 189189 207207 223839 243243 261261
279279 297297 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
 
First 40 odd zumkeller numbers not multiple of 5
runtime 0.373 s
81081 153153 171171 189189 207207 223839 243243 261261
3492 SumOfDivs 8918 CountOfDivs 18 isZk FALSE
279279 297297 351351 459459 513513 567567 621621 671517
14184 SumOfDivs 38610 CountOfDivs 24 isZk FALSE
729729 742203 783783 793611 812889 837837 891891 908523
58896 SumOfDivs 165230 CountOfDivs 30 isZk FALSE
960687 999999 1024947 1054053 1072071 1073709 1095633 1108107
236448 SumOfDivs 673218 CountOfDivs 36 isZk FALSE
1145529 1162161 1198197 1224531 1270269 1307691 1324323 1378377
954432 SumOfDivs 2737358 CountOfDivs 42 isZk FALSE
2549700 SumOfDivs 7994714 CountOfDivs 54 isZk FALSE
10884600 SumOfDivs 36560160 CountOfDivs 72 isZk FALSE
 
runtime 0.145 s
Real time: 0.659 s CPU share: 99.05 %
3492 SumOfDivs 8918 CountOfDivs 18 isZk FALSE
14184 SumOfDivs 38610 CountOfDivs 24 isZk FALSE
58896 SumOfDivs 165230 CountOfDivs 30 isZk FALSE
236448 SumOfDivs 673218 CountOfDivs 36 isZk FALSE
954432 SumOfDivs 2737358 CountOfDivs 42 isZk FALSE
2549700 SumOfDivs 7994714 CountOfDivs 54 isZk FALSE
10884600 SumOfDivs 36560160 CountOfDivs 72 isZk FALSE
runtime 0.079 s
Real time: 0.446 s CPU share: 99.21 %
</pre>
Tested odd zumkeller numbers up to 500,000,000 aka 5*10^8 , which takes 14h <BR>
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/292828"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information