Tonelli–Shanks algorithm

この記事ではTonelli–Shanks algorithmの原理の「気持ち」を説明します。
間違った説明はしていないつもりですが、（故意に）不正確な記述があるかもしれません。
ループ不変量に着目する考え方はmod_sqrtについて.md · GitHubを参考にしました。
（追記：上記サイトの更に参考元である英語版wikipediaにも同様の記述がありました）
アルゴリズムの説明のみであり、実装については一切説明していません。

Tonelli-Shanks algorithmとは次のような問題を解決するアルゴリズムである。
　問題：整数x$x$と奇素数p$p$、整数n$n$が与えられた時、g2≡xmodpn$g^2\equiv x\mathrm{mod}{p}^n$をみたす整数gを高速に求めよ

より一般に次が解ける。
（p$p$が奇素数ならば、(Z/pnZ)∗$(\mathbb{Z}/p^n\mathbb{Z}{)}^{\ast }$は位数(p−1)pn−1$(p-1)p^{n-1}$の巡回群になることに注意せよ）
　問題：巡回群(G,+)$(G,+)$とその元x∈G$x\in G$が与えられた時、g+g=x$g+g=x$を満たすg∈G$g\in G$を高速に求めよ

以下、群の演算は全て加法的に書いているため、実際に乗法群(Z/pZ)∗$(\mathbb{Z}/p\mathbb{Z}{)}^{\ast }$に対して用いる際は注意せよ

まずは簡単のためG$G$の位数が2ベキの場合を考える
　問題：位数2s$2^s$ (s>0)の巡回群(G,+)$(G,+)$とx∈G$x\in G$が与えられる。g+g=x$g+g=x$となるg∈G$g\in G$を求めよ。

Z/2sZ$\mathbb{Z}/2^s\mathbb{Z}$の代表元を0～2^s-1とする。
　補題：t∈G$t\in G$が与えられた時、tの末尾に0が何bit続くか求めることができる。
　証明：0になるまで2倍して（＝左ビットシフトして）、その回数をsから引くことで分かる。
＝＝＝＝お気持ちパート＝＝＝＝
この補題から「xの下位bitを見ながら順にgを構成する」という方針を思いつくが、これはうまくいかない。
例えばs=6、x=44で考える。（x=0b101100）
まずxの1bit目が0なのでgの0bit目は0。
次にxの2bit目が1なのでgの1bit目は1。
そしてxの3bit目は……と考えたいが、2bit目を0にしないと3bit目は見れない。しかし、2bit目だけを変える元（2^s-4）を見つけるのは難しい。困った……。
＝＝＝＝お気持ちパートおわり＝＝＝＝

そこで発想を転換して「誤差項errorを用意して、ans+ans=x+errorが常に保たれている状態のまま、errorを0にする」という方法を考える。
等式を満たす初期値としてans=error=xが取れる。
errorを0にするには下位bitから消していけば良い。
まず奇数c∈G$c\in G$を勝手にとる。奇数か偶数かの判定は上述の補題により行う。Gの元の半分は奇数なので、ランダムに選ぶことで十分すばやくそのようなものが取れる。
errorのk bit目を消すためにcの2^k倍をerrorに加えるが、この時、ansにもcの2^(k-1)倍を加えることで等式ans+ans=x+errorを保ったままにすることができる。
（errorの初期値はxなので最下位bitは0でり、k>0となることに留意する）
先ほど同様s=6、x=44を例としよう。cは奇数から勝手にとって良いので例えば13とする。
c=001101　←適当に選ぶ
error=101100　←初期値はx
ans=101100　←初期値はx
まずerrorの2bit目を消したいのでcを2^2倍してerror加える。同時にansにcを2^1倍して加える
c=001101
error=100000
ans=000110
次にerrorの5bit目を消したいのでcを2^5倍してerror加える。同時にansにcを2^4倍して加える
c=001101
error=000000
ans=010110
error=0となったのでこうしてans+ans=xとなるansが求められた。
（当たり前だがもう1つの解は最上位bitを反転させることで得られるのでcの2^(s-1)倍を加えれば良い）

それから位数が奇数の場合は簡単。
　補題：位数が奇数Nである巡回群Gの元g∈G$g\in G$に対しては、x=N+12g$x=\frac{N+1}{2}g$とおくとx+x=g$x+x=g$を満たす

一般の場合は、中国剰余定理を使うことで上の2つの場合に帰着できる。
　問題：位数N=q2s$N=q{2}^s$ (qは奇数,s>0)の巡回群(G,+)$(G,+)$とx∈G$x\in G$が与えられる。g+g=x$g+g=x$となるg∈G$g\in G$を求めよ。

先ほどと同様に等式ans+ans=x+errorを満たすようにansとerrorを更新していくことを考える。
初期値はans=q+12x$ans=\frac{q+1}{2}x$、error=qx$error=qx$とする。
＝＝＝お気持ちパート＝＝＝
中国剰余定理からG≅Z/qZ×Z/2sZ$G\cong \mathbb{Z}/q\mathbb{Z}\times \mathbb{Z}/2^s\mathbb{Z}$であるので、Gの元を0≦a≦q−1$0\leqq a\leqq q-1$と0≦b≦2s−1$0\leqq b\leqq 2^s-1$により(a,b)の形で表す。
x=(x1,x2)$x=(x_1,x_2)$とする。
ans=q+12x=(q+12x1,q+12x2)$ans=\frac{q+1}{2}x=(\frac{q+1}{2}{x}_1,\frac{q+1}{2}{x}_2)$、error=qx=(0,qx2)$error=qx=(0,q{x}_2)$である。
この時点で第1成分は既にうまく調整できているので、第2成分だけを更新していきたい。
つまりcとしては(0,奇数)となるようなものを取ってきたい。
ここで、「yがGの平方元⇔第2成分が偶数⇔N/2倍すると0」であることが容易にわかるので、cは非平方元を取ってからq倍して第1成分を潰せば良い。
＝＝＝お気持ちパートおわり＝＝＝
N2c′≠0$\frac{N}{2}{c}^{\prime }\ne 0$となるc'をとり、c=qc′$c=q{c}^{\prime }$とする。
あとは位数2ベキの場合と同様にしてcを使ってerrorの第2成分を下位bitから消していけばよい。