メモ（数論16）: モジュラー平方根
「チラ裏」は、きちんとまとまった記事ではなく、断片的なメモです。誤字脱字・間違いがあるかもしれません。
- 2023-07-09 Shanks のアルゴリズム 導入
- 2023-07-10 Tonelli vs. Shanks w とは AZC である
- 2023-07-11 Tonelli vs. Shanks（その2） w は 1 になっちゃうよ？
2023-07-09 Shanks のアルゴリズム 導入
ある数 a のモジュラー平方根とは x2 ≡ a (mod p) を満たす x のこと――要するに「2乗して p で割ると a 余るような整数」。ここで p は 3 以上の素数とする。
〔例〕 p = 41 としよう: mod 41 の世界では √10 = 16 が成り立つ。 41 × 6 = 246 に注意すると、確かに 162 = 256 = 41 × 6 + 10 ≡ 10 は 41 で割ると 10 余る数だ。この場合、平方してしまうので符号の違いは問題にならず、平方根は ±16 のどっちでもいい。
平方根が正しいかどうかの検算は上の例のように簡単にできるが、検算以前に、どうやって平方根を求めればいいのだろうか…。方法はいろいろある。Tonelli のアルゴリズムは原理が分かりやすく、実用性も高い。特に Shanks による実装は RESSOL* とも呼ばれ、よく知られている。RESSOL は、本質的には Tonelli のアルゴリズムと同じであり、しばしば二人の名前を併記して Tonelli–Shanks のアルゴリズムと呼ばれる。
* Residue（割り算の余り）の世界における Solution（解）なので Res + Sol ということらしい。
Tonelli 版と比べると、Shanks 版（RESSOL）の方が、計算量が節約されて実装上の効率がアップしている。半面 Shanks 版の仕組みは、ベタの Tonelli 版と比べると、少し分かりにくい。いきなり Shanks 版を紹介されると、天下り的に感じられるかもしれない。
このページの目標は三つある。第一に Tonelli 版と Shanks 版の関係を考え、何がどう改善されているのか明らかにすること。二つのバージョンの関係は「分かっている人」にとっては明白かもしれないが、一見したところ、それほど明らかではない。多くの文献でどちらか一つのバージョンが紹介されているものの、両者の関係を記述している資料が全くない。この隙間を埋めたい。
第二に「Tonelli 版には興味がない。モダンな Shanks 版だけ知りたい」という要望に対しても、明快な道筋を示すこと。これらのアルゴリズムに関連して、ほとんど全ての文献において、二重の指数が絡む計算が当たり前のように使われている――慣れている人にとっては何でもない処理だが、ほんの少し説明を追加・工夫するだけで、分かりやすさが格段に向上すると思われる。
第三に、抽象的な観点へのなだらかな橋渡しを試みること。最初は群論的な用語・概念を極力避けて説明を完結させるが、その後、多少は群論的な観点も検討したい…。
このメモは §64 から始まるが、このセクション番号は内部リンクの便宜上のもので、あまり連続性はない。冒頭の §1 の辺りでは
14 + 24 + 34 + … + 104
のような「累乗の和の計算」という、ほとんど無関係の話題を紹介していた。もともとベルヌーイ数の話だったのに、話がそれて、モジュラー平方根の話題になってしまった…。
§64. 上で例に挙げた mod 41 の世界の √10 は（つまり x2 ≡ 10 の解は）、次のようにすれば簡単に求めることができる。
例えば、a6 の平方根は、指数を半分にした a3。事実 a3 を平方すると:
(a3)2 = (a × a × a)2 = (a × a × a)(a × a × a) = a6
平方してしまうのだから、符号は反対でもいい。つまり a6 の平方根は ±a3 といえる。平方根を求めたい数 A が「何かの○乗」の形で表されていて指数○が偶数なら、このように、○を半分にすることで、A の平方根を表すことができる。一方、平方根を求めたい数が 1, 4, 9, 16 のように平方数のとき、その平方根がそれぞれ ±1, ±2, ±3, ±4 などであることは、言うまでもない。例えば
a20 ≡ 9 (mod p)
が成り立つとして、その両辺の平方根を考えると、一応、次が成り立つ:
±(a10) ≡ ±3 (mod p)
この場合の ± は「複号同順」とは限らないけど、とにかく一般に平方根は 2 種類ある。左辺の符号が + の場合だけを考えると:
a10 ≡ ±3 (mod p)
右辺の ± は、実際には + か − のどちらか一方だけが正しい: a10 を p で割った余りは一種類に定まるからだ。 +3 かもしれないし −3 ≡ p − 3 かもしれないが、通常その両方ということはない。
さて、話の大前提として、x2 ≡ a の形の式は a の値によって、解があるときと・ないときがある。普通の整数の世界でも x2 = 16 や x2 = 25 には解があるが、x2 = 17 や x2 = 24 には解がない。同様のことが x2 ≡ a (mod p) においてもいえる。これについては mod 3 の場合や mod 7 の場合について、別の場所で具体的に検討している…。 x2 ≡ a に解がない場合、無理やり a の平方根を求めようとしても、答えが求まるわけがない。一般論として、x2 ≡ a に解があるのか・ないのか事前に分かった方が都合がいい。
解の有無の判定法はいろいろあるが、a を (p−1)/2 乗したとき ≡ +1 になれば解があり、≡ −1 になれば解がない。これを Euler の基準という: p = 41 つまり mod 41 の世界では (p−1)/2 = 20 なので、20乗して ≡ +1 なら解がある。a = 10 は、この条件を満たす:
1020 ≡ +1 『あ』
「両辺の平方根を考える」ことを「開く」と略すことにする。『あ』を開くと:
1010 ≡ ±1 『い』
この場合の ± は、どちらか一方だけが題意に適する。実際に確かめてみると + が題意に適し、1010 ≡ +1 が成り立つ。それをさらに開くと…
105 ≡ ±1 『う』
実際に計算してみると、再び + が題意に適し、105 ≡ +1 が成り立つ。
〔注〕 例えば 105 = 100000 = 41 × 2439 + 1 なので 41 の倍数の違いを無視すると 105 ≡ +1 となる。ここで 105 ≡ −1 ≡ 40 ではないことに注意。10万を 41 で割った余りなので 1 or 40 のいずれか。もちろん「余りが 2 種類ある」なんて変なことは起きない！
結局 +1 ≡ 105 なので、その両辺を 10 倍して 10 ≡ 106、それを開くと √10 ≡ ±103 となる:
√10 ≡ ±103 = ±1000 ≡ ±16 『え』
『え』の ± は「符号はどっちでもいい」。なぜ『い』『う』では片方だけの符号が題意に適して、『え』では符号がどちらでもいのか。『え』は、平方して 41 で割ると 10 余る数という意味なので、平方によって符号の違いが消滅する。『い』『う』は平方しないで左辺を直接 41 で割ると余りが何か？なので、符号の違いは消滅せず、どちらか一種類の符号だけが正しい。
以上を要約すると a のモジュラー平方根を求めるには、次の手順が成り立ちそうだ: Euler の基準 a(p−1)/2 ≡ +1 から始めて、指数が奇数になるまで次々と開き、a奇数 ≡ +1 の形になったら、両辺を a 倍して
a偶数 ≡ a
それを開いて ±a偶数/2 ≡ √a
…このアイデアは果たして正しいであろうか？
mod 41 における √10 に関する限り、確かにそれでうまくいく。もっともそれは、『い』『う』において、正しい符号の選択がたまたま + になってくれたから。開いたとき ≡ −1 になってしまうと、微妙に話が変わってくる。
§65. 引き続き mod 41 において √23 を考えてみたい。
2320 ≡ +1 『ああ』
なので、Euler の基準により √23 は存在する。『ああ』を開いて、正しい符号を選択すると:
2310 ≡ +1 『いい』
これをさらに開いて、正しい符号を選択すると:
235 ≡ −1 『うう』
これは『う』に似ているが、符号がマイナスになってしまった。もし仮に両辺を 23 倍すると 236 ≡ −23、それを開くと:
√−23 ≡ ±233 ≡ ±31
23 の平方根を求めたかったのに −23 の平方根が求まってしまった。これは都合が悪い！
不都合の原因は『うう』の −1 なので、これを何とかして +1 に戻したい。単純思考で『うう』の両辺を −1 倍してみたら、どうだろう？
−(235) ≡ +1
両辺を 23 倍すると −(236) ≡ 23 だが、その左辺は (−1) × 236 だ。上の式を無理やり開くと:
±√−1 × 233 ≡ √23 『ううう』
ますます変な式になってしまった…。
『ううう』は正しい式には違いない。どうにかして mod 41 において −1 の平方根が ±9 であることを突き止められるなら（実際 92 = 41 × 2 − 1 ≡ −1 である）、『ううう』の左辺は ±9 × 233 ≡ ±33 となるが、これは 23 の正しい平方根。では、どうすれば −1 の平方根を突き止められるのか。再び Euler の基準が助けとなる: 「平方根が存在するような数」を 20 乗すれば ≡ +1 になり、「平方根が存在しない数」を 20 乗すれば ≡ −1 になるのだった。第三補充法則によれば mod 41 では 3 の平方根は存在しないので、次が成り立つ:
320 ≡ −1 『ええ』
これを開けば −1 の平方根が ±(310) であることが分かる。さらに良い考えとして、『うう』と『ええ』の左辺同士・右辺同士を掛け算すると:
235 × 320 ≡ +1
この両辺を 23 倍すれば
236 × 320 ≡ 23
それを開けば √23 ≡ ±(233 × 310) となる！
けれど「第三補充法則」なんて知らない場合、あるいは知っていてもその法則が使えない場合（つまり 3 に平方根がある場合）、どうすればいいのだろうか…。これは理論的には非常に難しい問題なのだが、実用上の観点では、ばかばかしいほど簡単な解法がある。Euler の基準で ≡ −1 になる数が見つかるまで、でたらめに選んだ数を次々と (p−1)/2 乗すればいい。mod p において、0 以外の種類の数は、確率 1/2 で平方根を持たない。ランダムな数を幾つか試せば、すぐに
z(p−1)/2 ≡ −1 『お』
を満たす z が見つかるだろう。このギャンブルが10回連続で失敗する確率は:
(1/2)10 = 1/1024
つまり0.1%以下。『お』を満たす z が見つかるまで、でたらめに選んだ z を試す――“解法”と呼ぶには、あまりに無責任なようだけど、実用上それでうまくいくので、以下では「必要なら『お』の性質を持つ z がいつでも利用可能」と認めよう: mod 41 の例では『ええ』のように z = 3 が利用可能。『お』の性質さえ持てば、他の数（例えば z = 7）でも構わない。この性質を持つ z は平方非剰余、略して非剰余と呼ばれ、Tonelli–Shanks のアルゴリズムにおいて（より一般的に、この種のいろいろなアルゴリズムにおいて）重要な役割を果たす。
§66. 整理すると: a の平方根が存在するなら a(p−1)/2 ≡ +1 が成り立つので、a の指数が奇数になるまで次々と開いていけばいい。開いたとき ≡ +1 になればいいけど、≡ −1 になってしまったら、『お』を掛け算することで ≡ +1 に戻してやる。開くたびに a の指数は半分になる（つまり 2 で割り算される）。無限回は 2 で割れないので、遅かれ早かれ a の指数は奇数になる。そうなったとき、両辺を a 倍して開けば、機械的に a の平方根が求まる。
mod 41 で √5 を求めたいとしよう。520 ≡ +1 なので、事実この平方根は存在する。それを開くと:
510 ≡ −1
不都合が起きるけど『ええ』を使って補正すると:
510 × 320 ≡ +1
それを開くと:
55 × 310 ≡ −1
再び不都合が起きるけど再び『ええ』を使って補正すると:
55 × 310 × 320 ≡ +1
両辺を 5 倍すると:
56 × 310 × 320 ≡ 5
これを開けば √5 ≡ ±(53 × 35 × 310) ≡ ±28
この方法がうまくいくのは、「偶数乗」の積 ≡ a を開くのは簡単だから。しかも a奇数 が現れた時点で、左辺にある「補正」の係数は z偶数 だから。つまり次の形になっている:
a奇数 × z偶数 × z偶数 × … ≡ +1
その両辺を a 倍すれば:
a偶数 × z偶数 × z偶数 × … ≡ a
となって簡単に開くことができる。
なぜ a奇数 が現れたとき z の指数は偶数だけと言い切れるか？ この計算は (p−1)/2 乗から始まって、次々に指数が半分になるわけだが、偶数 p−1 が 2 で何回割り切れるか？には限度がある: p−1 を 2 で e 回割ると奇数になるが、それより少ない回数 2 で割っても商はまだ偶数だとしよう: 例えば p = 41 のとき p−1 = 40 は 2 で 3 回まで割り切れる（つまり 23 = 8 で割り切れる）が、2 で 4 回は割り切れない（つまり 24 = 16 では割り切れない）。計算の初期設定 a(p−1)/2 ≡ +1 は Euler の基準により保証されているので、どんなに早く補正が必要になるとしても、それは
a(p−1)/4 ≡ −1
の時点、つまり p−1 を 2 で 2 回割った時点である。その補正結果は:
a(p−1)/4 × z(p−1)/2 ≡ +1
つまり z の指数は a の指数の 2 倍。もっと後で――例えば a(p−1)/8 が現れたときに――z(p−1)/2 による補正がかかるとすれば、そのとき z の指数は a の指数の 4 倍。要するに z の指数は、少なくとも a の指数の 2 倍、場合によっては 4 倍・8 倍・16 倍…なので z の指数は（a の指数と比べて）少なくとも 1 回多く 2 で割り切れる。そのため次々と開いてとうとう a の指数が奇数になった瞬間、z の指数はまだ（少なくとももう1回）2 で割り切れるので、偶数。
以上が元祖 Tonelli のアルゴリズムの原理。処理の流れは、別に難しくない。次回、このアルゴリズムを一応定式化した上で、そこに「節約の余地」があることを観察したい――実装を工夫して計算量を減らしたのが Shanks のアルゴリズムに当たる。（続く）
2023-07-10 Tonelli vs. Shanks w とは AZC である
Shanks のアルゴリズムは、それ自体として直接理解することも可能で、その方が分かりやすいかもしれない（この観点は後述）。けれど最初は Tonelli のアルゴリズムとの関係について、明らかにしたい。このメモでは、簡単な数値例によって、その入り口の部分を紹介する。
§67. 前回に引き続き mod 41 の例。 a が平方根を持てば:
a20 ≡ +1
が保証されている。この出発点を「第1ステップ」と呼ぶことにしよう。指数は p−1 = 40 を 2 で 1 回割ったもの（つまり 40 の半分）。
この式を開いたとき、次のどちらかが起きる:
☆ a10 ≡ +1
または
★ a10 ≡ −1
これを第2ステップと呼ぼう。指数は p−1 = 40 を 2 で 2 回割ったもの（つまり 40 の 4 分の 1）。a を実際に 10 乗して ≡ ±1 のどちらになるか確かめる必要がある。もし☆の +1 なら何もしなくていいが、★の −1 なら両辺に z20 を掛け算して値を ≡ +1 に戻す必要がある。つまり:
☆ a10 ≡ +1
または
★ a10 × z20 ≡ +1
第2ステップの（必要に応じて補正を加えた後の）式をさらに開くと、次のどれかが起きる:
☆☆ a5 ≡ +1 もしくは
☆★ a5 ≡ −1
または
★☆ a5 × z10 ≡ +1 もしくは
★★ a5 × z10 ≡ −1
これを第3ステップと呼ぼう。a の指数は p−1 = 40 を 2 で 3 回割ったもの（つまり 40 の 8 分の 1）。実際に a の 5 乗を含む計算をして ≡ ±1 のどちらになるか確かめる必要がある。もし☆☆または★☆の +1 なら何もしなくていいが、☆★または★★の −1 なら両辺に z20 を掛けて補正する必要がある。つまり:
☆☆ a5 ≡ +1 もしくは
☆★ a5 × z20 ≡ +1
または
★☆ a5 × z10 ≡ +1 もしくは
★★ a5 × z10 × z20 ≡ +1
これで a奇数 × z偶数 ≡ +1 になったので、その両辺を a 倍して開けば a の平方根が求まる。☆☆の場合、実質、補正は何も必要ないが、形式的に「0乗補正」 z0 = 1 が掛け算されていると解釈してもいい: z0 も z偶数 には違いない。
原理はシンプルで分かりやすいけど、実際の計算上では、もう少し効率的に整理することができる。まず a5 やら a10 やら a20 やらを、毎回 a1 から始めて計算するのは冗長だろう。最初に1回だけ A ≡ a5 を求めておけば、a10 と a20 はそれぞれ A2 と A4 になる。その方が計算も楽だし、見通しもいい。同様に、最初に1回だけ Z ≡ z5 を求めておけば、補正に関連する z10 と z20 はそれぞれ Z2 と Z4 になる。
実際、第1ステップについては、こう書くことができる:
(AZC)4 ≡ +1
ただしこの段階では C = 0, ZC = 1 なので、上の左辺は (A)4 ≡ +1 つまり (a5)4 ≡ +1 と同じ意味。
〔注〕 Tonelli バージョンの (AZC)B 表記については§43以下で具体例（§45で一般の場合）を解説しているが、ここではそれと無関係に、一応最初から説明する。上の a20 → a10 → a5 から分かるように A の指数は各ステップで半減: 丸かっこの外側の 4 → 2 → 1 に当たる。一方、補正が必要なときには毎回 z20 ≡ Z4 が掛け算される。これは指数 C が増えることを意味するが、丸かっこの外側が 4 → 2 → 1 なので、内側にある C の指数に足し算される数（足し算の結果である C 自体の値ではない）は 1 → 2 → 4 と倍増していく。
同様に、第2ステップの結論については、こう書くことができる:
(AZC)2 ≡ +1
ここでは補正がなければ C = 0 で ZC は無いのと同じことだが、★のケースでは補正が加わり C = 2, ZC = Z2 ≡ z10。その場合、上の式は次の意味を持つ:
(Az10)2 ≡ (a5z10)2 ≡ a10z20 ≡ +1
最後に、第3ステップに結論については:
(AZC)1 ≡ +1
第3ステップで新たに補正が加わらなければ C の値は、第2ステップと同じ。新たに補正が加わるなら、左辺は z20 倍、つまり Z4 倍されるのだから、丸かっこ内に Z4 が掛け算され、指数を一つの C にまとめるなら、C の値は 4 増加する。例えば、次のように:
古い C の値は 2、新しい C の値は 6:
(AZ2 × Z4)1 = (AZ6)1 ≡ +1
いずれにしても Z の指数 C は偶数（0 を含む）、A の指数は奇数（これは a奇数 を意味する）であり、両辺を a 倍して開けば a の平方根が得られる。
§68. 上記の計算には、さらに工夫の余地がある。丸かっこ内の AZC の値を細かく計算する限りに、w ≡ AZC と置いて w の値をまとめて更新していくことができる。もちろん第1ステップの初期値は C = 0, w = A だ。第2ステップにおいて、もし補正が必要ないなら C の値は変化しないので w の値も変化しない。そのような場合、第2ステップを考える必要すらなく、1ステップ飛ばして直ちに第3ステップに進むことができる。一方、第2ステップで補正が必要な場合…。上記の計算過程を検討すると、その場合、古い w を Z2 倍（つまり z10 倍）したものを新しい w とすればいい。
〔注〕 補正には、両辺を z(p−1)/2 ≡ −1 倍つまり z20 ≡ Z4 倍する必要がある。この場合、丸かっこの外に「2乗」があるので、丸かっこ内の値つまり w が Z2 倍されれば、左辺は全体として (Z2)2 倍され、目的が達成される。
同様に、第3ステップで補正が必要なら、古い w を Z4 倍（つまり z20 倍）したものを新しい w とすればいい。この場合、丸かっこの外には「1乗」しかないので、z20 ≡ Z4 倍の効果を得るためには、そのまんま Z4 を掛ける必要がある。
すなわち、第2ステップ・第3ステップで補正が必要なら、w の値がそれぞれ Z2 倍・Z4 倍される（補正が必要なければ、ステップ自体を省略して次のステップに進んでいい）。その際、倍率の Z2 や Z4 を毎回一から計算するのも面倒なので、例えば y = Z とでも置いておいて、次のようにすれば処理がシンプルになるだろう:
第2ステップの補正は y2 倍（古い w の y2 倍を新しい w にする） ♪
このステップでの y2 倍とは Z2 倍に他ならない
この補正の後で、現在の y = Z の 2 乗をあらためて y とする（新しい y は Z2 に当たる）
第3ステップの補正は y2 倍（古い w の y2 倍を新しい w にする） ♪
このステップでの y2 倍とは Z4 倍に他ならない
この二つの ♪ は同一の処理なので、単純な反復計算になる。
要するに、y = Z から始めて、各ステップで古い y の 2 乗をあらためて y とする。補正が必要なステップでは、古い w の y2 倍（この y は古い y）をあらためて w とする。それだけのことで (AZC)B の丸かっこ内の値 w = AZC が正しく求まる（B は丸かっこの外側の指数で、B = 4, 2, 1 のように、各ステップで半減。指数 C については、明示的に計算する必要すらない）。
Z ≡ z5 は mod 41 の場合の例だった。一般には、法 p から 1 を引いた偶数を 2 で割れるだけ割って（e 回割れるが e+1 回は割れないとしよう）、次の形にする:
p−1 = 2eq ここで q は奇数
この場合も y = Z (≡ zq) から始めて、各ステップで y の値を次々と 2 乗し、補正が必要なら、古い w の y2 倍を新しい w とすることができる。
この部分の計算に関する限り y = Z の代わりに初めから 2 乗して Y ≡ Z2 とでも置けば、補正処理は y2 倍の代わりに単に Y 倍になる。けれど、アルゴリズム全体としては、y2 と y の 2 種類の数を使い分けた方が便利なのだ（後述）。それに各ステップで古い y の 2 乗を新しい y とするのだから、古い y を単に y、新しい y を大文字の Y とすると Y = y2 であり、「y2 倍の代わりに単に Y 倍の方がシンプルでは？」という発想も、うまくやればアルゴリズムに組み込むことができそう…。
各ステップで y や w の値が更新され、「古い y と新しい y」「古い w と新しい w」というコンセプトが生じること（同じ変数名が新旧で別の値を持つこと）は、Shanks のアルゴリズムが若干分かりにくい原因となり得る。でも「同じ変数名の値が次々と更新されていく」というのはごく普通のループ処理で、本質的に難しいことは何もない: 各ステップで古い y の 2 乗が新しい y（この値を便宜上 Y と呼ぼう）になり、古い w の Y 倍（古い y から見ると y2 倍）が新しい w になるというだけ…。
「1ステップずつ」進めるなら話は簡単だが、Shanks のアルゴリズムの特徴として「必要ないステップはスキップできる」。例えば、第2ステップが必要なければ、第1ステップから第3ステップに直行できる。関連して二つの問題が生じる。
第一に「必要ないステップとは何か。必要なステップだけをたどるとして、どのステップからどのステップに直行できるのか」。第二に「必要なステップだけをたどる場合、y や w の新旧の値はどうなるか」。全ステップを真面目にやれば、毎回、古い y の 2 乗が新しい y になるだけだが、例えばステップを1個飛ばすと、古い y の 4 乗が新しい y になるし、ステップを2個飛ばすと古い y の 8 乗が新しい y になる。ステップごとの単純処理に比べると、処理内容が動的に変わり、ちょぴり複雑度が上がる――だけど必要ないステップを省略できるというのは、アルゴリズムの効率の上でとっても良いこと。ほんの少しの工夫（複雑化）で、場合によっては2倍・4倍といった高速化が得られるなら、ぜひやるべきだろう！
これは「実装上の工夫」であり、本質的には Shanks のアルゴリズムと Tonelli のアルゴリズムと同じ。どちらか一方を研究すれば十分ともいえる。でも Shanks のアルゴリズムは――実用上、効率が良いというだけでなく――数論的にも面白い観点を含んでいて、それ自体としても研究する価値がある。次回は「ステップを飛ばす」部分に取り組みたい。
2023-07-11 Tonelli vs. Shanks（その2） w は 1 になっちゃうよ？
具体例だとかえって話が見えにくいので、各ステップで何が起きているのか表にまとめてみる。
§69. mod が p = 97 の場合。 p−1 = 96 が 2 で 5 回も割れる（つまり 25 = 32 で割り切れる）: p = 2eq + 1 の形に当てはめると、指数 e = 5 で奇数 q = 3。この世界で a が平方根を持つ数（平方剰余）なら、Euler の基準から a(p−1)/2 ≡ +1。ここで (p−1)/2 = 2eq/2 = 2e−1q である。具体的数値で言うと p−1 = 96 の半分つまり 48 が 24 × 3 = 16 × 3 に等しい（当たり前）。従って Euler の基準を次のように表現できる:
a(p−1)/2 = a2e−1q = aq⋅2e−1 = (aq)2e−1 が ≡ +1
aq ≡ A と置くと A2e−1 ≡ +1
ついでに 2e−1 を β とすると Aβ ≡ +1 スッキリ♪
p = 97 の場合、A ≡ a3, β = 16, A16 ≡ +1 となる。
| ステップ | w ≡ AZC | B | 要補正なら w は何倍に？ |
|---|---|---|---|
| 1 | AZC | 16 | （補正不要: C = 0） |
| 2 | AZC | 8 | Z2 （C が2増える） |
| 3 | AZC | 4 | Z4 （C が4増える） |
| 4 | AZC | 2 | Z8 （C が8増える） |
| 5 | AZC | 1 | Z16 （C が16増える） |
各ステップで w ≡ AZC の B 乗、つまり wB は ≡ +1 になってほしい（B の初期値は β）。ステップ1では（C = 0 なので ZC = Z0 = 1 は無いのと同じで）この希望は Aβ ≡ +1 を意味し、上記にように既に成り立っている。数値例では (AZC)16 = A16 ≡ +1。それを開いたとき、A8 ≡ +1 になってくれれば一番楽なのだが、確率半々で A16 ≡ −1 になるかもしれず、そうなったら右辺を +1 に戻すために補正が必要なのだった（前回参照）。
補正とは、小文字の z を非剰余として z(p−1)/2 ≡ −1 を両辺に掛けることだった。上記の議論の小文字の a と大文字の A をそれぞれ小文字の z と大文字の Z に置き換えれば、要補正時に掛け算する数は Zβ ≡ −1 に他ならない。数値例では Z16。
1回開くごとに、w の指数 B は半減する。B の初期値は 2 の累乗なので、半々にしていけば、いつかは 1 になり w1 ≡ +1 が達成される。これは…
(AZC)1 = AZC ≡ +1
つまり aqZC ≡ +1
…を意味し q は奇数、C は偶数なので、両辺を a 倍すれば a偶数Z偶数 ≡ a となり、それを開けば a の平方根が得られる。それが Tonelli のアルゴリズムであった。
さて各ステップで補正が必要になった場合、上記数値例で Z16 (≡ −1) を掛け算すればいいのだから、ステップ2が (AZC)8 ≡ −1 になってしまったとすれば、丸かっこ内の値つまり w に Z2 を掛けてやればいい。確かに:
補正前の (AZC)8 から見て
(AZC × Z2)8 = (AZC)8 × (Z2)8 = (AZC)8 × Z16 は Z16 倍
同様に、もしステップ3で (AZC)4 ≡ −1 になっちゃった場合、丸かっこ内に Z4 を掛ければ全体が Z16 ≡ −1 倍されて、補正の目的が達成される。ステップ4で (AZC)2 ≡ −1 になっちゃった場合、丸かっこ内に Z8 を掛ければいい。以下同様。要するに、補正のとき w に掛け算する数は、ステップ2の Z2 から始まって、Z4, Z8, Z16 のように、各ステップごとに平方される――あるステップでのこの倍率を Y とすれば、次のステップでの倍率は Y2 である。次の例のように:
ステップ3での倍率 Y = Z4 から見ると
ステップ4での倍率は Y2 = (Z4)2 = Z8
この計算は1ステップずつ丁寧にやってもいいのだが、実際に補正が必要になるのは ≡ −1 が起きたときだけなので、≡ +1 になってくれたら、何もしないスルーできる。1ステップずつなら、あるステップの補正の倍率を Y として次のステップの補正の倍率は Y2 だが、もし1ステップ飛ばして一気に2ステップ進めるなら、2ステップ前から見て Y4 になる: だって「あるステップ」の次の（飛ばした）ステップで本来 Y2 になって、そのまた次のステップでさらに平方されるのだから、「あるステップ」から見れば (Y2)2 = Y4 じゃん。表の数値例で、ステップ2からステップ4に直行できるとすると、倍率は Z2 から Z8 になるが、後者は前者の 4 乗に他ならない: (Z2)4 = Z8。ステップを2個以上飛ばす場合も同様で、一般に一気に n ステップ、ジャンプできるならジャンプ前の倍率を Y として、ジャンプ後の倍率は Y2n に。
前回の末尾で発生した二つの疑問の答えは、今や明白。第一に「飛ばしていい（必要ない）ステップ」とは、右辺が ≡ +1 になるようなステップ。第二に、補正の倍率の変化は、1ステップごとの計算なら毎回、Y から Y21 = Y2 になるが、一気に2ステップ進めるなら Y22 = Y4 になり、3ステップ進めるなら Y23 = Y8 になり、一般に n ステップ進めるなら Y2n となる。この n に当たる数をアルゴリズムの実装上どうやって求めるかは考えどころだが、処理内容のコンセプトは単純明快。
§70. 上記の数値例は、具体例と言っても A と Z が不定の変数。本当に具体的な数値の例で、もう一度、内容を確認してみる。mod 97 において a = 31 の平方根を求めることにする。この場合、非剰余として z = 5 を選択できる。
A = aq = 313 ≡ 12
Z = zq = 53 ≡ 28
初期値 C = 0, B = 16
ステップ1 w ≡ AZC = A = 12
wB = 1216 ≡ +1
a が平方根を持つ場合、これが +1 になることは Euler の基準によって保証されているので、このステップでは決して補正は必要にならない。言い換えれば、このステップでは w の値は変化しない。
ステップ2 B が半減して 8 に。
wB = 128 ≡ −1
マイナスになってしまったので、補正を発動しよう。既に説明したように――＜表1＞からも一目瞭然だが―― w を Z2 倍すればいい:
この倍率を Y = Z2 = 282 ≡ 8 とすると
新しい w = 古い w × 8 = 12 × 8 = 96 ≡ −1
一応、検算してみると wB = (−1)8 ≡ +1 計画通り +1 に戻ってくれた！
ステップ3 B が半減して 4 に。現在の w の値は −1 なので、結果はすぐ分かる:
wB = (−1)4 ≡ +1
+1 になるので、何もする必要がない！
ステップ4 B が半減して 2 に。現在の w の値は −1 なので、ステップ3と同様に:
wB = (−1)2 ≡ +1
何もする必要がない！ うーん、こいつは楽だ♪
ステップ5 B が半減して 1 に。現在の w の値は −1 なので:
wB = (−1)1 ≡ −1
−1 になってしまった。＜表1＞によると Z16 倍の補正をすればいい。ステップ2で設定した Y から見ると、補正の倍率は 8 乗だ（だって Z16 は Z2 の 8 乗じゃん）:
新しい Y = (古い Y)8 = 88 ≡ −1
［注: アルゴリズムの説明のためにこのように書いているが、実際には Z16 = z(p−1)/2 ≡ −1 は Euler の基準から明白。］
新しい w = 古い w × Y = (−1) × (−1) ≡ +1
最終的に B = 1 になっているのだから、上の結論は wB = (AZC)B = (AZC)1 = AZC ≡ +1 を意味する。もし偶数 C の値が分かるなら AZC ≡ +1 つまり aqZC ≡ +1 の両辺を a 倍して開くことで、いつものように a の平方根を導くことができる。それが Tonelli のアルゴリズムだが、ここで考えている Shanks のアルゴリズムでは C の値を求めず AZC 全体の値 w だけを求めている。すると…
最終的な結論: w = 1 のとき wB = w1 ≡ +1
1 の 1 乗が 1 なんてことは最初から分かり切っている！ この式は何の情報も与えてくれない！ …という感じがする。その通り、これだけでは Shanks のアルゴリズムは目的を達成できない――平行して、もう一つの計算が必要なのだ。その「もう一つの計算」が上記の計算の「ついで」のようにできるので、トータルでは、それでもむしろ効率が良いのだが、その点については次回に検討しよう。今は Tonelli の方法で、結論を出しておく。計算内容と＜表1＞から分かるように、ステップ2の補正のとき、初期値 0 の C に 2 がプラスされる。同様にステップ5の補正のとき C に 16 が補正される。よって Tonelli 風の表記では:
(AZC)1 = AZ2+16 = AZ18 《お》
《お》の値が上記の w ≡ 1 と同じであることを確認しておこう:
AZ18 = 12 × 2818 ≡ +1 《おお》
［注: 確認の意味で記しているだけで、実際には計算するまでもなく、各ステップの結末が必ず ≡ +1 になるように、必要な補正がなされている。］
《おお》の両辺に a = 31 を掛け算:
aAZ18 ≡ a つまり a⋅a3Z18 ≡ a4Z18 ≡ a
それを開いて √a = ±(a2Z9) = ±(312 × 289) ≡ ±82 ≡ ∓15
このことから √31 = ±15 という答えを得る。
検算 (±15)2 = 225 = 194 + 31 = 97 × 2 + 31 ≡ 31 (mod 97)
確かに ±15 は 31 の平方根。
Tonelli の方法でも、別に大して難しくはないけど《お》以下――AZC ≡ +1 を a 倍して a の平方根を導く部分――は、それなりにゴチャゴチャ計算する必要がある。この末尾の部分について、《お》より前の各ステップの「ついで」に平行処理してしまうのが、Shanks の方法の大きな工夫といえるだろう。
