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

[hermann at stamm-wilbrandt.de]

On 7/25/2023 new largest known prime p =1 (mod 4) was proven and 
published:

-----  ------------------------------- -------- ----- ---- 
--------------
  rank  description                     digits   who   year comment
-----  ------------------------------- -------- ----- ---- 
--------------
...
     7d Phi(3,-465859^1048576)          11887192 L4561 2023 Generalized 
unique
...
    10  10223*2^31172165+1               9383761 SB12  2016
...

hermann at 7600x:~$ gp -q
? p=polcyclo(3,-465859^1048576);
? #digits(p)
11887192
? p%4
1
?

I determined sqrt(-1) (mod p) for that prime as well, in 6.7days with 
patched LLR software:
https://github.com/Hermann-SW/11887192-digit-prime#motivation

For some reason the computations from sum of squares to sqrtm1 and back 
with libgmpxx became faster for that bigger prime ...

hermann at 7600x:~/RSA_numbers_factored/c++$ 
./sqrtm1.11887192_digit.largest_known_1mod4_prime
a = y^(-1) (mod p) [powm]; a *= x; a %= p
0.651703s
[M,V] = halfgcdii(sqrtm1, p)
0.189161s
[x,y] = [V[2], M[2,1]]
0s
done, all asserts OK
hermann at 7600x:~/RSA_numbers_factored/c++$


Last, but not least, the 9.4million digit prime of previous email is 
largest Colbert number.
I computed sqrtm1 and sum of squares for the 5 smaller Colbert numbers 
as well.
The stored data is usable with PARI/GP, Python as well as C++ with 
libgmpxx:
https://github.com/Hermann-SW/Colbert_numbers#readme

hermann at 7600x:~/Colbert_numbers$ make
sed "s/C =//;\
      y/[]/{}/;\
      s/\([0-9a-fx][0-9a-fx]*\)/mpz_class\(\"\1\"\)/g;" Colbert.py > 
Colbert.h
g++ validate.cc -lgmp -lgmpxx -O3 -Wall -pedantic -Wextra -o validate

...
time -f %E\\n  ./validate
6 entries of the form [k,n,s,x,y], with p=k*2^n+1, s^2%p==p-1 and 
p==x^2+y^2
   5359*2^5054502+1 (1521561-digit prime)
  33661*2^7031232+1 (2116617-digit prime)
  28433*2^7830457+1 (2357207-digit prime)
  27653*2^9167433+1 (2759677-digit prime)
19249*2^13018586+1 (3918990-digit prime)
10223*2^31172165+1 (9383761-digit prime)
done, all asserts OK
0:03.90
...

Regards,

Hermann.

On 2023-08-06 16:01, hermann at stamm-wilbrandt.de wrote:
> 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.
> _______________________________________________
> gmp-discuss mailing list
> gmp-discuss at gmplib.org
> https://gmplib.org/mailman/listinfo/gmp-discuss