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

Tue Aug 22 02:00:31 CEST 2023

```Hermann

Can you please bear with a few very simple questions?

How did you show that this prime had only 1 representation as a 1 mod 4
number? If it has 2 or more, we instantly know it is composite and we
can quickly obtain the factors.

Thanks for explaining this a bit more for us.

- Randall

On 8/20/23 10:54, hermann at stamm-wilbrandt.de wrote:
> 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:
>
> 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
> _______________________________________________
> gmp-discuss mailing list
> gmp-discuss at gmplib.org
> https://gmplib.org/mailman/listinfo/gmp-discuss
```