from collections import Counter
from math import prod

primes = (3, 5, 7, 11, 13, 17, 19, 23, 29, 31)

def factor(m):
    n = m
    f = []
    for d in [2, 3, 5]:
        while n % d == 0:
            f.append(d)
            n //= d
    inc = (4, 2, 4, 2, 4, 6, 2, 6)
    i, d = 0, 7
    while d * d <= n:
        while n % d == 0:
            n //= d
            f.append(d)

        d += inc[i]
        i = (i + 1) % 8
    if n > 1:
        f.append(n)
    return Counter(f)  


def legendre(x, p):
    r = 1 if ((p - 1) * (x - 1)) % 8 == 0 else -1
    q = p % x
    for e in range(1, x):
        if mod_pow(e, 2, x) == q:
            return r
    return -r


def fst_non_residue(p):
    q = (p - 1) // 2
    for _c in primes:
        if mod_pow(_c, q, p) == -1:
            return _c
    return -1

def mod_pow(base, exponent, modulus):
    if modulus == 1:
        return 0
    result = 1
    base = base % modulus
    while exponent > 0:
        if exponent % 2 == 1:
            result = (result * base) % modulus
        exponent >>= 1
        base = (base * base) % modulus
    return result

def hermite(p):
    if p % 4 == 3:
        return None
    elif p == 2:
        return (1, 1)
    elif p % 8 == 5:
        c = 2
    elif p % 24 == 17:
        c = 3
    else:
        c = fst_non_residue(p)
    x0 = mod_pow(c, (p - 1) // 4, p)
    a, b = p, x0
    while 1:
        if b**2 < p:
            break
        a, b = b, a % b
    return b, a % b

def sum_two_squares(n):    
    f = factor(n)
    if all(v % 2 == 0 for v in f.values()):
        isqr = prod(k**(v // 2) for k, v in f.items())
        return [(isqr, 0)]
    if any(e % 2 and p % 4 == 3 for p, e in f.items()):
        return None
    # g, h = f.keys(), prod(k**v for k, v in f.items() if v % 2 == 0 and k % 4 == 3)
    g = []
    for k, v in f.items():
        g.extend([v] * k)
    print(g)
    ans = set([(1, 0)])
    for p in g:
        _ans = set()
        x, y = hermite(p)
        print(f"hermite {p}", x, y)
        for a, b in ans:
            _ans.add((a * x - b * y, a * y + b * x))
            _ans.add((a * x + b * y, a * y - b * x))
        ans = _ans
    return [(abs(x), abs(y)) for x, y in ans]
    
p = 4994669
x, y = hermite(p)
print(x, y, x**2 + y**2 == p)

n = 250
r = sum_two_squares(n)
print(r)