Project Euphoria

移行しました -> https://project-euphoria.netlify.app/

Tonelli-Shanksのアルゴリズムについて

はじめに

前回がCTFポエムだったので今回は真面目な事書きます。
今回は平方剰余の問題を解く際に登場するTonelli-Shanksのアルゴリズムについて勉強したのでそれをまとめました。といってもアルゴリズムの解説はだいたい英語版Wikipediaを訳したのに加筆修正を加えただけなので英語が読める方は該当記事を読んだ方が早いかもしれません。
また、CTFのCrypto問題の勉強の一環でこの記事を書いているのでもしかしたら競技プログラミングをやっている方と流儀が異なるかもしれません(特に言語、ここでの実装はPythonです)、そこはご了承ください。
前提知識が多少面倒で平方剰余に関する知識として"ルジャンドル記号"と"オイラー規準"が登場します、と言っても式と意味はそこまで難しくないので適当にググる程度で問題ないと思います(実は筆者もオイラー規準や平方剰余の相互法則の証明は眺めた程度で真面目に追ったり自力でしたりしていない)
それと一番大事なお断りですが、もし誤字脱字、誤植、論理の飛躍、そもそも根本的に違う箇所がありましたらコメントかTwitter(@Xornet_Euphoria)までよろしくお願いします(DMは解放していませんのでお手数ですがリプでお願いします)

解きたい問題

今回解きたい問題は次の通りです

r^2 \equiv a \mod p となる整数rを求めよ、但し、pは奇素数とする

簡単に言えば法の下での平方根を求めよ、ということになります。また、この問題はaの値によって解の存在の有無が変わってきますがここでは解ける場合を与えられたとします(ちなみにその判定については↓のオイラー規準で紹介します)

簡単に解ける場合

ここでpの値によって簡単に解ける場合があることを示します。pは奇素数なので整数kを用いて、p = 4k+1, 4k+3の2通りに表すことが出来ます。この内、後者(p = 4k+3)の場合はフェルマーの小定理を利用することで簡単に解くことが出来ます。
フェルマーの小定理からr^{p-1} = r^{4k+2} = (r^{2k+1})^2 \equiv 1 \mod pとなります。
ここでr^2 \equiv a \mod pであることから(a^kr)^2 \equiv 1 \mod pとなります。
したがってa^kr \equiv \pm 1 \mod pとなることから r \equiv \pm a^{-k} \mod pと簡単に解けます。
ただ、これをp = 4k+1の場合に適用しようとしてもrが上手く残ってくれる形にならないので使えません、そこでTonelli-Shanksのアルゴリズムが使われます。

オイラー規準

前述の通り、この形の問題には解が存在する場合と解が存在しない場合があります。例えばp=7の場合を考え、\mathbb{Z}/7\mathbb{Z}の0以外の要素を2乗してみると1,2,4しか現れません(実際に計算してみてください)。したがって、a = 1,2,4の場合には解が存在することになり、a = 3,5,6の場合には存在しないことになります。解が存在する時はapの平方剰余であるという言い方をします。
これを次のような記号を使って表します。

  1.  \left(\frac a p\right) = 1: (apの平方剰余である)
  2.  \left(\frac a p\right) = -1: (apの平方剰余でない)

この記号を用いて次のオイラー規準という関係が存在します。本記事ではこれを多用します。

\left( \frac a p \right) \equiv a^{\frac {p-1} 2} \mod p

これの証明やなぜ計算結果が\pm1になるのかということに付いてはここでは触れません、原始根を利用した証明等があるので興味がある人は調べるか自分で証明してみてください。
また、p:未満の正の整数の内、半分はpの平方剰余で残りは平方剰余ではないという性質もあります。実際に先程挙げたp=7の例でも平方剰余なものとそうでないものは3つずつになっています。
Tonelli-Shanksのアルゴリズムでは、平方剰余でない方を1つ取ってくる場面が存在します。この性質よりp未満の正の整数をランダムに取ってきたものが平方剰余でない可能性は1/2であるため、数のランダムな取得とオイラー規準による確認を何度も繰り返せば現実的な選択回数の中で平方剰余でないものを手に入れることができます。

アルゴリズムのお気持ち

このアルゴリズムはループを使用します、ということはやっぱり不変条件が登場し、それが満たされているようにループを繰り返しながら解を求めていく形になります。
その不変条件をいきなり提示しても意味が分からないので実際のアルゴリズムを多少追いながらこの不変条件と共に"""お気持ち"""を説明します。

まずpは奇素数なのでp-1は偶数になります、したがって別の奇数Qと正の整数Sを用いてp-1 = Q2^Sと書くことができます。
p=13を例にとると、p-1 = 12 = 2^2 \times 3となることからQ=3, S=2になります。
ではこのQを用いてR :\equiv n^{\frac{Q+1}2}となるような数Rを用意します。アルゴリズム中ではこのR:をループ毎に変えていくことで解であるrを目指すことになります。
これを両辺2乗するとR^2 \equiv a^{Q+1} = a \times a^Q \mod pとなります。Rの2乗がnとある数(ここではa^Q)の積と合同という形になったことから対象である問題の形に近くなりました。そこでR^2 \equiv at \mod pのようにt \equiv a^Q \mod pを用いて表します。
仮に t \equiv 1 \mod pであった場合、R^2 \equiv a \mod pとなるのでRがそのまま解になります。勿論最初から選んだこのようなQがそれを満たすことは珍しいので実際はループを繰り返しながらこのようなtを目指すということになります。
そしてこのR^2 \equiv at \mod pという関係がアルゴリズムで使われるループ不変条件の1つになります。この条件を満たすようにループ毎にRを上手いこと選んでそれに対応するこのようなtを目指します。

続いてこのように選んだtが実は整数Mを用いて1の2^{M-1}乗根である事を示します(指数部が普通の整数Mではなく、M-1のようにしている理由は直ぐわかります)。実はこのMSになります。実際に計算してみると次のようになります。
t^{2^{M-1}} = t^{2^{S-1}} \equiv (a^Q)^{2^{S-1}} = a^{\frac {p-1} 2} \equiv 1 \mod p

ここではapの平方剰余である場合を扱っているため、最後の変形ではオイラー規準によって1になります。
ここで2つ目の不変条件が登場します、もう既に勘の良い読者の方なら気付いているかもしれませんが(これが言いたかっただけです)、t12^{M-1}乗根である、になります。

1つ目の条件をループが回っても満たすにはどうすれば良いでしょうか?ループ毎に変えていくのはR,tであるのでRが定数倍された場合、tはその定数の2乗倍になります。式にすると、ループによってtt'に、RR'になったとして、ある整数bを用いてR' = Rb, t' = tb^2であればR'^2 = R^2b^2 \equiv atb^2 = at' \mod pとなりわかりやすい形にすると、R'^2 \equiv at' \mod pであるためループの前後で条件が保存されていることがわかります。
以上の議論から1つ目の条件を満たすためには単にある係数bを用意してt,Rにそれぞれ掛ければ良いことになります(tには2乗したものを掛ける)。この係数を2つ目の条件を満たすように選びます。
なお、この新たなtを得る段階でt=1となった場合は目的のRが得られたことになるのでそこでアルゴリズムは終了となります。

2つ目の条件ですがt=1となるようなMがどのようになるかを考えます。t12^{M-1}乗根であるという条件からM=1であれば確実にt=1となります。したがってループ毎にMを小さくしていき、この値を目指すことにします。
またこれは完全に私のお気持ちですが、tはその条件から1である可能性があります。但し、12^{M-1}乗根であるということは2^{M-1}乗個の選択肢があるため、tがそのうちの1つである可能性はそこまで高くありません。そこでMを減らしながらtを更新していくことで1が得られる可能性を高くしていくというお気持ちが見えます(結局Mを減らしていけば最終的に1になるので可能性もクソも無いと言えば無いんですが)。
アルゴリズム中ではMt^{2^{M-1}}を満たすもっとも小さいものを探します、目的となるMを一旦iと置くと0 \lt i \lt Mの範囲で見つかった場合にMiを代入して値を更新します。これはMの初期値であるS\log p以下なので1から順に探っていくだけでも十分小さい計算量で出来ますが、二分探索を行うアルゴリズムがより高速らしいです(面倒なので実装は前者の1から探索する方を採用しています)。
これによってt^{2^i} \equiv 1 \mod pであるが、t^{2^{i-1}} \not \equiv 1 \mod pであるような整数i (0 \lt i \lt M)が見つかったとします。ここで重要なのがt^{2^{i-1}} \equiv -1 \mod pが満たされているということです。これは次のようにt^{2^{M-1}} \equiv 1 \mod p平方根を両辺でとることで示されます。

\sqrt{t^{2^{M-1}}} = \left(t^{2^{M-1}}\right)^\frac 1 2 = t^{2^{M-2}} \equiv \pm 1 \mod p

これより、もしMをデクリメントした際の計算結果が1で無かった場合は-1になります。したがってt^{2^{i}} \equiv 1 \mod p \land t^{2^{i-1}} \not \equiv 1 \mod p \Rightarrow t^{2^{i-1}} \equiv -1 \mod pとなります。
これで次のループでパラメータを更新するための、bを求める準備が整いました。

ここで新たなパラメータとしてcが登場します、新しく出たので初期値を使って性質を説明します。初期値はc := z^Qとなります、但しzpの"平方剰余では無い"数です。オイラー規準のところで述べたようにこのようなzは現実的な計算量で取得することが(期待)できます。そしてこのように定義されたc-12^{M-1}乗根であることを示します(初期値なのでM=Sであることを用います)。

c^{2^{M-1}} = \left(z^Q\right)^{2^{S-1}} = z^{\frac{p-1}2} \equiv - 1\mod p

最後の変形はもちろんオイラー規準によるものです。実はこのc-12^{M-1}乗根であるというのもループ不変条件の1つになります。よってループ毎に以上の3性質を満たしているt,M,R,cを用意して、次のループでそれを更新してもそれらが満たされていれば、(数学的帰納法と同様に)同一のプロセスを更に次のループでも適用できることになります。
初期値で満たしていることを既に確認しているのは初期値でこれらを満たしていれば後はループが1つ進んでも満たされていることを示すだけでアルゴリズムの健全性が示されるからです(停止性は別で示します)。

このi, cの2つを用いて係数bb := c^{2^{M-i-1}}のように定義します。これによって既に3つのループ不変条件を満たしているパラメータt, M, R, cを次のように更新します('を付けた変数が更新後の値になる)。

M' = i, \ t' = tb^2, \ R' = Rb, \ c' = b^2

この更新によって1つ目の性質(R^2 \equiv at \mod p)が満たされていることは先程示した通りです(式を満たすように係数bを掛けているので)。したがって2つ目の性質(t12^{M-1}乗根)と出たばかりの3つ目の性質(c-12^{M-1}乗根)を示します。

tの方から示します。t' = tb^22^{M'-1} = 2^{i-1}乗すると

t'^{2^{M'-1}} = \left(tb^2\right)^{2^{i-1}} = t^{2^{i-1}}b^{2^{i}}になります。
ここでb^{2^i} = \left(c^{2^{M-i-1}}\right)^{2^i} = c^{2^{M-1}} \equiv -1 \mod pになります。また、it^{2^{i-1}} \equiv -1 \mod pとなるように選んでいることからこの2つを用いてt'^{2^{M'-1}} = t^{2^{i-1}}b^{2^i} = (-1)^2 = 1 \mod pとなり、2つ目の性質はループが回っても保存されることが示されました。

同様にcについても示します。c' = b^22^{M'-1} = 2^{i-1}乗すると

c'^{2^{M'-1}} = \left( b ^2 \right)^{2^{i-1}} = \left( c^{2^{M-i}} \right)^{2^{i-1}} = c^{2^{M-1}} \equiv -1 \mod pになることから3つ目の性質もループが回った後でも保存されることが示されました。

このようにして3性質を満たすようにしながらMを減らし、t1へ近づけるというのがこのアルゴリズムの"""お気持ち"""だと勝手に思っています。ちなみにMはループ毎に必ず減少することからループ回数に対して狭義の単調減少になります。よって必ずM=1に到達することは示されている為このアルゴリズムが停止することも示されます。

アルゴリズム

お気持ちがある程度分かったところでアルゴリズムを提示します。といっても上記である程度書いたのでシンプルに箇条書きにする程度です、ついでに例とPythonによる実装も置いておきます。

入出力

  • 入力: 奇素数p, 整数a
  • 出力: r^2 \equiv a \mod pとなるような整数r、但しオイラー規準を用いてこのようなrが存在するのを確認したものとする

手順

  1. p-1を2で割り続け、p-1 = Q2^SとなるようなQ,Sを入手する
  2. pの平方剰余ではない値zを取ってくる
  3. 各変数に対し、初期値として次のように値を設定する
    M = S, c = z^Q, t = a^Q, R = a^{\frac{Q+1}2}
  4. 次のループを繰り返す
    1. t=0ならr=0を返してアルゴリズムを終了する(a = 0の場合のみ)
    2. t=1ならr=Rを返してアルゴリズムを終了する
    3. t^{2^{i}} \equiv 1 \mod pとなるようなi \ (0 \lt i \lt M)の中で最も小さいものを取得する
    4. b = c^{2^{M-i-1}}と係数を設定しM, c, t, Rを次のように更新する(プログラミングの流儀にならって代入であることを強調する為に\leftarrowを使っています) M \leftarrow i, c \leftarrow b^2, t \leftarrow tb^2, R \leftarrow Rb

英語版Wikipedia同様p=41, a = 5の場合を例にとります。

p-1 = 40 = 5 \times 2^3より、Q = 5, S = 3になります。
pの平方剰余でない数zとしてz=3を取ってきます(3^{\frac{41-1}2} \equiv 40 = \equiv -1 \mod 41である)。
初期値として
M=S=3, c = z^Q = 3^5 \equiv 38 \mod 41, t = a^Q =5^5 \equiv 9 \mod 41, R = a^{\frac{Q+1}2} =5^{\frac{5+1}2} \equiv 2 \mod 41
を設定しループを開始します。
まず、tは0でも1でも無いのでt^{2^{i}} \equiv 1 \mod pとなるようなiを探します。i=1から順に探っていくとこの様な最小のii=2であることがわかります。
係数b = c^{2^{M-i-1}} = 38^{2^{3-2-1}} = 38 \mod 41を得ることが出来たので次のようにM,c,t,Rを更新します。
M \leftarrow i = 2, c \leftarrow b^2 = 38^2 \equiv 9 \mod 41, t \leftarrow tb^2 = 9 \times 38^2 \equiv 40 \mod 41, R \leftarrow Rb = 2 \times 38 \equiv 35 \mod 41

次のループに入ります、tは0でも1でも無いので再びiを探すとi=1となります。
係数b = c^{2^{M-i-1}} = 9^{2^{2-1-1}} \equiv 9を得ることが出来たので同様に値を更新します。
M \leftarrow 1, c \leftarrow 9^2 \equiv 40 \mod 41, t \leftarrow 40 \times 40 \equiv 1 \mod 41, R \leftarrow 35 \times 9 \equiv 28 \mod 41

次のループに入ります、t=1なのでこの時のR=28が解の1つとなります。
当然-Rも解になるので-28 \equiv 13 \mod 41も解になり最終的に求める解はr = 13, 28の2つになります。
(実際に計算してみると13^2 = 169 \equiv 5 \mod 41, 28^2 = 784 \equiv 5 \mod 41となっていることがわかる)

実装

Pythonで実装します、tonelli_shanks(5, 41)(28, 13)になるはずです。

from Crypto.Util.number import isPrime


def neg(x, n):
    return -x % n


def is_quadratic_residue(a, p):
    if a % p == 0:
        return True

    return legendre_symbol(a, p) == 1


def legendre_symbol(a, p):
    if not isPrime(p) or p == 2:
        raise ValueError("p must be a odd prime number")

    return pow(a, (p-1) // 2, p)


def get_q_s(p):
    q = p - 1
    s = 0
    while q % 2 == 0:
        q //= 2
        s += 1

    return (q, s)


def get_nonresidue(p):
    ret = 2
    while is_quadratic_residue(ret, p):
        ret += 1

    return ret


def tonelli_shanks(a, p):
    if not is_quadratic_residue(a, p):
        return ()

    if a == 0:
        return 0

    q, s = get_q_s(p)
    z = get_nonresidue(p)
    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, neg(r, p))

        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)
        b_pow = pow(b, 2, p)
        m, c, t, r = i, b_pow, t * b_pow % p, r * b % p


if __name__ == '__main__':
    # (28, 13)
    # 28 ** 2 = 784 \equiv 5 \mod 41
    # 13 ** 2 = 169 \equiv 5 \mod 41
    print(tonelli_shanks(5, 41))

結び(筆者近況を含む)

英語版Wikipediaによればこのアルゴリズム巡回群に一般化できるらしく、更にk乗根を求めるようにも一般化出来るらしいです、凄い
今回の記事ではそこまで追いませんが純粋に興味があるのでそのうち調べてみようと思います。

今回この問題を扱ったきっかけですがCryptoHackという暗号向けの常設CTFがあり、その中の数学ジャンルの問題に平方剰余の問題があったことでした。
問題文中でTonelli-Shanksのアルゴリズムについて言及されておりそこで調べていたら、アルゴリズムが何故動くかについても気になったのでこのようなクソ長お気持ち表明文が生成されたという経緯です。
このCryptoHackですが思ったより面白く、難易度も高すぎることは無いので興味のある方は是非挑戦してみてください。記事執筆以外にもこれをきっかけに暗号用のライブラリを整備し直す等、私がCTFに復帰するきっかけにもなっています。

ところで今回のアルゴリズム競技プログラミングでも使われているようで、調べていたらそちらでの実装例が幾つか見つかりました。CTFチームメイトの競プロ勢にも意見を募ったら過去の実装例が返ってきました。そういうわけでもしかしたらこの記事も競技プログラミング勢に捕捉されるかもしれませんが私はCTF勢ですのでもし作法や流儀が異なっていたら申し訳ないです(計算量の議論はしていない上に実装コードもC++では無くPythonであるし高速化も図っていない、試しに2048bitの素数pに設定して問題解いたら2sは確実に超えた)。 競プロ勢の皆さんもこれを機にCTFしませんか?面白いですよ。

そういえば今回この記事を書くにあたって一度HackMDにノートを作ったんですが最近HackMDの使い勝手の良さが気に入ったので技術系文書をそちらに移すことを検討しています。
具体的には編集がリアルタイムで反映される点や、同一プロジェクト間のリンクが張りやすい点が気に入りました。Scrapboxも結構好きだったんですが、Markdownが書けないのが玉に瑕なので最近は個人プロジェクトはHackMDを使っています(それに個人プロジェクトなら公開非公開設定も無料で楽にできる)。
実際にここの全技術記事を移すかどうかと今後どうするかは別にして今後もノート感覚で書いて公開できるぐらい書いたらTwitterやらに投下する方針にしようかなと思っています、その中でデカ目の記事になったら清書してはてブロに投下するかもしれません。
ちなみに今回実際に書いたノートはこちらです

hackmd.io

以下は筆者の近況になります
3月に全くCTF関連のことをせず、4月は進路と労働で精神を破壊されたので最近はあまり元気では無かったのですが、こうやってインプットに対するアウトプットが出来るぐらいにはなったので良かったです。CTFチームも自分の進路やメンバーが皆忙しそうにしていたりであまり出れていないので次はこれを再起動出来たら良いなと思っています。
以前、何かの記事で進路を物理にしたと書きましたが、結局最後まで悩み続け、自分の人生を投じる程魅力を感じなかったことから情報系の進路を目指すことにしました、暗号理論を勉強し始めているのはそれの一環でもあります(研究したい分野が暗号理論なので)。
そういうわけで今後も暗号系の記事が増えるんじゃないかなと勝手に思っています、本当はPwnについても色々やって書きたいことがありますが一旦進路が落ち着くまでは控えめにする予定です。

ここまで読んで頂きありがとうございました。冬場にpwnしていたはずが、前述の通りCrypto熱が再燃したので次回もそんな感じの記事になるんじゃないかと思います。
最後に参考文献を載せておきますのでそちらにも目を通してみてください。

参考文献