Jordan-Pólya numbers: Difference between revisions

Content added Content deleted
(→‎{{header|Free Pascal}}: Convert using Uint64. misinterpretate 2^53 with 1E53 tse tse tse...)
m (→‎{{header|Free Pascal}}: more comments and commatize)
Line 743: Line 743:
const
const
MaxIdx = 3800;//7279 < 2^62
MaxIdx = 3800;//7279 < 2^62
maxFac = 21;//21!> 2^63
maxFac = 25;//21!> 2^63
type
type
tnum = Uint64;
tnum = Uint64;
Line 749: Line 749:
tFac_mul = packed record
tFac_mul = packed record
fm_num : tnum;
fm_num : tnum;
fm_pow : tpow; //which Factorials where used for fm_num
fm_pow : tpow;
fm_high_idx : word;//memorize highest Factorial Index
fm_high_idx : word;
fm_high_pow : word;// and it's power
fm_high_pow : word;
end;
end;
tpFac_mul = ^tFac_mul;
tpFac_mul = ^tFac_mul;
tFacMulPow = array of tFac_mul;
tFacMulPow = array of tFac_mul;
tFactorial = array[0..maxFac-2] of tnum;

var
var
Factorial: array[0..maxFac-2] of tnum;
FacMulPowGes : tFacMulPow;
FacMulPowGes : tFacMulPow;
Factorial: tFactorial;
LastSearchFor :tFac_mul;
LastSearchFor :tFac_mul;
dblLimit : tnum;
dblLimit : tnum;


function CommatizeUint64(num:Uint64):AnsiString;
procedure QuickSort(var AI: tFacMulPow; ALo, AHi: Int32);
var
var
fromIdx,toIdx :Int32;
Tmp :tFac_mul;
Begin
Pivot : tnum;
str(num,result);
Lo, Hi : Int32;
fromIdx := length(result);
begin
Lo := ALo;
toIdx := fromIdx-1;
Hi := AHi;
if toIdx < 3 then
exit;
Pivot := AI[(Lo + Hi) div 2].fm_num;

toIdx := 4*(toIdx DIV 3)+toIdx MOD 3 +1 ;
setlength(result,toIdx);
repeat
repeat
result[toIdx] := result[FromIdx];
while AI[Lo].fm_num < Pivot do
result[toIdx-1] := result[FromIdx-1];
Inc(Lo);
result[toIdx-2] := result[FromIdx-2];
while AI[Hi].fm_num > Pivot do
Dec(Hi);
result[toIdx-3] := ',';
dec(toIdx,4);
if Lo <= Hi then
begin
dec(FromIdx,3);
until FromIdx<=3;
Tmp := AI[Lo];
AI[Lo] := AI[Hi];
AI[Hi] := Tmp;
Inc(Lo);
Dec(Hi);
end;
until Lo > Hi;
if Hi > ALo then
QuickSort(AI, ALo, Hi) ;
if Lo < AHi then
QuickSort(AI, Lo, AHi) ;
end;
end;


Line 795: Line 789:
fac,
fac,
num : tNum;
num : tNum;
i,j,pow : integer;
FacIdx,pow : integer;
begin
begin
num := fm.fm_num;
num := fm.fm_num;
FacIdx := fm.fm_high_idx;
write(num:20);
write(CommatizeUint64(num):25,' = ');
i := fm.fm_high_idx;

repeat
repeat
j := 0;
pow := 0;
fac := Factorial[i];
fac := Factorial[FacIdx];
while (num>=fac) AND (num mod Fac = 0) do
while (num>=fac) AND (num mod Fac = 0) do
Begin
Begin
num := num DIV Fac;
num := num DIV Fac;
inc(j);
inc(pow);
end;
end;
if j = 0 then
if pow = 0 then
write(' 1')
write(' 1')
else
else
if j = 1 then
if pow = 1 then
write(' ',i+2,'!')
write(' ',FacIdx+2,'!')
else
else
write(' (',i+2,'!)^',j);
write(' (',FacIdx+2,'!)^',pow);
if num = 1 then
if num = 1 then
BREAK;
BREAK;
dec(i);
repeat
dec(FacIdx);
while (i>=0 ) AND Not(i in fm.fm_pow) do
until(FacIdx<0) OR (FacIdx in fm.fm_pow);
dec(i);
until i < 0;
until FacIdx < 0;
writeln;
writeln;

end;
end;


procedure Out_I_th(i: integer);
procedure Out_I_th(i: integer);
begin
begin
if i < 0 then
write(i:8,': ');
write(i:8,' too small');
if i <= High(FacMulPowGes) then
if i <= High(FacMulPowGes) then
begin
write(i:6,'-th : ');
Out_MulFac(i,FacMulPowGes[i-1])
Out_MulFac(i,FacMulPowGes[i-1])
end
else
else
writeln('Too big');
writeln('Too big');
end;

procedure Out_First_N(n: integer);
var
s,fmt : AnsiString;
i,tmp : integer;
Begin
if n<1 then
EXIT;
writeln('The first ',n,' Jordan-Polia numbers');
s := '';
If n > Length(FacMulPowGes) then
n := Length(FacMulPowGes);
dec(n);
tmp := length(CommatizeUint64(FacMulPowGes[n].fm_num))+1;
fmt := '%'+IntToStr(tmp)+'s';
tmp := 72 DIV tmp;
For i := 0 to n do
Begin
s += Format(fmt,[CommatizeUint64(FacMulPowGes[i].fm_num)]);
if (i+1) mod tmp = 0 then
Begin
writeln(s);
s := '';
end;
end;
if s <>'' then
writeln(s);
writeln;
end;
end;


Line 850: Line 879:
inc(idx);
inc(idx);
end;
end;
Fillchar(LastSearchFor,SizeOf(LastSearchFor),#0);
Fillchar(LastSearchFor,SizeOf(LastSearchFor),#0);
LastSearchFor.FM_NUM := 0;
LastSearchFor.FM_NUM := 0;
dblLimit := 1 shl 53;// 1 shl 6;
// dblLimit := 1 shl 53;
dblLimit := 1 shl 5;
end;
end;


procedure ResetSearch;
procedure ResetSearch;
Begin
Begin
setlength(FacMulPowGes,0);
setlength(FacMulPowGes,0);
end;
end;


Line 863: Line 893:
//generating the first entry with (2!)^n
//generating the first entry with (2!)^n
var
var
res_p : tpFac_mul;
Fac_mul :tFac_mul;
Fac_mul :tFac_mul;
facPow,Fac : tnum;
facPow,Fac : tnum;
i,MaxIdx : integer;
i,MaxPowOfFac : integer;
begin
begin
fac := Factorial[idx];
fac := Factorial[idx];
MaxIDx := trunc(ln(dblLimit)/ln(Fac))+1;
MaxPowOfFac := trunc(ln(dblLimit)/ln(Fac))+1;
setlength(res,MaxIdx);
setlength(res,MaxPowOfFac);

fillchar(Fac_Mul,SizeOf(Fac_Mul),#0);
with Fac_Mul do
with Fac_Mul do
begin
begin
Line 878: Line 907:
fm_high_idx := 0;
fm_high_idx := 0;
end;
end;

res_p := @res[0];
res_p^ := Fac_Mul;
res[0] := Fac_Mul;
facPow := 1;
facPow := 1;
For i := 1 to MaxIdx-1 do
i := 1;
begin
repeat
facPow *= Fac;
facPow *= Fac;
if facPow >dblLimit then
BREAK;
with Fac_Mul do
with Fac_Mul do
begin
begin
Line 889: Line 920:
fm_high_pow := i;
fm_high_pow := i;
end;
end;
inc(res_p);
res[i] := Fac_Mul;
res_p^ := Fac_Mul;
inc(i);
until i = MaxPowOfFac;
end;
setlength(res,i);
end;
end;


procedure DelDoublettes(var FMP:tFacMulPow);
procedure DelDoublettes(var FMP:tFacMulPow);
//throw out doublettes,
//throw out doublettes,
//the one with highest power in the highest n! survives
var
var
pI,pJ : tpFac_mul;
pNext,pCurrent : tpFac_mul;
i, j,idx : integer;
i, len,idx : integer;
begin
begin
j := 0;
len := 0;
pJ := @FMP[0];
pCurrent := @FMP[0];
pI := pJ;
pNext := pCurrent;
For i := 0 to High(FMP)-1 do
For i := 0 to High(FMP)-1 do
begin
begin
inc(pI);
inc(pNext);
// don't increment pCurrent if equal
if pJ^.fm_num = pI^.fm_num then
// pCurrent gets or stays the highest Value in n!^high_pow
if pCurrent^.fm_num = pNext^.fm_num then
Begin
Begin
idx := pJ^.fm_high_idx;
idx := pCurrent^.fm_high_idx;
if idx < pI^.fm_high_idx then
if idx < pNext^.fm_high_idx then
pJ^ := pI^
pCurrent^ := pNext^
else
else
if idx = pI^.fm_high_idx then
if idx = pNext^.fm_high_idx then
if pJ^.fm_high_pow < pI^.fm_high_pow then
if pCurrent^.fm_high_pow < pNext^.fm_high_pow then
pJ^ := pI^;
pCurrent^ := pNext^;
end
end
else
else
begin
begin
inc(j);
inc(len);
inc(pJ);
inc(pCurrent);
pJ^ := pI^;
pCurrent^ := pNext^;
end;
end;
end;
end;
setlength(FMP,j);
setlength(FMP,len);
end;

procedure QuickSort(var AI: tFacMulPow; ALo, AHi: Int32);
var
Tmp :tFac_mul;
Pivot : tnum;
Lo, Hi : Int32;
begin
Lo := ALo;
Hi := AHi;
Pivot := AI[(Lo + Hi) div 2].fm_num;
repeat
while AI[Lo].fm_num < Pivot do
Inc(Lo);
while AI[Hi].fm_num > Pivot do
Dec(Hi);
if Lo <= Hi then
begin
Tmp := AI[Lo];
AI[Lo] := AI[Hi];
AI[Hi] := Tmp;
Inc(Lo);
Dec(Hi);
end;
until Lo > Hi;
if Hi > ALo then
QuickSort(AI, ALo, Hi) ;
if Lo < AHi then
QuickSort(AI, Lo, AHi) ;
end;
end;


procedure InsertFacMulPow(var res:tFacMulPow;Facidx:integer);
function InsertFacMulPow(var res:tFacMulPow;Facidx:integer):boolean;
var
var
Fac,newFac,limit : tnum;
Fac,FacPow,NewNum,limit : tnum;
l_res,l_NewMaxPow,idx,i,j : Integer;
l_res,l_NewMaxPow,idx,i,j : Integer;
begin
begin
fac := Factorial[Facidx];
fac := Factorial[Facidx];

if length(res)= 0 then
Begin
GenerateFirst(Facidx,res);
EXIT;
end;

if fac>dblLimit then
if fac>dblLimit then
EXIT;
EXIT(false);
l_NewMaxPow := trunc(ln(dblLimit)/ln(Fac))+1;
l_res := length(res);


if length(res)> 0 then
//calc new length, reduces allocation of big memory chunks
j := 0;
begin
l_NewMaxPow := trunc(ln(dblLimit)/ln(Fac))+1;
idx := High(res);
l_res := length(res);
For i := 1 to l_NewMaxPow do
//calc new length, reduces allocation of big memory chunks
Begin
//first original length + length of the new to insert
limit := dblLimit DIV fac;
j := l_res+l_NewMaxPow;
if limit < 1 then
//find the maximal needed elements which stay below dbllimit
BREAK;
// for every Fac^i
repeat
dec(idx);
idx := High(res);
until res[idx].fm_num<=limit;
FacPow := Fac;
For i := 1 to l_NewMaxPow do
inc(j,idx);
Begin
fac *=Factorial[Facidx];
limit := dblLimit DIV FacPow;
end;
if limit < 1 then
j += l_res+l_NewMaxPow+2;
BREAK;
setlength(res,j);
//search for the right position
repeat
dec(idx);
until res[idx].fm_num<=limit;
inc(j,idx);
FacPow *= fac;
end;
j += 2;
setlength(res,j);


fac := Factorial[Facidx];
idx := l_res;
idx := l_res;
FacPow := fac;
For j := 0 to l_NewMaxPow-1 do
For j := 1 to l_NewMaxPow do
begin
For i := 0 to l_res-1 do
begin
begin
NewFac := res[i].fm_num*Fac;
For i := 0 to l_res do
if NewFac>dblLimit then
begin
Break;
NewNum := res[i].fm_num*FacPow;
if NewNum>dblLimit then
res[idx]:= res[i];
with res[idx] do
Break;
Begin
res[idx]:= res[i];
fm_num := NewFac;
with res[idx] do
include(fm_pow,Facidx);
Begin
fm_high_idx := Facidx;
fm_num := NewNum;
fm_high_pow := j;
include(fm_pow,Facidx);
end;
fm_high_idx := Facidx;
inc(idx);
fm_high_pow := j;
end;
inc(idx);
end;
FacPow *= fac;
end;
end;
setlength(res,idx);
fac *= Factorial[Facidx];
QuickSort(res,Low(res),High(res));
end;
setlength(res,idx);
DelDoublettes(res);
end
QuickSort(res,Low(res),High(res));
else
DelDoublettes(res);
GenerateFirst(Facidx,res);
Exit(true);
end;
end;


Line 993: Line 1,058:
BEGIn
BEGIn
InitFirst;
InitFirst;

repeat
repeat
ResetSearch;
ResetSearch;
i := 0;
i := 0;
repeat
repeat
if Factorial[i] < dblLimit then
if Not(InsertFacMulPow(FacMulPowGes,i)) then
InsertFacMulPow(FacMulPowGes,i)
BREAK;
else
break;
inc(i);
inc(i);
until i > High(Factorial);
until i > High(Factorial);
//check if MaxIdx is found
if (Length(FacMulPowGes) > MaxIdx) then
if (Length(FacMulPowGes) > MaxIdx) then
begin
begin
Line 1,008: Line 1,073:
Begin
Begin
LastSearchFor := FacMulPowGes[MaxIdx-1];
LastSearchFor := FacMulPowGes[MaxIdx-1];
//the next factorial is to big, so search is done
if LastSearchFor.fm_num < Factorial[i] then
if LastSearchFor.fm_num < Factorial[i] then
break;
break;
end
end
else
else
Break;
Break;
end;
end;
if dblLimit> HIGH(tNUm) DIV 256 then
if dblLimit> HIGH(tNUm) DIV 256 then
BREAK;
BREAK;
dblLimit *= 256;
dblLimit *= 256;
until false;
until false;

write('Found ',length(FacMulPowGes),' Jordan-Polia numbers ');
write('Found ',length(FacMulPowGes),' Jordan-Polia numbers ');
writeln('up to ',dblLimit);
writeln('up to ',CommatizeUint64(dblLimit));
writeln;
writeln;


Out_First_N(50);
writeln('The first 50 Jordan-Polia numbers');
For i := 1 to 50 do
Begin
write(FacMulPowGes[i-1].fm_num:5);
if i mod 10 = 0 then
writeln;
end;
writeln;


write('The last < 1E8 ');
write('The last < 1E8 ');
Line 1,040: Line 1,100:
writeln;
writeln;


Out_I_th(0);
Out_I_th(1);
Out_I_th(100);
Out_I_th(100);
Out_I_th(800);
Out_I_th(800);
Line 1,050: Line 1,110:
{{out|@home}}
{{out|@home}}
<pre>
<pre>
Found 3876 Jordan-Polia numbers up to 9007199254740992
Found 3876 Jordan-Polia numbers up to 9,007,199,254,740,992


The first 50 Jordan-Polia numbers
The first 50 Jordan-Polia numbers
1 2 4 6 8 12 16 24 32 36
1 2 4 6 8 12 16 24 32 36 48 64
48 64 72 96 120 128 144 192 216 240
72 96 120 128 144 192 216 240 256 288 384 432
256 288 384 432 480 512 576 720 768 864
480 512 576 720 768 864 960 1,024 1,152 1,296 1,440 1,536
1,728 1,920 2,048 2,304 2,592 2,880 3,072 3,456 3,840 4,096 4,320 4,608
960 1024 1152 1296 1440 1536 1728 1920 2048 2304
5,040 5,184
2592 2880 3072 3456 3840 4096 4320 4608 5040 5184

The last < 1E8 99532800 (6!)^2 4! (2!)^3


0: 1 1
The last < 1E8 99,532,800 = (6!)^2 4! (2!)^3
100: 92160 6! (2!)^7
800: 18345885696 (4!)^7 (2!)^2
1050: 139345920000 8! (5!)^3 2!
1800: 9784472371200 (6!)^2 (4!)^2 (2!)^15
2800: 439378587648000 14! 7!
3800: 7213895789838336 (4!)^8 (2!)^16


1-th : 1 = 1
real 0m0,002s user 0m0,002s sys 0m0,000s</pre>
100-th : 92,160 = 6! (2!)^7
800-th : 18,345,885,696 = (4!)^7 (2!)^2
1050-th : 139,345,920,000 = 8! (5!)^3 2!
1800-th : 9,784,472,371,200 = (6!)^2 (4!)^2 (2!)^15
2800-th : 439,378,587,648,000 = 14! 7!
3800-th : 7,213,895,789,838,336 = (4!)^8 (2!)^16
real 0m0,004s user 0m0,004s sys 0m0,000s</pre>


=={{header|Phix}}==
=={{header|Phix}}==
Retrieved from "https://rosettacode.org/wiki/Jordan-Pólya_numbers"