mpz_set_str() question

[hermann at stamm-wilbrandt.de]

On 2024-02-29 16:58, Joerg Arndt wrote:
> Isn't it possible to do the Miller-Rabin thing with some
> nonresidues and observe whether you get sqrt(-1) at
> the final squarings?  Once you start with a primitive root
> (which is not terrible unlikely even with just random start
> values) you *will* catch sqrt(-1).
> 
> Best regards,   jj
> 
The patched LLR 4.0.5 software did that, with quadratic non-residue 
616139.
But as you can see, >39 million square(+multiply) steps were needed.
Even done with 16 threads forced on chiplet0 of AMD 7950X CPU that 
computation
did take 8.55 days ...

? p=polcyclo(3,-516693^1048576);
? kronecker(616139,p)
-1
? #binary((p-1)/4)
%8 = 39801737
?


Today I learned from mersenneforum.org @Neptune how to determine
sqrt(-1) (mod p) for p=Phi(3,_) primes instantaneously(!).

The whole single threaded computation of sum of two squares does take 
only 347ms now — quite a drop from 8.55 days ...

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

hermann at 7950x:~$ gp -q
? {
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
%19 = 1
? s^2%p==p-1
%20 = 1
? #digits(p)
%21 = 11981518
?


Regading thread topic, the new top10.sum_of_2_squares.primes repo
linked above has sum of two squares as well as sqrt(-1) (mod p)
stored for top 10 largest known primes p=1 (mod 4).
So these very large numbers can be used as input for mpz_set_str().
Currently only PARI/GP validate code, C++ with libgmpxx will be
committed to repo soon.

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


Regards,

Hermann.