def solve(a, p):
    r = pow(a, (p - 1) // 2, p)
    if (r + 1) % p == 0:
        return None
    print("ok", r)
    if p % 3 == 4:
        return pow(a, (p + 1) // 4, p)
    b = bin(p - 1)[2:]
    s = 0
    n = len(b)
    while b[n - 1 - s] == '0':
        s += 1
    q = (p - 1) >> s
    print(s, q)
    z = 2
    while True:
        res = pow(z, (p - 1) // 2, p)
        print("res :", res)
        if (res + 1) % p == 0:
            break
        z += 1
    print(z)
    m = s
    c = pow(z, q, p)
    t = pow(a, q, p)
    r = pow(a, (q + 1) // 2, p)
    print(m, c, t, r)
    while t > 1:
        i = 0
        while pow(t, 1 << i, p) != 1:
            i += 1
        print("i:", i)
        b = pow(c, 1 << (m - i - 1), p)
        m = i
        c = pow(b, 2, p)
        t = (t * c) % p
        r = (r * b) % p
    print("root :", r)
    print(pow(r, 2, p), a)
    return r     
    
r = solve(51546465468, 643793627263446114641670828109) 

Embed on website

To embed this program on your website, copy the following code and paste it into your website's HTML: