square roots modulo a prime

概要

ある奇素数 \(p\) と \(p\) の平方剰余(Quadratic residue) となる \(y\) について \(x^2 \equiv y\) (mod \(p\)) となる \(x\) を求めるアルゴリズム。

具体的な計算

\(p-1 = 2^mQ\) とする。(\(m \ge 0\), \(Q\) は奇数)

まず平方非剰余となる \(z\) を一つ見つける。 (平方非剰余となる数は1からp-1の間に半分は存在する)
この \(z\) は オイラーの規準(Euler’s criterion) より \(z^\frac{p-1}{2} \equiv z^{2^{m-1}Q} \equiv -1\) を満たす。

\(t \equiv y^Q\) とする。この \(t\) はオイラーの規準 より \(t^{2^{m-1}} \equiv 1\) を満たす。

\(i = m-2\) から 0 まで \(i\) を小さくしながら \(t^{2^i} \equiv -1\) を満たす \(i\) を見つけるたびに \(t\) に \(z^{2^{m-1-i}Q}\) を掛ける。
(この時、掛けたあとの \(t\) は \(t^{2^i} \equiv 1\) を満たす)

最終的に \(t\) は \(\displaystyle t \equiv y^Q z^{2^{j_1}Q} z^{2^{j_2}Q} \cdots z^{2^{j_k}Q} \equiv 1\) \((1 \le j_1 < j_2 < ... < j_k \le m-1)\) の形になるため、
\(\displaystyle x \equiv (ty)^\frac{1}{2} \equiv y^\frac{Q+1}{2} z^{2^{j_1-1}Q} z^{2^{j_2-1}Q} \cdots z^{2^{j_k-1}Q}\) を計算することで解が得られる。

実装

import random
random.seed()
def square_root(n, p):
    n %= p
    if pow(n, (p-1)>>1, p) != 1:
        return -1
    q = p-1; m = 0
    while q & 1 == 0:
        q >>= 1
        m += 1
    z = random.randint(1, p-1)
    while pow(z, (p-1)>>1, p) == 1:
        z = random.randint(1, p-1)
    c = pow(z, q, p)
    t = pow(n, q, p)
    r = pow(n, (q+1)>>1, p)
    if t == 0:
        return 0
    m -= 2
    while t != 1:
        while pow(t, 2**m, p) == 1:
            c = c * c % p
            m -= 1
        r = r * c % p
        c = c * c % p
        t = t * c % p
        m -= 1
    return r
print(square_root(67, 73))
# => "33" or "40"
print(square_root(4, 37))
# => "2" or "35"
print(square_root(6, 37))
# => "-1"
# (12603242261242278713399181313 - 1) = 2^66 * 3 * 17 * 83 * 40351
print(square_root(32189, 12603242261242278713399181313))
# => "4979547606364809208986559660" or "7623694654877469504412621653"

参考

戻る