Skip to content →

Tonelli-Shanks のアルゴリズム


整数 \(a\) と奇素数 \(p\) が与えられたとき、\(x^2 = a \bmod p\) を満たす \(x\) を \(O(\log^2 p)\) で求める Tonelli-Shanks のアルゴリズムを解説する。計算量は Cipollaのアルゴリズム \(O(\log p)\) に劣るが、2乗以外の冪根では Tonelli-Shanks が優位になる。

掛け算は難しいので、足し算に直す。原始根を一つ固定して、その原始根で表したときの冪指数部分に着目して \(x^2=a \bmod p\) を書き直す。つまり、\(\left((\mathbb{Z}/p\mathbb{Z})\setminus\{0\},\times)\cong(\mathbb{Z}/(p-1)\mathbb{Z},+\right)\) の右辺で考える(左辺を \((\mathbb{Z}/p\mathbb{Z})^\times\) と略記する)。ただし原始根を実際に求めると計算量が悪化するので、右辺で表示したときの情報を左辺での操作から引き出していくことになる。左辺でできる軽い操作と、右辺での対応は次の通り。ただし、\(x\) のべき指数部分を \(\mathrm{ind}(x)\) と置いた。

\((\mathbb{Z}/p\mathbb{Z})^\times\) \(\mathbb{Z}/(p-1)\mathbb{Z}\) 計算量
\(ab\) \(\mathrm{ind}(a)+\mathrm{ind}(b)\) \(O(1)\)
\(ab^{-1}\) \(\mathrm{ind}(a)-\mathrm{ind}(b)\) \(O(\log(b))\)
\(a^k\) \(k\cdot\mathrm{ind}(a)\) \(O(\log(k))\)
\(a^{2^k}\) \(2^k\cdot\mathrm{ind}(a)\) \(O(k)\)
\(1\) \(0\) \(O(1)\)

\(x^2=a\) の解の存在判定はオイラーの規準でできる。以下、解が存在すると仮定する。

中国剰余定理により \(\mathbb{Z}/(p-1)\mathbb{Z}\cong \mathbb{Z}/q\mathbb{Z}\times\mathbb{Z}/2^Q\mathbb{Z}\) (\(q\) : 奇数)と分解できる(+は省略した)。\(a\) を中国剰余定理で分解した右辺で見たとき、\((a_1,a_2)\) だとする。

\(x\) の指数の第一近似として \(\mathrm{ind}(x)=\frac{q+1}{2}(a_1,a_2)\) を選ぶ。\((\mathbb{Z}/p\mathbb{Z})^\times\) に翻訳すると、\(a^{(q+1)/2}\) を \(x\) の第一近似として選ぶということである。\(2\mathrm{ind}(x)=(a_1,(q+1)a_2)\) であるから第1成分は \(a\) と等しくなる。コードにすると次のようになる。

q = p - 1
Q = 0
while q % 2 == 0:
    Q += 1
    q //= 2
# p - 1 = q * 2^Q
x = pow(a, (q + 1) // 2, p)

あとは第2成分を揃えれば良い。\(2\mathrm{ind}(x)-\mathrm{ind}(a)\) の第2成分 \(w\) を誤差項と呼ぶことにする。\(w\) を \(0\) にしたい。

誤差項 \(w\) を二進法表示したときの最下位ビットは \(w\) をいくつ左シフトすると \(\bmod 2^Q\) で 0 になるかを調べることで分かる。この最下位ビットを順に 0 にしたい。\((\mathbb{Z}/p\mathbb{Z})^\times\) に翻訳すると \(k\) ビット左シフトすることは \(2^k\) 乗することに対応し、bit が \(0\) になることは \(1\) になることに対応する。右シフトの難しさと左シフトの簡単さの非対称性がこのアルゴリズムの肝だと思う。

第2成分を2進法表示したときの一桁目が 1 である数 \(b\) をランダムに持ってくる(ヒット率50%)。\(q\) は奇数だから \(b\) を \(q\) 倍すると第1成分が 0 になり、第2成分を2進法表示したときの一桁目は 1 のまま。よって、誤差項の最下位ビットが 0 になるように、\(bq\) を左シフトした数を \(w\) から引き算すれば誤差項の最下位ビットが 0 になる。

以上の操作を誤差項が 0 になるまで繰り返せば良い。誤差項の桁数が \(O(\log p)\) であり、誤差項の最下位ビットを見つけるのに \(O(\log p)\) かかるから全体で \(O(\log^2p)\) で求まることが分かった。

コード(Python)は次のようになる:

import random

# a^{-1} mod p を返す。
def inv(a, p):
    return pow(a, p - 2, p)

# a = x^2 mod p なる x を返す。
# そのような x が存在しなければ -1 を返す。
def tonelli_shanks(a, p) :
    a %= p
    if a == 0:
        return 0
    if p == 2:
        return a

    if pow(a, (p - 1) // 2, p) != 1:
        return -1

    b = 1
    while pow(b, (p - 1) // 2, p) == 1:
        b = random.randint(1, p - 1)

    q = p - 1
    Q = 0
    while q % 2 == 0:
        Q += 1
        q //= 2
    # p - 1 = q * 2^Q
    x = pow(a, (q + 1) // 2, p)
    b = pow(b, q, p)

    shift = 2
    while pow(x, 2, p) != a:
        error = inv(a, p) * pow(x, 2, p) % p # 誤差項
        if pow(error, 1 << (Q - shift), p) != 1:
            x = x * b % p
        b = pow(b, 2, p) % p
        shift += 1
    return x
この見方は sugarknri さんに教えてもらいました。工夫するとほかの冪根にも応用できます。より一般の場合について yukicoder に問題として投稿しているのでよければ覗いてみてください。

Published in アルゴリズム

Comments

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

CAPTCHA