Prime factors of Fermat numbers F0..F30

Fermat.gp

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__;

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