efficient transformations between "sum of squares" and "sqrt(-1) (mod p)" for 1M-digit prime

[hermann@stamm-wilbrandt.de]

I learned on this mailing list about "halfgcd()" allowing to efficiently 
compute sum of squares given "sqrtm1 = sqrt(-1) (mod p)", and for 
computing "sqrtm1" from sum of squares as "x*y^(-1) (mod p)".


I implemented those computations for

C++ with libgmpxx and PARI
https://github.com/Hermann-SW/RSA_numbers_factored/blob/main/c%2B%2B/sqrtm1.smallest_known_1million_digit_prime.cc

Python with cypari
https://github.com/Hermann-SW/RSA_numbers_factored/blob/main/python/sqrtm1.smallest_known_1million_digit_prime.py

PARI/GP
https://github.com/Hermann-SW/RSA_numbers_factored/blob/main/pari/sqrtm1.smallest_known_1million_digit_prime.gp

for smallest known 1,000,000-digit prime number 
10^999999+308267*10^292000+1.


.cc/.py/.gp runtimes for "x*y^(-1) (mod p)" are (GP slower than the 
other languages):
249/268/309 ms

.cc/.py/.gp runtimes for "halfgcd()" are (big overhead for cypari call 
of "halfgcd()"):
267/445/269 ms

Same numbers in diagram:
https://stamm-wilbrandt.de/images/xy_sqrtm1__xy.transformation_times.png


Outputs, with single process running and showing near maximal boost CPU 
frequency of 4.8GHz of i7-11850H:


c++$ perf stat ./sqrtm1.smallest_known_1million_digit_prime 2>err
a = y^(-1) (mod p) [powm]; a *= x; a %= p
0.249437s
[M,V] = halfgcdii(sqrtm1, p)
0.266697s
[x,y] = [V[2], M[2,1]]
0s
done
c++$ grep GHz err
     2,805,719,697      cycles:u                  #    4.701 GHz
c++$


python$ perf stat python sqrtm1.smallest_known_1million_digit_prime.py 
2>err
1000000 -digit prime p ( 3321925  bits)
PARI stack size set to 1000000000 bytes, maximum size set to 1000001536
[M,V] = pari.halfgcd(sqrtm1, p)
444.5ms
[x,y] = [mpz(V[1]), mpz(M[1,0])]
0.5ms
sqrtm1 = x * pow(y, -1, p) % p
268.4ms
done, all asserts OK
python$ grep GHz err
     4,903,646,857      cycles:u                  #    4.606 GHz
python$


pari$ perf stat gp-2.15 -q < 
sqrtm1.smallest_known_1million_digit_prime.gp 2>err
1000000-digit prime p (3321925 bits)
[M,V] = halfgcd(sqrtm1, p)
  ***   last result computed in 269 ms.
[x,y] = [V[2], M[2,1]]
  ***   last result computed in 0 ms.
sqrtm1 = lift(Mod(x/y, p))
  ***   last result computed in 432 ms.
sqrtm1 = lift(x*Mod(1/y, p))
  ***   last result computed in 311 ms.
sqrtm1 = lift(Mod(x, p)/y)
  ***   last result computed in 309 ms.
done, all asserts OK
pari$ grep GHz err
     7,549,438,197      cycles:u                  #    4.703 GHz
pari$


Regards,

Hermann.