halfgcd() computes in only 347ms sum of two squares for largest known 11,981,518 digit prime =1 (mod 4)

[hermann@stamm-wilbrandt.de]

I learned on this forum on fast "halfgcd()", faster than "gcd()" on 
gaussian integers.

Yesterday Bill showed how to compute sum of two squares for 1133 decimal 
digit
not yet factored composite part of F12:
https://pari.math.u-bordeaux.fr/archives/pari-users-2402/msg00084.html

mersenneforum.org @Neptune was impressed by that and questioned why I 
computed
sum of two squares for largest known 11,981,518 decimal digit prime =1 
(mod 4)
which is Phi(3,-516693^1048576) in 8.55 days. Based on Bill's solution 
with "gcd()"
replaced with "halfgcd()" he did determine sum of two squares for that 
prime really
fast. Like Bill used the fact that 2^2^11 is sqrt(-1) (mod F12) with 
F12=2^2^12+1,
for Phi(3,_) primes sqrt(-1) (mod p) can be computed immediately. And 
with that
sum of two squares:

? {
b=516693;
e=1048576;
p=polcyclo(3,-b^e);
s=b^(e*3/2);
[M,V]=halfgcd(s,p);
[x,y]=[V[2],M[2,1]];
}
? ##
  ***   last result: cpu time 347 ms, real time 347 ms.
? x^2+y^2==p
1
? s^2%p==p-1
1
?

For this very large prime "halfgcd()" is 6.6× faster than "gcd()":

? #digits(p)
11981518
? [M,V]=halfgcd(s,p);
? ##
  ***   last result: cpu time 190 ms, real time 190 ms.
? gi=gcd(s+I,p);
? ##
  ***   last result: cpu time 1,247 ms, real time 1,263 ms.
? norm(gi)==p
1
?


top10.sum_of_2_squares.primes repo runtimes section
https://github.com/Hermann-SW/top10.sum_of_2_squares.primes?tab=readme-ov-file#runtimes

shows that only for the three Proth primes real computation is needed
(with 16 threads forced on chiplet0 of AMD 7950X CPU, with LLR 4.0.5 
software):

rank	description	        digits	        who	year	runtime	comment
7e	Phi(3,-516693^1048576)	11,981,518	L4561	2023	347ms
8	Phi(3,-465859^1048576)	11,887,192	L4561	2023	358ms
11	10223*2^31172165+1	9,383,761	SB12	2016	6:33:01h
15	1963736^1048576+1	6,598,776	L4245	2022	—	x^2+1
16	1951734^1048576+1	6,595,985	L5583	2022	—	x^2+1
17	202705*2^21320516+1	6,418,121	L5181	2021	4:09:02h
19	1059094^1048576+1	6,317,602	L4720	2018	—	x^2+1
21	919444^1048576+1	6,253,210	L4286	2017	—	x^2+1
22	81*2^20498148+1	        6,170,560	L4965	2023	—	x^2+1
23	7*2^20267500+1	        6,101,127	L4965	2022	1:59:20h


It is really cool to have PARI/GP, and especially "halfgcd()".


Regards,

Hermann.