Cornacchia's algorithm for Gaussian Integers paper

Cornacchia.py

#Complete walkthrough: https://leetarxiv.substack.com/p/computation-of-discrete-logarithms
import sympy as sp

p = sp.Integer("5213619424271520371687014113170182341777563603680354416779")

# Step 1: find root of -2 mod p
r = sp.sqrt_mod(-2, p, all_roots=False)   # choose one root
if r is None:
    raise ValueError("No solution exists.")

# Step 2: Cornacchia
a0, a1 = p, r
while a1 > sp.sqrt(p):
    a0, a1 = a1, a0 % a1

T = a1

# Step 3: compute V
rhs = (p - T*T) // 2
if rhs < 0 or not int(sp.sqrt(rhs))**2 == rhs:
    raise ValueError("No integer V exists; try the other root.")

V = int(sp.sqrt(rhs))

print("T =", T)
print("V =", V)
print("Check:", (T*T + 2*V*V) % p)