Tonelli-Shanksのアルゴリズム

# 0. はじめに
はじめまして高校二年のもちもちと言います。
自分はセキュリティキャンプ全国大会'21の暗号のゼミの修了生で、またSECCON BEGINNERS(ctf4b)ではcrypto(暗号)の作問をしています。
Tonelli-Shanksのアルゴリズムとは$p$を奇素数として、法$p$においての$a$の平方根、つまり
```math
x^2 \equiv a\ mod\ p
```
となるような$x$を求めるアルゴリズムです。お気持ちとしてはとある後ほど述べる不変条件を保ったまま、近似解を実際の解に近づけていきます

# 1. 前提知識
まずオイラーの基準というものを紹介する。
奇素数$p$と互いに素である$a$に対して
```math
(\frac{a}{p}) = a^{(p-1)/2} \equiv
\begin{cases}
1 & \text{if } a \equiv x^2 \\
-1 & otherwise.
\end{cases}
```
が成り立つ。またこの時の$(\frac{a}{p})$をルジャンドル記号と呼ぶ。
これを一般化したヤコビ記号というものもあるのだがここでは割愛する。

https://ja.wikipedia.org/wiki/%E3%83%A4%E3%82%B3%E3%83%93%E8%A8%98%E5%8F%B7

証明
pは素数なのでラグランジュの定理より$x^2 \equiv a\ mod\ p$は各$a$に対して最大2個の解をもつ。つまりこれは0のほかに$p$を法とする異なる平方剰余が$(p-1)/2$個存在していることを意味している、
またこれは実際に
```math
(p-x)^2 \equiv p^2-2px+x^2 \equiv x^2\ mod\ p
```
となっている。
ここでフェルマーの小定理より互いに素である$a,p$に対し
```math
a^{p-1} \equiv 1\ mod\ p
```
であることを利用し
```math
(a^{(p-1)/2}+1)(a^{(p-1)/2}-1) \equiv 0
```
が言える。すると$a \equiv x^2$である時
```math
a^{(p-1)/2} \equiv x^{p-1} \equiv 1\ mod\ p
```
であるため、先ほどのフェルマーの小定理の式が満たされることとなる。

※ラグランジュの定理
pが素数、$x ∈ Z/pZ$、$f(x)$が整数係数の多項式である時、以下のいずれかが成り立つ
- $f(x)$のすべての係数が$p$で割り切れる
- $f(x) ≡ 0\ mod\ p$は、高々$deg\ f(x)$個の解を持つ($deg$は次数)
(pが素数でないなら、$deg\ f(x)$個よりも多い解を持つ可能性あり)

証明
$g(x) ∈ (Z/pZ)[x]$を、$f(x)$の係数を$mod\ p$で取ることで得られる多項式とする。
- (1)$f(k)$が$p$で割り切れるのは、$g(k) = 0$のとき、またその逆も真である
- (2)$g(x)$は、高々 $deg\ g(x)$個の根を持つ
まず$deg\ g(x)<=deg\ f(x)$なのは明らかである。ここで$g(k)$は$f(k)$を$mod\ p$でとったものであることから(1)は成立する。
また$p$は素数のため$Z/pZ$は有限整域であり、これは体である。さらに体上の非ゼロ多項式は多項式における除算アルゴリズムによりその次数以下の根しか持たない。
$f(k_1)\equiv f(k_2) \equiv 0$の時$k_1 \equiv k_2$である。またこうでない時$k_1 ≢ k_2$となる。
よって(1)による不一致の解の数は$g(x)$の根の数と同じであり、これは(2)により高々$deg\ g(x)$であり、これは高々$deg\ f(x)$となる。

https://en.wikipedia.org/wiki/Lagrange%27s_theorem_(number_theory)

# 2. pが8で割って余りが3,7の時
平方根は
```math
x=a^{(p+1)/4}
```
となる。
証明
先ほどのフェルマーの小定理を利用して
```math
x^2=a^{(p+1)/2}=a*a^{(p-1)/2}=a
```
# 3. pが8で割って余りが5の時
```math
x=a^{(p+3)/8}\ mod\ p
```
とする。この時$x^2≢a$でないなら
```math
x=x*2^{(p-1)/4}\ mod\ p
```
とする。
証明
先ほどと同様に
```math
x^2=a^{(p+3)/4}=a*a^{(p-1)/4}=a\ mod\ p
```
```math
x^2=x*2^{(p-1)/2}=x\ mod\ p
```
# 4. pが8で割って余りが1の時-お気持ち編
$p$は奇素数より
```math
p-1=Q*2^S
```
とおける。ここで
```math
R=a^{(Q+1)/2}
```
とすると、
```math
R^2=a*a^Q
```
となるため、二乗の値が目的の$a$と$a^Q$との籍で表現できた。
そこで
```math
t=a^Q
```
というふうに置いてやると
```math
R^2=at\ mod\ p
```
となり、$t \equiv 1$となった時点で$R$が目的の$x$となったことになる。
実はこの時に選んだ$t$が$M=S$とし、$1$の$2^{M-1}$乗根であることを証明する。
```math
t^{2^{M-1}}=t^{2^{S-1}}=a^{Q{2^{S-1}}}=a^{(p-1)/2}=1\ mod\ p
```
ここで出てきた二つの式
```math
R^2=at\ mod\ p
```
```math
t^{2^{S-1}}=1\ mod\ p
```
が冒頭で述べた不変条件となる。
一つ目の式より$R$が定数倍、つまり
```math
R'=Rb,t'=tb^2
```
である時
```math
R'^2=R^2b^2=atb^2=at'
```
となりこの条件を満たしている。
二つ目の条件だが、アルゴリズムでは$M$が$t^{2^{M-1}}$を満たす最も小さいものを探す。これによって得られた$M$を$i$とおくと$0<i<M$で$i$が見つかった場合に$M<-i$をすることでより小さい$i$を見つけられる。
```math
\sqrt{t^{2^{M-1}}} \equiv t^{2^{M-2}} \equiv ±1
```
であるから$t^{2^{i}} \equiv 1$であるが$t^{2^{i-1}} ≢ 1$の時、$t^{2^{i-1}} \equiv -1$となることがいえる。
ここで$p$の平方非剰余$z$を用い、$c=Z^Q$というパラメーターを用意する。
この時
```math
c^{2^{M-1}} \equiv {(z^Q)}^{2^{S-1}} \equiv z^{(p-1)/2} \equiv -1
```
であることから$c$は$-1$の$2^{M-1}$乗根となり、実はこれも3つ目の不変条件として含まれることとなる。このような$t,M,R,c$を初期値で用意してやると、次のループでもこれが満たされてることが証明されてることから、以後数学的帰納法と同様にこのプロセスを続けられる。
$b$を$2^{M-i-1}$とし、$M'=i,t'=tb^2,R'=Rb,c'=b^2$とする。
一つ目の条件は仮定より明らか、また
```math
t'^{2^{M'-1}}=t^{2^{i-1}}b^{2^i}=(-1)*c^{2^{{M-i-1}^{2^i}}}=(-1)*c^{2^{M-1}}=(-1)^2=1
```
より二つ目も示せた。さらに
```math
c'^{2^{M'-1}}=(b^2)^{2^{i-1}}=c^{{2^{M-1}}^{2^{i-1}}}=c^{2^{M-1}}=-1
```
であることから三つめも示された。
このようにしながら$M$を単調減少させつつ$t$を1に近づけるのがこのアルゴリズムの目標である。

# 5. pが8で割って余りが1の時-アルゴリズム編

- 1
$p-1=Q2^S$となるような$Q,S$を入手する
- 2
$p$の非平方剰余$z$を用意する
- 3
$M=S,c=z^Q,t=a^Q,R=a^{(Q+1)/2}$とする。
- 4
次のループを繰り返す。
    - 1. $t=0$なら$return\ r=0$
    - 2. $t=1$なら$return r=R$
    - 3. $t^{2^i} \equiv 1$となる$0<i<M$を満たす最小の$i$をとってくる
    - 4.$b=c^{2^{M-i-c}}$とし$M \gets i,c \gets b^2,t \gets tb^2,R \gets Rb$とこれらを更新する

# 6. 実装

```py
from Crypto.Util.number import *


def quadratic_residue(a, p):
    a %= p
    if not isPrime(p) or p == 2:
        return -1
    if p % 8 == 3 or p % 8 == 7:
        return pow(a, (p + 1) // 4, p)
    if p % 8 == 5:
        x = pow(a, (p + 3) // 4, p)
        if (x ^ 2) % p == a:
            return x
        else:
            x *= 2 ^ ((p - 1) // 4)
            return x % p
    # p%8==1,tonelli_shanks
    q = p - 1
    s = 0
    while q % 2 == 0:
        s += 1
        q //= 2
    z = 2
    # legendre_symbol
    while pow(z, (p - 1) // 2, p) == 1:
        z += 1
    m, c, t, r = s, pow(z, q, p), pow(a, q, p), pow(a, (q + 1) // 2, p)
    while True:
        if t == 1:
            return r
        i = m
        for j in range(1, m):
            if pow(t, pow(2, j), p) == 1:
                i = j
                break
        b = pow(c, pow(2, m - i - 1, p), p)
        m, c, t, r = i, (b * b) % p, (t * b * b) % p, (r * b) % p


a = 8479994658316772151941616510097127087554541274812435112009425778595495359700244470400642403747058566807127814165396640215844192327900454116257979487432016769329970767046735091249898678088061634796559556704959846424131820416048436501387617211770124292793308079214153179977624440438616958575058361193975686620046439877308339989295604537867493683872778843921771307305602776398786978353866231661453376056771972069776398999013769588936194859344941268223184197231368887060609212875507518936172060702209557124430477137421847130682601666968691651447236917018634902407704797328509461854842432015009878011354022108661461024768
p = 30531851861994333252675935111487950694414332763909083514133769861350960895076504687261369815735742549428789138300843082086550059082835141454526618160634109969195486322015775943030060449557090064811940139431735209185996454739163555910726493597222646855506445602953689527405362207926990442391705014604777038685880527537489845359101552442292804398472642356609304810680731556542002301547846635101455995732584071355903010856718680732337369128498655255277003643669031694516851390505923416710601212618443109844041514942401969629158975457079026906304328749039997262960301209158175920051890620947063936347307238412281568760161
print(quadratic_residue(a, p))
print(
    p
    - 2362339307683048638327773298580489298932137505520500388338271052053734747862351779647314176817953359071871560041125289919247146074907151612762640868199621186559522068338032600991311882224016021222672243139362180461232646732465848840425458257930887856583379600967761738596782877851318489355679822813155123045705285112099448146426755110160002515592418850432103641815811071548456284263507805589445073657565381850521367969675699760755310784623577076440037747681760302434924932113640061738777601194622244192758024180853916244427254065441962557282572849162772740798989647948645207349737457445440405057156897508368531939120
)
```
追記
たまに壊れていることが確認できました。
ひとまずただしい実装はこちらの方のものを参考にしてください。

https://qiita.com/taiyaki8926/items/6720a1c19e9c9afbfd12#%E5%85%A8%E4%BD%93%E3%81%AE%E3%82%BD%E3%83%BC%E3%82%B9%E3%82%B3%E3%83%BC%E3%83%89

参考

yukicoderのものは二乗根ではなく一般根に関する議論が載っています

https://yukicoder.me/problems/no/981

https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm


https://xornet.hatenablog.com/entry/2020/05/11/112508