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

[Joerg Arndt <arndt@jjj.de>]

Btw. if you happen to know the decomposition p = A^2 + B^2
then sqrt(-1) = +- A/B ( = -+ B/A )

Here is why, setting i = sqrt(-1):
p = A^2 + B^2 = A^2 - (i^2) * B^2 == 0
so (i^2) = A^2 / B^2, so i = +- A/B

This is especially nice if B = 1 or 2,
1^(-1) = 1 (of course) and 2^(-1) = (p+1)/2.

Best regards,   jj


* hermann@stamm-wilbrandt.de <hermann@stamm-wilbrandt.de> [Aug 03. 2023 16:10]:
> On 2023-06-20 02:02, hermann@stamm-wilbrandt.de wrote:
> > 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 just finished big project to determine "sqrt(-1) (mod p)" for largest
> known 9,383,761-digit prime p =1 (mod 4):
> https://github.com/Hermann-SW/9383761-digit-prime#readme
> 
> I stopped the 75.4 days expected runtime run after 9d 20.9h, because I
> learned about llr tool based on gwnum library, and computed "sqrt(-1) (mod
> p)" in 13.2h plus 10s post processing. That is 137x faster than with my
> libgmpxx based code.
> 
> gwnum library is x86 specific, but allows for huge speedups using AVX512
> instructions, and multi-core processing a single square or square+mult (mod
> p). And 31,172,177 such operations had to be done ...
> 
> New PARI/GP script with that sqrtm1 constant:
> https://github.com/Hermann-SW/RSA_numbers_factored/blob/main/pari/sqrtm1.9383761_digit.largest_known_1mod4_prime.gp
> 
> Nice, computing sum of squares for that huge prime from sqrtm1 in 2.9s only,
> and computing sqrtm1 given x and y in 4.2s:
> 
> hermann@7600x:~/RSA_numbers_factored/pari$ gp -q
> ? \r sqrtm1.9383761_digit.largest_known_1mod4_prime
> 9383761-digit prime p (31172179 bits)
> [M,V] = halfgcd(sqrtm1, p)
>   ***   last result computed in 2,854 ms.
> [x,y] = [V[2], M[2,1]]
>   ***   last result computed in 1 ms.
> sqrtm1 = lift(Mod(x/y, p))
>   ***   last result computed in 5,929 ms.
> sqrtm1 = lift(x*Mod(1/y, p))
>   ***   last result computed in 4,531 ms.
> sqrtm1 = lift(Mod(x, p)/y)
>   ***   last result computed in 4,203 ms.
> done, all asserts OK
> ?
> 
> 
> Reards,
> 
> Hermann.