晴耕雨読

working in the fields on fine days and reading books on rainy days

Tonelli-Shanks Algorithm

Tonelli-Shanks Algorithm は平方剰余の根を求めるアルゴリズムです。 以下にPythonでの実装例を載せます。

def legendre(a, p):
    return pow(a, (p - 1) // 2, p)

def tonelli_shanks(a, p):
    if legendre(a, p) != 1:
        raise Exception("not a square (mod p)")

    # Step 1. By factoring out powers of 2, find q and s such that p - 1 = q 2^s with Q odd
    q = p - 1
    s = 0
    while q % 2 == 0:
        q >>= 1
        s += 1

    # Step 2. Search for a z in Z/pZ which is a quadratic non-residue
    for z in range(2, p):
        if legendre(z, p) == p - 1:
            break

    # Step 3.
    m = s
    c = pow(z, q, p) # quadratic non residue
    t = pow(n, q, p) # quadratic residue
    r = pow(n, (q + 1) // 2, p)

    # Step 4.
    t2 = 0
    while True:
        if t == 0: return 0
        if t == 1: return r
        t2 = (t * t) % p
        for i in range(1, m):
            if t2 % p == 1:
                break
            t2 = (t2 * t2) % p

        b = pow(c, 1 << (m - i - 1), p)
        m = i
        c = (b * b) % p
        t = (t * c) % p
        r = (r * b) % p


if __name__ == '__main__':
    from sympy.ntheory.residue_ntheory import sqrt_mod

    ttest = [
        (10, 13), (56, 101), (1030, 10009), (44402, 100049),
	    (665820697, 1000000009), (881398088036, 1000000000039),
        (41660815127637347468140745042827704103445750172002, 10**50 + 577)
    ]

    for n, p in ttest:
        r = tonelli_shanks(n, p)
        roots = [r, p-r]
        r2 = sqrt_mod(n, p)
        roots2 = [r2, p-r2]
        assert (roots[0] * roots[0] - n) % p == 0
        assert (roots[1] * roots[1] - n) % p == 0
        assert min(roots, roots2)
        assert max(roots, roots2)
        print("n = %d p = %d" % (n, p), end=' ')
        print("roots : %d %d" % (r, p - r))

参考文献