Fast libgmpxx and libpari computations on 9,383,761-digit biggest known prime p that is =1 (mod 4)

[hermann at stamm-wilbrandt.de]

I determined "sqrt(-1) (mod p)" for that prime p, rank 10 of largest 
known primes list
https://t5k.org/primes/lists/all.txt

  rank  description                     digits   who   year comment
-----  ------------------------------- -------- ----- ---- 
--------------
...
    10  10223*2^31172165+1               9383761 SB12  2016
...

in 13.2h with llr tool with 24 threads:
https://github.com/Hermann-SW/9383761-digit-prime#fast-sqrt-1-mod-p-for-9383761-digit-prime-p-1-mod-4


 From that I determined unique sum of squares of p=x^2+y^2.
The 4,691,881- and 4,691,880-digit numbers x and y are defined in C++ 
code
https://github.com/Hermann-SW/RSA_numbers_factored/blob/main/c%2B%2B/sqrtm1.9383761_digit.largest_known_1mod4_prime.cc

by just mpz_class without issues:
...
     mpz_class x("223757 ... 534644");
     mpz_class y("236151 ... 476249");
...

That demo code starts with x,y and p and computes sqrtm1 from that in 
only 4.23s (i7-11850H CPU).
Then it uses libpari "halfgcdii()" function to compute x and y from just 
sqrtm1 and p in 3.72s.
Computing x,y from sqrtm1 is possible with gaussian integer gcd, but 
that is orders of magnitude slower than using "halfgcdii()".
All intermediate results are verified with asserts.

$ f=sqrtm1.9383761_digit.largest_known_1mod4_prime
$ g++ $f.cc -lgmp -lgmpxx -O3 -o $f -lpari -DPARI
$ ./$f
a = y^(-1) (mod p) [powm]; a *= x; a %= p
4.22922s
[M,V] = halfgcdii(sqrtm1, p)
3.71779s
[x,y] = [V[2], M[2,1]]
1e-06s
done
$


Nice that such fast computations for more than 31million bit numbers are 
possible with libgmpxx.

Regards,

Hermann.