def legendre(a, p):
    return pow(a, (p - 1) // 2, p)
 
def tonelli(n, p):
    assert legendre(n, p) == 1, "not a square (mod p)"
    q = p - 1
    s = 0
    while q % 2 == 0:
        q //= 2
        s += 1
    if s == 1:
        return pow(n, (p + 1) // 4, p)
    for z in range(2, p):
        if p - 1 == legendre(z, p):
            break
    c = pow(z, q, p)
    r = pow(n, (q + 1) // 2, p)
    t = pow(n, q, p)
    m = s
    t2 = 0
    while (t - 1) % p != 0:
        t2 = (t * t) % p
        for i in range(1, m):
            if (t2 - 1) % p == 0:
                break
            t2 = (t2 * t2) % p
        b = pow(c, 1 << (m - i - 1), p)
        r = (r * b) % p
        c = (b * b) % p
        t = (t * c) % p
        m = i
    return r

H = 381704527450191606347421195235742637659723827441243208291869156144963
factors = [1504073, 20492753, 59833457464970183, 467795120187583723534280000348743236593] # all of them are coprime

# The Tonelli–Shanks algorithm is used in modular arithmetic to solve for 
# r in a congruence of the form r^2 ≡ n (mod p), where p is a prime number.
# Note: it can be applied only if " Euler's criterion" is fullfiled or "Legendre symbol" = 1
if __name__ == '__main__':
    ttest = [(H,factors[0]),(H,factors[1]), (H,factors[2]), (H,factors[3])]
    for n, p in ttest:
        print("=====================================================")
        print("Legendre symbol: ",  legendre(n, p))
        r = tonelli(n, p)
        assert (r * r - n) % p == 0
        print("n = %d p = %d" % (n, p))
        print("\t  roots : %d %d" % (r, p - r))


# Results:
# p0 = 1504073: 
# 399666 
# 1104407 

# p1 = 20492753:
# 7111848
# 13380905 
  
# p2 = 59833457464970183:
# 34240854883018057
# 25592602581952126 

# p3 = 467795120187583723534280000348743236593:
# 308269479959806774875048102517512730884 
# 159525640227776948659231897831230505709