メモ（数論16）: mod p の平方根（Shanks 版）
「チラ裏」は、きちんとまとまった記事ではなく、断片的なメモです。誤字脱字・間違いがあるかもしれません。
- 2023-07-09 Shanks のアルゴリズム 導入
- 2023-07-10 Tonelli vs. Shanks w とは AZC である
- 2023-07-11 Tonelli vs. Shanks（その2） w は 1 になっちゃうよ？
- 2023-07-12 Tonelli vs. Shanks（その3） 正義の変数セーラー v
- 2023-07-16 Tonelli vs. Shanks（その4） y のスルー
- 2023-08-05 Tonelli vs. Shanks（その5） スルーの実装
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 になる。
〔注〕 小文字の z と大文字の Z の違いを強調するため、今回に限って、大文字の Z を太字にしておく。
実際、まず第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 と倍増していく: つまり、各ステップにおいてもし補正が必要なら、ステップに応じて、現在の ZC の Z1 倍 → Z2 倍 → Z4 倍が補正後の新しい ZC となる（実際には、第1ステップでは補正は不要）。
同様に、第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 になってくれれば一番楽なのだが、確率半々で A8 ≡ −1 になるかもしれず、そうなったら右辺を +1 に戻すために補正が必要なのだった（前回参照）。
補正とは、小文字の z を非剰余として z(p−1)/2 ≡ −1 を両辺に掛けることだった。上記の議論の小文字の a と大文字の A をそれぞれ小文字の z と大文字の Z に置き換えれば、要補正時に掛け算する数は Zβ ≡ −1 に他ならない。数値例では Z16 ≡ −1。
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 倍
* ステップ2が始まった時点では C = 0, ZC ≡ 1 なので、実質的にはこの式は A8 ≡ −1 を意味する。
同様に、もしステップ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個以上飛ばす場合も同様で、一気に s ステップのスルーが可能なら、ジャンプ前の倍率を Y として、スルー後の倍率は Y2s に。
前回の末尾で発生した二つの疑問の答えは、今や明白。第一に「スルーていい（必要ない）ステップ」とは、右辺が ≡ +1 になるようなステップ。第二に、補正の倍率の変化は、1ステップごとの計算なら毎回、Y から Y21 = Y2 になるが、一気に2ステップ進めるなら Y22 = Y4 になり、3ステップ進めるなら Y23 = Y8 になり、一般に s ステップ進めるなら Y2s となる。この s に当たる数をアルゴリズムの実装上どうやって求めるかは考えどころだが、処理内容のコンセプトは単純だ。
§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 の方法の大きな工夫といえるだろう。
2023-07-12 Tonelli vs. Shanks（その3） 正義の変数セーラー v
元祖 Tonelli 版では (AZC)1 = AZC ≡ +1 の両辺を a 倍することで、機械的に平方根が求まる。丸かっこ内の AZC を単に w と置く Shanks 版では、この最終ステップは (w)1 = w ≡ +1 に当たる: 要するに w ≡ 1 の 1 乗が 1 という当たり前の式。当たり前過ぎて、何の情報も与えてくれないように思える。
鍵となるのは、次の事実。「w ≡ +1 の両辺を a 倍した aw ≡ a は、何かの平方で表される」。だって平方根を持つ a を考えてるんだもん、x2 ≡ a を満たす x が存在するはず…
§71. w ≡ 1 になる最終ステップにおいて: aw ≡ a ≡ x2 を満たす x が求める平方根であり、Tonelli のアルゴリズムでは、
aw ≡ a(AZC) ≡ a(aqZC) ≡ aq+1ZC ≡ a
を開くことで、この x を決定する。明示的に書けば:
x ≡ ±a(q+1)/2ZC/2
q+1 と C は偶数なので、やれば機械的にできる計算だけど、ゴチャゴチャしてちょっと面倒くさい。これを単に、こう書いたらどうだろう？
aw ≡ v2 『か』
最終ステップでは w ≡ 1 なので、上の式は a ≡ v2 と同じことで、この v は求める x に他ならない。
なぜ、わざわざ文字を変えて x ではなく v としたのか？ まぁ別に変数名 x を使い回してもいいんだけど、要は最終ステップ（w ≡ 1 になる）以外でも『か』の関係を考えることができるってこと…。w ≡ AZC は最終ステップ以外では ≡ 1 とは限らないけど、とにかくその w と a の積を v の平方として表現する。どうやって？
第1ステップでは、初期値 C = 0 なので w ≡ AZC ≡ A ≡ aq であることを思い出そう。つまり:
aw ≡ aA ≡ a(aq) = aq+1 『き』
ここで q は奇数、q+1 は偶数なので…
v ≡ a(q+1)/2 『く』
…は『か』の関係を満たす。最初のステップはこれでOK。
〔注〕 なぜなら、そのとき v2 ≡ [a(q+1)/2]2 = aq+1 だが、『き』の通り、これは ≡ aw である。
初期状態（第1ステップ）における『か』の関係について、a は与えられた値だし（この値の平方根を求めたい）、w の初期値 A ≡ aq も分かっているし、v の初期値『く』も簡単に計算できるので、われわれは最初から『か』という鍵を持っている！
一般には w ≡ 1 ではないので、『か』が分かったからといって v は a の平方根ではないけど（v は aw の平方根である）、w の値は、補正が必要な各ステップにおいてだんだん変化して、最終ステップでは w ≡ 1 になるわけである。ということは、各ステップにおいて『か』の関係が維持されるように v の値をメモしておけば、w ≡ 1 となる最終ステップにおいては、自動的に a の平方根 v が確定する…。
初期値さえ正しければ、各ステップで『か』の関係を維持することは、簡単にできる。補正が不要のステップでは w の値は変化しないのだから、何もしなくても（つまり v の値をいじらずに）『か』の関係は維持される; 補正が必要なステップでは w の値は Z2 倍されたり、Z4 倍されたり、Z8 倍されたりするわけだが（前回参照）、例えば w が Z8 倍されるなら v を Z4 倍してやれば、『か』の関係が保たれる。その場合、仮に y = Z4, y2 = Z8 と置くとして:
aw ≡ v2
…の w を y2 倍して v を y 倍すると:
a(wy2) ≡ (vy)2 『け』
『け』は、既に成り立ってる aw ≡ v2 の両辺を y2 倍することに当たるので、もちろん成り立つ。つまり『け』において wy2 をあらためて w とし、vy をあらためて v とするなら、関係『か』が維持される（適切な y の値は状況によって変わるが、y の値が何であっても、この事実は不変）。しかも y2（上記の例では Z8 に当たる）の具体的な値を決定するとき、どうせ途中計算で y（上記の例では Z4 に当たる）の値も決定されるだろうから、w の値に補正をかけるついでに v の値に同様の補正をかけることは、安上がりな処理だ。
各ステップで『か』の関係を維持することにより、最終ステップに達したとき、求めたい平方根が v として（正確に言えば ±v だが）自動的に求まる！
§72. §70の具体例（mod 97 における a = 31 の平方根）について、この v という変数の値がどうなるか、具体的にトレースしてみよう。最初に§70の例を簡単にまとめておくと:
q = 3: 出発点は A ≡ aq = 313 ≡ 12 と Z ≡ zq = 53 ≡ 28
w ≡ A ≡ 12 からスタート 各ステップで wB ≡ +1 が維持されるように必要な補正を行う
第2ステップで w8 ≡ −1 になってしまうので w を Z2 倍して w ≡ +1 に戻す
第3・第4ステップでは補正不要
第5ステップでは w1 ≡ −1 になってしまうので w を Z16 倍して w ≡ +1 に戻す
これだけでは 11 ≡ +1 という分かり切った結論が出るだけだが、そこに v の計算を絡めると…
まず v の初期値は『く』の通り: v ≡ a(q+1)/2 = 312 ≡ 88
〔検〕 実際には必要な計算じゃないけど、aw ≡ v2 が成り立ってるってことを検証しておく:
左辺 31 × 12 ≡ 81 右辺 882 ≡ (−9)2 ≡ 81 (mod 97) ちゃんと一致している
さて、第2ステップで w を Z2 倍するのだから、それに対応して v を Z1 倍しよう:
新しい v ≡ 88 × 28 ≡ 39
〔検〕 新しい aw ≡ 31 × (−1) ≡ −31
新しい v2 ≡ 392 ≡ 66 ≡ −31 一致
第3・第4ステップでは w の値は変わらないので v の値も変わる必要ない。最後に、第5ステップで w を Z16 倍するんだから、それに対応して v を Z8 倍しよう。Z8 ≡ 22 という値は、どうせ Z16 の途中計算で既に求まってるはずなんで、それを再利用できるだろう:
新しい v ≡ 39 × 22 ≡ 82 ≡ −15
〔検〕 新しい aw ≡ 31 × 1 ≡ 31
新しい v2 ≡ (−15)2 ≡ 31
w ≡ 1 なので、この一致は a ≡ v2 を意味する
a = 31 の平方根は、最終ステップの v を使って ±v に当たる。つまり ±15。 w の計算のついでに、平行して v を計算しておくことで、§70と同じ答えが自然に求まった！
Shanks のアルゴリズムは、元祖 Tonelli のアルゴリズム（機械的な反復で極めて単純だが、少々無駄がある）と同じ内容を、より効率的に計算できるように工夫したもの。いきなり提示されると天下り的にも思えるが、Tonelli 版との比較で言えば:
Shanks の w とは Tonelli の AZC のことである （前々回）
最終ステップで w ≡ +1 のとき 11 ≡ 1 になるが、それだけでは役立たない （前回）
平行して aw ≡ v2 の v を計算すれば万事うまくいく （今回）
→ 愛と正義の美少女変数・セーラー v （意味不明）
今回は、この v の部分のコンセプトと、一つの具体例を提示した。次回以降、この内容をきちんと整理し、アルゴリズムの一般形を得たい（細部には、さらなる改善の余地もある）。「この計算で合ってる」という表面的な話だけでなく、内容的に何がどうなってるのか、群論や離散対数の観点からの考察も試みたい。
2023-07-16 Tonelli vs. Shanks（その4） y のスルー
前回までで分かったこととして、Shanks のアルゴリズムの原理はシンプル。まず v と w の初期値を求める。次に、各ステップにおいて、もし wB が ≡ −1 になったら v を y 倍、w を y2 倍する。この y とは初期値では zq のことだが、各ステップで y の値も次々と平方されていく。w が 1 になったときの v が求める平方根。
本質は、ただそれだけ！
この意味が一目で分かるように、今回は Shanks 版を表にまとめてみる。その上で、実装を完成させるため「スルーありの y 計算」に取り組みたい。
§73. 下の＜表2＞が＜表1＞の Tonelli 版と同じ意味であることは、前回までに順を追って説明した。要点として w は Tonelli の AZC の値。 wB ≡ −1 のとき、Tonelli 版では補正として Z2, Z4 などが掛け算される。この掛け算される値は下記 Shanks 版の y2 に他ならない（なぜなら y の初期値 zq とは、大文字の Z のこと）。 v とは何か――というと、Tonelli 版の「最終ステップの後で必要になる計算」を先回りして、ついでに（効率よく）実行している。 v の具体的な意味は「aw の平方根」。だから w ≡ 1 になれば v は「a の平方根」。
q とは p−1 を 2 で割れるだけ割ったとき出てくる奇数の商: p = q⋅2e + 1 と書ける（この場合 q = 3, e = 5）。そして B の値は 2e−1 から始まって次々半減する。 y の値はステップごとに 2 乗されるので、前のステップの y2 をあらためて y として（表の右斜め上の数をコピペ）、その平方を新しい y2 とすればいい。どの計算結果も p = 97 以上になったら p で割った余りで置き換えて適当に小さくしておく（大きいままでもいいのだが、桁が多いと後々計算が面倒）。 96 ≡ −1 (mod 97) なので、必要に応じて 96 の代わりに −1 を使う（その方がシンプルで分かりやすい）。
| ステップ | 各ステップの処理内容 | B | y | y2 |
|---|---|---|---|---|
| 1 | w ≡ aq と v ≡ a(q+1)/2 を求める | 16 | y ≡ zq を求める | |
| 2 | wB を計算する ・もし wB ≡ +1 なら何もしない ・もし wB ≡ −1 なら v を y 倍し w を y2 倍 w ≡ 1 になったときの v が a の平方根 | 8 | 28 | 8 |
| 3 | 4 | 8 | 64 | |
| 4 | 2 | 64 | 22 | |
| 5 | 1 | 22 | −1 | |
【例1】 §72の例を再び取り上げると: a = 31 の平方根を求める場合、まず次の値を決定する:
w ≡ 313 ≡ 12 ← Tonelli 版の A に当たる
v ≡ 31(3+1)/2 ≡ 312 ≡ 88
ステップ2において 128 ≡ −1 なので補正を発動:
現在の v ≡ 88 の 28 倍を新しい v とする: v ≡ 88 × 28 ≡ 39
現在の w ≡ 12 の 8 倍を新しい w とする: w ≡ 12 × 8 ≡ −1
これで現在の w は −1 なので、ステップ3の (−1)4 ≡ +1 とステップ4の (−1)2 ≡ +1 においては、何もする必要ない。
ステップ5において (−1)1 ≡ −1 なので補正を発動:
現在の v ≡ 39 の 22 倍を新しい v とする: v ≡ 39 × 22 ≡ 82
現在の w ≡ −1 の −1 倍を新しい w とする: w ≡ (−1) × (−1) ≡ 1
めでたく w ≡ 1 になったので、このときの v ≡ 82 が a の平方根。符号はどちらでもいいので +82 or −82 となる（言い換えれば −15 or +15）。
Tonelli バージョンと本質的に同じ内容とはいえ、アルゴリズムとして洗練されている。Tonelli は19世紀、Shanks は20世紀なので、多少の進化があって当然だろう。
【例2】 念のため、もう一つ例を挙げておく。同じ mod 97 において a = 43 の平方根を求めよう。これは Niven らの教科書の例題である――この文献は、教科書の中では一番親切で初心者向けのように思われるが、長々と「文章で」計算過程が説明されていて、少々見通しが悪い。＜表2＞を使って、整理してみよう…。まず次の値を決定:
w ≡ 433 ≡ 64
v ≡ 43(3+1) ≡ 432 ≡ 6
ステップ2では 648 ≡ +1 なので、何もしない。
ステップ3において 644 ≡ −1 なので補正を発動:
現在の v ≡ 6 の 8 倍を新しい v とする: v ≡ 6 × 8 ≡ 48
現在の w ≡ 64 の 64 倍を新しい w とする: w ≡ 64 × 64 ≡ 22
ステップ4において 222 ≡ −1 なので補正を発動:
現在の v ≡ 48 の 64 倍を新しい v とする: v ≡ 48 × 64 ≡ 65
現在の w ≡ 22 の 22 倍を新しい w とする: w ≡ 22 × 22 ≡ −1
ステップ5において (−1)1 ≡ −1 なので補正を発動:
現在の v ≡ 65 の 22 倍を新しい v とする: v ≡ 65 × 22 ≡ 72
現在の w ≡ −1 の (−1) 倍を新しい w とする: w ≡ (−1) × (−1) ≡ 1
めでたく w ≡ 1 になったので、このときの v が a の一つの平方根。すなわち ±72 言い換えれば ±25 が √43。
〔検算〕 252 = 625 = 582 + 43 = 97 × 6 + 73 ≡ 43 ちゃんと平方すると a になる。
§74. Shanks のアルゴリズムは、上の形式で問題なく動作する。けど、本来の完成形としては「ステップをスルーできる部分」をきっちり実装する必要がある。
スルーできるステップも含めて、全ステップを処理するなら、あるステップで一定の値 y（それを Y としよう）が使われた後、次のステップではその古い Y の平方が新しい y になる。もし2ステップまとめて処理できるなら、Y の平方の平方（つまり 4 乗）が新しい y に。3ステップまとめてなら Y の平方の平方の平方（つまり 8 乗）が新しい y に。4ステップまとめてなら 16 乗。等々。一般に、まとめて s ステップ進める場合:
新しい y ≡ Y の 2s 乗
ここで s = 1 の場合は、単に「すぐ次のステップ」での処理（古い y が平方される）。他方、「ステップ2」（y の初期値 Y に対応する）自体で補正が必要なら、1ステップも進まないで y の初期値 Y がそのまま使われる（ステップ2を出発点に、ステップ2自身が要補正）。これは s = 0 に当たる:
新しい y ≡ Y の 20 乗（つまり 1 乗） 要するに、新しい y = Y = 古い y
…ここまでは割と簡単。だけど、この s って数をどこから引っ張ってきましょう？
この点がほんのちょっとややこしい。
次のように B を 2m として表記してみる: ステップ1での m の値は e−1、この例では e = 5 なので m = 4。なぜならステップ1の B は 2e−1 だから。
〔注〕 ステップ1とは Euler の基準 a(p−1)/2 ≡ +1 に他ならない。 p = q⋅2e + 1 つまり
p−1 = q⋅2e
と仮定しているので:
(p−1)/2 = q⋅2e/2 = q⋅2e−1
ゆえに a(p−1)/2 = aq⋅2e−1 = (aq)2e−1 ≡ w2e−1
→ なぜなら初期状態では w ≡ aq
ステップ1では w の指数 2e−1 が B で、指数の指数 e−1 が m。
ただし y の初期値 28 はステップ2で使われるのだから、ステップ2における m = 3 つまり e−2 が「y の初期値に対応する m の初期値」。
| ステップ | m | B = 2m | y | y2 |
|---|---|---|---|---|
| 1 | 4 | 16 | y ≡ zq を求める | |
| 2 | 3 | 8 | 28 | 8 |
| 3 | 2 | 4 | 8 | 64 |
| 4 | 1 | 2 | 64 | 22 |
| 5 | 0 | 1 | 22 | −1 |
ステップ2以下で wB ≡ w2m が最初に ≡ −1 になる場所を探す――「補正が必要になる最初のステップ」だ。例えば、それがステップ2の m = 3 だとすると、それは m の初期値自身であり、初期値から「1ステップも進めない」ことになる。つまり s = 0 であり、この場合、単純に y の初期値（Y とする）をそのまま y として使うことになる:
補正に必要な y ≡ Y ≡ Y1 ≡ Y20
また例えば、ステップ3の m = 2 で最初の補正が必要になるとすると、これは（出発点のステップ2を基準に）「1ステップ進める」ことなので s = 1 であり:
補正に必要な y ≡ Y2 ≡ Y21
さらにステップ4の m = 1 で最初の補正が必要になるとすると、「2ステップ進める」ことなので s = 2 であり:
補正に必要な y ≡ Y4 ≡ Y22
一般に s ステップ進めるのなら:
補正に必要な y ≡ Y2s
補正が済んでしまえば、そのステップの結末では wB ≡ +1 が回復する。そのときの m を M として y を Y として、次の要補正ステップでの m の値を使うと、上記と同様 M から m まで s = M−m ステップ進めるのだから:
次の補正に必要な y ≡ Y2s
この考え方は、コンセプト的にはだいたい正しそうだ。もっとも、技術的には問題もある。第一、w2m ≡ −1 になる場所（そのような m）を探す…って言っても、もし w ≡ 1 になってたら
w20, w21, w22, w23, …
はいずれも 1 の累乗なので ≡ +1 であり、どこまで行っても ≡ −1 の場所が見つからない。捜索が終わる保証がない…。その他にも若干の短所がある。けどとりあえず、第一の欠点については、捜索を始める前に既に w ≡ 1 になっていないか確かめれば（もし w ≡ 1 ならそのときの v が答えなので、それ以上、何もする必要ない）、回避可能だろう。
従って、「スルー可能なステップはスルー」という最適化は、一応次のように実装できる。
（暫定的実装） ステップ2 の m = e−2（この数値を M とする）が出発点。既に w ≡ 1 になっていないか確認。なっていたら処理は終わっている。なっていなければ:
w20, w21, w22, w23, …
の中で ≡ −1 になるものを探す。それを w2m ≡ −1 としよう。すると s = M−m ステップ進めることができるので、現在の y を Y として:
新しい y ≡ Y2s
v にこの y を掛け算、w に y2 を掛け算すると、その結果の新しい w の値は w2m ≡ +1 を満たす（補正完了）。ここで w ≡ 1 になっていないか再確認。なってたら終了。なってなければ、現在の m をあらためて M として:
w20, w21, w22, w23, …
の中で ≡ −1 になるものを探す。それを w2m ≡ −1 としよう。すると s = M−m ステップ進めることができるので、現在の y を Y として:
新しい y ≡ Y2s
v にこの y を掛け算、w に y2 を掛け算すると…（以下同様）。
これを w ≡ 1 になるまで反復すればいい！
具体例で見てみましょう。
【例1（スルーあり）】 mod 97 における a = 31 の平方根（§73前半）を再考。初期値は v ≡ 88, w ≡ 12, y ≡ 28 だった。m の初期値は M = e−2 = 3 である。w ≡ 1 ではないので、捜索を開始しよう:
w20 = w1 = 12 （≡ −1 ではない）
w21 = w2 = 122 ≡ 47 （≡ −1 ではない）
w22 = w4 = 472 ≡ 75 （≡ −1 ではない）
w23 = w8 = 752 ≡ 96 （≡ −1 である！）
よって（初期値の後に発見される）最初の m は 3 で s = M−m = 3−3 = 0。従って y の値として、現在の y の値 28 を基準に、次の値を使わねばならない:
新しい y ≡ 2820 ≡ 281 ≡ 28
大騒ぎした割には、要するに「現在の y をそのまま使え」という結論が出た。そこで y ≡ 28 と y2 ≡ 8☆ を使って現在の v, w をそれぞれ補正すると:
新しい v ≡ 88 × 28 ≡ 39 新しい w ≡ 12 × 8 ≡ −1
ここまでは§73と同様に進んでいる。補正に使われた 28 と 8☆ も、＜表2＞あるいは＜表3＞のステップ2の欄と一致。新しい w ≡ −1 は、まだ ≡ 1 になっていない。現在の m = 3 を M として、次の捜索に進もう:
w20 = w1 = −1 （≡ −1 である！）
よって今度は m = 0 であり s = M−m = 3−0 = 3。うれしいことに一気に 3 ステップ進める。つまり現在の y の値 28 を基準に:
新しい y ≡ 2823 ≡ 288 ≡ 22
この y ≡ 22 と y2 ≡ −1☆ を使って、現在の v, w をそれぞれ補正すると:
新しい v ≡ 39 × 22 ≡ 82 新しい w ≡ (−1) × (−1) ≡ 1
w ≡ 1 になってので終了。このときの v すなわち 82 が a = 31 の平方根。§73の結果と一致！
2回目の補正で使われた y ≡ 22, y2 ≡ −1☆ は、当然ながら、＜表3＞のステップ5の欄と一致。
【例2（スルーあり）】 mod 97 における a = 43 の平方根（§73後半）を再考。初期値は v ≡ 6, w ≡ 64, y ≡ 28 で、m の初期値は M = e−2 = 3。 w ≡ 1 じゃないので、捜索開始:
w20 = w1 ≡ 64 （≡ −1 ではない）
w21 = w2 ≡ 642 ≡ 22 （≡ −1 ではない）
w22 = w4 ≡ 222 ≡ 96 （≡ −1 である！）
よって最初の m = 2 で s = M−m = 3−2 = 1。つまり 1 ステップ進めてステップ3へ。 y の値として、現在の y の値 28 を基準に、次の値を使わねばならない:
新しい y ≡ 2821 ≡ 282 ≡ 8
そこで y ≡ 8 と y2 ≡ 64★ を使って現在の v, w をそれぞれ補正すると:
新しい v ≡ 6 × 8 ≡ 48 新しい w ≡ 64 × 64 ≡ 22
ここまでは§73と同様に進んでいる。補正に使われた 8 と 64★ も、＜表3＞のステップ3の欄の値と一致。新しい w ≡ 22 は、まだ ≡ 1 ではない。現在の m = 2 を M として、次の捜索へ:
w20 = w1 ≡ 22 （≡ −1 ではない）
w21 = w2 ≡ 222 ≡ 96 （≡ −1 である！）
よって今度は m = 1 であり s = M−m = 2−1 = 1。つまり 1 ステップだけ進める。現在の y の値 8 を基準に:
新しい y ≡ 821 ≡ 82 ≡ 64
この y ≡ 64 と y2 ≡ 22★ を使って、現在の v, w をそれぞれ補正すると:
新しい v ≡ 48 × 64 ≡ 65 新しい w ≡ 22 × 22 ≡ −1
まだ w ≡ 1 になってないので、現在の m = 1 を M として、次の捜索へ:
w20 = w1 = −1 （≡ −1 である！）
よって今度は m = 0 であり s = M−m = 1−0 = 1。従って現在の y の値 64 を基準に:
新しい y ≡ 6421 ≡ 642 ≡ 22
この y ≡ 22 と y2 ≡ −1★ を使って、現在の v, w をそれぞれ補正すると:
新しい v ≡ 65 × 22 ≡ 72 新しい w ≡ (−1) × (−1) ≡ 1
a の平方根 72 が得られた。§73の結果と一致！
補正に使われる y, y2 の値★は、順に、＜表3＞のステップ3・ステップ4・ステップ5の欄と一致している。
上記「スルーありの y の値の更新法」は、いろいろな文献にある標準の方法と微妙に違うが、ある意味むしろ効率が良く、教科書のやり方より理解しやすい。ただこの実装だと m の初期値が e−2 なので、もし e が 2 未満だと、指数が負になってしまい具合が悪い。実用上は e が 2 未満なら、そもそもこんなアルゴリズム、必要ないんだけどね…
おとなしく教科書通りの形を検討しとくのも後々役立つので、次回は文献にある標準の方法（今回のと実質同じだが）を紹介したい。
2023-08-05 Tonelli vs. Shanks（その5） スルーの実装
p = 97 のとき、mod p での a = 31 の平方根の決定は、次の合同式を軸としていた（§72）。
aq(Z2)(Z16) ≡ 1 『さ』
ここで q = 3 は p−1 を 2 で割れるだけ割ったとき出てくる「奇数部分」。 Z ≡ zq は、任意に選択した非剰余 z の q 乗だが、具体的に z ≡ 5, Z ≡ 28 を使うことができる。『さ』の両辺を a 倍して、両辺の平方根を考える:
aq+1(Z2)(Z16) ≡ a
つまり a(q+1)/2(Z1)(Z8) ≡ √a 『し』
『さ』の左辺は w = aq から始めて、w の値をまず Z2 倍し、次に Z16 倍したものに他ならない。この処理のついでに『し』も計算できる。v = a(q+1)/2 から始めて、v の値をまず Z1 倍し、次に Z8 倍すれば、最終的な v の値（つまり『し』の左辺）は a の平方根。
w が Z2 倍されるとき v は Z1 倍される。
w が Z16 倍されるとき v は Z8 倍される。
一般に w の値が Z2n 倍されるとき、v の値は Z2n/2 = Z2n−1 倍される（前者の倍率から見て、後者の倍率では Z の指数が半分: 2 を n 個掛け合わせた指数 2n の半分は、2 を 1 個少なく掛け合わせた 2n−1）。言い換えると、v が Y ≡ Z2n−1 倍されるなら、w は Y2 倍される。というのも…
Y ≡ Z2n−1 なら Y2 ≡ (Z2n−1)2 ≡ Z(2n−1)⋅2
この右端の指数 (2n−1)⋅2 は、2 を n−1 個掛け、さらに 2 を 1 個掛けるのだから、トータルでは 2 を n 個掛けることに当たり:
Y ≡ Z2n−1 なら Y2 ≡ Z(2n−1)⋅2 ≡ Z2n
v, w の補正の進め方 v が Y 倍されるとき w は Y2 倍される
Y も Y2 も Z2□ の形を持つが、前者の方が □ に当たる数が 1 小さい
§75. §73【例2】でいえば a = 43 の平方根について:
aq(Z2)(Z4)(Z8) ≡ 1 『ささ』
a(q−1)/2(Z1)(Z2)(Z4) ≡ √a 『しし』
『さ』にせよ、『ささ』にせよ、要は w ≡ aq から始めて Z2, Z4, Z8, Z16 のうち必要な因子を順に掛け算。『し』にせよ、『しし』にせよ、v ≡ a(q+1)/2 から始めて、上記の因子に対応してそれぞれ Z1, Z2, Z4, Z8 のどれかを順に掛け算。
Z2, Z4, Z8 などの因子（補正係数）がケースバイケースで必要になるけど、これらの値は Z ≡ zq から始めて次々に平方するだけなので、単純計算で求まる。比較的面倒なのは、具体的な各ケースにおいて Z2, Z4, Z8 などのうち「どの因子が必要か」（掛け算しなければならないか）。直接的には、現在の w の値を 8 乗、あるいは 4 乗、あるいは 2 乗…などしたときに ≡ 1 になれば何もする必要がないが ≡ −1 になれば適切な因子を掛け算する「補正」が必要なのだった。
具体的に、ベタで書くと…
素数 p = q⋅2e + 1 について（q: 奇数）、mod p での平方剰余 a の平方根を求めたいとき、初期値 w ≡ aq, v ≡ a(q+1)/2, y ≡ zq ≡ Z から始めて…
- w2e−2 が ≡ 1 なら何もせず ≡ −1 なら v を y 倍、w を y2 倍。
- 上記 y2 の値をあらためて y とする。w2e−3 が ≡ 1 なら何もせず ≡ −1 なら補正: v を y 倍、w を y2 倍。
- 上記 y2 の値をあらためて y とする。w2e−4 が ≡ 1 なら何もせず ≡ −1 なら補正: v を y 倍、w を y2 倍。
- 上記 y2 の値をあらためて y とする。以下同様。（y の値は各ステップで平方される。）
これを w ≡ 1 になるまで続ければ、そのときの v の値が求める平方根の一つ。 w の指数 2e−2, 2e−3, 2e−4, … は、だんだん減っていき、しかも
w2e−j ただし j = 2, 3, 4, …
の値は、各 j に対して最初から ≡ 1 であるか、さもなければ ≡ −1。後者の場合、補正後には ≡ 1。従って、最悪でも j = e になったとき、
w2e−e = w20 = w1 = w
…が（もともと、あるいは補正後において） ≡ 1 になる。これは w1 = w ≡ 1 を意味するので、「w ≡ 1 になるまで続ける」というループは、必ずいつかは完了してストップする！
このベタなやり方で、理論的には何も問題ないけど、実装上 w2e−j ≡ 1 なら何もしない（補正は必要ない）のだから、そのようなステップはスルーして、補正が必要なステップだけを扱った方が効率的かも…
この「不要ステップはスルー」の実装法は、ざっと3通り考えられる――どれでも実質同じ意味。第一の方法は、上記「ベタ」に基づく（前回紹介したもの）。2番目の方法は、その亜種。「どのステップで処理が必要になるか？」という目印になるステップ番号について、1 ずつずらしただけ（全部 1 ずつずれるので、ステップの間隔は元のまま）。3番目の方法は、一般の文献や教科書に出てるやり方で、2番目の方法とほとんど同じ。
§76. スルーの実装その1。
ステップ2（ベタな方法で w2e−2 が ≡ 1 か −1 か判定する部分）が出発点: 2 の指数の初期値を m とすると m = e−2。
w20 (=w), w21 (=w2), w22 (=w4), …
を次々と計算して w2j ≡ −1 を満たす j を見つける。そのような j は、もしかすると m 自身かもしれない――その場合、ステップ2で早くも補正が必要になるのだから、補正係数 y は y の初期値 Z のまま（つまり v が y 倍され、w が y2 倍される）。でも、もしかすると j は m−1 という可能性もある――その場合、ステップ2はスルーでき、ステップ3（w2e−3 ≡ −1）で最初の補正が必要になるのだから、v の補正係数は、y の初期値から見て y2 に変化している（別の言い方をすると: y2 をあらためて y とした上で、その新しい y を使って v が y 倍され、w が y2 倍される）。
あるいはまた、もしかすると j = m−2 という可能性もある――その場合、ステップ2も3もスルーでき、ステップ4（w2e−4 ≡ −1）で最初の補正が必要になる（まとめて3ステップ処理できる）。v の補正係数は、y の初期値から見て y4 に。
一般に j = m−s なら、v の補正係数は、y の初期値から見て y2s = y2m−j。
いずれにせよ最初の補正が終わったら、現在の j をあらためて m として、再び
w20 (=w), w21 (=w2), w22 (=w4), …
を次々と計算し w2j ≡ −1 を満たす j を見つける。進め方は上と同じで、古い y の m−j 乗を新しい y とし、その y を使って v を y 倍 w を y2 倍。
この処理は、必ずいつかは w ≡ 1 になって終了する。なぜなら、もし w ≡ 1 でないなら、次々に平方される数…
w20 (=w), w21 (=w2), w22 (=w4), … 『す』
…が全部 ≡ 1 ってことは、あり得ない（特に左端の w20 = w は ≡ 1 じゃない）。けれど、この数たちの右端には ≡ 1 になるものがある。というのも、出発点となる w2m = w2e−2 についていえば、そのまた平方
(w2e−2)2 = w2e−1 ≡ (aq)2e−1 = aq⋅2e−1 = a(p−1)/2 【注1】
が ≡ 1 であることは、Euler の基準によって保証されている。他方、出発点以外のステップでは w2j ≡ −1 が発生し、それを補正した後では w2j ≡ 1 が成立。『す』の右の方には ≡ 1 になる数があり、一度そうなれば、そのまた平方も ≡ 1 なので、『す』の数たちは、ある場所から右は全部 ≡ 1。ところが『す』の左端には ≡ 1 でない数もあるので、『す』のどこかには「そこから右は全部 ≡ 1 だけど、その手前は ≡ 1 ではない」という境目の数がある。この最初に ≡ 1 になる手前の ≡ 1 ではない数（それを X とする）は、「その右の数（つまり X2）が ≡ 1 になるものの、X 自身は ≡ 1 ではない」という性質を持つ。つまり X ≡ −1 でなければならない。よって、w ≡ 1 になっていない限り、『す』のどこかに ≡ −1 の数があり、従って w2j ≡ −1 を満たす j が、必ずどこかにある！
しかも、出発点となる w2m は、それ自体が ≡ −1 かもしれないが、補正を済ませてしまえば w2m ≡ 1 になるので、それ以降 w2j ≡ −1 を満たす j は、m より小さい。それに対して補正を行い w2m ≡ 1 を確立してから、再び w2j ≡ −1 を満たす j を捜索すれば、その j はやはり m より小さい。従って j の値（補正後はその j をあらためて m とする）は、補正をするたびに、だんだん小さくなる。他方 j = 0 なら補正後の w20 ≡ 1 は w ≡ 1 と同じ意味なので、j の減少はどこかで止まる（w ≡ 1 になればそこで終了だが、j はだんだん小さくなり、j = 0 のときは w ≡ 1 になるしかないので、どっちにしても、この処理は終了する）。
【注1】 この式の右端の等号の理由は、次の通り。 p = q⋅2e+1 つまり p−1 = q⋅2e と仮定しているので p−1 乗と q⋅2e 乗は同じ意味:
aq⋅2e = ap−1
その両辺の指数を 2 で割ると（2e/2 = 2e−1 なので）:
aq⋅2e−1 = a(p−1)/2
仮定により a は平方剰余なので、この最後の右辺は Euler の基準により ≡ +1。
§77. スルーの実装その2。「実装その1」で具合が悪いのは、補正の目印になる m の初期値が e−2 という点。例えば
p = 11 = 5⋅21 + 1 ⇒ q = 5, e = 1
のように、e が 2 より小さいと、指数の初期値 e−2 が負になってしまう。実用上そのケースは別扱いすればいいとはいえ、理論上の整理としては、やはり全部のケースを統一的に扱えるアルゴリズムの方がいい。そこで指数の初期値を 1 大きくして、m = e−2 の代わりに m = e−1 という実装も考えられる。目印の数を全部 1 ずつ大きくすると、結果的に y の設定は「実装その1」と全く同じに…。目印を 1 増やすためには、w2j ≡ −1 を満たす j を捜索する代わりに、それより 1 大きい指数 k（それは j+1 に等しい）を捜索する必要がある。具体的に:
w20, w21, w22, w23, …
の中で w2j ≡ −1 になる場所を探す「実装その1」の代わりに、同じ数たちの中で w2k = 1 になる最小の k を探せばいい。左端から見て、初めて w2k = 1 になる場所は、w2j ≡ −1 になる場所の右隣なので、この実装の k は「その1」の j よりちょうど 1 大きい数。
【例】 mod 97 における a ≡ 31 の平方根。97 = 3⋅25 + 1 なので q = 3, e = 5。 v, w の初期値は、それぞれ v ≡ a(q+1)/2 ≡ 88, w ≡ aq ≡ 12。非剰余 z = 5 を選ぶと y の初期値は 5q ≡ 28。
実装その1では、2 の指数の初期値を m = e−2 = 3 として、w2j ≡ −1 を満たす j を探すと j = 3。新しい y は、古い y を使って y2m−j = y1 = y ≡ 28 と表される（古い y の値と同じ）。それを使って v, w を補正すると:
新しい v ≡ 88 × 28 ≡ 39 新しい w ≡ 12 × 282 ≡ −1
現在の j = 3 をあらためて m とする。新しい w ≡ −1 を使って w2j ≡ −1 を満たす j を探すと j = 0。今度の新しい y は現在の y の値を使って y23−0 = 288 ≡ 22。それを使って v, w を補正すると:
新しい v ≡ 39 × 22 ≡ 82 新しい w ≡ −1 × 222 ≡ 1
w ≡ 1 になったので、現在の v ≡ 82 は a の平方根。
実装その2では、2 の指数の初期値を m = e−1 = 4 として、w2k ≡ 1 を満たす最小の自然数 k を探すと k = 4。新しい y は、古い y を使って y2m−k = y1 = y ≡ 28 と表される（古い y の値と同じ）。 v, w の補正は、実装その1と同じ。現在の k = 4 をあらためて m とする。新しい w ≡ −1 を使って w2k ≡ 1 を満たす最小の自然数 k を探すと k = 1。今度の新しい y は現在の y の値を使って y24−1 = 288 ≡ 22。以下同じ。
「その1」でも「その2」でも、基準となる 2 の指数 m と、新しい 2 の指数（「その1」では j、「その2」では k）の差を使って、次々と y の値を更新していく。「その2」において、これらの指数は「その1」における指数より全部 1 ずつ大きいので、「差」という意味ではどちらでも同じで、y の値の設定に違いはない（結果として v, w の値の更新も全く同じ）。
§78. スルーの実装その3。文献に掲載されているスタンダードな方法。「実装その2」とほとんど同じだが、2 の指数の初期値を m = e−1 ではなく、さらに 1 大きく r = e とする。この場合も w2k ≡ 1 を満たす最小の自然数 k を探し、古い y を使って y2r−k という値を考える。r の値は、最初だけ「その2」と比べて 1 大きい（「その2」では e−1 だが、「その3」では e である）。その結果、最初に使われる y2r−k という値は、「その2」における y2m−k（この値を Y とする）と比べ、y の指数が 2 倍になっていて、Y2 に当たる:
Y = y2m−k とすると、初期値において r は m より 1 大きいので…
y2r−k = y2m+1−k = y2m−k+1 = y(2m−k)⋅2 = (y2m−k)2 = Y2
v の補正係数ではなく、その平方（すなわち w の補正係数）が得られる。初期値以外では「その2」の k と全く同じ値（w2k ≡ 1 を満たす最小の自然数）をあらためて r と置く: つまり 2r−k という値は、最初の補正以外では「その2」と同一。結局、「その2」と比べ、最初だけ補正係数が1回余分に平方されていて、それ以降は「その2」と同様に（同じ指数乗で）補正係数が更新されていく。だから、このようにして次々と得られる補正係数たち――更新のたびに古い y を使って y2r−k で表される――は、「その2」の各ステップでの y2 に等しい――言い換えれば、w の補正係数に当たる。
各ステップでの v の補正係数については、y の指数が半分になればいいので、y2r−k−1 を使えばいい:
v の補正係数を Y = y2r−k−1 とすると
w の補正係数は Y2 = (y2r−k−1)2 = y(2r−k−1)⋅2 = y2r−k
この Y, Y2 はそれぞれ「その2」の y, y2 に等しい
「実装その3」は、文字の使い方が少し違うだけ。各ステップにおける補正（の結果としての v, w の値の変化）は、「その1」「その2」と完全に一致する。
