Skip to content

Instantly share code, notes, and snippets.

@Hermann-SW
Last active June 29, 2024 19:50
Show Gist options
  • Save Hermann-SW/df18355362f5cce70da6a2d5f20a22ca to your computer and use it in GitHub Desktop.
Save Hermann-SW/df18355362f5cce70da6a2d5f20a22ca to your computer and use it in GitHub Desktop.
Prime factors of Fermat numbers F0..F30
F__={[\\ needs parisizemax >= 600M; from http://www.prothsearch.com/fermat.html
[3], \\ F0
[5],
[17],
[257],
[65537],
[641, 6700417],
[274177, 67280421310721],
[59649589127497217, 5704689200685129054721],
[1238926361552897, "P62"],
[2424833, 7455602825647884208337395736200454918783366342657, "P99"],
[45592577, 6487031809, 4659775785220018543264560743076778192897, "P252"], \\ F10
[319489, 974849, 167988556341760475137, 3560841906445833920513, "P564"],
[114689, 26017793, 63766529, 190274191361, 1256132134125569, 568630647535356955169033410940867804839360742060818433, "C1133"],
[2710954639361, 2663848877152141313, 3603109844542291969, 319546020820551643220672513, "C2391"],
[116928085873074369829035993834596371340386703423373313, "C4880"],
[1214251009, 2327042503868417, 168768817029516972383024127016961, "C9808"],
[825753601, 188981757975021318420037633, "C19694"],
[31065037602817, 7751061099802522589358967058392886922693580423169, "C39395"],
[13631489, 81274690703860512587777, "C78884"],
[70525124609, 646730219521, 37590055514133754286524446080499713, "C157770"],
["C315653"], \\ F20
[4485296422913, "C631294"],
[64658705994591851009055774868504577, "C1262577"],
[167772161, "C2525215"],
["C5050446"],
[25991531462657, 204393464266227713, 2170072644496392193, "C10100842"],
[76861124116481, "C20201768"],
[151413703311361, 231292694251438081, "C40403531"],
[1766730974551267606529, "C80807103"],
[2405286912458753, "C161614233"],
[640126220763137, 1095981164658689, "C323228467"] \\ F30
]};
assert(b)=if(!(b),error());
isGint(a)=type(real(a))=="t_INT"&&type(imag(a))=="t_INT";
posGint(a)=if(real(a)<0,posGint(-a),if(imag(a)<0,conj(a),a));
bigGint(a)=if(real(a)>imag(a),a,imag(a)+real(a)*I);
sqrtGint(p)=assert(p%4==1&&isprime(p));my(s=factor(p,I)[3,1]);assert(norml2(s)==p);posGint(s);
divGint(a,b)=my(c=a/b);if(!isGint(c),c=a/conj(b));assert(isGint(c));posGint(c);
sqrtGintF(n)={
assert(n>0&&n<=30);
my(SF=2^2^(n-1)+I, pf=F__[1+n]);
foreach(pf[1..#pf-1],p, SF=divGint(SF,sqrtGint(p)));
SF;
}
sqrtGintFallKnown(n)={
my(S=[sqrtGintF(n)], pf=F__[1+n]);
foreach(pf[1..#pf-1],p,
f=sqrtGint(p); S=concat([[posGint(s*f),posGint(s*conj(f))]|s<-S])
);
vecsort([bigGint(s)|s<-S]);
}
sumOf3distinctSquares(n)=assert(n>=3);X=2^2^(n-1);[(2*X+1)\3,(2*X-2)\3,(X+2)\3];
\\ (26) in https://mathworld.wolfram.com/FermatNumber.html
isProbablePrimeCofactor(n)={
assert(n>11&&n<=30);
my(pf=F[1+n], F=vecprod(pf[1..#pf-1]), C=pf[#pf], Fn=2^2^n+1;);
my(A=Mod(3,Fn)^(Fn-1), B=Mod(3,Fn)^(F-1));
lift(A-B)%C != 0;
}
{
assert(1+30==#F__);
i=0;
foreach(F__,v,s=v[#v];f=2^2^i+1;
if(type(s)=="t_STR",
p=vecprod(v[1..#v-1]);
assert(f%p==0);
if(Vec(s)[1]=="P",
assert(#digits(f\p)==eval(strjoin(Vec(s)[2..#s])))
);
F__[1+i][#v]=f\p;
);
assert(vecprod(v)==f);
i+=1
);
n=0;
foreach(F__,v, \\ Lucas 1878
if(n>1, \\ print(n);
foreach(v,f,
assert(valuation(f-1,2)>=n+2)
)
);
n+=1;
);
n=0;
foreach(F__,v, \\ all F5..F30 prime factors have sum of squares
foreach(v[1..(n>4)*#v-(n>11)],p, assert(p%4==1&&isprime(p)));
n+=1;
)
}
F__;
@Hermann-SW
Copy link
Author

Hermann-SW commented May 6, 2024

On slow Intel Celeron J4125 CPU Fermat.gp read and initialization time is only 3.89× slower than that of AMD 7950X CPU:

hermann@j4125:~$ gp -q
? F=readvec("Fermat.gp");F=F[#F];
? ##
  ***   last result: cpu time 20,122 ms, real time 12,912 ms.
? version()
[2, 15, 4]
?

The other commands from previous comment complete instantaneously on J4125.

Not bad, J4125 time reported for F30 is only 4.73× slower than that for AMD 7950X CPU.
Passmark benchmark single thread ratio is 4,288 / 1165 = 3.68:
https://www.cpubenchmark.net/compare/5031vs3667/AMD-Ryzen-9-7950X-vs-Intel-Celeron-J4125

hermann@j4125:~$ gp -q < Fermat_comp.gp 
F12 has 6 prime factors; determining composite sum of squares: 1ms
F13 has 4 prime factors; determining composite sum of squares: 9ms
F14 has 1 prime factors; determining composite sum of squares: 27ms
F15 has 3 prime factors; determining composite sum of squares: 2ms
F16 has 2 prime factors; determining composite sum of squares: 1ms
F17 has 2 prime factors; determining composite sum of squares: 22ms
F18 has 2 prime factors; determining composite sum of squares: 2ms
F19 has 3 prime factors; determining composite sum of squares: 5ms
F20 has 0 prime factors; determining composite sum of squares: 3ms
F21 has 1 prime factors; determining composite sum of squares: 16ms
F22 has 1 prime factors; determining composite sum of squares: 41ms
F23 has 1 prime factors; determining composite sum of squares: 76ms
F24 has 0 prime factors; determining composite sum of squares: 87ms
F25 has 3 prime factors; determining composite sum of squares: 462ms
F26 has 1 prime factors; determining composite sum of squares: 879ms
F27 has 2 prime factors; determining composite sum of squares: 1962ms
F28 has 1 prime factors; determining composite sum of squares: 3915ms
F29 has 1 prime factors; determining composite sum of squares: 8107ms
F30 has 2 prime factors; determining composite sum of squares: 18194ms

real(SF12) = 200632848085394229198405077309776409669556160755822894920478194045891524675173877582799789843512719390209285348887171584058267613825062519170949236869832740299611688879431491248560122275125138227835639875304442149679485916420376715785002453587853905329008047468218821526665318251417289791164787502264540469658007753188396466487968753988674615092615847790001421479841641921279595503860736218792224235350272376658369292603790019796500735806899786991660195728966759044116399240680328117271881207382080232786405040556863376322477213246700048245459183343930058344600346916
imag(SF12) = 11512882899820054257144225772505994511430981968359355559240636997087397239461885404688940982112272498773691260355731224763278685518244745544198267923163368736091123701779226072209279679342867029500044275233215203437226071842172804234583591297137729569486761340213325710137879698831126615998659706343950808674850862574868322314902443424081205544133789500128645355501388833990928089030944977862262874243179626287736961093227838096073086612878632276868708056678373714902078426666851025890207418013027573248367464970951431311736356210867866665430397629513384884406535591
hermann@j4125:~$ 

@Hermann-SW
Copy link
Author

Hermann-SW commented May 6, 2024

Some more CPUs for runtime comparison:

GP read&init [s] F30 comp [s]
AMD Ryzen 9 7950X 2.16.1 3.316 3.85
AMD Ryzen 7 7600X 2.15.4 4.585 3.914
Intel i7 11700K 2.16.1 5.111 4.581
Intel Xeon 6126 2.15.4 10.086 9.0
Intel Xeon e5-2680 2.15.4 11.291 9.695
Qualcomm Kryo 480 2.15.4 below 15.823 --+
Intel Celeron J4125 2.15.4 12.912 18.194 !
Intel Celeron J4105 2.15.5 15.456 19.787 !
Raspberry Pi5 (3GHz) 2.15.5 18.299 40.677 !
Raspberry Pi5 Cortex-A76 2.15.5 22.537 50.867 !
Qualcomm Kryo 480 2.15.4 26.205 above --+

On my Nokia XR20 under termux running "Fermat_comp.gp" to determine "F30 comp" time kills termux because of memory issues. I created stripped down script F30.gp that allows to determine that time (with gp -q -s 700000000 < F30.gp) on systems with memory problem (it is unclear why termux terminates when running while Pi5 does not — both have 4GB RAM). I verified that F30.gp reports same time as Fermat_comp.gp on J4105 CPU.

Under termux on Kryo 480 I had to start gp -q -s 700000000 to determine read&init time with:

F=readvec("Fermat.gp");F=F[#F];
##

Kryo 480 has two rows in above table because "read&init" is extremely slow, while "F30 comp" is quite fast.

@Hermann-SW
Copy link
Author

Hermann-SW commented May 13, 2024

For n≥3 one can use sumOf3distinctSquares():
https://math.stackexchange.com/a/941933

hermann@7950x:~$ gp -q Fermat.gp
? sumOf3distinctSquares(3)
[11, 10, 6]
? norml2(sumOf3distinctSquares(3))
257
? s=sumOf3distinctSquares(30);
? ##
  ***   last result: cpu time 151 ms, real time 199 ms.
? #s
3
? norml2(s)==2^2^30+1
1
? ##
  ***   last result: cpu time 5,109 ms, real time 5,462 ms.
? s[1]>s[2]&&s[2]>s[3]&&s[3]>0
1
? 

@Hermann-SW
Copy link
Author

Hermann-SW commented Jun 29, 2024

New isProbablePrimeCofactor(n) based on
(26) in https://mathworld.wolfram.com/FermatNumber.html
gets quite slow for bigger n, even on 32GB RAM AMD 7950X fast CPU:

? gettime();for(n=12,30,print("F",n,": ",isProbablePrimeCofactor(n)," ",gettime()))
F12: 1 15
F13: 1 73
F14: 1 414
F15: 1 2057
F16: 1 10677
F17: 1 58440
F18: 1 298709

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment