-
-
Save Hermann-SW/df18355362f5cce70da6a2d5f20a22ca to your computer and use it in GitHub Desktop.
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__; |
https://mathworld.wolfram.com/FermatNumber.html
In 1878, Lucas increased the exponent of 2 by one, showing that factors of Fermat numbers must be of the form 2^(n+2)*L+1 for n>1.
This is verified above in lines 50-58 for prime as well as composite factors.
Fermat factor excess over exponent of 2 guaranteed by Lucas:
? {
n=0;
foreach(F__,v, \\ Lucas 1878
if(n>1, \\ print(n);
foreach(v,f,
print1(valuation(f-1,2)-(n+2)," ");
)
);
n+=1;
print();
)
}
0
3
10
0 0
0 0
0 0
1 1
5 0 0
0 2 0 1
0 0 1 0 0
0 2 2 0 0 1 0
1 2 4 4 1
0 0
4 0 0 2
1 2 1
0 0 5
0 3 0
0 0 1 2
1048554
0 0
0 0
0 0
16777190
2 0 0 1
1 1
1 0 0
6 6
0 0
0 1 0
?
From "Prime factors k*2^n+1 of larger Fermat numbers Fm" than this gist F0..F30:
http://www.prothsearch.com/fermat.html#Larger
m k n Year Discoverer
...
11075 171369935 11077 29 Jun 2020 G. B. Gostin
...
50078 7619 50081 06 Mar 2002 G. Axelsson, Jobling, Woltman & Gallot
...
100518 79425 100520 22 Oct 2014 P. Samidoost, Fougeron & Penné
...
303088 3 303093 13 Jan 1998 J. Young
...
Proving k*2^n+1 as factor of Fm can be easily done with PARI/GP modular arithmetic:
? lift(Mod(2,171369935*2^11077+1)^(2^11075)+1)
%5 = 0
? ##
*** last result: cpu time 162 ms, real time 163 ms.
? lift(Mod(2,7619*2^50081+1)^(2^50078)+1)
%6 = 0
? ##
*** last result: cpu time 5,928 ms, real time 5,928 ms.
? lift(Mod(2,79425*2^100520+1)^(2^100518)+1)
%7 = 0
? ##
*** last result: cpu time 30,638 ms, real time 30,642 ms.
? lift(Mod(2,3*2^303093+1)^(2^303088)+1)
%8 = 0
? ##
*** last result: cpu time 6min, 55,824 ms, real time 6min, 55,847 ms.
?
New sqrtGauss(p) and sqrtGaussF(n) allow to determine sum of squares for the composite factors of F12..F30 easily ...
hermann@7950x:~$ cat Fermat_comp.gp
F=readvec("Fermat.gp");F=F[#F];
{
sqrt(Mod(-1,641));gettime(); \\ warmup
n=0;
foreach(F,v,
if(n>11,
print1("F",n," has ",#v-1," prime factors; ");
SF=sqrtGintF(n);
assert(norml2(SF)==(2^2^n+1)/vecprod(v[1..#v-1]));
print("determining composite sum of squares: ",gettime(),"ms");
);
n+=1;
);
SF12=sqrtGintF(12);
print("\nreal(SF12) = ",real(SF12));
print( "imag(SF12) = ",imag(SF12));
}
hermann@7950x:~$
... and fast:
hermann@7950x:~$ gp -q < Fermat_comp.gp
F12 has 6 prime factors; determining composite sum of squares: 0ms
F13 has 4 prime factors; determining composite sum of squares: 4ms
F14 has 1 prime factors; determining composite sum of squares: 7ms
F15 has 3 prime factors; determining composite sum of squares: 0ms
F16 has 2 prime factors; determining composite sum of squares: 0ms
F17 has 2 prime factors; determining composite sum of squares: 2ms
F18 has 2 prime factors; determining composite sum of squares: 0ms
F19 has 3 prime factors; determining composite sum of squares: 1ms
F20 has 0 prime factors; determining composite sum of squares: 1ms
F21 has 1 prime factors; determining composite sum of squares: 3ms
F22 has 1 prime factors; determining composite sum of squares: 11ms
F23 has 1 prime factors; determining composite sum of squares: 16ms
F24 has 0 prime factors; determining composite sum of squares: 19ms
F25 has 3 prime factors; determining composite sum of squares: 92ms
F26 has 1 prime factors; determining composite sum of squares: 180ms
F27 has 2 prime factors; determining composite sum of squares: 413ms
F28 has 1 prime factors; determining composite sum of squares: 837ms
F29 has 1 prime factors; determining composite sum of squares: 1710ms
F30 has 2 prime factors; determining composite sum of squares: 3850ms
real(SF12) = 200632848085394229198405077309776409669556160755822894920478194045891524675173877582799789843512719390209285348887171584058267613825062519170949236869832740299611688879431491248560122275125138227835639875304442149679485916420376715785002453587853905329008047468218821526665318251417289791164787502264540469658007753188396466487968753988674615092615847790001421479841641921279595503860736218792224235350272376658369292603790019796500735806899786991660195728966759044116399240680328117271881207382080232786405040556863376322477213246700048245459183343930058344600346916
imag(SF12) = 11512882899820054257144225772505994511430981968359355559240636997087397239461885404688940982112272498773691260355731224763278685518244745544198267923163368736091123701779226072209279679342867029500044275233215203437226071842172804234583591297137729569486761340213325710137879698831126615998659706343950808674850862574868322314902443424081205544133789500128645355501388833990928089030944977862262874243179626287736961093227838096073086612878632276868708056678373714902078426666851025890207418013027573248367464970951431311736356210867866665430397629513384884406535591
hermann@7950x:~$
hermann@7950x:~$ gp -q
? F=readvec("Fermat.gp");F=F[#F];
? ##
*** last result: cpu time 7,685 ms, real time 3,316 ms.
? version()
[2, 16, 1]
?
New sqrtGintFallKnown(n) determines all gaussian integer square roots of Fn.
- number of square roots for Fn is 2^#primefactors(Fn).
- gaussian integers have positive real and imaginary parts, with real bigger than imaginary
- the vector is sorted by real part of gaussian integers
- the L2-norm of each gaussian integer is Fn (so all are "sum of squares" representations of Fn)
? S=sqrtGintFallKnown(12);
? #S
64
? foreach(S,s,assert(real(s)>imag(s)&&imag(s)>0))
? o=0;foreach(S,s,assert(real(o)<real(s));o=s)
? foreach(S,s,assert(norml2(s)==2^2^12+1))
?
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:~$
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.
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
?
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
First command reads array of Fermat number factors into variable F.
It validates all data while reading, and replaces all prime constants (like "P252")
by their values. Composite numbers (like "C1133" for F12) are not replaced.
Second command outputs F0..F11 with their replaced factors.
Access to the array without replaced "P..." values (for F8 below) is possible as well (use F=F[1] instead of F=F[#F]):