トネリ・シャンクスのアルゴリズム

**トネリ・シャンクスのアルゴリズム** (英：**Tonelli-Shanks algorithm**、シャンクス自身は RESSOL アルゴリズムと呼んでいる) は、奇素数 p を法とする合同算術(剰余算、モジュラー算法、mod算) において、与えられた整数 n (平方剰余) について合同式 r^2 = n \pmod p の解(つまり n の平方根 r)を多項式時間(\Omicron({(\log_2 p)}^4) のオーダー)で求めるためのアルゴリズムである。

このアルゴリズムは合成数を法とする場合には使えない：合成数を法とする平方根を求めることは、素因数分解 と同等の計算問題である

。

このアルゴリズムと同等だがやや冗長なバージョンが、1891年にによって開発された

。

本稿で説明するバージョンは、1973年にによって独立して開発されたものであり、シャンクスは次のように釈明している。

私がこれらの歴史的文献を知るのが遅れたのは、ディクソンの 数論の歴史の第1巻を友人に貸してしまい、返ってこなかったためである。


ディクソン によれば、 トネリの(オリジナルの)アルゴリズムは、素数を法とした場合以外に、素数の累乗 p^\lambda を法とした場合でも x の平方根を求めることができる。

## 予備知識
整数全体 \mathbb{Z} は代数的には環(有理整数環)を成しており、素数 p の倍数の集合 p\mathbb{Z} は \mathbb{Z} の極大イデアルを成している。

\mathbb{Z}/p\mathbb{Z} は 0 から p-1 の整数で構成される集合であるが、これは有限体を成しており、\mathbf{F}_p と表記される。ただし、\mathbf{F}_p における加法と乗法は、\mathbb{Z} におけるpを法とする剰余算(モジュラー算法、mod算)と同じものになる。以降、\mathbf{F}_p における演算は mod算の形式で表記することにする(r^2 = n \pmod p など)。\mathbf{F}_p においては -1 = p-1 \pmod p である。

\mathbf{F}_p から 0 を除いた集合を {\mathbf{F}_p}^\times と表記することにすれば、{\mathbf{F}_p}^\times は \mathbf{F}_p の乗法について位数が p-1 の巡回群(pを法とする整数の乗法群)を成す。従って、{\mathbf{F}_p}^\times の任意の元 n について、

:n^{p-1} = 1 \pmod p

が成り立つ。これをフェルマーの小定理と呼ぶ。

以降、p は奇素数(2ではない素数)とする。p-1 は0でない偶数であるから \frac{p-1}{2} は 0 でない整数となる。{\mathbf{F}_p}^\times の任意の元 n について、(n^{\frac{p-1}{2}})^2 = n^{p-1} = 1 \pmod p  である。\mathbf{F}_p は体であるから、x^2 = 1 \pmod p の解は x = \pm 1 \pmod p のみである。従って、n^{\frac{p-1}{2}} = \pm 1 \pmod p  である。

オイラーの規準によれば、n が {\mathbf{F}_p}^\times で平方根を持つ(つまり、n が平方剰余である)のは、次の場合のみである。

:n^{\frac{p-1}{2}} = 1 \pmod p  

一方、{\mathbf{F}_p}^\times の元 z が {\mathbf{F}_p}^\times で平方根を持たない場合(つまり、z が平方非剰余である)、オイラーの規準によれば、次のようになる。

:z^{\frac{p-1}{2}} = -1 \pmod p  

{\mathbf{F}_p}^\times は位数 p-1 の巡回群であり、その正確に半数が平方非剰余であるので、このような z を見つけるのは難しくない。したがって、以降、このような平方非剰余に常にアクセスできると仮定する。群論的には、全ての平方非剰余は巡回群 {\mathbf{F}_p}^\times の生成元である。

## 基本となるアイデア
以降、n^{\frac{p-1}{2}} = 1 \pmod p  を満たす(つまり平方剰余である) {\mathbf{F}_p}^\times の元 n の {\mathbf{F}_p}^\times における平方根を求めるものとする。

p = 3 \pmod 4 となる素数に対しては、トネリ・シャンクスのアルゴリズムを用いなくても、r = \pm n^{\frac{p+1}{4}}\pmod p と置けば、

:r^2 = (n^{\frac{p+1}{4}})^2 = n^{\frac{p-1}{2}+1} = n^{\frac{p-1}{2}}\cdot n = n \pmod p

であるから、r が求める n の平方根である。

以降は、p = 1 \pmod 4 となる素数 p について考える。

p-1 は、2 で繰り返し割ることで、Q 2^Sと表すことができる。ここで Q は奇数、 S は非負整数とする。p = 1 \pmod 4 であるから p-1 は 4 の倍数である。従って、S の最小値は 2 となる。

:R = n^{\frac{Q+1}{2}} \pmod p

とすると、R^2 = n^{Q+1} = (n)(n^Q) \pmod p となる。t = n^Q \pmod p と置く。t = 1 \pmod p であれば、R は n の平方根となる。それ以外の場合、R と t は次式を満たす。

* R^2 = nt \pmod p；かつ
* t^{2^{S-1}} = n^{Q 2^{S-1}} = n^{\frac{p-1}{2}} = 1 \pmod p

これを群論的に考えると、t は巡回群 {\mathbf{F}_p}^\times の元であり、その位数は 2^{S-1} の約数ということである。この位数を 2^M (M は整数)と置けば、1 \le M \le S-1 である(M が 0 の場合、t の位数は 1 となるが t \ne 1 を仮定しているので、これはあり得ないことになる)。M の値は t を順次2乗することで得られる(最初に1になるまでに2乗した回数が M となる)。

(t^{2^{(M-1)}})^2 = t^{2^M} = 1 \pmod p であるから、前節で説明したように、t^{2^{(M-1)}} = \pm 1 \pmod p である。t^{2^{(M-1)}} = 1 \pmod p であれば、t の位数は 2^{(M-1)} となるので仮定に反する。従って t^{2^{(M-1)}} = -1 \pmod p でなければならない。
 
ここでもし、-1 の 2^M 乗根となるような {\mathbf{F}_p}^\times の元b が存在するとしよう。b^{2^M} = (b^2)^{2^{(M-1)}} = -1  \pmod p であるから t^{2^{(M-1)}} (b^2)^{2^{(M-1)}} = 1 \pmod p である。R、t  と b の関係は次のようになる。

* (Rb)^2 = nt \cdot b^2 = n(t \cdot b^2) \pmod p；かつ
* (t \cdot b^2)^{2^{(M-1)}} = 1 \pmod p。従って、t \cdot b^2 の位数は 2^{(M-1)} の約数である。

t \cdot b^2 の位数を 2^{M'} と置く。M' は明らかにM より小さい整数である。M' の値は tb^2 を順次2乗することで得られる(最初に1になるまでに2乗した回数が M' となる)。

Rb、t \cdot b^2、M' を新たに R、 t、M と書き換えて、以降同じ手順を M が 0 となるまで繰り返せば、t = 1 \pmod p となり、R^2 = nt = n \pmod p となるので、このときの R が求めていた n の平方根となる。

以上のアルゴリズムを群論的に要約すれば、各ステップでは、t の正確な位数を測定し、同じ位数の要素を乗算することで、t をより小さな位数の部分群の生成元に書き換えるということである。

では、-1 の 2^M 乗根となるような都合の良い b をどのように見つけるかが次の問題になるが、これを解決するためのトリックは、既知の平方非剰余である z を利用することである。上に示した z にオイラーの規準を適用すると、z^Q は -1 の 2^{S-1} 乗根であることが示される。したがって、(z^Q)^{2^{S-1}} = z^{Q \cdot 2^{S-1}} = -1 = b^{2^M} \pmod p であるから b = (z^Q)^{2^{S-1-M}} \pmod p と置けば良いことが分かる。

以上がこのアルゴリズムのアウトラインである。

## アルゴリズム
**入力**
:* p: &nbsp; p = 1 \pmod 4 を満たす素数
:* n: &nbsp; {\mathbf{F}_p}^\times の平方剰余

**出力**
:* r: &nbsp; r^2 = n \pmod p を満たす \mathbf{F}_p の元

**入力チェック**
# n = 0 \pmod p の場合、r = 0 を返して終了。
# n が平方剰余でない場合(オイラー規準でチェック)、エラーを返して終了。
# p が奇素数でない場合、エラーを返して終了。
# p = 3 \pmod 4 の場合、r = n^{\frac{p+1}{4}}\pmod p を返して終了。

**初期設定**
# p-1 を 2 で繰り返し割って、p-1 = Q 2^S となる 奇数 Q と 整数 S を求める。
# {\mathbf{F}_p}^\times の平方非剰余を1つ探し z と置く。(z^{\frac {p-1}{2}} = -1 \pmod p)
# ループ変数 t, \, M, \, R, \, c \,  初期値設定
#* M = S \pmod p 
#* t = n^Q \pmod p 
#** n は平方剰余であるので、t^{2^{(M-1)}} = (n^Q)^{2^{(S-1)}} = n^{Q 2^{(S-1)}} = n^{\frac{p-1}{2}} = 1 \pmod p である。
#* R = n^\frac{Q+1}{2} \pmod p  (R^2 = tn \pmod p)
#* c = (z^Q)^{2^{S-1}} \pmod p 
#** z は平方非剰余であるので、c^{2^{(M-1)}} = (z^Q)^{2^{(S-1)}} = z^{Q 2^{(S-1)}} = z^{\frac{p-1}{2}} = -1 \pmod p である。

**ループ不変条件**
ループの各反復処理の開始時には、ループ変数間で以下の関係式(ループ不変条件)が成立する。
:* t^{2^{(M-1)}} = 1 \pmod p
:* R^2 = tn \pmod p
:* c^{2^{(M-1)}} = -1 \pmod p

**ループ**
# t = 1 の場合、r = R を返して終了(ループ終了条件)。
#* ループ不変条件の R^2 = tn \pmod p から R^2 = n \pmod p 
# t を繰り返して2乗し、t^{2^{M'}} = 1 \pmod p となる最小の非負整数 M' を求める。
#* t^{2^{(M'-1)}} = -1 \pmod p
# b = c^{2^{(M-1-M')}} \pmod p とする。
#* ループ不変条件の c^{2^{(M-1)}} = -1 \pmod p から (b^2) = c^{2^{(M-M')}} \pmod p、従って (b^2)^{2^{M'-1}} = c^{2^{(M-1)}} = -1 \pmod p
#* t^{2^{(M'-1)}} (b^2)^{2^{(M'-1)}} = (-1)(-1) = 1\pmod p、従って (tb^2)^{2^{(M'-1)}} = 1\pmod p。 ループ不変条件の R^2 = tn \pmod p から (Rb)^2 = tb^2n \pmod p。
#ループ変数を更新。
#* M \leftarrow M'
#* t \leftarrow tb^2
#* R \leftarrow Rb
#* c \leftarrow b^2
:: 更新されたループ変数はループ不変条件を満たしている。ループ不変条件の t^{2^{(M-1)}} = 1 \pmod p から M' は更新前の M-1 以下であり、M は各反復で厳密に小さくなるため、アルゴリズムは停止することが保証される。

## 例
合同式 *r*2 ≡ 5 (mod 41) を解く。41 は要求どおり素数であり、41 ≡ 1 (mod 4) を満たす。5 はオイラーの基準により平方剰余である: 5^{\frac{41-1}{2}} = 5^{20} = 1 (前述と同様に、(\mathbb{Z}/41\mathbb{Z})^\times 内の演算は暗黙的に mod 41 である)。

# p-1 = 40 = 5 \cdot 2^3  なので、 Q \leftarrow 5、 S \leftarrow 3
# z の値を求める。
#* 2^{\frac{41-1}{2}} = 1 なので、2 はオイラーの定理により平方剰余である。
#* 3^{\frac{41-1}{2}} = 40 = -1 なので、3 は平方剰余ではない。set z \leftarrow 3
# Set
#*M \leftarrow S = 3
#*c \leftarrow z^Q = 3^5 = 38
#*t \leftarrow n^Q = 5^5 = 9
#*R \leftarrow n^{\frac{Q+1}{2}} = 5^{\frac{5+1}{2}} = 2
# Loop:
#* 最初の反復:
#** t \neq 1 なので、まだ終了していない
#** t^{2^1} = 40、t^{2^2} = 1 なので、i \leftarrow 2
#** b \leftarrow c^{2^{M-i-1}} = 38^{2^{3-2-1}} = 38
#** M \leftarrow i = 2
#** c \leftarrow b^2 = 38^2 = 9
#** t \leftarrow tb^2 = 9 \cdot 9 = 40
#** R \leftarrow Rb = 2 \cdot 38 = 35
#* 2回目の反復:
#** t \neq 1 なので、まだ完了していない
#** t^{2^1} = 1 なので i \leftarrow 1
#** b \leftarrow c^{2^{M-i-1}} = 9^{2^{2-1-1}} = 9
#** M \leftarrow i = 1
#** c \leftarrow b^2 = 9^2 = 40
#** t \leftarrow tb^2 = 40 \cdot 40 = 1
#** R \leftarrow Rb = 35 \cdot 9 = 28
#* 3回目の反復:
#** t = 1 となり、終了。 r = R = 28 を返す。

確かに、282 ≡ 5 (mod 41) であり、(-28)2 ≡ 132 ≡ 5 (mod 41) である。したがって、このアルゴリズムは、合同式の2つの解を正しく導出している。

## アルゴリズムの速度
トネリ・シャンクスのアルゴリズムは、(すべての可能な入力(平方剰余と平方非剰余)にわたって平均すると)

:2m+2k+\frac{S(S-1)}{4} +\frac{1}{2^{S-1}} - 9

回の剰余乗算を必要とする。ここで、m は p の2進表現の桁数、k は p の2進表現における1の数である。必要な平方非剰余 z を、ランダムに取った数 y が平方非剰余かどうかを調べることによって求める場合、(平均して) 2回のルジャンドル記号の計算が必要である
。
ルジャンドル記号の2 回の計算の平均は次のように説明される。y は、確率 \tfrac{\tfrac{p+1}{2}}{p} = \tfrac{1 + \tfrac{1}{p}}{2} で平方剰余になる。これは 1 より小さいが、\geq \tfrac{1}{2} である。したがって、平均すると y が平方剰余であるかどうかを 2 回確認する必要がある。

これは本質的に、法 p がランダムである場合、つまり S が p の2進表現の桁数に対してそれほど大きくない場合、トネリ・シャンクスのアルゴリズムが非常にうまく機能することを示している。前述のように、Cipollaのアルゴリズム は、S(S-1) > 8m + 20 の場合(そしてその場合のみ)、トネリ・シャンクスのアルゴリズムよりもうまく機能しする。
しかし、代わりにサザーランドのアルゴリズムを使用して \mathbb{F}_p^{\ast} の 2-シロー部分群における離散対数計算を実行する場合、S(S-1) を O(S\log S/\log\log S) によって漸近的に制限される式に置き換えることができる
。
明示的には、c^e\equiv n^Q となるような e を計算し、R\equiv c^{-e/2} n^{(Q+1)/2} が R^2\equiv n を満たすようにする(n は平方剰余なので、e は2の倍数であることに注意)。

このアルゴリズムでは、平方剰余でない数 z を求める必要がある。このような z を多項式時間で求める決定論的アルゴリズムは知られていない。しかし、一般化されたリーマン予想が正しい場合、平方非剰余 z 、

が存在し、上記の限度まですべての z をチェックし、多項式時間 内で適切な z を見つけることができる。ただし、これは最悪のシナリオであることに留意してしてもらいたい。一般に、z は上記のように平均 2 回の試行(つまり確率1/2)で見つかる。

## 利用
トネリ・シャンクスのアルゴリズムは(当然のことながら)素数を法とする平方根が必要となるあらゆる処理に使用できる。例えば、楕円曲線上の点を求めるのに使用できる。また、デジタル署名におけるラビン署名アルゴリズムや、素因数分解における二次ふるい法のふるい分けステップの計算においても有用である。

## 一般化
トネリ・シャンクスのアルゴリズムは、任意の巡回群((\mathbb{Z}/p\mathbb{Z})^\timesの代わりに)と任意の整数 k に対する k乗根、特に有限体の元のk乗根を求めることに一般化できる
。

同一の巡回群において多数の平方根を求める必要があり、Sがそれほど大きくない場合、2乗位数の元の平方根の表を事前に用意しておくことで、アルゴリズムを以下のように簡略化・高速化でる。
# *p* − 1 から 2 のべき乗を因数分解する。*Q* と *S* を次のように定義する。p-1 = Q2^S (*Q* は奇数)
# R \leftarrow n^{\frac{Q+1}{2}}, t\leftarrow n^Q \equiv R^2/n とする
# 表から b を求め、b^2 \equiv t  とし、R \equiv R/b とする
# *R* を返す。

### トネリのアルゴリズムは法 ''p<sup>λ</sup>'' でも適用可能
ディクソンの「数論」によれば、


A. トネリ は、x^{2}=c \pmod{ p^{\lambda}} の平方根について明確な公式を与えた。


ディクソンの文献には、x^{2}\bmod{p^{\lambda}} の平方根について次の公式が示されている。

:p=4\cdot7+1、またはs=2(この式ではsは2でなければならない)かつa=7であって29=2^{2}\cdot7+1となる場合
::x^{2}\bmod{p^{\lambda}}\equiv cの場合
:::x \bmod{p^{\lambda}}\equiv \pm (c^{a}+3)^{\beta}\cdot c^{(\beta+1)/2}、ただし\beta \equiv a\cdot p^{\lambda-1}となる場合

23^{2} \bmod{ 29^{3}}\equiv 529であること、また\beta = 7\cdot29^{2} のとき

:(529^{7} + 3)^{7\cdot 29^{2}}\cdot 529^{(7\cdot29^{2} + 1)/2}\bmod{ 29^{3}}\equiv 24366 \equiv -23

別の例を挙げると: 2333^{2} \bmod{ 29^{3}}\equiv 4142 そして

::(4142^{7} + 3)^{7 \cdot29^{2}}\cdot 4142^{(7\cdot29^{2} + 1)/2}\bmod{ 29^{3}}\equiv 2333

ディクソンはまた、次の式もトネリによるものであるとしている。

:X\bmod{p^{\lambda}}\equiv x^{p^{\lambda-1}}\cdot c^{(p^{\lambda}-2p^{\lambda-1}+1)/2} ただし、X^{2}\bmod{p^{\lambda}}\equiv c かつ x^{2}\bmod{p}\equiv c。

p=23 と p^{3} を法として用いると、計算は次のようになる。

:1115^{2}\bmod{ 23^{3}}=2191

まず、p を法とするモジュラー平方根を求める。これは、どちらか一方の根に対して、通常のトネリのアルゴリズムで求めることができる。

:1115^{2}\bmod{ 23}\equiv 6 となり、\sqrt{6}\bmod{23}\equiv 11 となる。

これに、トネリの式(上記参照)を適用する。

:11^{23^{2}}\cdot 2191^{(23^{3} - 2 \cdot23^{2} + 1)/2}\bmod {23^{3}} \equiv 1115

ディクソンの文献  は、トネリのアルゴリズムが p^{\lambda} を法として動作することを明確に示している。

## 脚注

## 参考文献
* 
* Daniel Shanks. Five Number Theoretic Algorithms. Proceedings of the Second Manitoba Conference on Numerical Mathematics. Pp.&nbsp;51–70. 1973.
* Alberto Tonelli, Bemerkung über die Auflösung quadratischer Congruenzen.  Nachrichten von der Königlichen Gesellschaft der Wissenschaften und der Georg-Augusts-Universität zu Göttingen. Pp.&nbsp;344–346. 1891.
* 
* Gagan Tara Nanda - Mathematics 115: The RESSOL Algorithm
* Gonzalo Tornaria