Talk:Factor-perfect numbers: Difference between revisions
m
→Constructing p1**n *p2xp3xp3...: mpz version
m (→Constructing p1**n *p2xp3xp3...: mpz version) |
|||
| (6 intermediate revisions by 3 users not shown) | |||
Line 76:
=== wp link deleted ===
I have deleted the "see also" link to [[wp:Enumerative_combinatorics|Wikipedia: Enumerative Combinatorics]] because that page contains nothing pertinent to this task. --[[User:Petelomax|Petelomax]] ([[User talk:Petelomax|talk]])
: I think rc has a problem with tasks where a schoolboy with a pencil could do better than the solutions presented. There are many of them and my time is limited, a loop over the set of integers testing each to see if it meets the tasks requirements is I think a dumb solution, so I added a link to a Wikipedia entry which I think describes an intelligent solution. To work it requires a generating function which this task supplies as F(n). As any schoolboy worth the effort of educating (better the rest are just sent out to play football) will tell you functions of this form depend only on the characterization of n's prime decomposition. So let n=P1<sup>A1</sup>*P2<sup>A2</sup> ... Pk<sup>Ak</sup>. 48 is 3*2<sup>4</sup> so let me consider those n such that the sum of A1..A5 is less than 6.
A1 A2 A3 A4 A5 n<sub>min</sub> F(n)
1 0 0 0 0 2 1
2 0 0 0 0 4 2
1 1 0 0 0 6 3
3 0 0 0 0 8 4
2 1 0 0 0 12 8
1 1 1 0 0 30 13
4 0 0 0 0 16 8
3 1 0 0 0 24 20
2 2 0 0 0 36 26
2 1 1 0 0 60 44
1 1 1 1 0 210 75
5 0 0 0 0 32 16
4 1 0 0 0 48 48
3 2 0 0 0 72 76
3 1 1 0 0 120 132
2 2 1 0 0 180 176
2 1 1 1 0 420 308
1 1 1 1 1 2310 541
n<sub>min</sub> is 2<sup>A1</sup>*3<sup>A2</sup>*5<sup>A3</sup>*7<sup>A4</sup>*11<sup>A5</sup>. 7*97<sup>4</sup> has the same characteristic as 3*2<sup>4</sup> so without further calculation I suggest F(619704967) is 48. Interestingly, the above table suggests that F(n) is odd only if all A are 1. I think I can prove this and that no such factorization meets the task requirements, which implies that n cannot be odd. See [[Talk:Untouchable numbers#Nice_rec]] for similar logic and [[9 billion names of God the integer]] for A<sup>1</sup>..A<sup>k</sup>.--[[User:Nigel Galloway|Nigel Galloway]] ([[User talk:Nigel Galloway|talk]]) 13:06, 23 March 2023 (UTC)
=== How about counting directed knots connecting in Hasse-Diagram ===
<pre>
48 = 2^4 * 3
48-24--12---6---3
| | | | | ^
| | | | | ^ factor 3 ( prime b )
| | | | | ^
16--8---4---2---1 < factor 2 (prime a)
^4 ^3 ^2 ^1 (pot of prime a )
Starting in knot 1 and than only to knots that increase the distance to 1
like and don't change side (prime-factor) more than once.
1->2->6->12->24->48
or 1->48, 1->12->48 , 1->2->4->8->16->48
or 1->6->24->48
but not 1->3->4->24->48
</pre>
But one advantage is, that one solution of a primfactorisation a^n*b^m .. gives always the same solution.
: Sounds to me like you had a ''really'' good Paddy's night, and that's the result of the hangover... --[[User:Petelomax|Petelomax]] ([[User talk:Petelomax|talk]]) 12:40, 20 March 2023 (UTC)
=== Constructing p1**n *p2xp3xp3... ===
Most of the Factor-perfect numbers are of a form 2*n *p2xp3xp3..x p_n<br>
Using bigint will probably find new ones.
<syntaxhighlight lang=pascal>
program PerfectFactor;
//https://mathworld.wolfram.com/OrderedFactorization.html
{$IFDef WINDOWS}
{$APPTYPE Console}
{$ENDIF}
{$IFDEF FPC }
//{$R+,I+,O+}
{$MODE Delphi}
{$OPTIMIZATION On,ALL}
{$ENDIF}
uses
SysUtils,gmp;
const
smlPrimes: array[0..5] of byte = (2,3,5,7,11,13);
type
tPrimes = array[0..15] of Uint64;
tarrMPInt = array of mpz_t;
ptarrMPInt = ^tarrMPInt;
tP_x = 0..6;
var
wheeldiff: array [0..(2-1)*(3-1)*(5-1)*(7-1)*(11-1)*(13-1)-1] of byte;
Primes: tPrimes;
powers: array[0..15] of Uint32;
ps : array[0..16383] of AnsiChar;
P : Array[tP_x] of tarrMPInt;
primLmtCheck : mpz_t;
CntPrimes : integer;
procedure InitWheelPrimes;
var
sieve : array of byte;
i,j,p,ws: integer;
begin
ws := 1;
For i := low(smlPrimes) to High(smlPrimes) do
ws *= smlPrimes[i];
setlength(sieve,ws+1);//implicit fillchar(sieve[0],ws+1,#0);
sieve[0] := 1;
For i := low(smlPrimes) to High(smlPrimes) do
Begin
p := smlPrimes[i];
j := p;
while j <= ws do
begin
sieve[j] := 1;
inc(j,p);
end;
end;
j := 0;
p := 1;
For i := p+1 to High(sieve) do
begin
if sieve[i] = 0 then
begin
wheeldiff[j] := i-p;
p := i;
inc(j);
end;
end;
wheeldiff[j] := 2;
end;
procedure Init_P1_P6(len:Int32);
//p1^x .. p1^x*p2*p3*p4*p5*p6
var
ofs,p1,p1_f2,tmp_mpz: mpz_t;
f1,f2,f3,f4,f5 : Uint64;
i : integer;
begin
For i := low(tP_x) to High(tP_x) do
setlength(P[i],len+1);
mpz_init_set_ui(P[0,1],1);
For i := 2 to len do
mpz_init_set_ui(P[0,i],i-1);
mpz_init_set_ui(P[1,0],1);
mpz_init_set_ui(ofs,1);
For i := 1 to len do
Begin
mpz_init_set(P[1,i],P[1,i-1]);//2*(P2[i-1])+p;
mpz_add(P[1,i],P[1,i],P[1,i]);
mpz_add(P[1,i],P[1,i],ofs);
mpz_add(ofs,ofs,ofs);
end;
mpz_init_set_ui(p1,1);
mpz_init_set_ui(P[2,0],3);
mpz_set_ui(ofs,1);
f1 := 4;
For i := 1 to len do
Begin
mpz_set(p1,P[1,i]);// p1 := P1[i];
mpz_init_set(P[2,i],p1);
mpz_mul_ui(P[2,i],P[2,i],f1);//f1*p1
mpz_add(P[2,i],P[2,i],ofs);//f1*p1+p;
mpz_add(ofs,ofs,p1);//p += p1;
f1+= 1;
end;
mpz_init(p1_f2);
mpz_init_set_ui(P[3,0],75);
mpz_set_ui(ofs,1);
f1 := 5;
f2 := 3;
For i := 1 to len do
Begin
mpz_set(p1,P[1,i]);//p1 := P2[i];
mpz_set(p1_f2,p1);
mpz_mul_ui(p1_f2,p1_f2,f2);//f2*p1
mpz_init_set(P[3,i],P[2,i]);
mpz_mul_ui(P[3,i],P[3,i],f1);//f1*P3[i]
mpz_add(P[3,i],P[3,i],p1_f2);//f1*P3[i]+f2*p1;
mpz_add(P[3,i],P[3,i],ofs);//f1*P3[i]+f2*p1+ofs;
f1 +=1;
f2 +=2;
mpz_add(ofs,ofs,p1);//ofs += p1;
end;
// 3 13 75 541 7* 75 + 1*13+1*3
// 8 44 308 2612 8*308 + 3*44+2*8
// 20 132 1076 10404 9*1076 +5*132+3*20
mpz_init(tmp_mpz);
mpz_set_ui(ofs,1);
f1 := 7;
f2 := 1;
f3 := 1;
For i := 1 to len do
Begin
mpz_set(tmp_mpz,p[1,i]);
mpz_mul_ui(tmp_mpz,tmp_mpz,f3);//f3*P2[i];
mpz_init_set(P[4,i],P[3,i]);
mpz_mul_ui(P[4,i],P[4,i],f1);//f1*P4[i]
mpz_add(P[4,i],P[4,i],tmp_mpz);//f1*P4[i]+f3*P2[i]
mpz_set(tmp_mpz,p[2,i]);
mpz_mul_ui(tmp_mpz,tmp_mpz,f2);//f2*P3[i]
mpz_add(P[4,i],P[4,i],tmp_mpz);//f1*P4[i]+f2*P3[i]+f3*P2[i];
f1 +=1;
f2 +=2;
f3 += 1;
end;
// 3 13 75 541 4683 8*541 +4* 75 +4*13+1*3
// 8 44 308 2612 25988 9*2612 +7*(308+44)+2*8
// 20 132 1076 10404 116180 10*10404+10*(1076+132)+3*20
// P6[i] :=f1*P5[i]+f2*(P4[i]+P3[i])+p*P2[i];
f1 := 8;
f2 := 4;
f3 := 1;
For i := 1 to len do
Begin
mpz_init_set(P[5,i],P[4,i]);
mpz_mul_ui(P[5,i],P[5,i],f1);
mpz_set(tmp_mpz,p[2,i]);
mpz_add(tmp_mpz,P[3,i],tmp_mpz);
mpz_mul_ui(tmp_mpz,tmp_mpz,f2);
mpz_add(P[5,i],P[5,i],tmp_mpz);
mpz_set(tmp_mpz,p[1,i]);
mpz_mul_ui(tmp_mpz,tmp_mpz,f3);
mpz_add(P[5,i],P[5,i],tmp_mpz);
f1 +=1;
f2 +=3;
f3 += 1;
end;
{ P5 P4 P3 P2
13;75;541;4683;47293 9* 4683 +9*541 +3*75 +4*13
44;308;2612;25988;296564 10*25988 +13*2612 +8*308 +6*44
132;1076;10404;116180;1469892 11*116180+17*10404+13*1076 +8*132
}
f5 := 9;
f4 := 9;
f3 := 3;
f2 := 4;
For i := 1 to len do
Begin
mpz_init_set(P[6,i],P[5,i]);
mpz_mul_ui(P[6,i],P[6,i],f5);
mpz_set(tmp_mpz,p[4,i]);
mpz_mul_ui(tmp_mpz,tmp_mpz,f4);
mpz_add(P[6,i],P[6,i],tmp_mpz);
mpz_set(tmp_mpz,p[3,i]);
mpz_mul_ui(tmp_mpz,tmp_mpz,f3);
mpz_add(P[6,i],P[6,i],tmp_mpz);
mpz_set(tmp_mpz,p[2,i]);
mpz_mul_ui(tmp_mpz,tmp_mpz,f2);
mpz_add(P[6,i],P[6,i],tmp_mpz);
f5 +=1;
f4 +=4;
f3 += 5;
f2 += 2;
end;
mpz_clear(ofs);
mpz_clear(p1);
mpz_clear(p1_f2);
mpz_clear(tmp_mpz);
end;
procedure OutMPZ(var tmp:mpz_t);
begin
mpz_get_str(@ps[0],10,tmp);
writeln(pChar(@ps[0]));
end;
procedure OutPow(const pr:tPrimes);
var
i,cnt :integer;
begin
cnt := CntPrimes-1;
IF cnt<0 then
EXIT;
For i := 0 to cnt-1 do
if powers[i] = 1 then
write(pr[i],'*')
else
write(pr[i],'^',powers[i],'*');
if powers[cnt] = 1 then
write(pr[cnt])
else
write(pr[cnt],'^',powers[cnt]);
writeln;
end;
function CalcNewLimit64(n:Uint64;cnt:Int32):Uint64;
//new limit p1*p2_p_cnt >= n ~ p1^cnt >= n lmt = n^(1/cnt)
begin
cnt := cnt-CntPrimes;
if cnt <= 1 then
cnt := 2;
result := trunc(exp(ln(n)/Cnt));
end;
procedure GetPrDcmp(n: UInt64);
var
p,lmt : Uint64;
FlipFlop : Int32;
begin
fillchar(powers,SizeOf(powers),#0);
lmt := trunc(sqrt(n));
CntPrimes := 0;
if Not(Odd(n)) then
Begin
Primes[CntPrimes] := 2;
repeat n := n SHR 1;inc(powers[CntPrimes]);until ODD(N);
inc(CntPrimes);
end;
if n Mod 3 = 0 then
Begin
Primes[CntPrimes] := 3;
repeat n := n div 3; inc(powers[CntPrimes]); until n mod 3 <> 0;
inc(CntPrimes);
end;
if n Mod 5 = 0 then
Begin
Primes[CntPrimes] := 5;
repeat n := n div 5; inc(powers[CntPrimes]); until n mod 5 <> 0;
inc(CntPrimes);
end;
if n Mod 7 = 0 then
Begin
Primes[CntPrimes] := 7;
repeat n := n div 7; inc(powers[CntPrimes]); until n mod 7 <> 0;
inc(CntPrimes);
end;
if n Mod 11 = 0 then
Begin
Primes[CntPrimes] := 11;
repeat n := n div 11; inc(powers[CntPrimes]); until n mod 11 <> 0;
inc(CntPrimes);
end;
if n Mod 13 = 0 then
Begin
Primes[CntPrimes] := 13;
repeat n := n div 13; inc(powers[CntPrimes]); until n mod 13 <> 0;
inc(CntPrimes);
end;
p := 1;
flipflop := High(WheelDiff)-1;
repeat
p := p + wheeldiff[flipflop];
flipflop -= 1;
if flipflop < 0 then
flipflop := High(wheelDiff)-1;
if p > lmt then
BREAK;
if n Mod p = 0 then
Begin
Primes[CntPrimes] := p;
repeat
n := n div p;
inc(powers[CntPrimes]);
until n mod p <> 0;
inc(CntPrimes);
end;
until false;
if (n > 1) then
begin
Primes[CntPrimes] := n;
powers[CntPrimes] := 1;
inc(CntPrimes);
end;
end;
function GetPrimeDecomp64(n: UInt64;p:Uint64;FlipFlop,cnt:integer):boolean;
var
lmt: NativeUInt;
begin
if n Mod 3 = 0 then
Begin
n := n div 3;If n mod 3 = 0 then EXiT(false);
Primes[CntPrimes] := 3;powers[CntPrimes] := 1;inc(CntPrimes);
if CntPrimes>cnt then EXIT(false);
end;
if n Mod 5 = 0 then
Begin
n := n div 5;If n mod 5 = 0 then EXiT(false);
Primes[CntPrimes] := 5;powers[CntPrimes] := 1;inc(CntPrimes);
if CntPrimes>cnt then EXIT(false);
end;
if n Mod 7 = 0 then
Begin
n := n div 7; If n mod 7 = 0 then EXiT(false);
Primes[CntPrimes] := 7;powers[CntPrimes] := 1;inc(CntPrimes);
if CntPrimes>cnt then EXIT(false);
end;
if n Mod 11= 0 then
Begin
n := n div 11; If n mod 11 = 0 then EXiT(false);
Primes[CntPrimes] := 11;powers[CntPrimes] := 1;inc(CntPrimes);
if CntPrimes>cnt then EXIT(false);
end;
if n Mod 13 = 0 then
Begin
n := n div 13;If n mod 13 = 0 then EXiT(false);
Primes[CntPrimes] := 13;powers[CntPrimes] := 1;inc(CntPrimes);
if CntPrimes>cnt then EXIT(false);
end;
result := true;
lmt := CalcNewLimit64(n,cnt);
while result do
Begin
p := p + wheeldiff[flipflop];
flipflop -= 1;
if flipflop < 0 then
flipflop := High(wheeldiff);
if p > lmt then
BREAK;
if n Mod p = 0 then
Begin
n := n div p;
If n mod p = 0 then
EXiT(false);
Primes[CntPrimes] := p;
powers[CntPrimes] := 1;
inc(CntPrimes);
if CntPrimes>cnt then
EXIT(false);
lmt := CalcNewLimit64(n,cnt);
end;
end;
if (n > 1) then
begin
Primes[CntPrimes] := n;powers[CntPrimes] := 1;inc(CntPrimes);
end;
result := CntPrimes = cnt;
end;
function CalcprimLmt(var lmtMpz,tmpMpz:mpz_t;cnt:integer):Uint64;
//new limit p1*p2*..*p_cnt >= n ~ p1^cnt >= n lmt = n^(1/cnt)
begin
mpz_set(lmtMpz,tmpMpz);
cnt := cnt-CntPrimes;
if cnt <= 1 then
cnt := 2;
mpz_root(lmtMpz,lmtMpz,cnt);
result := 0;
if mpz_fits_ulong_p(lmtMpz)<> 0 then
result := mpz_get_ui(lmtMpz);
if result > High(Uint32) then
result := High(Uint32);
end;
function GetPrimeDecomp(var tmpMpz : mpz_t;cnt:Integer):Boolean;
var
lmtMpz : mpz_t;
p,n: UInt64;
flipflop : integer;
begin
For FlipFlop := 0 to High(smlPrimes) do
Begin
p := smlPrimes[FlipFlop];
if mpz_divisible_ui_p(tmpMpz,p)<>0 then
begin
mpz_divexact_ui(tmpMpz,tmpMpz,p);
if mpz_divisible_ui_p(tmpMpz,p)<>0 then EXIT(false);
Primes[CntPrimes] := p;
powers[CntPrimes] := 1;
inc(CntPrimes);
end;
end;
mpz_init(lmtMpz);
n := CalcprimLmt(lmtMpz,tmpMpz,cnt);
p := 1;
flipflop := High(wheeldiff)-1;
while (p <= n) do
begin
p := p + wheeldiff[flipflop];
flipflop := flipflop - 1;
if flipflop < 0 then
flipflop := High(wheeldiff);
if mpz_divisible_ui_p(tmpMpz,p)<>0 then
begin
mpz_divexact_ui(tmpMpz,tmpMpz,p);
Primes[CntPrimes] := p;
powers[CntPrimes] := 1;
inc(CntPrimes);
if (CntPrimes>cnt) OR (mpz_divisible_ui_p(tmpMpz,p)<>0) then
Begin
mpz_clear(lmtMpz);
EXIT(FALSE);
end;
if CntPrimes = cnt-1 then
BREAK;
if mpz_fits_ulong_p(tmpMpz)=1 then
Begin
mpz_clear(lmtMpz);
n := mpz_get_ui(tmpMpz);
Exit(GetPrimeDecomp64(n,p,FlipFlop,cnt));
end;
n := CalcprimLmt(lmtMpz,tmpMpz,cnt);
end;
end;
if (n > 1) AND (mpz_probab_prime_p(tmpMpz,15) > 0) then
begin
Primes[CntPrimes] := mpz_get_ui(tmpMpz);
powers[CntPrimes] := 1;
inc(CntPrimes);
end;
result := CntPrimes = cnt;
end;
procedure CheckForFactorPerfect(p : ptarrMPInt;maxpow,maxcnt:integer);
var
numMpz,tmpMpz : mpz_t;
n : Uint64;
max2pow : integer;
begin
primes[0] := 2;
n := 0;
mpz_init(numMpz);
mpz_init(tmpMpz);
writeln('check count of prime = ',maxcnt);
For max2pow := 0 to maxpow-1 do
Begin
if ODD(max2pow)
then continue;
mpz_set(numMpz,p^[max2pow]);
mpz_set(tmpMpz,numMpz);
if mpz_divisible_2exp_p(tmpMpz,max2Pow) <> 0 then
if mpz_divisible_2exp_p(tmpMpz,max2Pow+1) = 0 then
Begin
powers[0] := max2Pow;
CntPrimes := 1;
mpz_tdiv_q_2exp(tmpMpz,tmpMpz, max2Pow);
if mpz_fits_ulong_p(tmpMpz)=1 then
begin
n := mpz_get_ui(tmpMpz);
if GetPrimeDecomp64(n,1,High(wheeldiff)-1,maxcnt) then
begin
OutPow(primes);//OutMPZ(numMpz);
end;
end
else
if GetPrimeDecomp(tmpMpz,maxcnt) then
begin
OutPow(primes);//OutMPZ(numMpz);
end;
end;
end;
mpz_clear(tmpMpz);
writeln;
end;
var
len,
i : integer;
begin
InitWheelPrimes;
len := 160;// -> 33000-> 10000 digits
primes[0] := 2;
mpz_init(primLmtCheck);
Init_P1_P6(len);
For i := low(tP_x) to High(tP_x) do
CheckForFactorPerfect(@P[i],len,i+1);
end.</syntaxhighlight>
{{out}}
<pre>check count of prime = 1
check count of prime = 2
2^4*3
2^8*5
2^12*7
2^20*11
2^24*13
2^32*17
2^36*19
2^44*23
2^56*29
2^60*31
2^72*37
2^80*41
2^84*43
2^92*47
2^104*53
2^116*59
2^120*61
2^132*67
2^140*71
2^144*73
2^156*79
check count of prime = 3
2^6*3*13
2^12*3*37
2^14*11*13
2^18*3*73
2^30*3*181
2^32*13*47
2^40*13*71
2^42*3*337
2^46*11*109
2^48*3*433
2^50*23*61
2^54*3*541
2^56*37*47
2^60*3*661
2^68*11*229
2^72*3*937
2^78*3*1093
2^80*11*313
2^82*23*157
2^86*37*107
2^92*13*347
2^104*59*97
2^108*3*2053
2^110*13*491
2^112*11*601
2^114*3*2281
2^118*13*563
2^120*3*2521
2^122*73*107
2^124*11*733
2^128*23*373
2^130*37*239
2^132*3*3037
2^134*11*853
2^136*13*743
2^138*3*3313
2^142*23*457
2^146*11*1009:
check count of prime = 4
2^34*3*47*193
2^40*3*103*137
2^44*7*29*271
2^52*3*13*2243
2^58*3*7*5657
2^64*3*167*313
2^82*3*109*971
2^86*5*7*10399
2^90*11*37*1019
2^98*13*59*691
2^102*29*73*281
2^104*13*29*1669
2^112*3*11*23629
2^136*3*5*91411
2^146*5*71*4751
check count of prime = 5
2^18*3*7*59*107
2^54*3*7*17*16829
2^72*3*29*53*3793
2^88*7*17*23*13619
2^116*7*139*211*521
2^124*7*19*173*6007
check count of prime = 6
2^28*3*5*23*71*883
2^48*5*11*17*67*3593
check count of prime = 7
2^156*3*23*53*71*233*154409</pre>
all 2^(2*(p2-1) * p2 where p2 is prime >2 are all Factor-perfect numbers
</pre>
--[[User:Horst|Horst]] ([[User talk:Horst|talk]]) 06:52, 12 April 2023 (UTC)
| |||