整数 \(a\) と奇素数 \(p\) が与えられたとき、\(x^2 = a \bmod p\) を満たす \(x\) を \(O((\log p)^2)\) で求める Tonelli-Shanks のアルゴリズムを解説する。計算量は Cipollaのアルゴリズム \(O(\log p)\) に劣るが、２乗以外の冪根では Tonelli-Shanks の方が速くなる（yukicoder981）。
原始根で考える
掛け算は難しいので、足し算に直す。原始根を一つ固定して、その原始根で表したときの冪指数部分に着目して \(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\) と略記する）。
例えば、\(3 \cdot 9 = 27 = 1 \pmod {13}\) と計算するのではなく、原始根 \(2\) を用いて、
\begin{align}
3 \cdot 9 &= 2^4 \cdot 2^8 \\
&= 2^{(4 + 8) \bmod {(13-1)}} \\
&= 2^{0} \\
&= 1 \\
\end{align}
のように計算したいのだ。
残念ながら、原始根を実際に求めると計算量が悪化するので、上記のような計算はできない。しかし、原始根で表示したときの情報を通常の \((\mathbb{Z}/p\mathbb{Z})^\times\) での操作から引き出していくことはできる。この発想が Tonelli-Shanks で最も重要だ。左辺でできる軽い操作と、原始根のべき指数の対応は次の通り。ただし、\(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)
誤差項を0に
あとは第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 p)^2)\) で求まることが分かった。
コード（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
Comments