メモ（数論15）: mod p の平方根

チラ裏　＞　数論 i　＞　メモ15

「チラ裏」は、きちんとまとまった記事ではなく、断片的なメモです。誤字脱字・間違いがあるかもしれません。

• 2023-04-20　「しげちー」の1円玉集め　圧死の危険…ｗ
• 2023-04-21　チョコレートの世界の平方根
• 2023-04-23　バニラの世界の平方根　習うより慣れろ
• 2023-04-26　春の世界の平方根　忘れられたイタリア人
• 2023-05-02　Tonelli のアルゴリズム　源流・1891年バージョン
• 2023-05-05　Atkin の平方根アルゴリズム　mod 8k+5
• 2023-05-10　Atkin の平方根アルゴリズム（その2）　幾何学的イメージ
• 2023-05-13　Müller のアルゴリズム　mod 16k+9
• 2023-05-17　Müller のアルゴリズムの簡単化に成功（速報）
• 2023-05-23　Müller のアルゴリズムの簡単化　2023年版
• 2023-05-27　Müller のアルゴリズムと無料平方根　25%のラッキー
• 2023-06-01　Dirichlet の美しい証明　x⁴ ≡ 2 が解を持つ条件
• 2023-06-24　相互法則の簡単な紹介　分かりやすく？！
• 2023-07-02　平方を 7 で割った余り・平方を p で割った余り　相互法則の例

2023-04-20　「しげちー」の1円玉集め　圧死の危険…ｗ

#数論　#べき和　#平方根 mod p

「しげちー」こと重清（しげきよ）クンには、「落ちてる1円玉を引き寄せ集める」という奇妙なスタンド能力がある。このスタンドには「ルール」があり、集められる枚数は、1日目には 1×1×1×1 = 1 枚、2日目には 2×2×2×2 = 16 枚、3日目には 3×3×3×3 = 81 枚…といった具合。漫画の設定とは違うが、それはそれ、これはこれ、ということで。

1円玉が16枚や81枚あっても、仕方ない？　表をよく見てもらおう。10日目の収穫は 10×10×10×10 つまり1万円。一見せこい1円玉集めだが、22日目には合計100万を突破するのだから、これは侮れないぞ。使おうとしても21枚以上は「受け取り拒否」されるかもしれないという小さな欠点と、そもそも「重くて持ち歩けない」という重大な欠点に、目をつぶればだがｗ…（1枚1グラムとして100万枚は1トン）

＜表1＞　「しげちー」の1円玉集め: その日の収穫と集めた合計枚数

	日	月	火	水	木	金	土
日付		1	2	3	4	5	6
収穫		1	16	81	256	625	1296
合計	(0)	1	17	98	354	979	2275
日付	7	8	9	10	11	12	13
収穫	2401	4096	6561	10000	14641	20736	28561
合計	4676	8772	15333	25333	39974	60710	89271
日付	14	15	16	17	18	19	20
収穫	38416	50625	65536	83521	104976	130321	160000
合計	127687	178312	243848	327369	432345	562666	722666

面白い現象: 1日目を月曜として日・水・土曜には合計枚数が必ず 7 の倍数になる。特に、水曜日には 49 の倍数になる。7 で割り切れるかどうかは分かりにくいが、電卓などで確かめてみよう！

〔例〕　10日目の 25333 = 7×7×11×47 は 7×7 = 49 の倍数。

このページは「べき和の秘密」の第2部に当たり、§31 から始まるけど、第1部とあんまり連続性はない。「べき和の秘密」の「1⁴+2⁴+3⁴+…+10⁴ のような和」その4・その5で、たまたま「しげちー現象」に気付いて、話がこんな方向に…。こーゆー何の役に立つのか分からんよーなことでも「自分で発見した数のパターン」には愛着が湧くものだ。せっかくだから、とことん追究するとしよう。途中でまた脱線するだろうけど…ｗ

§31.　n 日目の1円玉の合計枚数を Q(n) とする。
　　しげちー関数　Q(n) = 1⁴ + 2⁴ + 3⁴ + … + n⁴

「月・水・土に Q(n) が 7 の倍数になる」ってことを記号で書くと、こうなる:
　　n ≡ 0, 3, 6 (mod 7) ⇒ 7 | Q(n)
　　意味　日付 n を 7 で割った余りが 0, 3, 6 のどれかなら 7 は Q(n) を割り切る

1週間の日数 p が 7 とは限らないパラレルワールドを考えても、p が 7 以上の素数である限り、必ず週の最初（普通の世界で言えば「日曜」に当たる）、真ん中（水曜に当たる）、最後（土曜に当たる）には Q(n) が p で割り切れる。例えば p = 11 の場合:
　　n ≡ 0, 5, 10 (mod 11) ⇒ 11 | Q(n)
つまり日付が 11 の倍数のとき・11 の倍数より 5 大きいとき・11 の倍数より 10 大きいときに Q(n) は 11 の倍数になる。最後の「11 の倍数より 10 大きい」は、さらにあと 1 足せば 11 で割り切れるんだから「11 の倍数より 1 小さい」とも言い換えられる。

〔例〕　Q(10) = 25333 = 7×7×11×47 は 11 で割り切れる。

同様に p = 13 として:
　　n ≡ 0, 6, 12 (mod 13) ⇒ 13 | Q(n)

一般に p が 7 以上の素数なら:
　　n ≡ 0, (p − 1)/2, p − 1 (mod p) ⇒ p | Q(n)　　【ア】

言葉で説明すると…

基本しげちー理論　p を 7 以上の素数とする。しげちーが集めた1円玉の枚数は、p 日を周期に、各周期の 0 日目と (p − 1)/2 日目と p − 1 日目に p の倍数になる。この3種類の日を（p 日周期の）基本の日と名付けよう。

p の値によっては、基本の日以外の日にも Q(n) が p の倍数になることがある。例えば p = 17 として n ≡ 0, 8, 16 (mod 17) の日に Q(n) が 17 の倍数になるのは上記の通りだが、それ以外に n ≡ 2, 14 (mod 17) の日にも Q(n) は 17 の倍数。

〔例〕　Q(2) = 17, Q(14) = 127687 = 17 × 7511 は 17 の倍数。つまり 2 と 14 は、17日周期のスペシャルデー。次の周期（17日目～33日目）では 17 と 31 がスペシャルデーとなり Q(19) と Q(31) が 17 の倍数になる。

実は、しげちー関数は、次の値を持つ:

Q(n) = n(n + 1)(2n + 1)(3n² + 3n − 1)/30　　【イ】

見ての通り、この分子は 4 種類の数 n, (n + 1), (2n + 1), (3n² + 3n − 1) の積。もし n が p の倍数なら、もちろん【イ】の分子は p の倍数。 n が p の倍数より 1 小さければ、2番目の (n + 1) が p の倍数。同様に、n が p の倍数より 1 小さい数の半分なら、3番目の (2n + 1) が p の倍数。これらのどの場合でも、4個の数の積は、当然 p の倍数――つまり因子 p を持つ。分母は 2 × 3 × 5 なので、p が 7 以上の素数なら、分子に発生した因子 p が約分によって消えてしまうこともない。このことによって、基本しげちー理論は容易に説明がつく。

4番目の数 3n² + 3n − 1 が p の倍数になることがあれば、それ以外の日にも分子は p の倍数になり得る。言い換えると:
　　3n² + 3n − 1 ≡ 0 (mod p)　　【ウ】
に解があれば、p = 17 の場合のように、スペシャルデーが発生。

【ウ】の解は、形式的にはこーなる:
　　n ≡ ±√21 − 3/6 (mod p)
これが意味を持つためには、その世界に √21 という数が存在する必要がある。通常の数（実数）の世界には √21 = 4.58… という無理数が存在するが、mod p の世界には r² ≡ 21 (mod p) を満たす r が存在することも、しないこともある。もしそのような r が存在するなら、
　　n ≡ (±r − 3)/6 (mod p)　　【エ】
の場合がスペシャルデーで、そのときにも Q(n) は p の倍数。

p = 7 は非常に特殊なケースに当たり、次のように、基本の日とスペシャルデーが重なる。21 ≡ 0 (mod 7) の平方根は r ≡ 0 なので:
　　n ≡ (±r − 3)/6 ≡ −3/6 ≡ −3/−1 ≡ 3 (mod 7)
となり［mod 7 では 6 ≡ −1 である］、【ウ】の解 n ≡ 3 が週の真ん中にある基本の日 (p − 1)/2 ≡ (7 − 1)/2 と一致。つまり n ≡ 3 (mod 7) のとき、【イ】の分子の積の4因子のうち (2n + 1) と (3n² + 3n − 1) の両方が 7 の倍数。結果として、水曜の合計枚数 Q(n) は 7² = 49 の倍数！

§32.　11 以上の、好きな素数 p を考える（例えば p = 37）。果たして p 日周期のスペシャルデーは存在するだろうか？　細かい証明などは後回しにして、まずロジックの全体像を記す。

第三補充法則　p を 5 以上の素数とする。 p が12k±1型のとき、そしてそのときに限って、x² ≡ 3 (mod p) は解を持つ。

第三補充法則は、普通、平方剰余（じょうよ）の相互法則と Legendre（ルジャンドル）記号を使って証明されるが、実は直接証明も可能。この法則の意味は: もし p が 11 や 13 や 23 や 37 のように「12の倍数と1違い」の素数なら、mod p の世界には x² ≡ 3 を満たす x が存在する。例えば 5² ≡ 25 ≡ 3 (mod 11) なので、p = 11 の場合 x ≡ ±5 が解となる。他方、もし p が 17 や 19 や 29 や 31 のように「12の倍数と1違い」ではない素数なら、x² ≡ 3 (mod p) に解はない。別の言い方をすると、mod 11 の世界には √3 という数が存在するが、mod 17 の世界には √3 という数は存在しない。

「第七補充法則」　p を 11 以上の素数とする。 p ≡ ±1, ±3, ±9 (mod 28) のとき、そしてそのときに限って、x² ≡ 7 (mod p) は解を持つ。

p ≡ ±1, ±3, ±9 (mod 28) とは、p が「28k±1型」「28k±3型」「28k±9型」のどれかの素数であること。例えば、素数 37 は 28×1+9 なので「28k+9型」。

第三補充法則と同系統だが、少しややこしい。素数 p が「28の倍数と1違い・3違い・9違い」なら x² ≡ 7 (mod p) に解がある。例えば 9² ≡ 81 ≡ 7 (mod 37) なので、確かに mod 37 の世界には「平方すると 7 になる数」が存在している。この世界では √7 ≡ 9 というわけだ。

37 という数は「28の倍数と9違い」であると同時に「12の倍数と1違い」でもあるから、第三補充法則によって x² ≡ 3 (mod 37) も解を持つ。実際、37 × 6 = 222 は 37 の倍数なので 15² ≡ 225 ≡ 3 (mod 37)。

まとめると mod 37 の世界では √3 ≡ 15 で √7 ≡ 9。このことから、単純に考えて…
　　r ≡ √21 ≡ √3 × √7 ≡ 15 × 9 ≡ 135
…は r² ≡ 21 (mod 37) を満たすはず。前節【エ】によれば、この r（または −r）から 3 を引いて 6 で割った数が 3n² + 3n − 1 ≡ 0 の解を与え、すなわちスペシャルデーとなる。
　　r ≡ 135 とすると　(135 − 3)/6 ≡ 132/6 ≡ 22
　　r ≡ −135 とすると　(−135 − 3)/6 ≡ −138/6 ≡ −23 ≡ 14 (mod 37)

事実、Q(14) = 127687 = 37 × 3451 は 37 の倍数。そして 14 は「37 周期」の基本の日（0, 18, 36）に含まれないスペシャルデー。同様に、22日目は 37 の倍数になるスペシャルデー。
　　Q(22) = 1151403 = 37 × 31119
22日目は、しげちーが集めた1円玉の合計が100万円を突破する記念日でもある…。だが、重さが1トンを超えるので、その頃には圧死する危険が（笑）。ともあれ、この例は、次の事実を暗示している。

第三補充法則と第七補充法則を組み合わせて √3 と √7 が両方とも存在するような p を選択して、それら二つの平方根の積を r = √21 とすれば、【エ】によって p 日周期のスペシャルデーを決定できる。

この素直な考えは、間違ってはいない。p を 11 以上の素数として二つの法則を組み合わせると、p ≡ ±1, ±25, ±37 (mod 84) のとき、そしてそのときに限って mod p では √3 と √7 が両方存在し、それらの積は √21 だ。証明は簡単で、もし
　　x² ≡ 3 (mod p)
を満たす x と
　　y² ≡ 7 (mod p)
を満たす y が両方存在するのなら、左辺同士・右辺同士を掛け算して:
　　x²y² ≡ 21　つまり　(xy)² ≡ 21 (mod p)
――このように、積 xy は、平方すると確かに ≡ 21 になる。

p ≡ ±1, ±25, ±37 (mod 84) とは、p が「84k±1型」「84k±25型」「84k±37型」のどれかの素数であること。84 は 12 の倍数で 1, 25, 37 は12の倍数より 1 大きいので、これらのタイプの素数は第三補充法則の条件を満たす。84 は 28 の倍数でもあり ±25 ≡ ∓3, ±37 ≡ ±9 (mod 28) なので、第七補充法則の条件も満たす。

しかしこれは真実の半分に過ぎない！　√3 が存在せず √7 も存在しない世界においても、√21 は存在し、スペシャルデーが発生する。実例として p = 17 の場合を既に観察した: 17日周期で、2日目と14日目はスペシャルデーになり、それは mod 17 の世界に「平方すると 21 になる数」が存在していることと関係しているのだった。事実 2² ≡ 4 ≡ 21 (mod 17)。だが、17 は 12 の倍数と 1 違いではないので、第三補充法則から、この世界には「平方すると 3 になる数」つまり √3 は存在しないし、同様に、第七補充法則から √7 も存在しない。
　　√3 が実在せず √7 も実在しない世界で
　　√3 × √7 に当たる数が実在する
…というこの性質は、あたかも「虚数空間を通ってワープ！」のような、面白い感じがする。

スペシャルデーの特徴付け　p を 11 以上の素数とする。
　　(i)　素直なスペシャルデー　p ≡ ±1, ±25, ±37 (mod 84) のとき √3 と √7 は両方とも存在し、当然 √21 も存在して、スペシャルデーが生じる。
　　(ii)　ワープ・スペシャルデー　p ≡ ±5, ±17, ±41 (mod 84) のとき √3 と √7 は両方とも存在しないが、その場合も √21 は存在し、スペシャルデーが生じる。
　　(iii)　それ以外の場合、√3 と √7 のどちらか一方だけが存在し、√21 は存在しないため、スペシャルデーは発生しない。

(i) については既述。 (iii) の性質も簡単に証明できる: まず √3 だけが存在するケースについて…
　　x² ≡ 3 (mod p)　　‥‥①
を満たす x が存在するが、
　　y² ≡ 7 (mod p)　　‥‥②
を満たす y が存在しないとして、それにもかかわらず、もしも
　　z² ≡ 21 (mod p)　　‥‥③
を満たす z が存在するとしたら、矛盾が生じる。なぜなら③の左辺・右辺を①の左辺・右辺でそれぞれ割り算すると:
　　z²/x² ≡ 7　つまり　(z/x)² ≡ 7
ここで y ≡ z/x と置くと、「②を満たす y は存在しない」という前提に反する。√7 だけが存在するケースについても、①と②の右辺の値を入れ替えれば、同様。通常の数の世界で「平方数と非平方数の積が、平方数になるわけない」というのは直観的に明らかだが、(iii) の性質はそれと似ている。

通常の平方数は、素因数分解したとき、各因数の指数が偶数。例えば 9 = 3², 144 = 12² = (2²×3)² = 2⁴×3²。このタイプの数を2個、掛け合わせた積は、素因数分解するとやはり各因数の指数が偶数なので、再び平方数になる。一方、非平方数は、少なくとも一つの因数の指数が奇数乗（1乗も含む）。そのような数と平方数の積を素因数分解すると、やはり少なくとも一つの指数が奇数になるので、結果は非平方数。

(ii) はそれほど直観的ではなく、「①と②がどちらも解なしなら、③は解あり」という意味。通常の世界でも、「非平方数と非平方数の積が平方数」という現象は「起こり得ること」だが（例: 2 × 18 = 36）、「必ず起きること」ではない。mod p の世界では、この現象が必ず起きる。

原始根 g を使うと、平方剰余は g の偶数乗、平方非剰余は g の奇数乗。このことから: (i) 平方剰余同士の積は再び g の偶数乗で平方剰余。 (ii) 平方非剰余と平方非剰余の積は g の奇数乗と g の奇数乗の積なので、g の偶数乗になって、平方剰余。 (iii) 平方剰余と非剰余の積は g の奇数乗で、平方非剰余。

以上、大ざっぱな説明だが、一応スペシャルデーの発生条件（つまり √21 の存在条件）を明らかにした。でも √21 が存在するかどうか判定することと、具体的に √21 を求めることは、全然別の問題。後者については、検討の余地が多い。（続く）

2023-04-21　チョコレートの世界の平方根

#数論　#平方根 mod p

ある数 a について、mod p の世界で r² ≡ a を満たす r が存在するとき、r を a の平方根と呼ぶ（要するに、平方すると a になる数）。ここで p は 3 以上の素数。

〔例〕　p = 17 とする。5² = 25 ≡ 8 (mod 17) だから（意味: 25 も 8 も 17 で割った余りは同じ）、r² ≡ 8 (mod 17) を満たす r が、少なくとも一つは存在する: r ≡ 5。つまり mod 17 の世界では 5 は 8 の平方根。

平方根の計算は、p が「4 の倍数より 1 小さい素数」（例: 3, 7, 11, 19）なら簡単にできる。

次の四つの表現は、どれも本質的に同じ意味。　①「2乗して 47 で割ると 21 余る」ような整数を求める。　②「2乗すると 47 の倍数より 21 大きくなる」ような整数を求める。　③ mod 47 で √21 を求める。　④ x² ≡ 21 (mod 47) の解を求める。

「しげちー」のスペシャルデーを解明する鍵となるのは mod p での 21 の平方根だが、21 に限らず、一般的に平方根の求め方を考えてみたい。

§33.　普通の数（実数）の場合、正の数は 2 種類の平方根を持ち、0 は 1 種類の平方根（0 自身）だけを持ち、負の数は（実数の範囲では）平方根を持たない。一方、mod p の場合、0 が 1 種類の平方根を持つところは普通と同じだが、それ以外の数では平方根を持つ場合と持たない場合があり、どっちになるかについて（正の数・負の数のような）簡単な区別はできない。

〔例1〕　x² ≡ 3 (mod 13) には解 x ≡ ±4 があるが、x² ≡ 5 (mod 13) には解がない。mod 13 では、0 を別にすると 1, 3, 4, 9, 10, 12 が平方根を持ち、2, 5, 6, 7, 8, 11 が平方根を持たない。1, 4, 9 のように普通の意味での平方数のときは、もちろん平方根を持つが、それ以外の場合、平方根を持つか持たないかは、単純なパターンでは説明がつかない。

平方根の有無について、理論上重要なのは「Euler（オイラー）の基準」。0 以外の種類の数 a について、それを (p − 1)/2 乗すると ≡ ±1 になるのだが、この符号がプラスなら a は平方根を持ち、マイナスなら平方根を持たない（3 以上の素数 p を mod として）。

Euler の基準
　　(i)　a^(p−1)/2 ≡ 1 (mod p) なら a の平方根は存在する。
　　(ii)　a^(p−1)/2 ≡ −1 (mod p) なら a の平方根は存在しない。

〔補足〕　mod p において a の平方根が存在する場合・しない場合について、それぞれ a を「平方剰余（じょうよ）」「平方非剰余」ともいう（後者については、略して単に「非剰余」とも）。

〔例2〕　例1で見たように x² ≡ 3 (mod 13) には解がある。つまり mod 13 では 3 の平方根が存在する（言い換えれば 3 は平方剰余）。これは Euler の基準で a = 3, p = 13 のケース:
　　3^(13−1)/2 ≡ 3⁶ = 729 ≡ 1 (mod 13)
――確かに Euler の基準で (i) になる。一方、x² ≡ 5 (mod 13) には解がないのだった。Euler の基準によると:
　　5^(13−1)/2 ≡ 5⁶ = 15625 ≡ 12 ≡ −1 (mod 13)
――確かに (ii) になる。

Euler の基準を実際に計算する場合、素直に 3⁶ = 729 や 5⁶ = 15625 のような大きな指数計算をする必要はない。次の例のように「指数を小さく分けて」少しずつ計算した方が、はるかに楽。
　　3⁶ = (3³)² = 27² ≡ 1² = 1　　なぜなら 27 ≡ 1 (mod 13)
　　5⁶ = (5²)³ = 25³ ≡ (−1)³ = −1　　なぜなら 25 ≡ −1 (mod 13)

a の平方根が存在する場合、Euler の基準 (i) と下記〔補足〕から、こうなる。
　　a = 1 × a ≡ a^(p−1)/2 × a = a^(p+1)/2　　‥‥①
もし右辺の指数 (p+1)/2 が偶数なら、その指数を半分にした次の数を r とすると:
　　r ≡ a^(p+1)/4　従って
　　r² = r × r ≡ a^(p+1)/4 × a^(p+1)/4 = a^(p+1)/2　　‥‥②
②の右辺は①により ≡ a なので r² ≡ a となり mod p での a の平方根 r が求まる！

〔補足〕　指数法則により a^e × a^f = a^e+f が成り立つ。例えば:
　　a² × a³ = (a×a) × (a×a×a) = a×a×a×a×a = a⁵　　［← e = 2, f = 3, e+f = 5］
②の場合 e = f = (p+1)/4 なので:
　　e+f = (p+1)/4 + (p+1)/4 = 2(p+1)/4 = (p+1)/2
①の場合 a = a¹ だから e = (p−1)/2, f = 1 なので:
　　e+f = (p−1)/2 + 1 = (p−1)/2 + 2/2 = (p−1+2)/2 = (p+1)/2

p がチョコレート素数のときの平方根　a が mod p の世界で平方根を持つとする。もし (p+1)/2 が偶数なら:
　　a^(p+1)/4
は、a の一つの平方根。

(p+1)/2 が偶数なら (p+1)/2 = 2k と書ける（k: 整数）。両辺を 2 倍して p+1 = 4k つまり p = 4k−1。つまり (p+1)/2 が偶数という条件は、素数 p が「4k−1型」（4の倍数より 1 小さい）ということ（「4k+3型」ともいえる）。このタイプの素数をチョコレート素数という。

p がチョコレート素数なら、mod p での a の平方根 r の求め方は簡単で、単に a を (p+1)/4 乗したものを r とすればいい。r が a の平方根なら r² ≡ a だが、もちろん (−r)² ≡ r² なので、−r も a の平方根となる。r ≡ 0 の場合（つまり a ≡ 0 の場合）を別にすれば、r と −r は別の種類の数であり、a は（平方根を持つのなら）2種類の平方根を持つ。

なぜ r と −r は別の種類の数と言い切れるのか。もしも r ≡ −r なら両辺に r を足すと 2r ≡ 0、両辺を 2 で割ると r ≡ 0。つまり r ≡ −r と仮定すると r ≡ 0 になってしまう。「r ≡ 0 の場合を別にすれば」という話の前提と矛盾するので、その仮定は無理。（ただし合同式の両辺を 2 で割るためには、mod が奇数でなければならない。例外的ケースとして p = 2 なら 1 ≡ −1 (mod 2) が成り立つ。でも p = 2 なら (p+1)/2 が偶数という大前提が満たされていないので、この例外は上の話と関係ない。）

§34.　mod p での a の平方根を求めたいとする。 y = (p+1)/4 とすると a^y と −a^y がその答え―― p がチョコレート素数ということと、そもそも平方根が存在することが条件だが…。p がチョコレートでない場合（例えば p = 17）、y = 18/4 = 4.5 となって a^4.5 は（単純な整数乗としては）計算できないので、この方法は使えない。では p がチョコレート素数だとして、もし平方根が存在しない場合に無理やりこの式を使うと、何が起きるのだろうか？

mod 11 において √2 は存在しない。事実、Euler の基準によると:
　　2^(11−1)/2 = 2⁵ = 32 ≡ −1 (mod 11)
あるいは（もっと手っ取り早い判定法として）、3 以上の素数 p は、8k±1型のとき、そしてそのときに限って √2 (mod p) を持つ（第二補充法則）。存在しないのに、無理やり上の方法で a = 2 の平方根を求めようとすると…
　　y = (11+1)/4 = 3　　a^y = 2³ = 8
…となって、何やら ±8 という「答え」が出る（ちなみに mod 11 では 8 ≡ −3, −8 ≡ 3 なので ±3 と言っても同じ意味）。でも 2 の平方根はそもそも存在しないのだから、これは「間違った答え」のはず。検証のため、平方してみよう（平方根なら、平方すれば 2 になるはず）:
　　(±8)² = 64 = (11 × 6) − 2 ≡ −2 (mod 11)
　　あるいは同じことだが　(±3)² = 9 ≡ −2 (mod 11)
平方すると 2 ではなく、符号が逆の −2 になってしまう。言い換えると √−2 が求まった。mod 11 では −2 ≡ 9 なので、そして 9 は 3 の平方なので、√−2 ≡ √9 ≡ ±3 (mod 11) という関係自体は、当たり前といえば当たり前。

この例からも予想可能だが、a の平方根がないのに無理やり a^y を計算すると、a の平方根ではなく −a の平方根が求まってしまう。それもそのはずで、平方根が存在しないのなら Euler の基準から…
　　−a = −1 × a ≡ a^(p−1)/2 × a = a^(p+1)/2
…となり、その場合 a^y = a^(p+1)/4 の平方、すなわち a^(p+1)/2 は、≡ −a となる。

p がチョコレート素数のときの √a
　　y = (p+1)/4 を指数として、平方根を求めたい a を y 乗する: r ≡ a^y (mod p)
　　もし r² ≡ a なら ±r が a の平方根。もし r² ≡ −a なら a の平方根は存在しない。

√a が存在することが（Euler の基準などにより）事前に分かっている場合には、自動的に r² ≡ a になるので、理論上、そのことを確認する必要はない（実用上は、検算の意味で確かめてもいいが）。一方、√a が存在するかしないか分からない状態で、やみくもに r ≡ a^y を計算した場合、r² が ≡ a になるか必ず確認しなければならない。さもないと、約50%の確率で、実際には解 √a は存在しないのに √−a を解としてしまう羽目になる！

ついでに分かることとして、p がチョコレートの場合、a が平方剰余なら −a は非剰余。a が非剰余なら −a は平方剰余。±a のどちらか一方は（そして一方だけが）平方剰余になる（a ≡ 0 の場合は a = −a なので例外）。

§35.　上記の方法では p がチョコレート素数の場合しか扱えないが、素数の約半分はチョコレートなので、これで問題の50%は片付く。実例として、p 日周期の「しげちー」スペシャルデーを分析してみよう。p が「84k±1」「84k±25」「84k±37」型の素数なら、素直なスペシャルデーが発生し、p が「84k±5」「84k±17」「84k±41」型の素数なら、ワープ・スペシャルデーが発生する（§32）。これらの場合、√21 が存在すること自体は（補充法則から）分かっているのだが、具体的にどーやって √21 を求めるか。

84k+1型、84k+25型などは、4で割って1余るタイプ（バニラ素数）なので、上の計算法を使えない。一方:
　　84k−1型素数: 83, 167, 251, …
　　84k−25型素数: 59, 227, 311, …
　　84k−37型素数: 47, 131, 383, …
などはチョコレート素数。

具体例として p = 47 で試してみる。つまり mod 47 の世界で √21 を求めよう。 y = (47+1)/4 = 48/4 = 12 なので 21¹² (mod 47) がその答えだが、どういう数値になるのだろうか。
　　21² = 441 = 470 − 29 ≡ −29 ≡ 18 (mod 47)
　　だから　21⁴ = (21²)² ≡ 18² = 324 ≡ −5　　← 47×7 = 329
　　同様に　21⁸ = (21⁴)² ≡ (−5)² = 25 ≡ −22
　　結局　21¹² = 21⁸ × 21⁴ ≡ (−22) × (−5) ≡ 110 ≡ 16　　← 47×2 = 94
　　すなわち　√21 ≡ ±16 (mod 47)

〔検算〕　16² = 256 ≡ 21 (mod 47)　　← 47×5 = 235

実際のスペシャルデーを求めるには、3 を引いて 6 で割ればいいのだった。 16 − 3 = 13 は直接 6 では割れないが、この場合 13 ≡ 60 (mod 47) なので 13/6 ≡ 60/6 ≡ 10。だから10日目は、47日周期のスペシャルデー。事実、
　　Q(10) = 25333 = 7² × 11 × 47
は、47 の倍数。10日目の1円玉の合計枚数 25333 が49の倍数になること（水曜日の特殊ケース）、11の倍数になること（11日周期の基本日）は既に分かっていたが、47日周期のスペシャルデーでもある。Q(10) = 25333 の全部の素因数について、それぞれの正体がはっきりした。

47日周期のもう一つのスペシャルデーは (−16 − 3)/6 ≡ −19/6 ≡ −66/6 ≡ −11 ≡ 36 (mod 47)。［実は具体的に計算するまでもなく、一つのスペシャルデー s が分かれば、もう一つのスペシャルデーは (p−1) − s になる。］

チョコレート型のスペシャルデー　11 以上のチョコレート素数 p に対して Q(n) がスペシャルデーを持つ場合、それは次の日付。
　　(±21^y − 3)/6 (mod p)　ただし y = (p+1)/4
具体的には p ≡ −1, −25, −37; −5, −15, −41 (mod 84) の場合に、このタイプのスペシャルデーが生じる。

§36.　mod 47 の世界には √3 と √7 も存在するので（第三・第七補充法則）、それらを求めて掛け合わせることでも √21 が得られる。その方法だと平方根の計算を2回しなければいけないので、かえって手間が増えると予想されるが、物は試し、まぁやってみよう。

まず √3 を求める。3 を (47+1)/4 乗つまり 12 乗すればいいのだった。
　　3³ = 27 ≡ −20, 3⁶ ≡ (−20)² ≡ 400 ≡ −23　　← 47×9=423
　　3¹² ≡ (−23)² = 529 ≡ 12　← 47×11=517
　　検算　12² = 144 ≡ 3　← 47×3=141

次に √7 を求める。
　　7² = 49 ≡ 2, 7⁶ ≡ 2³ = 8, 7¹² ≡ 8² = 64 ≡ 17
　　検算　17² = 289 ≡ 7　　← 47×6=282

従って √21 ≡ 12 × 17 = 204 ≡ 16　［← 47×4=188］。符号が逆の −16 も平方根。後は前節と同じ。

ある数の累乗を何かで割った余り、つまり a^e (mod p) の形の計算（べき剰余）は筆算では少々面倒だが、やってること自体は単純。機械的にできる。

mod がバニラ素数の場合は、はるかに奥が深い。バニラ素数でも、8k+5型に限っては、上記と似た方法が成り立つ。（続く）

2023-04-23　バニラの世界の平方根　習うより慣れろ

#数論　#平方根 mod p

普通の数の世界で x² = 21 の解 x = ±√21 が ±4.58… であることは一応（？）簡単に計算できるが、mod p の世界で √21 に当たる数を決定することは、少々面倒くさい。だから面白い！

もっとも p が「4の倍数より1小さい素数」の場合（例: p = 3, 7, 11, 19）――言い換えると p−1 を 2 で割ると奇数になる場合――には、前回見たように、簡単なショートカットがある。

p−1 を 4 で割ると奇数になる場合（例: p = 5, 13, 29, 37）にも、似た方法が使える。今回は、これを考えてみたい。

p−1 を 4 で割った結果が偶数の場合（このパターンが一番難しい。例: p = 17, 41, 73）に取り組むための、準備ともなるだろう。

§37.　狙ってやったわけではないが、節の番号が 37。この素数 p = 37 は「1 を引いて 4 で割ると奇数になる」という性質を持つ。
　　(37 − 1)/4 = 36/4 = 9

奇数は偶数より 1 大きいので、k を整数として 2k+1 と書ける（例えば、奇数 9 は k = 4 の場合）。このことから、「1 を引いて 4 で割ると奇数になる」という性質は:
　　(p − 1)/4 = 2k + 1
両辺を 4 倍して:
　　p − 1 = 8k + 4　つまり　p = 8k + 5
要するに p が「8 の倍数より 5 大きい素数」言い換えれば「8 で割って 5 余る素数」の場合に当たる。このような素数の愛称は「秋素数」。奇数の素数は「バニラ」と「チョコ」に分かれ、バニラ素数は「春素数」と「秋素数」に分かれる――それぞれ 8 で割って 1 余る場合と 5 余る場合だ。前回、p がチョコレート素数の場合は全面解決したので、今回 p がバニラ素数の場合のうち「秋素数」を片付ければ、トータルで問題は75%解決したことになる。

mod 37 において x² ≡ 21 の解 x = ±√21 を求めたいとしよう。「そもそも解があるかないか」が問題だが、37 という周期では、しげちーの「素直なスペシャルデー」が存在する（§32）: すなわち 37 は 12 の倍数と 1 違いなので、この世界には √3 が存在し、さらに 28 の倍数とのずれが ±1, ±3, ±9 のどれかなので（具体的には +9）、この世界には √7 が存在する。必然的に、それらの積 √21 も存在する。従って Euler の基準から、次のような関係が成り立つ。
　　3¹⁸ ≡ 1, 7¹⁸ ≡ 1, 21¹⁸ ≡ 1 (mod 37)
指数の 18 は (p − 1)/2 に当たる。

p がチョコレート素数であれば、この指数 (p − 1)/2 が奇数なので、例えば 3¹⁸ ≡ 1 に当たる式の両辺を 3 倍して 3¹⁹ ≡ 3 として、そこから直ちに √3 を確定できるところだが（§33）、p がバニラ素数だと、指数 (p − 1)/2 が偶数なので、その手は通用しない。だがしかし！

指数をさらに半分にして (p − 1)/4 が奇数になるとするなら、次のように、そこを突破口として、チョコレート素数の場合と同様の処理が成り立つ。今回のテーマ8k+5型素数（愛称: 秋素数）は、まさにこのパターン。例えばの話、Euler の基準から
　　7¹⁸ ≡ 1 (mod 37)　　【カ】
は分かってるわけで、もし指数を半分にした…
　　7⁹ ≡ 1 (mod 37)　　【キ】
…が成り立つなら、その両辺に 7¹ を掛け算して…
　　7¹⁰ ≡ 7 (mod 37)　　【ク】
…そのことから直ちに次の結論に至る:
　　7⁵ ≡ √7 (mod 37)　　【ケ】

7¹⁰ が ≡ 7 なら（【ク】）、7⁵ は「平方すると ≡ 7 になる数」。事実 (7⁵)² = 7¹⁰。より正確に言うと、平方される数の符号は ± どっちでもいい:
　　7⁵ ≡ ±√7 (mod 37)

【ク】の段階で「偶数乗」になるので、容易に両辺の平方根を考えることができるのだ。さて、上の例はたまたま正しいが、一般には【カ】から【キ】に進める保証がない。というのは、ある数 a について a¹⁸ ≡ 1 つまり (a⁹)² ≡ 1 のとき、X ≡ a⁹ と置くと X² ≡ 1。従って X は平方すると 1 になる数: X ≡ ±1 つまり a⁹ ≡ ±1 となり、【キ】に当たる式の右辺が +1 or −1 のどっちになるか、分からない。【キ】に当たる部分の符号がマイナスになってしまうと、次の例のような状況になる。
　　3¹⁸ ≡ 1 (mod 37)　　【ガ】
　　3⁹ ≡ −1 (mod 37)　　【ギ】
その両辺に 3¹ を掛け算して…
　　3¹⁰ ≡ −3 (mod 37)　　【グ】
従って:
　　3⁵ ≡ √−3 (mod 37)　　【ゲ】
平方すると 3 になる数を求めたいのに、符号が逆転して「平方すると −3 になる数」が求まってしまった。【ケ】【ゲ】どちらの流れになるかはケースバイケース。後者になっても補正する方法があるので、心配は要らない。補正法を研究する前に、上の実例を具体的に計算してみよう。

まず【ケ】の 7⁵ を求める。
　　7² = 49 ≡ 12 (mod 37)
　　7⁴ = 12² = 144 = 148 − 4 ≡ −4 (mod 37)　　← 37×4=148
　　7⁵ = 7⁴ × 7¹ = (−4) × 7 = −28 ≡ 9 (mod 37)　符号は逆でもいい（±9）

確かに (±9)² = 81 = 74 + 7 ≡ 7 (mod 37) なので［37×2=74］、7 の平方根は ±9。

次に【ゲ】の 3⁵ を求める。
　　3² = 9
　　3⁴ = 9² = 81 = 74 + 7 ≡ 7 (mod 37)　　← 37×2=74
　　3⁵ = 3⁴ × 3¹ = 7 × 3 = 21 ≡ −16 (mod 37)

(−16)² = 256 = 259 − 3 ≡ −3 (mod 37) なので［37×7=259］、−3 の平方根が求まってしまった。−3 ≡ 34 なので、34 の平方根と言ってもいい。求めたいのは「平方すると 3 になる数」であって「平方すると −3 になる数」ではない。両者は似ているが、別の存在。「−16 で符号が逆なら +16 にすりゃいいじゃん？」と言いたいところだが（それで済めば楽なのだが）、±16 のどっちも平方すれば 256 ってことに変わりはないので、それでは問題は解決しない。

§38.　a の平方根を求めたいとき r² ≡ a になる r が一つでも求まれば、答えは ±r で、この部分の符号はどっちでもいい。しかし r² ≡ a を求めたいのに r² ≡ −a になる r が求まってしまった場合、±r のどちらを平方しても結果は −a であり、求めたいものと違う。この場合、√a が欲しいのに √−a が求まってしまったのだから、要するに √−1 倍してやればいい:
　　√−a × √−1 ≡ √a
問題は、この処理に必要な −1 の平方根をどーやって求めるか？　Euler の基準によれば、z が平方非剰余なら
　　z^(p−1)/2 ≡ −1 (mod p)　ゆえに　z^(p−1)/4 ≡ √−1 (mod p)
われわれの p は 8k+5 型なので p−1 は 4 で割り切れ、この計算を実行できる。でも、その原料となる平方非剰余 z をどーやって手に入れればいいのだろうか…。z は平方非剰余でありさえすれば何でもいいのだから、
　　z にいろんな数をランダムに当てはめ、Euler の基準を使って、見つかるまで探しまくる
　　つまり z^(p−1)/2 ≡ −1 になる z を試行錯誤で見つける
…というのも、一つの手ではある。でも、もっといい考えがある。p が8k+5型（秋素数）なら必ず 2 は平方非剰余（第二補充法則）。だから z = 2 を使うのが手っ取り早い！

8k+5型の素数 p を mod として、2^(p−1)/4 は −1 の平方根の一つ。

〔例〕　p = 37 のとき。この場合 6² = 36 ≡ −1 なので −1 の平方根が ±6 であることは最初から見え見えだが、まぁ説明のための例として、上の手順通りにやってみよう。 (37−1)/4 = 9 なので 2⁹ を求めればいい。
　　2² = 4, 2⁴ = 4² = 16, 2⁸ = 16² = 256 = −3 (mod 37)　　← 37×7=259
　　2⁹ = 2⁸ × 2¹ ≡ (−3) × 2 ≡ −6　符号は逆でもいい（±6）

従って √3 ≡ √−3 × √−1 ≡ 16 × 6 = 96 ≡ 22。符号は反対でもいいので ∓22 ≡ ±15。めでたく √−3 から √3 が求まった！

〔検算〕　(±15)² = 225 = 222 + 3 ≡ 3 (mod 37)　　← 37×6=222
あるいは同じことだが　(∓22)² = 484 = 481 + 3 ≡ 3 (mod 37)　　← 37×13=481
平方するとちゃんと ≡ 3 になってくれる。

この √3 ≡ 15 を、前節で求めてある √7 ≡ 9 と掛け算すれば:
　　√21 ≡ 15 × 9 = 135 = 148 − 13 ≡ −13 (mod 37)　　← 37×4=148
　　結果の符号は逆でもいい（±13）

〔検算〕　(±13)² = 169 = 185 − 16 ≡ −16 ≡ 21 (mod 37)　　← 37×5=185
確かに 13 は 21 の平方根だ！　（169 = 148 + 21 ≡ 21 とした方が少し楽だが、結論には変わりない。）

8k+5型素数 p を mod とする平方根の求め方　a の平方根 r が存在するなら a^(p−1)/4 × a ≡ +a or −a　　〔注1〕
　その左辺は = a^(p+3)/4 なので a^(p+3)/8 ≡ √a or √−a　　〔注2〕
　　(i)　前者の場合 a の平方根が既に求まっている。
　　(ii)　後者の場合 −a の平方根が得られる。それを √−1 倍すれば a の平方根。ここで √−1 ≡ 2^(p−1)/4
　いずれの場合も、最終的な平方根の符号は ± どちらもいい。

〔注1〕　Euler の基準により a^(p−1)/2 ≡ 1。ここで X ≡ a^(p−1)/4 と置くと X² ≡ a^(p−1)/2 ≡ 1 なので X ≡ +1 or −1 になる（§37参照）。

〔注2〕　a^(p−1)/4 × a = a^(p−1)/4 × a^4/4 = a^(p−1+4)/4 = a^(p+3)/4
その平方根を s とすると s² ≡ a^(p+3)/4。 s ≡ a^(p+3)/8 はこの条件を満たす:
　　s² = [a^(p+3)/8]² ≡ a^(p+3)/4　　←　指数 (p+3)/8 は割り切れて正の整数
s の前に付く符号はプラスでもマイナスでもいいが、s² が +a になるか −a になるかはコントロールできない。後者の場合、補正が必要。

(ii) の場合、√−1 ≡ 4^(p+3)/8 ÷ 2 が成り立つ。

なぜなら (ii) から 2⋅√−1 ≡ 2¹⋅2^(p−1)/4 = 2^(p+3)/4 = 2^2⋅(p+3)/8 = (2²)^(p+3)/8 = 4^(p+3)/8。両辺を 2 で割ると上の式。

この書き方を使うと (ii) の場合:
　　a^(p+3)/8 × √−1 ≡ a^(p+3)/8 × 4^(p+3)/8 ÷ 2
　　= (a × 4)^(p+3)/8 ÷ 2 = (4a)^(p+3)/8 ÷ 2

p = 8k+5 の形なので、指数 (p+3)/8 は (8k+5+3)/8 = (8k+8)/8 = k+1 に等しい。こう整理できる。

Pocklington の公式（1917年^†）　p を8k+5型素数、a を mod p の平方剰余とする。
　　(i)　a^(p−1)/4 ≡ +1　なら　√a ≡ a^k+1
　　(ii)　a^(p−1)/4 ≡ −1　なら　√a ≡ (4a)^k+1 / 2

〔例〕　p = 37 = 8⋅4 + 5 の場合 k = 4。 a = 7 について、【キ】で見たように 7^(37−1)/4 ≡ +1 なので √7 ≡ 7⁴⁺¹ ≡ 9（前記計算と一致）。 a = 3 について、【ギ】で見たように 3^(37−1)/4 ≡ −1 なので √3 ≡ (4⋅3)⁴⁺¹ / 2 ≡ 22（前記計算と一致）。平方根の符号は逆でもいい（直前の例でいえば −22 ≡ 15 でもいい）。

†　H. C. Pocklington: The Direct Solution of the Quadratic and Cubic Binomial Congruences with Prime Moduli, Proceedings of the Cambridge Philosophical Society, Vol. 19, pp. 57–59.
https://archive.org/details/proceedingsofcam1920191721camb (Cf. [6], p. 222)
改良版ともいえる Atkin のアルゴリズムがある。

〔付記〕　もし a の平方根が存在しないなら a^(p−1)/2 ≡ −1 なので、a^(p−1)/4 ≡ √−1。 −1 の平方根は2種類あるが、そのどちらになるかは a に応じて a^(p−1)/4 ≡ a^(8k+5−1)/4 ≡ a^2k+1 によって定まる。その値を ≡ i とすると ai ≡ a^2k+2 つまり a^k+1 ≡ √ai。例えば p = 13, k = 1 の場合: mod p の非剰余 a ≡ 5 に対して i ≡ a^2k+1 ≡ 5³ = 130 − 5 ≡ −5 は −1 の一つの平方根（その平方 25 は ≡ −1）。そのとき a^k+1 ≡ 5² = 25 ≡ −1 は ai ≡ −25 ≡ +1 の一つの平方根。一方、非剰余 a ≡ 7 に対して、a² ≡ −3 なので i ≡ a³ ≡ −21 ≡ +5 も −1 の一つの平方根。そのとき a^k+1 ≡ 7² = 49 ≡ −3 は ai ≡ 35 ≡ 9 の一つの平方根。このように、a の平方根が存在しないのに無理やり (i) の計算を行うと a の平方根ではなく ai の平方根が求まる。

§39.　説明のため √3 と √7 を別々に求めて掛け算したが、実際には √21 を直接求めた方が簡明。 a = 21, p = 37 として…
　　a^(p+3)/8 = 21^(37+3)/8 = 21⁵ ≡ (−16)⁵ (mod 37)
…を計算すると:
　　(−16)² = 256 ≡ −3 (mod 37)　　← 37×7=259
　　(−16)⁴ = (−3)² = 9
　　(−16)⁵ = (−16)⁴ × (−16)¹ = 9 × (−16) = −144 ≡ 4 (mod 37)　　← 37×4=148
得られた 4 は 4² = 16 ≡ −21 (mod 37) なので −21 の平方根。−1 の平方根 6 を使って補正すると:
　　√−21 × √−1 ≡ 4 × 6 = 24 ≡ −13 (mod 37)
従って ±13 は 21 の平方根。前節の結果と一致。

しげちーのスペシャルデーは、21 の平方根から 3 を引いて 6 で割ったもの。13 − 3 = 10 は直接 6 で割れないが、10 ≡ 47 ≡ 84 (mod 37)。 84 を 6 で割れば答えは 14。ちなみに、もう一つのスペシャルデーは (37 − 1) − 14 = 22 になる。結局、集めた1円玉の合計枚数は n ≡ 0, 18, 36 (mod 37) なら 37 の倍数で（基本の日）、n ≡ 14, 22 (mod 37) のときも 37 の倍数（スペシャルデー）。

基本の日の例:　Q(18) = 432345 = 37 × 11685, Q(36) = 12948594 = 37 × 349962, Q(37) = 14822755 = 37 × 400615
スペシャルデーの例:　Q(14) = 127687 = 37 × 3451, Q(22) = 1151403 = 37 × 31119

これで p が「8 で割って 3, 5, or 7 余る素数」の場合、mod p での平方根の計算の仕方が分かった！　つまり 8 で割って余りが 3 or 7 ならチョコレート素数だし（解決済み）、余りが 5 なら今回の手が使える。

残るは p を 8 で割って 1 余るとき。その場合、mod p での平方根の計算はもっと難しい…。（続く）

2023-04-26　春の世界の平方根　忘れられたイタリア人

#数論　#平方根 mod p

前々回は「p−1 を 2 で割ると奇数」の場合。前回は「p−1 を 4 で割ると奇数」の場合。次の自然な一歩は、「p−1 を 8 で割ると奇数」になる場合だろう。

§40.　p を 3 以上の素数とする。mod p において x² ≡ a に解があり x² ≡ z に解がないとすると、Euler の基準から:
　　a^(p−1)/2 ≡ +1　　【サ】
　　z^(p−1)/2 ≡ −1　　【シ】
もし (p − 1)/2 が奇数なら、【サ】の両辺を a 倍すると
　　a の偶数乗 ≡ +a
の形になるので、その両辺を 1⁄2 乗すれば、直ちに √a が求まる。一方、(p − 1)/2 が偶数だとしても、もし (p − 1)/4 が奇数なら、【サ】から、次のどちらかが成り立つ:
　　a^(p−1)/4 ≡ +1　　【ス】
　　a^(p−1)/4 ≡ −1　　【ズ】
もし【ス】なら、その両辺を a 倍すれば
　　a の偶数乗 ≡ +a
の形になるので、最初と同様に、直ちに √a が求まる。もし【ズ】なら、その左辺・右辺にそれぞれ【シ】の左辺・右辺を掛けると:
　　a^(p−1)/4 × z^(p−1)/2 ≡ +1　　【セ】
仮定によって (p−1)/4 は奇数、(p−1)/2 は偶数なので、上の式の両辺を a 倍すれば:
　　a の偶数乗 × z の偶数乗 ≡ +a　つまり　(az) の偶数乗 ≡ +a
この場合も、容易に √a が求まる。

以上が前回までの内容だが、(p − 1)/4 がまだ偶数になってしまう場合のために、アルゴリズム拡張したい。その場合でも【ス】あるいは【セ】自体は成り立つが、それらの左辺は a の偶数乗を含むため、両辺を a 倍すると a の奇数乗が発生してまい、a の平方根を求めることができない。それでもとにかく【ス】【セ】を一つにまとめて、こう書くことにしよう:
　　a^(p−1)/4 × z^S(p−1)/2 ≡ +1　　【ソ】
ここで S は、【ス】【ズ】の場合分けに応じてオンオフされるスイッチで、【ス】のケースでは S = 0 とし、【ズ】のケースでは S = 1 とする。前者の場合、S = 0 なので z^S(p−1)/2 = z⁰ = 1 という因子は無いのと同じことで、【ソ】は【ス】と同じ意味。後者の場合、S = 1 なので z^S(p−1)/2 = z^(p−1)/2 であり、【ソ】は【セ】と同じ意味。

仮定によって (p−1)/2 も (p−1)/4 も偶数なので、【ソ】の両辺の自然な平方根を考えることができる（a の平方根ではない）――【ソ】左辺の自然な平方根とは、各指数を半分にしたもの。【ソ】右辺 1 の平方根は、左辺の自然な平方根に応じて、自動的に ≡ +1 または ≡ −1 になる。つまり、次のどちらかが成り立つ（【ス】or【ズ】を導いた前回の発想と同様）:
　　a^(p−1)/8 × z^S(p−1)/4 ≡ +1　　【タ】
　　a^(p−1)/8 × z^S(p−1)/4 ≡ −1　　【ダ】

ここで (p − 1)/8 がようやく奇数になってくれたとすると、【タ】左辺はついに a の奇数乗を含むので、両辺を a 倍することで √a が求まるパターンに。【ダ】になってしまった場合は、上記【ズ】のときと同様、両辺に【シ】を掛けると:
　　a^(p−1)/8 × z^S(p−1)/4 × z^(p−1)/2 ≡ +1　　【チ】
これも a の奇数乗 × z の偶数乗に整理されるので、両辺を a 倍すれば (az) の偶数乗 ≡ a の形になり、問題なく √a を引き出せる。

(p − 1)/8 がしつこくまだ偶数の場合に備えて【タ】と【チ】を一つにまとめておくと:
　　a^(p−1)/8 × z^S(p−1)/4 × z^T(p−1)/2 ≡ +1
　　つまり　a^(p−1)/8 × z^{[S(p−1)/4 + T(p−1)/2]} ≡ +1　　【ツ】
ここでスイッチ T は、【タ】の場合には T = 0 でその項はないのと同じこと（そのとき【ツ】は【タ】自身と同じ意味）、【ダ】の場合には T = 1 で、そのとき【ツ】は【チ】と同じ意味。もしも (p−1)/8 が偶数なら、【ソ】のときと同様に【ツ】の自然な平方根を考えることができ、同様に進めると、いつかは a の指数が奇数になる（p−1 は有限の数なので、2 で有限回しか割れないので）。

§41.　具体例で確かめるため、美しさや効率にこだわらず、上記の（途中までの）内容を機械的に公式化しておく。

Tonelli の公式（とりあえずの動作確認用）　mod p で √a が存在する場合…
　　(i)　p−1 を 4 で割ると奇数になるなら、次式から √a を求めることができる。
　　a^(p−1)/4 × z^S(p−1)/2 ≡ +1　　〔注3〕
　　すなわち　a^(p+3)/8 × z^S(p−1)/4 ≡ √a　　【テ】
ここで z は mod p の任意の平方非剰余（p−1 を 4 で割ると奇数になるなら z = 2 を選択可能）。 a^(p−1)/4 ≡ −1 なら S = 1、 a^(p−1)/4 ≡ +1 なら S = 0（ないのと同じ）。
　　(ii)　p−1 を 8 で割ると奇数になるなら、(i) の S 設定と次式から √a を求めることができる。
　　a^(p−1)/8 × z^{[S(p−1)/4 + T(p−1)/2]} ≡ +1　　〔注4〕
　　すなわち　a^(p+7)/16 × z^{[S(p−1)/8 + T(p−1)/4]} ≡ √a　　【ト】
ここで a^(p−1)/8 × z^S(p−1)/4 ≡ −1 なら T = 1 で、そうでなければ T = 0（ないのと同じ）。

〔注3〕　(p−1)/4 が奇数なら (p−1)/4 + 1 = (p−1)/4 + 4/4 = (p+3)/4 は偶数。〔注3〕の式の両辺を a 倍すると:
　　a^(p+3)/4 × z^S(p−1)/2 ≡ a　［なぜなら a^(p−1)/4 × a¹ = a^{(p−1)/4 + 4/4}］
その両辺の自然な平方根を考えると「すなわち」の式になる。

〔注4〕　(p−1)/8 が奇数なら (p−1)/8 + 1 = (p−1)/8 + 8/8 = (p+7)/8 は偶数。〔注4〕の式の両辺を a 倍すると:
　　a^(p+7)/8 × z^{[S(p−1)/4 + T(p−1)/2]} ≡ a　［なぜなら a^(p−1)/8 × a = a^{(p−1)/8 + 8/8}］
その両辺の自然な平方根を考えると「すなわち」の式になる。

例えば mod を p = 41 とすると (p−1)/8 が奇数 5 なので (ii) を適用できる。次の4種類の数（平方剰余である）の平方根を求めてみる: a ≡ ±18, a ≡ ±20。第三補充法則から非剰余 z = 3 を選択でき、そのとき 3⁵ ≡ −3, 3¹⁰ ≡ 9。

〔補足〕　3⁴ ≡ 81 ≡ 41 × 2 − 1 ≡ −1, 3⁵ ≡ 3⁴ × 3 ≡ −1 × 3 ≡ −3, 3¹⁰ ≡ (3⁵)² ≡ (−3)² ≡ 9。ちなみにもう一度平方すると 3²⁰ ≡ (3¹⁰)² ≡ 9² ≡ 81 ≡ 41×2 − 1 ≡ −1。これは 3^(p−1)/2 なので Euler の基準からも 3 は非剰余。

a ≡ 18 と a ≡ 23 ≡ −18 は、a⁵ ≡ (±18)⁵ ≡ ±1 なので〔注5〕、どちらも a¹⁰ ≡ +1 を満たし S = 0。まず a ≡ 18 のとき a^(p−1)/8 ≡ 18^(41−1)/8 ≡ 18⁵ ≡ +1 なので T = 0。【ト】（左辺後半は z⁰ ≡ 1）を使い、左辺と右辺を逆順に書くと:
　　√18 ≡ 18^(41+7)/16 = 18³ ≡ 10
実際 10² ≡ 18。一方 a ≡ −18 のとき a^(p−1)/8 ≡ (−18)⁵ ≡ −1 なので T = 1。【ト】から:
　　√23 ≡ (−18)^(41+7)/16 × 3^(41−1)/4 ≡ −10 × 3¹⁰ ≡ −10 × 9 ≡ −8
実際 (−8)² ≡ 23。平方される数の符号は逆でもいい（8² も ≡ 23）。

次に a ≡ 20 と a ≡ 21 ≡ −20 は、a⁵ ≡ (±20)⁵ ≡ ±32 ≡ ∓9 なので、どちらも a¹⁰ ≡ 81 ≡ −1 を満たし S = 1。 a ≡ 20 のとき a^(p−1)/8 × z^(p−1)/4 ≡ a⁵ × 3¹⁰ ≡ −9 × 9 ≡ +1 なので T = 0。【ト】から:
　　√20 ≡ 20^(41+7)/16 × 3^(41−1)/8 ≡ 20³ × 3⁵ ≡ 26
実際 26² ≡ 20。最後に a ≡ 21 ≡ −20 のとき a^(p−1)/8 × z^(p−1)/4 ≡ a⁵ × 3¹⁰ ≡ −1 なので T = 1。【ト】から:
　　√21 ≡ 21^(41+7)/16 × 3^{[(41−1)/8 + (41−1)/4]} ≡ 21³ × 3¹⁵ ≡ 12
実際 12² ≡ 21。

具体的な数値例からも、Tonelli の方法はうまくいきそうだという感触が得られた。

〔注5〕　指数計算（べき剰余）の例。18² ≡ 324 ≡ 328 − 4 ≡ −4　［←41×8=328］
18³ ≡ 18² × 18 ≡ (−4) × 18 ≡ −72 ≡ −82 + 10 ≡ 10　［←41×(−2)=−82］
18⁴ ≡ (18²)² ≡ (−4)² ≡ 16
18⁵ ≡ 18⁴ × 18 ≡ 16 × 18 ≡ 288 ≡ 1　［←41×7=287］
この場合、18⁵ については、こうした方が少し速い: 18⁵ ≡ 18² × 18³ ≡ −4 × 10 ≡ −40 ≡ 1

しげちーの文脈では、√21 − 3 ≡ 12 − 3 ≡ 9 ≡ 132 なので、その 1⁄6 の 22 がスペシャルデー。40 − 22 ≡ 18 も。要するに n ≡ 0, 20, 40 のとき Q(n) は 41 で割り切れ（基本の日）、n ≡ 18, 22 のときも Q(n) は 41 で割り切れる。事実 Q(18) = 432345 = 41 × 10545。 Q(22) = 1151403 = 41 × 28083。 1個41円の物をまとめ買いしたいなら、これらの日に買い物すれば、集めた1円玉を無駄なく使い切れる！

Alberto Tonelli（アルベルト・トネッリ）は、1849年か1850年 [1] のクリスマスの日に、イタリア中部トスカーナ地方の地中海に近い町 Lucca（ルッカ）で生まれた。地元トスカーナの Pisa（ピサの斜塔がある所）の大学を卒業後、ドイツの Göttingen に留学。この時代、ドイツの数論・代数学は世界一ィィィだったのだ。帰国後、1877年、シチリア島の Palermo 大学で解析学を教え、1879年にローマ大学に移り、1904年から1919年まで、ローマ大学の総長を務めた [2]。

1891年の Bemerkung über Auflösung quadratischer Congruenzen [3] （二次剰余の解法についての注記）は、ローマで書かれた。簡潔な内容ながら、mod p での平方根を求める史上初の確率的多項式時間アルゴリズムの発見・記述であり（少なくとも離散対数的アプローチにおいて）、数学史に残る論文だ。Tonelli はイタリア人だが、この論文をドイツ語で書き、ドイツの専門誌で発表した（ちなみに当時のドイツ語では Kongruenz は Congruenz と表記されていたようだ――Dirichlet もこの表記を使っている）。今で言う「Tonelli–Shanks のアルゴリズム」の Tonelli という部分は、このパイオニア論文を指す。

Lucca の町の人は、いろいろと地元の文化に誇りを持っているようだが、Lucca 出身の Alberto Tonelli の存在は忘れられているようだ。Gauß もできなかった直接解法を数学史上、初めて確立したのだから、もっと評価されてもいいのに…。Lucca の名物は、Buccellato（ブッチェッラート）という、レーズンとアニシードの入った輪型のパンだそうだ…。

Tonelli にとってある意味、不運な展開として、1891年の画期的論文は、1973年に Daniel Shanks が独立に発見した同様のアルゴリズム、いわゆる RESSOL によって「上書き」されてしまった。Shanks 版の方が実装上、少し効率がいいので、旧バージョンはもう要らないと…。ここではアルゴリズムの効率より原理・内容の理解を優先し、歴史的順序に従って、まずはシンプルな Tonelli 版を研究したい。

1900年前後は、この問題が大きく進歩した時期だった。1907年には、同じイタリアの Michele Cipolla が非常に美しい別のアルゴリズムを発見している。その他、同時期に Lucas 数列的なアルゴリズムも発見されている。実装上の諸問題は、21世紀に入っても完全には解決していない。楕円曲線の Ed25519 で有名な Bernstein は「p−1 が 2 の高いべきで割り切れるケース」を annoying（厄介）と表現し、2001年に Faster square roots in annoying finite fields という preprint [10] を公開している。大別して「離散対数」と「体（たい）の理論」の2系列のアルゴリズムが共存しているものの、どちらもヘルパーとなる非剰余を試行錯誤で見つける必要がある――実用上は試行錯誤で問題ないものの、「うまくいくまで、ひたすら試行錯誤」という曖昧さが理論上のネックとなってしまい、高速性が保証された決定論的アルゴリズムが得られていない。Bernstein の方法は、繰り返し必要になるデータを事前計算して「キャッシュ」しておくことで高速化を図るもので、p が固定されている場合、特に有効なようだ（p が固定されているなら、事前に非剰余も用意しておける）。一般の場合においても、条件に応じてキャッシュの大きさを適宜調整することで、Tonelli–Shanks のアルゴリズムは「Tonelli–Shanks–Bernstein のアルゴリズム」に発展するかもしれない。少なくとも理論上では、p−1 が 2 で何度も割り切れる「厄介」なケースでは、この方法より Cipolla 系のアルゴリズムの方が高速なようだ。

Gauß 自身の方法 methodus exclusionis (DA 319–322) は、p が小さい場合に実用性があり、歴史的には長年にわたって活用されてきた。時代遅れで非効率のようだが、遊びとして試してみると、得られるものもあるかも…。半面、Gauß が「直接的解法は面倒で実用にならない」と断言して間接的解法を推奨したことは（DA 152）、「Gauß がそう言うならそうなんだろう」という先入観を生じさせ、直接的解法の研究が遅れる一因になった可能性がある。

〔注〕　Alberto Tonelli は Leonida Tonelli とは別人（ウィキペディア、フランス語版ではこの二人が混同されている）。

2023-05-02　Tonelli のアルゴリズム　源流・1891年バージョン

#数論　#平方根 mod p

平方根 mod p の計算法の、現在の事実標準は Shanks の RESSOL（1973年）。それを研究する前に、1891年の元祖 Tonelli バージョンを紹介したい。

あえてシンプルな Tonelli 版を考えることで、問題の性質がはっきりして、Shanks 版（Tonelli 版をわずかに改善）の工夫も味わえる。しかも Tonelli 版なら、抽象代数的概念を使わず、古典的な言葉で分かりやすく説明できる。

何の役に立つのか？　現在のネットの通信では、楕円曲線というものが重要な基礎となっている。この曲線は、
　　y² ≡ x の3次式
のような形を持つため、曲線上の点の x 座標が与えられて y 座標を求めたいとき、
　　y² ≡ 何か（3次式の値）
を解く必要があり、「何か」の平方根が欲しい。数論として面白いだけでなく、現実の通信でも活躍するアルゴリズムといえよう。とはいえ、これは娯楽の数論なので、ワクワクできれば現実的に役立つかどうかは、どうでもいい！

この問題は古典的な数論の範囲を超えていて、2020年代の現在でも研究が続いている。既存の文献は、多少の予備知識がないと理解しにくい。ここでは mod の世界での足し算・掛け算、その世界で平方根があるかないか――といった、ごく単純な概念だけを使う。分かりやすさの代償として、群論的・現代的なエレガントで広い観点が得られないけど、それは後のお楽しみ…？

§42.　前回まで、3回にわたって p が 2q+1 型、4q+1 型、8q+1 型の場合を個別に扱った（q: 奇数）。 8q+1 ともなると、式が個別的過ぎて、かえって真意が分かりにくい。視点をもう少しだけズームアウトして、全体像を眺めてみたい。

3 以上の任意の素数 p = 2^eq+1 を考える。ここで e は p に応じて定まる正の整数、2^e は 2, 4, 8, 16, … などの数、q は正の奇数。

〔補足〕　p は（3 以上の素数だから 2 では割り切れずに）奇数。だから p−1 は偶数。つまり p−1 は必ず 2 で割り切れる。割り算した商がまた偶数なら、もう一度 2 で割り切れる。その商も偶数なら、さらにもう一度 2 で割り切れる…。2 で無限回は割り切れないので、いつかは商が奇数 q となる。つまり、素数 p が何であれ、p−1 は
　　p−1 = 2の累乗 × 奇数 = 2^eq
の形で表される。両辺に 1 を足すと:
　　p = 2^eq + 1
p が決まれば、その p に対して e と q の値は確定的。e は「p−1 が 2 で何回割り切れるか？」を表す。

〔例〕　11 = 2¹⋅5+1（このとき e = 1, q = 5）　13 = 2²⋅3+1（このとき e = 2, q = 3）
17 = 2⁴⋅1+1（このとき e = 4, q = 1）　19 = 2¹⋅9+1（このとき e = 1, q = 9）

【1】　Euler の基準から x² ≡ a (mod p) が解を持てば a^(p−1)/2 = +1 だが、もし
　　a^奇数 ≡ +1 (mod p)　　【ナ】
の形を作れれば、問題が解決することは、既に見た通り。なぜなら、その場合、両辺を a 倍すれば
　　a^偶数 ≡ a (mod p)　　【ニ】
となり、その自然な平方根を考えることができる。ここでは p−1 = 2^eq と考えているので、【ナ】の形が（Euler の基準の結果として）生じるのは、
　　(p−1)/2 = 2^e−1q　が奇数の場合
つまり e = 1 の場合（そのとき 2^e−1q = 2⁰q = 1q = q = 奇数）。

【2】　たとえ (p−1)/2 = 2^e−1q が奇数でないとしても、もしも次のような変形が成り立つなら、【ナ】の形にできる。

〔例1〕　p = 37, a = 7 のとき（§37【カ】参照）、p−1 = 2²⋅9 で (p−1)/2 = 2¹⋅9 = 18:
　　7¹⁸ ≡ +1 (mod 37)
　　「両辺の平方根から」　7⁹ ≡ +1 (mod 37)

〔例2〕　p = 41, a = 18 のとき（§41）、p−1 = 2³⋅5 で (p−1)/2 = 2²⋅5 = 20:
　　18²⁰ ≡ +1 (mod 41)
　　「両辺の平方根から」　18¹⁰ ≡ +1 (mod 41)
　　再び「両辺の平方根から」　18⁵ ≡ +1 (mod 41)

指数が奇数になるまで何度でも「両辺の平方根」を考え、指数を半々にしていけば、いつかは必ず奇数乗になる。 Tonelli のアイデアの出発点ともいえる単純明快なロジックだが、実際には、そんなに簡単でもない: 一般に X を「何か」として X² ≡ +1 のとき、両辺の平方根を考えると、X ≡ +1 になるかもしれないし X ≡ −1 になるかもしれない。〔例1〕で X ≡ 7⁹ と置けば、7¹⁸ の式は X² ≡ +1 なので X ≡ +1 or X ≡ −1 だが、たまたま前者が正しいので、7⁹ ≡ 1 という結論も正しい。〔例2〕も、2回連続で、平方根の計算がたまたま正しい。一般には「たまたま正しい」という保証がないので、この部分を積極的にコントロールする必要がある。
　　X^偶数 ≡ +1
…の形のとき、両辺の平方根を取った結果として、
　　X^{偶数の半分} ≡ +1
　　X^{偶数の半分} ≡ −1
のどちらが生じても、適切に対応しよう。前者の場合、何もしなくていいが、後者の場合、両辺を −1 倍して、+1 に戻す「補正」が必要。単に両辺を −1 倍するのでは、うまくいかない。

〔例3〕　p = 37, a = 3 のとき（§37【ガ】参照）:
　　3¹⁸ ≡ +1 (mod 37)
　　3⁹ ≡ −1 (mod 37)
ここでもし両辺を単純に −1 倍すると:
　　−3⁹ ≡ +1 (mod 37)
両辺に a = 3 を掛けると:
　　−3¹⁰ ≡ +3 (mod 37)
　　つまり　3¹⁰ × (−1) ≡ +3 (mod 37)
両辺の平方根を考えると:
　　3⁵ × √−1 ≡ ±√3 (mod 37)
3⁵ は普通に計算できるので √3 が求まりかけてはいるものの、邪魔くさい √−1 が入り込み、話がややこしくなっている…。これを容認すると、一般のケースでは √−1 のそのまた平方根の平方根の…みたいなゴチャゴチャが発生してしまう。

同じ −1 倍でも、Euler の基準の非剰余の場合 z^(p−1)/2 ≡ −1 を使うのが鍵となる！

〔例3・改善版〕　p = 37, a = 3 のとき、仮に z = 2 が非剰余であることが分かっているとしよう。すると Euler の基準から:
　　2¹⁸ ≡ −1 (mod 37)
この情報があれば、
　　3⁹ ≡ −1 (mod 37)
が発生したとき、両辺を単純に −1 倍する代わりに、左辺同士・右辺同士を掛け算して、こうできる:
　　3⁹ × 2¹⁸ ≡ +1 (mod 37)
正しい符号で a の奇数乗が得られたので、両辺に a = 3 を掛けると:
　　3¹⁰ × 2¹⁸ ≡ +3 (mod 37)
両辺の平方根を考えると:
　　3⁵ × 2⁹ ≡ ±√3 (mod 37)
3⁵ も 2⁹ もただの整数なので、後は単純計算。

【3】　計算途中で平方根を取ったとき、たとえ右辺が何度 −1 になってしまっても、それを検出してそのたびに z^(p−1)/2 ≡ −1 を掛け算すれば…
　　何かの偶数乗 ≡ −1　　← 符号が不適切
…という状況を、すぐ次の形に戻せる:
　　何かの偶数乗 × z^(p−1)/2 ≡ +1　　← 適切な符号を回復

Tonelli のアルゴリズム　a^(p−1)/2 ≡ +1 からスタート。 a の奇数乗が出るまで、繰り返し両辺の平方根を取り、指数をどんどん半分にする。途中計算で、右辺が +1 なら何もしないが −1 になったら、z^(p−1)/2 ≡ −1 を掛け算して、符号を回復させる。
　イメージ的には「自動追尾」: ターゲットが何度180度のUターンをしても、ロックオンを続けて逃さない。

〔例4〕　p = 41, a = 21 すなわち mod 41 の世界で √21 を求める。非剰余 z = 3 つまり 3²⁰ ≡ −1 を使って自動追尾:
　　21²⁰ ≡ +1　　← Euler の基準による出発点
　　21¹⁰ ≡ −1　　← ターゲットの反転を検出しました！
　　21¹⁰ × 3²⁰ ≡ +1　　← 自動追尾・発動しました！
　　21⁵ × 3¹⁰ ≡ −1　　← 再びターゲットの反転を検出！
　　21⁵ × 3¹⁰ × 3²⁰ ≡ +1　　← 自動追尾・発動！
a の奇数乗が得られたので、両辺を a = 21 倍して:
　　21⁶ × 3¹⁰ × 3²⁰ ≡ +21
左辺の指数は全部偶数なので、自然な平方根を考えられる:
　　21³ × 3⁵ × 3¹⁰ ≡ ±√21
あとは単純計算。なかなかクールなアルゴリズムだ。

§43.　以上のことを定式化したい。一般論の前に、ミニチュアの例を使って、道筋を確認しておこう。

【1】　〔例4〕の p = 41, a = 21, z = 3 のケースにおいて…
　　21⁵ × 3¹⁰ × 3²⁰ ≡ +1
…の左辺は a^(p−1)/8 × z^(p−1)/4 × z^(p−1)/2 に当たるが、次のように解釈した方が見通しがいい:
　　a^q × z^2q × z^4q　ただし p = 8q + 1 で q = 5

〔補足〕　(p−1)/2, (p−1)/4, (p−1)/8 をそれぞれ 4q, 2q, 1q と再解釈している。一般論として、例えば (p−1)/2 は 2^e−1q に当たる。

【2】　ミニチュア（本番前の予行演習）。 p = 2⁴q + 1 = 16q + 1 を素数、q を奇数とすると、出発点は…
　　a^(p−1)/2 ≡ +1 (mod p)
これも次のように解釈した方が見通しがいい:
　　a^8q ≡ +1　　なぜなら p−1 = 16q, (p−1)/2 = 16q/2 = 8q
この両辺の平方根を考えたとき a^4q ≡ +1 になってくれれば一番楽だが、a^4q ≡ −1 になってしまっても、
　　z^(p−1)/2 ≡ z^8q ≡ −1　　ここで z は mod p の非剰余
を使って補正すればいい。補正あり・なし両方のケースをまとめて書くと（§40【ソ】参照）:
　　a^4q × z^8qS ≡ +1　　【ヌ】
　　ただし a^4q = +1 or −1 に応じて S = 0 or 1
【ヌ】の自然な平方根を考えると:
　　a^2q × z^4qS ≡ +1 or −1
+1 or −1 に応じて符号を補正をすると（§40【ツ】参照）:
　　a^2q × z^4qS × z^8qT ≡ +1　　【ネ】
同様に【ネ】の自然な平方根を考えて:
　　a^q × z^2qS × z^4qT ≡ +1 or −1　　【ノ】
　　つまり　a^q × z^2qS × z^4qT × z^8qU ≡ +1　　【ハ】
ここで U は、もちろん【ノ】が +1 か −1 かに応じて 0 or 1 の値を持つ。

とうとう【ハ】では a の指数が奇数 q になるので、後は一本道で √a が求まる。【ハ】の両辺を a 倍すると:
　　a^q+1 × z^2qS × z^4qT × z^8qU ≡ a　　【ヒ】
この左辺で q+1 は偶数だし（なぜなら q は奇数）、それ以外の指数は明らかに偶数なので、両辺の自然な平方根として ±√a が得られる。

【3】　注目点として、出発点となる指数も、補正で使われる z の指数も、どちらも初期値は (p−1)/2 = 8q なので、何度も半々にしていくと、そのうち必ず奇数になる。けれど、補正で使われる z^8q は、a^8q の指数を1回以上、半分にした後で掛け算される（【ヌ】の時点）――補正が必要になる可能性があるのは、平方根を考えた後、つまり指数を半分にした後だから。結果として: z^8q の指数も次第に半減するものの、a 側の指数が繰り返し半減してとうとう奇数 q になった瞬間において、z の肩に乗る指数は、最小でも 2q であり、こっちはまだ偶数。そのため【ヒ】の各指数は偶数であり、【ヒ】の自然な平方根を考えることができる。

この微妙なタイミングの差（z 側の指数の方が半減が一歩遅い）によって、Tonelli の戦略は成功する。

さて【ハ】の左辺は、指数法則により、こう整理される:
　　a^q × z^{(2qS + 4qT + 8qU)} ≡ a^q × z^{q(2S + 4T + 8U)} ≡ a^q × (z^q)^{(2S + 4T + 8U)}
簡潔化のため大文字の A, Z を使って A ≡ a^q, Z ≡ z^q とすると、上の式はこうなる:
　　A × Z^(2S+4T+8U) ≡ A × Z^2(S+2T+4U)

この速記法を使うと:

出発点は　A⁸ ≡ +1　　← 何かの8乗
　　【ヌ】は　A⁴ Z^8S ≡ [AZ^(2S)]⁴ ≡ +1　　← 何かの4乗
　　【ネ】は　A² Z^4S+8T ≡ [AZ^(2S+4T)]² ≡ +1　　← 何かの2乗
　　【ハ】は　AZ^2S+4T+8U ≡ [AZ^(2S+4T+8U)]¹ ≡ +1　　← 何かの1乗

結局、素数 p = 16q+1 = 2⁴q+1 を mod として:
　　AZ^2S+4T+8U ≡ AZ^2(S+2T+4U) ≡ +1

【4】　次の節で、この 2⁴ = 16 という数を 2^e に一般化する。結論を先取りすると、任意の素数 p = 2^eq+1 を mod として、次が成り立つ:
　　AZ^2c ≡ +1　ただし c = S + 2T + 4U + 8V + … + 2^e−2○　　【フ】
ここで ○ は、S, T, U, … などのスイッチ（0 or 1）の連鎖における「最後の変数名」。例えば、スイッチが S, T, U, V, W まであるとすると ○ は W。

〔補足〕　Tonelli の原論文では、このようなぎこちないネーミングを避けて、上記 S, T, U, …, ○ をそれぞれ ε₀, ε₁, ε₂, …, ε_s−2 としている（s が e に当たる）。その場合、われわれの c に当たる整数は、次のように明快に表現される:
　　ε₀ + ε₁2 + … + ε_s−22^s−2

S, T, U, V, …, ○ はそれぞれ 0 or 1 なので、c については、e−1 桁の二進数とも解釈可能。つまり最下位桁（2⁰ の位。ビット0）が S。 2桁目（2¹ の位。ビット1）が T。 3桁目（2² の位。ビット2）が U。等々となり、最上位桁（2^e−2 の位。ビット e−2）が ○。ただし上位桁が 1 とは限らない。Z の指数 2c は、その 2 倍なので、ちょうど e 桁の二進数に当たる（最下位桁は 0。上位桁は 1 とは限らない）。「自動追尾」のイメージからすると、2c は「ターゲットが反転して自分も反転する」という動作を記録した「姿勢制御コマンド」履歴に当たる。

〔参考〕　AZ^(2c) ≡ 1 を満たす 2c を決定することは、Z^2c ≡ 1/A を 2c について解くことと同等（1/A ≡ m として未知数を x = 2c と書けば Z^x ≡ m）。あるいは同じことだが A ≡ Z^−2c を解くことと同等（1/Z ≡ n とすれば A ≡ n^x）。離散対数問題の特別な場合に当たる。特別性は、Euler の基準から A ないし A⁻¹ の位数が 2^e−1 の約数（従って 2 のべき）になること、そして Z ないし Z⁻¹ の位数が 2^e に等しいこと。【3】の速記は「S, T, U, … として 2c を1ビットずつ決定できる」ことを示唆する。

§44.　証明の本番！　【フ】が成り立つことを確かめる。

変数名に使う文字が足りなくなってきたので、以下ではギリシャ文字 β（ベータ）、γ（ガンマ）、δ（デルタ）… を使う。

素数 p = 2^eq+1 を mod として（q は奇数）、a が平方剰余、z が非剰余であるとき、r² ≡ a の解 r を求めたい。e が 1 以下のときは自明なので、e を 2 以上の整数とする。出発点となるのは、Euler の基準:
　　a^(p−1)/2 ≡ +1, z^(p−1)/2 ≡ −1
つまり:
　　a の 2^e−1q 乗 ≡ +1, z の 2^e−1q 乗 ≡ −1
毎回 a の 1 乗から上がっていくのも面倒なので、A = a^q と Z = z^q という速記を使い、上記をこう書く:
　　A の 2^e−1 乗 ≡ +1, Z の 2^e−1 乗 ≡ −1
第1ステップは、前者の自然な平方根の計算:
　　a^(p−1)/4 ≡ ±1　つまり a の 2^e−2q 乗
　　速記を使うと　A の 2^e−2 乗 ≡ ±1
この数が ≡ +1 になれば問題ないが、もし −1 になるなら、両辺に Z の 2^e−1 乗 ≡ −1 を掛け算して、符号を補正する必要がある。上記の平方根が +1 or −1 のどちらになるか応じて S = 0 or 1 とすれば、結局、左辺に Z の 2^e−1S 乗を掛け算することに相当する（S = 0 なら Z⁰ = 1 の掛け算であり、実質何もしないのと同じ）。次の関係が成り立つ:
　　A の 2^e−2 乗 × Z の 2^e−1S 乗 ≡ +1
A や Z の指数が「そのまた指数」を含んでる。つまり、こーゆーこと:

A^(2e−2) × Z^(2e−1S) ≡ +1

読みやすくするため β = 2^e−2 と置くと、2^e−1 = 2(2^e−2) = 2β なので 2^e−1S = 2βS となり、上の式をこう書ける:
　　A^β × Z^2βS ≡ +1　言い換えれば　[AZ^(2S)]^β ≡ +1

もし e = 2 なら β = 2²⁻² = 1 なので、上の式では A¹ = a^q が現れて、これは a の奇数乗だから √a を求めることができる。一方、もし e が 3 以上のときは、β = 2^e−2 = 2, 4, 8, 16, etc. なので、第2ステップとして、もう一度、両辺の自然な平方根を求め、指数を半減させる必要がある:
　　[AZ^2S]^β/2 ≡ ±1
この数が +1 or −1 のどちらになるかに応じて T = 0 or 1 として、第1ステップ同様、左辺に Z の 2^e−1T 乗を掛け算すると――言い換えると Z^2βT を掛け算すると――、右辺を +1 に保つことができる:
　　[AZ^2S]^β/2 × Z^2βT ≡ +1　つまり　[AZ^2S × Z^4T]^β/2 ≡ +1
　　指数法則から　[AZ^(2S+4T)]^β/2 ≡ +1
今 γ = β/2 = 2^e−2/2 = 2^e−3 と置くと:
　　[AZ^(2S+4T)]^γ ≡ +1
2^e−1 = 2²(2^e−3) = 4γ ということに注意しておく。

もし e = 3 なら γ = 2³⁻³ = 1 なので、上の式では A¹ = a^q が現れて、これは a の奇数乗だから以下同文。一方、もし e が 4 以上のときは、γ = 2^e−3 = 2, 4, 8, 16, etc. なので、第3ステップとして、もう一度、両辺の自然な平方根を求め、指数を半減させる必要がある:
　　[AZ^(2S+4T)]^γ/2 ≡ ±1
この数が +1 or −1 のどちらになるかに応じて U = 0 or 1 として、第2・第3ステップ同様、左辺に Z の 2^e−1U 乗を掛け算すると――言い換えると Z^4γU を掛け算すると――、右辺を +1 に保つことができる:
　　[AZ^(2S+4T)]^γ/2 × Z^4γU ≡ +1　つまり　[AZ^(2S+4T) × Z^8U]^γ/2 ≡ +1
　　指数法則から　[AZ^(2S+4T+8U)]^γ/2 ≡ +1
今 δ = γ/2 = 2^e−3/2 = 2^e−4 と置くと:
　　[AZ^(2S+4T+8U)]^δ ≡ +1
2^e−1 = 2³(2^e−4) = 8δ ということに注意しておく。

この調子で進めると β = 2^e−2, γ = 2^e−3, δ = 2^e−4, … は、どんどん半減するので、いつかは角かっこの外側の指数が 1 = 2^e−e になって、こうなる:
　　[AZ^{(2S+4T+8U+16V+…)}]¹ ≡ +1　　【ヘ】
問題は、これが起きた瞬間、丸かっこ内の指数
　　2S + 4T + 8U + 16V + …
の末尾の項は何か？　丸かっこの中身は――
　　e = 2 なら 2S　　← 最後の係数は 2¹
　　e = 3 なら 2S + 4T　　← 最後の係数は 2²
　　e = 4 なら 2S + 4T + 8U　　← 最後の係数は 2³
　　………
　　一般の e に対して 2S + 4T + 8U + 16V + … + 2^e−1○　　← 最後の係数は 2^e−1
――であるから、【ヘ】の一般形は:
　　AZ^(2c) ≡ +1　ただし 2c = 2S + 4T + 8U + 16V + … + 2^e−1○
この「ただし」の後ろの式の両辺を 2 で割ると、【フ】が得られる。Tonelli の公式が証明された！

§45.　例題。素数 p = 673 = 2⁵⋅21+1 を mod として（e = 5, q = 21）、a = 21 は平方剰余、z = 5 は非剰余であることが分かっているとしよう。この条件で r² ≡ 21 の解 √21 を求める。

【1】　Tonelli 自身の方法を額面通り実行する。指数 (673−1)/4 = 168 からスタート。
　　21¹⁶⁸ ≡ −1 (mod 673)
なので、5^(673−1)/2 ≡ 5³³⁶ ≡ −1 を使って補正:
　　21¹⁶⁸⋅5³³⁶ ≡ +1　両辺の平方根から
　　21⁸⁴⋅5¹⁶⁸ ≡ −1　なので
　　21⁸⁴⋅5¹⁶⁸⋅5³³⁶ ≡ +1
　　つまり　21⁸⁴⋅5^168+336 ≡ +1
　　21⁸⁴⋅5⁵⁰⁴ ≡ +1　両辺の平方根から
　　21⁴²⋅5²⁵² ≡ −1　なので
　　21⁴²⋅5²⁵²⋅5³³⁶ ≡ +1
　　つまり　21⁴²⋅5^252+336 ≡ +1
　　21⁴²⋅5⁵⁸⁸ ≡ +1　両辺の平方根から
　　21²¹⋅5²⁹⁴ ≡ −1　なので
　　21²¹⋅5²⁹⁴⋅5³³⁶ ≡ +1
　　つまり　21²¹⋅5^294+336 ≡ +1
　　21²¹⋅5⁶³⁰ ≡ +1　ゆえに
　　21²²⋅5⁶³⁰ ≡ 21
　　√21 ≡ ±21¹¹⋅5³¹⁵ ≡ ±201
　　検算　201² = 40401 ≡ 21 (mod 673)

【2】　現実的には、最初に A ≡ a^q ≡ 21²¹ ≡ 464 ≡ −209★ を求めて、次のようにした方が明快かもしれない:
　　21⁴² ≡ (−209)² ≡ 609 ≡ −64 (≡ A²)
　　21⁸⁴ ≡ (−64)² ≡ 58 (≡ A⁴)
　　21¹⁶⁸ ≡ 58² ≡ −1 (≡ A⁸)
同様に:
　　Z ≡ z^q ≡ 5²¹ ≡ 118
　　5⁴² ≡ −209 (≡ Z²)
ここから★以下を再利用できる:
　　5⁸⁴ ≡ −64 (≡ Z⁴)
　　5¹⁶⁸ ≡ 58 (≡ Z⁸)
　　5³³⁶ ≡ −1 (≡ Z¹⁶)

21¹⁶⁸ ≡ −1 は計算済みなので:
　　21¹⁶⁸⋅5³³⁶ ≡ +1　　←　A⁸Z¹⁶ ≡ +1
事前計算から:
　　21⁸⁴⋅5¹⁶⁸ ≡ (58)(58) ≡ 58² ≡ −1　　←　A⁴Z⁸ ≡ −1
　　21⁸⁴⋅5¹⁶⁸⋅5³³⁶ ≡ +1　　←　A⁴Z⁸⁺¹⁶ ≡ +1
事前計算から:
　　21⁴²⋅5⁸⁴⋅5¹⁶⁸ ≡ (−64)(−64)(58) ≡ (58)(58) ≡ −1　　←　A²Z⁴ ≡ −1
　　21⁴²⋅5⁸⁴⋅5¹⁶⁸⋅5³³⁶ ≡ +1　　←　A²Z⁴⁺⁸⁺¹⁶ ≡ +1
事前計算から:
　　21²¹⋅5⁴²⋅5⁸⁴⋅5¹⁶⁸ ≡ (−209)(−209)(−64)(58) ≡ (−64)(−64)(58) ≡ (58)(58) ≡ −1　　←　A¹Z² ≡ −1
　　21²¹⋅5⁴²⋅5⁸⁴⋅5¹⁶⁸⋅5³³⁶ ≡ +1　　←　A¹Z^2+4+8+16 ≡ +1
両辺を 21 倍して:
　　21²²⋅5⁴²⋅5⁸⁴⋅5¹⁶⁸⋅5³³⁶ ≡ 21　　←　a^(q+1)Z⁴⁺⁸⁺¹⁶ ≡ a
ゆえに:
　　√21 ≡ ±21¹¹⋅5²¹⋅5⁴²⋅5⁸⁴⋅5¹⁶⁸ ≡ ±21¹¹(118)(−209)(−64)(58) ≡ ±201　　←　a^(q+1)/2Z²⁺⁴⁺⁸

何か面白いことが起きていて、両辺の平方根を取ると、毎回 −1 が出てくる。平方根を求めようとしている a = 21 が、たまたま q = 21 に等しい。それと関係あるのかも…？

【3】　前節 §44 は [5] の Algorithm 2.3.8 に基づく。2c を C と置いて、整理すると（β, γ, δ, … などの 2 の累乗を B で表す）:

C = 2c = 0 とループ変数 j = 2 が初期値。
　　B = 2^e−j が 2^e−e = 1 になるまで j = 2, 3, 4, … と増やし B を半減させる。各ステップにおいて、
　　　　(AZ^C)^B を計算して、もし ≡ −1 になったら、
　　　　C の値を更新: 古い C の値に 2^j−1 を足したものをあらためて C とする。

上記核心部をアルゴリズム全体に埋め込むと、こうなる。

Tonelli のアルゴリズム　法 p = 2^eq + 1 と平方剰余 a と非剰余 z が与えられたとき（q は奇数）:
定数 A ≡ a^q, Z ≡ z^q を求め、変数 C の初期値を 0 とせよ。
　　j = 2, 3, 4, …, e に対して:
　　　　B = 2^e−j として
　　　　もし (AZ^C)^B が ≡ −1 になったら、C の値を更新:
　　　　→　古い C の値に 2^j−1 を足したものをあらためて C とする（★）
　　　　もし (AZ^C)^B が ≡ +1 になったら、何もしない（C の値は元のまま）
ループを抜けたとき AZ^C ≡ +1 なので（なぜなら最後の B は 2⁰ = 1）:
　　aAZ^C ≡ a　つまり　a^q+1⋅Z^C ≡ a
　　→　指数 q+1 と指数 C は偶数なので √a が自然に求まる。

§43 のミニチュアや、§44 を検討すれば、そうなってること・それでうまくいく理由を具体的に確認できるだろう。すなわち (AZ^C)^B が ≡ −1 になってしまった場合には（★）によって、丸かっこ内の旧 Z が Z^2j−1 倍される。丸かっこの外に B 乗つまり 2^e−j 乗があるので、(AZ^C)^B 全体としては、値が (Z^2j−1)^2e−j 倍つまり Z^2e−1 倍される。Z ≡ z^q なので、この倍率は ≡ (z^q)^2e−1 = z^(2e−1q) = z^(2eq)/2 = z^(p−1)/2 ≡ −1。結局、C の数値の更新により、「悪い」値 (AZ^C)^B ≡ −1 は −1 倍されて (AZ^C)^B ≡ +1 という「良い」関係が回復する（mod p において）。

〔注〕　文献では 2 の肩に乗る e を 1 小さくして、B = 2^e−j の代わりに B = 2^(e−1)−j のように表記していることがある。その場合、引く数も 1 小さくして、j = 2, 3, 4, … e の代わりに j = 1, 2, 3, … e−1 とする。どっちの書き方でも、内容は同じ。後者は「ループ変数が 1 から始まり e 未満」というところが、整然として気持ちがいいし、補正の数値についても 2^j−1 の代わりにスッキリ 2^j と書ける。半面、B の表現については、前者の 2^e−j の方が、後者の 2^e−1−j よりシンプル。ところで 2^e−1q は (p−1)/2 だが、a^(p−1)/2 ≡ +1 は Euler の基準から分かり切っている。一方 2^e−2q は (p−1)/4 だが、a^(p−1)/4 ≡ ±1 は、符号がどっちになるか分からない。そこから計算が始まる。だから 2 の指数 e−j が j = 2 から始まるのは、素直な表記だろう（エレガントかどうかはともかく）。

例題では e = 5 なので、j = 2, 3, 4, 5 と動き、B = 2³, 2², 2¹, 2⁰ と動く。つまり B = 8, 4, 2, 1。最初に:
　　[(−209)(118)⁽⁰⁾]⁸ ≡ (−209)⁸ ≡ −1　なので C = 0 から C = 2 に更新
　　[(−209)(118)⁽²⁾]⁴ ≡ [(−209)(−209)]⁴ ≡ −1　なので C = 2+4 = 6 に更新
　　[(−209)(118)⁽²⁺⁴⁾]² ≡ [(−209)(−209)(−64)]² ≡ −1　なので C = 2+4+8 = 14 に更新
　　[(−209)(118)^(2+4+8)]¹ ≡ [(−209)(−209)(−64)(58)]¹ ≡ −1　なので C = 2+4+8+16 = 30 に更新
すなわち:
　　(−209)(118)^(2+4+8+16) ≡ (−209)(118)³⁰ ≡ +1
両辺を 21 倍して:
　　21⋅21²¹⋅118³⁰ ≡ 21
　　21²²⋅118³⁰ ≡ 21
　　√21 ≡ ±21¹¹⋅118¹⁵ ≡ ±21¹¹⋅118^1+2+4+8 (mod 673)
この先は【2】と同じ。

このように整理すると、元祖 Tonelli バージョンには潜在的に「無駄」があることが感じられる。上の例では、第1ステップで 8 乗が計算されるが、8 乗を求める途中で 2 乗や「4 乗」も求まる。第2ステップで、今度は「4 乗」が計算される。もし「姿勢制御」がなければ（約50%の確率）、角かっこ内は一定なので、第1ステップで「4 乗」を求め、第2ステップで同じ数の「4 乗」を再計算するのは、無駄な二度手間。Shanks 版では、この点が改善されている。

手計算の時代の Tonelli にとって、「既に4乗が計算してあれば、それを再利用する」というのは「言うまでもない当たり前のこと」だったに違いない。Tonelli のアルゴリズムに「無駄」があるというより、「常識なので、いちいち明記してないだけ」とも考えられる。

【4】　しげちーの文脈では、√21 ≡ 201 (mod 673) ということから (201 − 3)/6 ≡ 33 がスペシャルデーになる。つまり:
　　Q(33) = 1⁴ + 2⁴ + 3⁴ + … + 33⁴ = 8432017
は、673 で割り切れる。同じことは 672 − 33 = 639 にも当てはまる。一方 n ≡ 0, 336, 672 (mod 673) は基本の日。1個673円の物をまとめ買いしたいなら、33日目がちょうど無駄なく1円玉を使い切るチャンス。843万枚の8トンの1円玉を動かす方法があって、受け取り拒否されなければだが…ｗ

Tonelli は論文の末尾で「これらの式は、単に理論的に成り立つだけでなく、実用的な計算法だ」という趣旨のコメントをして、実用性を強調している。他方、このアルゴリズムは、理論的にも、幾つかの深い問題と関連している。特に重大なのは「非剰余 z を一つ選択する」という部分。

1 から p−1 までの p−1 種類の数の中には、平方剰余と非剰余が50%ずつ含まれている。でたらめに選んだとして、10回連続で非剰余にならない確率は1000分の1以下（およそ 1/2¹⁰ = 1/1024）、20回連続で非剰余にならない確率は100万分の1以下。事実上、必ず非剰余が見つかる。けど、これは Las Vegas 型のギャンブル。極めて低い確率とはいえ、アンラッキーが連続して「1日ルーレットを回し続けても、なぜか一度も非剰余が出ず、平方根の計算を開始できない」という可能性がある。

「非剰余が簡単に見つかる」ことを前提とするなら、平方根の計算は「易しい」。その意味では、「非剰余が簡単に見つかる」ことと「平方根を簡単に計算できる」ことは、等価なのかもしれない。この「等価性」は、いささか逆説的だ: p−1 種類（この値は10億くらいかもしれない）の要素の中に、非剰余は50%（この値は5億くらいかもしれない）も含まれているが、平方根は 2 個しかないのだから、前者は見つけやすく、後者は見つけにくいはず――そう感じるのは錯覚で、本質的には、どちらも同じ難易度なのだろうか？

〔追記〕　n が求める平方根かどうかは 2 乗するだけで確かめられるが、n が非剰余かどうか Euler の基準によって確かめようとすると (p−1)/2 乗する必要がある。「非剰余の存在密度は大きいけれど、個々の n について、平方根かどうかの判定に比べ、非剰余かどうかの判定は、計算量が大きい」ということによって、上記の等価性をある程度は「正当化」できる。もっとも、ブルートフォースで平方根を探すことは一般には実用的でないが、非剰余を全数検索で見つけることは実用的; 非剰余のブルートフォース探索は「最悪でもすぐ終わる」ように思える。ところが最悪の場合の計算時間について、理論上の評価は、非常に弱い。「現実に経験・観測されることが、まだ理論的に説明できない」というのが、上記の奇妙さの本質だろう――「現実には…と感じていることが実は錯覚で、理論的には間違いだった」という可能性を（まだ）排除できない。（2023年7月4日）

〔参考文献〕
　　[1] Edizione Nazionale Mathematica Italiana - Alberto Tonelli
　　　　https://matematicaitaliana.sns.it/autori/1273/
　　[2] Alberto Tonelli - Wikipedia
　　　　https://it.wikipedia.org/wiki/Alberto_Tonelli
原論文:
　　[3] Tonelli, Bemerkung über Auflösung quadratischer Congruenzen
　　　　https://gdz.sub.uni-goettingen.de/dms/resolveppn/?PPN=GDZPPN002525739
　　　　簡易正誤表: 345ページ6行目 c の指数 α2^s−1 を α2^s−2 に、7行目 ε₁ = 1 を ε₀ = 1 に。
Shanks 版ではなく Tonelli 版を紹介しているのは:
　　[4] Bach & Shallit, Algorithmic Number Theory, vol. 1
　　[5] Crandall & Pomerance, Prime Numbers: A Computational Perspective (2nd ed.), §2.3.2
　　[6] Dickson, History of the Theory of Numbers, vol. 1, p. 215
　　　　https://archive.org/details/historyoftheoryo01dick
Shanks 版を紹介しているのは:
　　[7] Niven & Zuckerman & Montgomery, An Introduction to the Theory of Numbers (5th ed.)
　　[8] Cohen, A Course in Computational Algebraic Number Theory, §1.5
　　　Cf. Knuth: The Art of Computer Programming, Vol. 2, §4.6.2, Exercise 15
　　[9] Wikipedia, Tonelli–Shanks algorithm
　　　　https://wikiless.esmailelbob.xyz/wiki/Tonelli%E2%80%93Shanks_algorithm?lang=en
　　　　https://wikiless.northboot.xyz/wiki/Tonelli%E2%80%93Shanks_algorithm?lang=en
　　　　https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
Tonelli 版 Shanks 版その他の分析・改善:
　　[10] Bernstein, Faster square roots in annoying finite fields
　　　　https://cr.yp.to/papers.html#sqroot
離散対数問題の観点から:
　　[11] Galbraith, Mathematics of Public Key Cryptography, chap. 2, §2.9; chap. 13
　　　　https://www.math.auckland.ac.nz/~sgal018/crypto-book/crypto-book.html
　　　次の文献は参考になる:
　　　[22] Victor Shoup, A Computational Introduction to Number Theory and Algebra (Version 2), §12.5; §11.2
　　　　https://www.shoup.net/ntb
追加の資料:
　　[12] Eric Bach and Klaus Huber, Note on Taking Square-Roots Modulo N
　　　　https://minds.wisconsin.edu/bitstream/handle/1793/11024/file_1.pdf?sequence=1&isAllowed=y
　　[13] Rajeev Kumar, An algorithm for finding square root modulo p
　　　　https://arxiv.org/abs/2008.11814
　　[14] Palash Sarkar, Computing square roots faster than the Tonelli-Shanks/Bernstein algorithm
　　　　https://www.aimsciences.org/article/doi/10.3934/amc.2022007
　　　　https://eprint.iacr.org/2020/1407
　　[15] Tsz-Wo Sze, On Taking Square Roots without Quadratic Nonresidues over Finite Fields
　　　　https://arxiv.org/abs/0812.2591

[4][11] は視野が広い。[5][7] は比較的平明。[8] はプログラミング志向。[13] だけ見ると怪しげだが、このアイデアが [14] のヒントになったらしい。今回のメモは [3][5] をベースにした。

2023-05-05　Atkin の平方根アルゴリズム　mod 8k+5

#数論　#平方根 mod p

8k+5型素数 p を mod とする平方根について、比較的最近（1992年）の Atkin のアルゴリズムを紹介したい。a が平方剰余のとき y ≡ (2a)^k, i ≡ (2a)y² と置くと √a ≡ ay(i − 1) という、手品のような計算。

§46.　p が8k+5型の素数のとき、mod p の平方根については、既に紹介したように古典的解法もあるし、Pocklington の公式もあるが（§38）、そのどちらでも a^(p−1)/4 の値を調べて ≡ +1 か −1 かによって、場合分けが必要。A. O. L. Atkin の方法は、このチェックと場合分けを無くしたもの。1992年に報告された。k = 0 の場合（p = 5）はトリビアルなので、以下 k ≥ 1 とする。

p = 8k+5 なら k = (p−5)/8 だが、このとき
　　a^(p−1)/4 = a¹⋅a^(p−5)/4 = a[a^(p−5)/8]² = a[a^k]²　　【マ】
は ≡ ±1 のどちらになるか分からない。【マ】を出発点にすると、場合分けが必要―― a は平方剰余なので a^(p−1)/2 ≡ +1 となり、その平方根 a^(p−1)/4 は ≡ +1 にも −1 にもなり得る。一方 Atkin のアルゴリズムでは、【マ】の a を 2a で置き換えて、その式の値を i とする:
　　i ≡ (2a)^(p−1)/4 = (2a)[(2a)^k]²　　【ミ】

2 は非剰余（第2補充法則より）、従って 2a も非剰余なので、今度は (2a)^(p−1)/2 ≡ −1 となり、その平方根 (2a)^(p−1)/4 は「平方すると −1 になる数」。この数自体にも ±√−1 という符号の曖昧さがあるものの、符号を気にせず【ミ】の値が i だ、と定義してしまおう。すると i² ≡ −1 となり、以下の巧妙な扱いが成り立つ。まず:
　　(i − 1)² = i² − 2i + 1² = −2i　　【ム】
当たり前の式 a² = a/2(2a) と上記【ミ】を使うと、「何かの平方 ≡ a の簡単な式」を作れる:
　　a²[(2a)^k]² = a/2(2a)[(2a)^k]² ≡ a/2× i
さらに、この式と【ム】を辺々掛け合わせると:
　　a²[(2a)^k]²(i − 1)² ≡ a/2× i(−2i)　　ここで i(−2i) = i(−i)2 = −i²2 = −(−1)2 = 2
　　つまり　[a(2a)^k(i − 1)]² ≡ a　　【メ】
【メ】から a(2a)^k(i − 1) は a の平方根（の一つ）。簡潔化のため y ≡ (2a)^k と置き、【ミ】と今の結論をまとめると:

Atkin の平方根アルゴリズム（1992年）　素数 p = 8k+5 を mod として a を平方剰余とする。このとき
　　y ≡ (2a)^k　言い換えれば　y ≡ (2a)^(p−5)/8
と置くと…
　　i ≡ (2a)y²　すなわち　i ≡ (2a)^(p−1)/4
…は −1 の平方根の一つ。 y, i を使って、こう書ける:
　　√a ≡ ay(i − 1)

この仕組みについては、次回、あらためて整理・検討する。

【例題１】　37 は8k+5型素数（k = 4）。 mod 37 で 21 は平方剰余（§32）。Atkin の方法で √21 を求める（a = 21, 2a = 42）。
　　y ≡ 42⁴ ≡ 5⁴ ≡ 25² ≡ (−12)² ≡ 144 ≡ −4　　← 37×4=148
　　i ≡ 42(−4)² ≡ 5⋅16 ≡ 80 ≡ 6　　← 37×2=74
　　　検算　6² ≡ 36 ≡ −1　［確かに 6 は −1 の平方根］
　　√21 ≡ ay(i − 1) ≡ 21(−4)(6 − 1) ≡ (−16)(−4)5 ≡ 320 ≡ −50 ≡ −13　　← 37×10=370
　　　検算　(±13)² ≡ 169 ≡ 21　　← 37×4=148　［§38と一致］

注意　−6 も −1 の平方根だが、だからといって i を −i ≡ −6 に置き換えると、【ム】が (−i − 1)² = −1 + 2i + 1 = (+2i) となって【メ】の右辺が −a になってしまう（実際、例題の途中計算の 21(−4)(6 − 1) を 21(−4)(−6 − 1) ≡ −4 に置き換えると (−4)² ≡ −21）。この計算では i と −i は違う意味を持つ。ただし、2種類の数 ±i のどちらがプラスでどちらがマイナスなのかは、a によって決まり、見掛け上の + − とは必ずしも一致しない（a の値によっては +i ≡ −6 で −i ≡ +6 ということもあり得る）。

【例題2】　101 も8k+5型素数（k = 12）。 mod 101 でも 21 は平方剰余。
　　42² ≡ 47, 42³ ≡ 42⋅47 ≡ −46, 42⁶ ≡ (−46)² ≡ −5
　　従って y ≡ 42¹² ≡ 25
　　i ≡ 42⋅25² ≡ 42⋅625 ≡ 42⋅19 ≡ −10　（検算: 平方すると ≡ −1）
　　√21 ≡ 21⋅25(−11) ≡ 20(−11) ≡ −18　（検算: 平方すると ≡ 324 ≡ 21）

【例題3】　109 も8k+5型素数（k = 13）。 mod 109 でも 21 は平方剰余。
　　42² ≡ 20, 42⁴ ≡ 400 ≡ −36, 42⁸ ≡ (−36)² ≡ −12
　　従って 42¹² ≡ (−12)(−36) ≡ −4, y = 42¹³ ≡ (−4)42 ≡ 50
　　i ≡ 42⋅50² ≡ 42(−7) ≡ 33　（検算: 平方すると^* ≡ −1）
　　√21 ≡ 21⋅50⋅32 ≡ 28　（検算: 平方すると^** ≡ 21）

*　33² = 11²⋅3² = 121⋅9 ≡ 12⋅9 = 108 ≡ −1　［121=109+12≡12］
**　28² = 14²⋅2² = 196⋅4 = −22⋅4 ≡ −88 ≡ 21　［109×2=218, 196=218−22≡−22］

例題2では (√21 − 3)/6 ≡ −21/6 ≡ 282/6 ≡ 47 と 100 − 47 ≡ 53 (mod 101) はスペシャルデー。Q(47) = 48343448 と Q(53) = 87633963 は 101 で割り切れる。例題3では (√21 − 3)/6 ≡ 25/6 ≡ −84/6 ≡ −14 ≡ 95 と 108 − 95 ≡ 13 (mod 109) はスペシャルデー。Q(13) = 89271 と Q(95) = 1588572976 は 109 で割り切れる。

英国人の数論学者 A. Oliver L. Atkin (1925–2008) は、擬素数テストの研究の副産物としてこのアルゴリズムを発見した。アルゴリズムは Siguna Müller により16+9k型素数に拡張され、2013年、ルーマニアの Rotaru & Iftene [16] により、任意の4k+1型素数に一般化された。

[16] A Complete Generalization of Atkin’s Square Root Algorithm
　　https://profs.info.uaic.ro/~siftene/fi125(1)04.pdf

2023-05-10　Atkin の平方根アルゴリズム（その2）　幾何学的イメージ

#数論　#平方根 mod p

8k+5型の素数 p について、Atkin の平方根アルゴリズムによると、a の平方根 mod p が存在するなら、それは2種類あって ay(i − 1) と −ay(i − 1) ≡ ay(1 − i) に等しい。ここで y ≡ (2a)^k, i ≡ (2a)y²。

そうなる理由について、以下の明快な説明が成り立つ。本題の Tonelli–Shanks との関連では脱線になるが、Atkin のアルゴリズムも研究に値する。

§47.　Atkin のアルゴリズム（§46）について。便宜上、i − 1 ではなく、符号を逆にした 1 − i を考える。 −i については、単位円上の偏角 −90° の点とイメージできる。平方すると（つまり偏角を2倍にすると）、偏角 −180° の −1 になる。

1 − i は、画像の水色の正方形の右下の点に当たる。この点は、偏角 −45° なので、イメージ的には単位円上の点 A の方向――つまり √−i の方向――にある。ただし一辺の長さ 1 の正方形の対角線上にあるのだから、原点からの距離は √2 のはずで
　　1 − i = √2 × √−i = √−2i
となるだろう。

この解釈は、複素数の世界では、完全に正しい。実際 1 − i は、平方すると −2i になるのだから、−2i の平方根の一つ:
　　(1 − i)² = 1² − 2i + i² = 1 − 2i − 1 = −2i

上の計算は 8k+5 型素数 p を mod とする世界でも、そのまま成り立つ。ただし、この mod p の世界では −i や √−2i に当たる数は存在するものの、その「成分」ともいえる √2 や √−i に当たる数は存在しない。 √2 が存在しないこと、つまり x² ≡ 2 (mod p) に解がないことは、第二補充法則による。√−i つまり x² ≡ i の解が存在しないことは、あまり明らかではないが、もしも x² ≡ i を満たす x があるなら、同じ x は、両辺を平方した x⁴ ≡ −1 を満たすはず――だがそれは「四乗剰余の第一補充法則」に反する（別の説明: g を原始根として i ≡ g^(p−1)/4 ≡ g^(8k+5−1)/4 ≡ g^2k+1 は奇数乗なので、平方非剰余）。同じ理由から √i も存在しない。

〔例〕　mod 13 では 0 を別にして、次の5種類の平方数（平方剰余）が存在する: (±1)² ≡ 1, (±2)² ≡ 4, (±3)² ≡ 9 ≡ −4, (±4)² ≡ 3, (±5)² ≡ −1, (±6)² ≡ −3。よって平方すると −1 になる数 ±i は（どちらを +i とするかは別として） +5 と −5 ≡ +7 で、どちらも平方数ではない。例えば −i ≡ −5 とすると、mod 13 では「平方すると −5 になる数」はないのだから、√−i に当たる数もない。それでも、平方すると −2i ≡ (−2)(−5) ≡ 10 ≡ −3 になる数は存在するのだから、√−2i に当たる数も存在する。その一つは 1 − i ≡ 1 − (−5) ≡ +6 に他ならない。

今 mod p での任意の平方剰余 a が与えられたとき、y ≡ √i/2a としよう（根号内の分数 i/2a は i/2a の意味。以下同様）。この場合も、2 は非剰余なので 2a も非剰余。その逆数 1/2a も非剰余。他方、上記のように i も非剰余。従って、非剰余と非剰余の積 i/2a は平方剰余となり y ≡ √i/2a は存在する（具体的な値は a によって異なり、図には書き込めないが、とにかくどこかに存在）。

まとめると:
　　y ≡ √i/2a　そして　1 − i ≡ √−2i　　‥‥①
　　辺々掛けると　y(1 − i) ≡ √i/2a × √−2i　　‥‥②
　　整理すると　y(1 − i) ≡ √1/a ≡ 1 / √a　　〔注1〕
これは y(1 − i) が √a の逆数であることを意味する。両辺を a 倍すると:
　　ay(1 − i) ≡ a / √a ≡ √a　　‥‥③
a の一つの平方根が求まる！

もっと透明な説明:
　　y ≡ √i/2a　つまり　y² ≡ i/2a
　　同様に　(1 − i)² ≡ −2i　　‥‥④
　　だから　y²(1 − i)² ≡ i/2a × (−2i) ≡ 1/a　　〔注1〕
　　従って　[ay(1 − i)]² ≡ a²y²(1 − i)² ≡ a²/a ≡ a
このように ay(1 − i) の平方が ≡ a なのだから、ay(1 − i) は a の平方根。

〔注1〕　②の根号内の積は i/2a × (−2i) ≡ i × (−2i)/2a ≡ (−2)i²/2a ≡ (−1)(−1)/a ≡ 1/a

Atkin のアルゴリズムは（1 − i ではなく）符号が反対の i − 1 を使って表現されることが多いようだが、どうせ④で平方されてしまうので、どちらでも大勢に影響ない。

〔追記〕　複素数の世界では i と −i の違いや平方根の主値は、一応、明確に定まっている: ①に対応する式 1 − i = √−2i が成り立つ一方、左辺の符号を逆にすれば右辺の符号も逆になって i − 1 = −√−2i となる。一方 mod p の世界では −1 の2種類の平方根のどちらが +i でどちらが −i なのかは相対的な問題で、a の値に応じて、ある場合の i が別の場合には −i になる。複素数の世界に当てはめると「i と −i の名前が入れ替わることがある」。とはいえ、a の値によらない一般論としては、複素数の世界と同様にイメージした方が分かりやすいだろう。（2023年5月21日）

最後に y ≡ (2a)^k, i ≡ (2a)y² の正しさについて。 p = 8k+5 なので k = (p−5)/8、つまり y ≡ (2a)^(p−5)/8。だから:
　　(2a)y² = (2a)¹[(2a)^(p−5)/8]² ≡ (2a)^4/4(2a)^(p−5)/4 ≡ (2a)^(p−1)/4
Euler の基準から (2a)^(p−1)/2 ≡ −1 なので、その平方根 (2a)^(p−1)/4 は、確かに −1 の平方根 i。
　　i = (2a)^(p−1)/4
この両辺を 2a で割ると:
　　i/2a ≡ (2a)^(p−5)/4　　← (2a)^(p−1)/4 ÷ (2a) ≡ (2a)^(p−1)/4(2a)⁻¹
両辺の平方根を取ると:
　　√i/2a ≡ (2a)^(p−5)/8
この左辺を y とするのだから、確かに y ≡ (2a)^(p−5)/8 ≡ (2a)^k となる。

Atkin のアルゴリズムの動作原理　p を8k+5型素数とする。mod p において…
　　理論　−1 の平方根 i と i/(2a) の平方根 y が得られるなら ±ay(i−1) は平方剰余 a の平方根
　　実装　y ≡ (2a)^k と i ≡ (2a)y² は上の条件を満たす

y は i/(2a) の平方根。 y の符号は任意: i の計算では y は平方され、符号の違いは関係なくなるし、最終的な √a の計算では、y の符号の違いは ±ay(i−1) の複号に吸収される。

2023-05-13　Müller のアルゴリズム　mod 16k+9

#数論　#平方根 mod p

Tonelli のアルゴリズムは mod p での平方根を求めるシンプルで汎用的な方法。素数 p が8k+5型のときには、直接計算することもできるし（§38）、巧妙な Atkin のアルゴリズムを使うこともできる。

Atkin のアルゴリズムは「8k+1型素数」にも（部分的には）拡張可能。この拡張（Siguna Müller のアルゴリズム）は、意外と奥が深い。本題の Tonelli–Shanks から見るとますます脱線になるが、寄り道もまた楽し！

このメモでは p は常に 3 以上の素数を表し、a は mod p での平方剰余を表す（要するに a の平方根が存在すると仮定）。簡潔化のため「平方根を求める」こと（つまり「開平する」こと）を単に「開く」と表現する。

§48.　「8k+1型素数」は、もちろん 8 の倍数より 1 大きい。小さい順に 17, 41, 73, 89, 97, 113, … など。「8k+1型」の p と「8k+5型」の p の一つの大きな違いとして、後者を法（mod）とすると 2 は（従って 2a も）平方非剰余、前者を法とすると 2は（従って 2a も）平方剰余。これは、第二補充法則による。そうすると Euler の基準から:
　　p が8k+1型　⇒　(2a)^(p−1)/2 ≡ +1 (mod p)　両辺を開いて　(2a)^(p−1)/4 ≡ +1 or −1
　　p が8k+5型　⇒　(2a)^(p−1)/2 ≡ −1
後者では (2a) の累乗として自然に −1 を作れるが、前者では +1 or −1 のどっちに転ぶか運任せ。

Atkin のアルゴリズム（p = 8k+5 の場合）では、(p−1)/2 = 4k+2 が偶数なので、上の式を開いて (2a)^(p−1)/4 ≡ (2a)^4k+2 ≡ i を直ちに構成できる（i は −1 の平方根の一つ）。のみならず i/(2a) ≡ (2a)^(p−5)/4 の平方根 y ≡ (2a)^(p−5)/8 ≡ (2a)^k も自然に求まり、この y と 1−i から（または y と i−1 から） √a を計算できるのだった。

p = 8k+1 の場合にも、もしたまたま運よく (2a)^(p−1)/4 ≡ −1 (mod p) なら、(p−1)/4 = 2k は偶数なので、両辺を開いて (2a)^(p−1)/8 ≡ (2a)^k ≡ i を構成できる。この場合、計算上では i/(2a) ≡ (2a)^(p−9)/8 の平方根 y ≡ (2a)^(p−9)/16 も自然に求まる。

〔注〕　i ≡ (2a)^(p−1)/8 なら、両辺を 2a で割ると i/(2a) ≡ (2a)^(p−1)/8 ÷ (2a)¹ ≡ (2a)^{(p−1)/8−1} その指数は (p−1)/8 − 8/8 = (p−9)/8。

ただし y の指数 (p−9)/16 が自然数になるためには、もちろん p−9 が 16 で割り切れることが大前提。上記 8k+1型素数のうち、例えば p = 41 の場合、p−9 = 32 が 16 で割り切れるので問題ないけど、p = 17 の場合 p−9 = 8 が 16 で割り切れないので、上の方法を使えない。

説明の便宜上、文字を k から n に変え、8k+1型の素数を「8n+1」型と表現すると…
　　p = 17 = 8⋅2+1 （n = 2)　や　p = 97 = 8⋅12+1 （n = 12)　のように n が偶数なら、p−9 は 16 で割り切れず、
　　p = 41 = 8⋅5+1 （n = 5)　や　p = 73 = 8⋅9+1 （n = 9)　のように n が奇数なら、p−9 は 16 で割り切れる。

今、文字 k を再利用して n = 2k を偶数とすると（k: とある整数）、p = 8n+1 = 8(2k)+1 = 16k+1 なので p−9 = 16k−8 は 16 で割り切れない。一方、n = 2k+1 を奇数とすると、p = 8n+1 = 8(2k+1)+1 = 16k+9 なので p−9 = 16k は 16 で割り切れる。要するに「8n+1」型の素数を法として Atkin のアルゴリズムと同様のことができるとしたら、それは n が奇数の場合、要するに素数 p が「16k+9」型の場合に限られる――「8n+1」型素数は、n の偶奇により「16k+1」型と「16k+9」型に細分されるが、以下では法が「16k+9」型の素数の場合を扱う。

説明を後回しにして、全体像を一覧表にしてみる。

＜表2＞　Atkin のアルゴリズムとその拡張: 平方剰余 a の平方根

	p = 8k+5	p = 16k+9
2a は…	平方非剰余	平方剰余
		四乗非剰余	四乗剰余
Euler	(2a)4k+2 ≡ −1	(2a)8k+4 ≡ +1
		(2a)4k+2 ≡ −1	(2a)4k+2 ≡ +1
i	(2a)2k+1 ≡ (2a)y2		(2a)2k+1⋅z4k+2 ≡ (2a)y2
i の供給	無料（内部の非剰余から）		外部の非剰余 z から
y	(2a)k ≡ √i/2a		(2a)k⋅z2k+1 ≡ √i/2a
y2	i/(2a)
y2(1 − i)2	1/a　なぜなら (1 − i)2 = −2i
√a	ay(1 − i)　なぜならその平方 = a2/a

左の列の p = 8k+5 の場合が、通常の Atkin のアルゴリズム。中央の列の「2a が四乗非剰余」のケース（その意味は後述）も、通常の Atkin のアルゴリズムと実質同様に扱える。どちらの場合も i ≡ (2a)^奇数 の形にできれば、y ≡ (2a)^整数 となり（y は i/2a ≡ (2a)^偶数 の平方根）、そこから √a が求まる。そのためには、単に y ≡ (2a)^k, i ≡ (2a)y² とすればいい―― k は p の形から明白、a は与えられた数なので、全ては機械的な単純計算。

他方、右の列の「2a が四乗剰余の場合」は比較的複雑で、Tonelli のアルゴリズム同様、mod p の平方非剰余 z を外部から持ってくる必要がある（次回）。

§49.　具体例で試してみよう！

16k+9型素数の例として p = 89 を考える（k = 5）。これは「しげちー」のワープ・スペシャル素数: mod 89 では √3 も √7 も存在しないのに、それらの積のような √21 は存在する（§32）。けれど「存在する」というだけで、具体的な値は未知。実は a = 21 に対して…
　　(2a)^4k+2 ≡ 42²² ≡ −1 (mod p)　　【ヤ】
…なので、＜表2＞中央の「2a は四乗非剰余」というケースに当たり（表の「Euler」の欄の (2a)^4k+2 ≡ −1 に当たる）、普通の Atkin のアルゴリズムと全く同様に平方根を求めることができる。

といっても、第一印象として「そもそも四乗非剰余って何だよ」「だいたい 42²² って、42 を 22 個掛けることだろ。そんな超でかい数、面倒くさくて計算したくねーよ！」といったご不満もおありでしょう。

まず【ヤ】の計算だが、表の「四乗非剰余」（←意味はさておき）の欄を見ると i の計算法は分かっている: i ≡ (2a)y² だ。 a = 21 は与えられた条件なので 2a = 42。後は y さえ求めれば i が求まる（そして i は −1 の平方根なのだから、i² ≡ −1 を確認すれば【ヤ】を確かめたのと同じことになる）。必要な y の計算法も、表にある通り: y ≡ (2a)^k ≡ 42⁵ だ（この例では k = 5 なんで）。y と i についてのこれらの計算は、Atkin のアルゴリズムと全く同じ。

【1】　てなわけで、最初の一歩は y ≡ 42⁵ (mod 89) の計算。「22乗よりはマシだけど、42の5乗なんて面倒くさくて計算したくねーな」と感じる人もいるだろう。確かに、ばか正直に
　　42 × 42 × 42 × 42 × 42 を 89 で割った余り
を求めようとしたら、そりゃ面倒かも…。電卓などで計算するという方法も、もちろんある。けど、こつが分かれば、実は手計算（場合によっては暗算）でも大したことではない。このような指数計算をどーやってやるのか、以下では、その秘密の一部始終をお見せ致しましょう（なんか偉そうｗ）。「べき剰余」の計算の基本は、まず 2 乗、そのまた 2 乗（つまり最初の 4 乗）、そのまた 2 乗（最初の 8 乗）…等々。この例では、コンセプト的には、
　　42² ≡ A (mod 89), 42⁴ ≡ (42²)² ≡ A² ≡ B,
　　42⁵ ≡ 42 × 42⁴ = 42B
のように進めることができる。初めの一歩は 42² を 89 で割った余り A。
　　42² = (2⋅21)² = 2²⋅21² = 4⋅441
　　――ここで気を利かせて 89 × 5 = 445 を使うと 441 ≡ −4 (mod 89) なので:
　　42² ≡ 4⋅(−4) ≡ −16
　　だから　42⁴ ≡ (−16)² = 256
　　――ここで気を利かせて 89 × 3 = 267 を使うと 256 ≡ −11 (mod 89) なので:
　　42⁵ = 42⋅42⁴ ≡ 42⋅(−11) = −462　☆
　　――ここで気を利かせて 89 × (−5) = −445 を使うと −462 ≡ −17 (mod 89)
それが 42⁵ なので y ≡ 42⁵ ≡ −17 となる。「気を利かせて」というのは、要するに、その数に近い 89 の倍数を足すか引くかして、およそ −89/2 ～ +89/2 の範囲に「簡約」すること。そうすれば、その先は絶対値 45 未満の小さい数の計算で済む。概念的には 89 で割った余りだが、例えば 441 割る 89 について、5 余り −4 と考えている。素朴に 4 余り 85 としてもいいんだけど、85 より −4 の方が扱いやすいんで臨機応変に。商はどうでもよく、p で割った余りだけが重要。余りが p/2 以上の自然数 r のとき、商を 1 増やして余りは、負の数 r−p だとして構わない（商はどうでもいいので、好きに加減できる）。

【2】　次に i ≡ (2a)y² を求めよう。まず y² = (−17)² = 289 ≡ 22
　　［89 × 3 = 267 は上で計算済みなので、289 を 89 で割れば 289−267 = 22 余る］
だから i ≡ (2a)y ≡ 42⋅22。これを今までと同様、普通に計算してもいいのだが、上記☆の成り行きから 42⋅11 ≡ 17。で 42⋅22 はその2倍なので ≡ 34。場合によっては、再利用できる計算を再利用してもいいだろう（別にしなくてもいいけど）。

【3】　最後に i² ≡ −1 を確かめる:
　　34² ≡ (2⋅17)² ≡ 2²⋅17² ≡ 4⋅289 ≡ 4⋅22 = 88 ≡ −1 (mod 89)
　　［上記の通り 289 を 89 で割れば 289−267 = 22 余る］

このように順を追って進めれば、【ヤ】は、たわいもない。コンピューターを使って強引に
　　42²² = 5146 1730 8132 8524 0070 0537 6493 5345 7664
を計算しても構わないけど、いくらコンピューターでも1000乗、2000乗、あるいは100万乗をばか正直に計算するのは非効率。2乗の2乗のそのまた2乗の…という計算法なら、たとえ100万乗でも 2²⁰ 乗くらいなので、素早く決着がつく。

【4】　これで y ≡ −17 と i ≡ 34 が求まったので:
　　ay ≡ 21(−17) = −357 ≡ −1　　［← −357 は、89 の倍数 89×(−4) = −356 より 1 小さいので］
　　ゆえに　√21 ≡ ay(i − 1) ≡ (−1)(34 − 1) ≡ −33　符号は逆でも構わない

【5】　検算　(±33)² = 3²⋅11² ≡ 9⋅121 ≡ 9⋅32^* ≡ 288 ≡ 21^**
　　　［*121 を 89 で割ると 32 余る　　**288 を 89 で割ると 21 余る（89×3 = 267）］
(±33) の 2 乗は確かに ≡ 21。つまり mod 89 での 21 の平方根は ±33 ということが分かった。平方根の符号が ± どっちでもいいことは分かり切ってるので、複号を省略して √21 ≡ 33 としておこう。

「しげちー」の文脈では、(√21 − 3)/6 ≡ (33 − 3)/6 ≡ 5 はスペシャルデー。すなわち――
　　Q(5) = 1⁴ + 2⁴ + 3⁴ + 4⁴ + 5⁴ = 979
――は、p = 89 で割り切れる。 5 は11日周期の基本の日でもあるので、Q(5) は 11 の倍数。以上から Q(5) = 979 = 89 × 11 は完全に説明がつく。

〔付記〕　5個の数の4乗和 1⁴ + 2⁴ + … + 5⁴ が 89 で割り切れるということは、一見意味のない「ランダム」な現象のようだが、mod 89 での √21 という、それなりに高級な理論が背景にある。われわれは法 89 が「16k+9」型素数であることに着目して、21 の平方根を――Tonelli–Shanks のような標準的アルゴリズムではなく―― Atkin のアルゴリズムの拡張版という奇手によって求めた。「n 日目に n⁴ 枚ずつ1円玉を集めると、合計いくらになるか？」というくだらない問題との関連で言えば、結構なアクロバットだ！

§50.　四乗剰余・四乗非剰余という言葉の意味は…。前述のように「16k+9」型素数 p は「8n+1」型素数の一種（n が奇数の場合）なので、第二補充法則により、mod p では 2 は平方剰余。もし a も平方剰余なら、当然 2a も平方剰余: Euler の基準により (2a)^(p−1)/2 ≡ +1 となる。両辺を開くと:
　　(2a)^(p−1)/4 ≡ +1 or −1　　【ユ】
その右辺が +1 か −1 かに応じて 2a を四乗剰余・四乗非剰余と呼ぶ。表面的にはそれだけのことだけど、その真意は以下の通り。

2a が平方剰余なら、当然 r² ≡ 2a (mod p) を満たす r が存在する。では、この r という数（2a の平方根である）それ自体は、平方剰余だろうか、非剰余だろうか。つまり x² ≡ r を満たす x があるか、ないか？

r² ≡ 2a に留意すると、もし【ユ】の右辺が +1 なら:
　　(r²)^(p−1)/4 ≡ +1
　　つまり　(r)^2(p−1)/4 = r^(p−1)/2 ≡ +1
従って、この場合 Euler の基準から r 自体も平方剰余であり、x² ≡ r を満たす r が存在する。このとき、r² ≡ 2a なので (x²)² ≡ 2a となって、つまり x⁴ ≡ 2a (mod p) を満たす x が存在する。言葉で言うと「四乗して p で割ると、剰余（余り）が 2a になるような x」が存在する。このことを「2a は四乗剰余」と表現する（その場合、2a の四乗根は x である）。「2a ってのは、何かを四乗して p で割った余りだよ～ん」という意味。四乗した余りだから、四乗剰余。

逆にもし【ユ】の左辺が −1 なら:
　　(r²)^(p−1)/4 ≡ −1
　　つまり　(r)^2(p−1)/4 = r^(p−1)/2 ≡ −1
となって、Euler の基準から r は非剰余: x² ≡ r を満たす r は存在しない。従って x⁴ ≡ 2a を満たす x も存在しない。言葉で言うと「四乗して p で割ると余りが 2a になるような x は、存在しない」。このことを「2a は四乗非剰余」と表現する（その場合、2a の四乗根は存在しない）。「どんな数を四乗した余りにも、ならないよ～ん」ってことなんで、四乗非剰余。

法 p が16k+9型の素数のとき、2 は（あるいは 2a は）平方剰余。だからその平方根 r が存在するけど、「r のそのまた平方根」は、存在するかもしれないし、しないかもしれない。最初の数から見ると、平方根の平方根（つまり四乗根）があるか、ないか。これが四乗剰余 vs. 四乗非剰余の違いだ！

〔例〕　16k+9型の素数 p = 41 を mod とする。2 は平方剰余なので、当然 4 は四乗剰余（4 の平方根はもちろん 2 だが、その 2 が再び平方剰余＝平方根を持つ）。一方 8 = 2 × 4 は平方剰余と平方剰余の積で平方剰余だが、その平方根 7 は（7² = 49 ≡ 8）、平方非剰余。この事実については、(±1)², (±2)², …, (±20)² を一つずつ計算してもどれも ≡ 7 にならないことから、直接確かめることもできるし、Euler の基準 7²⁰ ≡ −1 から確かめることもできる。だが最も速いのは「第七補充法則」を使うことだろう: 奇素数 p = 41 は ±1, ±3, ±9 (mod 28) のどれでもないから、直ちに 7 は mod p の平方非剰余と結論される。まとめると: この例では 4 と 8 はどちらも平方剰余だが、そのうち 4 は四乗剰余（平方根が再び平方剰余になる）、8 は四乗非剰余（平方根が再び平方剰余にならない）。

ある数が四乗剰余かどうか？というのは、その数の平方根 r が平方剰余か非剰余か？という問題だから、Euler の基準を使って――つまり r^(p−1)/2 ≡ ±1 の符号によって――判定可能。われわれのアルゴリズムでは (r²)^(p−1)/2 ≡ (2a)^(p−1)/2 すなわち r^p−1 ≡ (2a)^(p−1)/2 なので、その両辺を開くと:
　　r^(p−1)/2 ≡ (2a)^(p−1)/4
r についての Euler の基準による判定は、(2a)^(p−1)/4 ≡ ±1 の符号の違いによる判定に当たる。p−1 = 16k+8 なので、指数 (p−1)/4 は 4k+2 に等しい。＜表2＞を再び見ていただくと、中央の列・右側の列において、四乗非剰余・四乗剰余の区別が Euler の欄の (2a)^4k+2 = (2a)^(p−1)/4 の符号に対応してることが分かるだろう。

さらに＜表2＞の左側の列と、中央の列を見比べると…。8k+5型素数に対する Atkin の平方根アルゴリズムは、2a が四乗非剰余なら、ほとんどそのまま16k+9型素数にも適用できる。問題は、16k+9型の素数を法として 2a が四乗剰余になる場合（右端の列）。それについては、次回に検討する予定。

上記の計算法は、あまり汎用性がない。法が16k+9型の素数というのは、ひどく限定的だし（最初の三つの16k+9型素数は 41, 73, 89――この型の素数は 0～500 の間に 9 個しかない）、Atkin 風の簡単な方法が使えるのは、2a が四乗非剰余の場合だけ。そんなこともあって、ほとんど紹介されることのないニッチなアルゴリズム。けれど、内容的には興味深いものを秘めている。

2023-05-17　Müller のアルゴリズムの簡単化に成功（速報）

#数論　#平方根 mod p　#独自研究

ルーマニアの Rotaru & Iftene (2013) によると、Siguna Müller のアルゴリズムは、他の複数の研究者によって改善された。われわれは当初 Rotaru の記述に沿って改善版を紹介し、まず「簡単な場合」を検討したが、後回しにした「面倒な場合」について、3日後の2023年5月16日、簡単化に成功した。

第一のポイントは 2a が四乗剰余のときの y の定義。Rotaru が使っている値（Müller のオリジナルの設定なのだろう）の z 倍を y とした方が、計算量を節約でき、式もすっきり（Atkin の公式と全く同じ形式になる）。第二のポイントとして、2a が四乗剰余（面倒なケース）だとしても、約25%の確率で、「無料」の（平方非剰余 z を必要としない）直接計算が成り立つ――理論的に興味深いことだが、a が八乗剰余なら、mod 16k+9 での平方根は無料。単に a^k+1 を計算すればいい。

小さな改善に過ぎないが、なかなか面白い！

以下のメモの内容については「Müller のアルゴリズムの簡単化」参照。

We have simplified and streamlined Siguna Müller’s algorithm (improved by Kong et al. as quoted by Rotaru & Iftene^*) for calculating a square root of a quadratic residue a mod p, where p is a prime of the form 16k+9. Rotaru’s Lines 8–10 have been cleaned up. [Rotaru’s version requires 1 exponentiation, 2 squares, 6 multiplications, whereas our version only requires 1 exponentiation, 1 square, 3 multiplications, and is more transparent.] Line 5a, added by us, shows an optional shortcut: even when we don’t have a square root of −1, if some additional calculation (inexpensive) is acceptable and if we’re lucky (~25%), then we don’t need an auxiliary quadratic non-residue, and can get a square root of a immediately. This shortcut becomes deterministic and “free” once the value 2^2k+1 is known.

* https://profs.info.uaic.ro/~siftene/fi125(1)04.pdf (see “Figure 2”)

b = Square Root of a, modulo p = 16k+9, by Siguna Müller (2004)
as quoted by Rotaru & Iftene (2013) + our simplification (May 16/17, 2023)

	Müller–Kong	Simplified	Comments
Input: a, a quadratic residue modulo p, where p is a prime of the form 16k+9, i.e. k = (p−9)/16 Output: b, a squre root of a
1	y := (2a)k		y is a square root of i/(2a), where i is a fourth root of 1; Rotaru’s α [alpha]
2	i := (2a)y2		Rotaru’s β. It’s convenient if i is a square root of −1, but it may be a squre root of 1
3	if i2 = −1	if i ≠ ±1	trivially reduce 1 square
4	then b := ay(1−i)		Atkin’s “free square root”: from y2 = i/(2a) and (1−i)2 = −2i, we have y2(1−i)2 = 1/a, so (ay(1−i))2 = a2/a = a. One may use (1−i) or (i−1), whichever is fine, as this is to be squared
5	else		// i = 1 or −1. From Lines 1 and 2, we see that i = (2a)2k+1, so here +1 = i2 = (2a)4k+2; hence 24k+2 and a4k+2 are both +1 or both −1 [having the same biquadratic character, because 2 and a are both quadratic residues]. We may be lucky if 24k+2 = a4k+2 = +1 [i.e. if both 2 and a are biquadratic residues]
5a		(optionally) if i = 22k+1 then b := ak+1 else	if we are lucky (~25%), we don’t need a quadratic non-residue z (Line 7) even when i = ±1. That is, if i = (2a)2k+1 is = 22k+1, then (by dividing both by 22k+1) we have a2k+1 = 1 [i.e. a is an octic residue modulo p], so a2k+2 = a, i.e. ak+1 = b. One could calculate b = ak+1 without knowing the value 22k+1 and see if b2 = a holds. This shortcut becomes deterministic and “free” once the value 22k+1 is known; it works ~50% if 22k+1 = +1 or −1, depending on the current value of i, and never works (at least in this context*1) if 22k+1 is a square root of −1 [i.e. 2 must be a biquadratic residue (24k+2 = +1) for this to work]
6	begin
7	generate z, a quadratic non-residue mod p		Rotaru’s d; given that p ≡ 1 (mod 4), one can use z = 2 if p ≡ 5 (mod 8), z = 3 if p = 2 (mod 3) or z = 5 if p ≢ ±1 (mod 10); z is harder to get if p ≡ 1 or 49 (mod 60)(mod 120)
8	y := yz2k	y := yz2k+1	here our y = (2a)kz2k+1 is z times bigger than Rotaru’s; this works better (see Lines 9, 10). [from Lines 1–3, the old y = (2a)k is a square root of ±1/(2a); since z2k+1 is a square root of i (where this new i is a square root of −1), the square of the new y is i/(2a)]
9	i := 2ay2z2	i := (2a)y2	by redefining y as in Line 8, we reduce 1 square and 1 multiplication here [i = (2a)2k+1⋅z4k+2, where (2a)2k+1 — the old i — is ±1 and z4k+2 is a square root of −1]*2,
10	b := ay(1−i)z	b := ay(1−i)	and we reduce 1 multiplication here. Moreover, our reasoning is transparent, always using the same formula of Atkin’s (Line 4), as the new i is a square root of −1 and the new y is a square root of i/(2a)
11	end		Lines 9, 10 are now identical to Lines 2, 4, no longer using z
12	return b		so not only faster, this is cleaner and simpler

*1 In fact, whatever the value 2^2k+1 is (±1 or not), if (2a)^2k+1 = 2^2k+1, then a^2k+1 = 1 and b = a^k+1. While this is nice, an Atkin-like method (Line 4) is faster if we already know the suitable i, a primitive fourth root of unity [May 19, 2023]. On the other hand, provided that we do know a primitive fourth root of unity already, a similar shortcut always works in Line 5 if a^4k+2 = +1 holds (~50%). These shortcuts are deterministic and “free” once we know the value 2^2k+1 and it is +1 or −1 i.e. if 2 is a biquadratic residue modulo p [May 25]. Moreover, when 2 is a biquadratic residue, a quadratic non-residue is never necessary to find √a. Indeed, in such a case we can immediately calculate √a itself or √−a. Thus knowing √−1 is enough, for which we only need a biquadratic non-residue ζ. It is sufficient but not necessary that ζ is a quadratic non-residue (ζ^8k+4 = −1); even if ζ^8k+4 = +1, it works if ζ^4k+2 = −1. How to extract √−1 from ζ^even ≡ −1 is obvious [May 30].

*2 We can reduce one more multiplication by setting h := z^2k+1, y := yh, i := ih² [May 19]. Let Y be the old y, ε be the old i. Y = √±1/2a, where the sign (+ or −) is the same as that of ε = ±1; so the new y (Line 8) is √±1/2a⋅z^2k+1 = Yh, from which the new i (Line 9) is 2a(±1/2a)(z^2k+1)² = ±z^4k+2 = εh² [May 22, 2023].

2023-05-23　Müller のアルゴリズムの簡単化　2023年版

#数論　#平方根 mod p　#独自研究

Siguna Müller の平方根アルゴリズムの簡単化に成功した。Rotaru & Iftene 版（2013年）では1回のべき・2回の平方・6回の掛け算が必要だが、われわれは1回のべき・1回の平方・3回の掛け算で同じ結果に至り、内容も透明になった。世界最先端（？）の最新成果を分かりやすく解説する。

ところで、通常は平方非剰余を必要とするケースでも、4回に1回の割合で「無料」の（平方非剰余を必要としない）平方根計算が成り立つ――四乗剰余・八乗剰余と関係していて、検討してみると面白そうだ。全体としては初等的でさまつなことだけど、これを入り口に視野が広がるかも…。

§51.　p を 3 以上の素数とする。x² ≡ a (mod p) が解を持つとき、a を平方剰余という。以下 a を平方剰余（ただし a ≢ 0）とする。平方すると ≡ a になる数を a の平方根と呼ぶ; b が平方根なら −b も平方根で b と −b は不合同（異なる種類の数）。

p が4k+3型のとき a の平方根を求めることは易しい; 1760年代、既に Lagrange [17] によって記述されている。

[17] https://fr.wikisource.org/wiki/Page:Joseph_Louis_de_Lagrange_-_%C5%92uvres,_Tome_2.djvu/501

p が4k+1型の場合は、それに比べ大変難しい。厳密に言うと、21世紀の現在でも、一般的・決定論的な多項式時間アルゴリズムが得られていない。Gauß (1777–1855) も、効率的な解法を発見できなかった（現在でもそうなのだから当然だが）。実用性の高い最初の確率的アルゴリズムを発見したのは、イタリアの Tonelli で、1891年。現在でも広く使われているのは、Tonelli の方法と本質的に同じ Shanks のアルゴリズムであり、しばしば二人の名前を併記して Tonelli–Shanks と呼ばれる。Tonelli のアルゴリズムについては既に記した。Shanks の方法についても、別系列の Cipolla のアルゴリズムについてもそのうち記す予定だが、どちらも有名なので、検索すればいろいろウェブ上で読むこともできるだろう。これらは汎用的なアプローチであり、以下の限定的ケースより実用性が高い。

いろいろな角度から問題を眺め、理解を深め、可能性を探るため、汎用的アルゴリズムに甘んじず、あえて限定的ケースを検討するのである。「p が8k+5型や16k+9型のとき、以下のように進める必要がある」という意味ではない。

4k+1型のうちでも、8k+5型に限っては、比較的簡単な直接的解法がある。Legendre も、それを知っていたらしい。1917年、Pocklington は簡潔な表現を記述（§38）。さらに、20世紀の Atkin は、別の研究の副産物として、8k+5型のケースについて巧妙な解法を発見。オーストリアの Siguna Müller（スィグゥナ・ムュラー）は、北米滞在中に Atkin の方法を16k+9型素数にまで拡張し、2004年に論文発表。この Müller のアルゴリズムは Kong らによって2006年に少し改良され、2013年にルーマニアの Rotaru & Iftene（以下 Rotaru と略）によって一般化されている。

*Siguna Müller (Mueller) はオーストリアの数学者。離散数学・暗号学の分野では、オーストリアで最初に教員資格を得た女性の一人だという。本国の Klagenfurt 大学で学位を得た後、カナダの Calgary 大学の研究員を経て、米国で生物医学の博士号を取得。2020年代の現在では、Covid-19 の原因となるウイルスに関連する研究でも知られる。数論出身だけあって、その論考は透明・厳密。「ワクチンが効果的だからこそ、耐性のある変異株がまん延する」といった基礎的な話から始まって、「分からないことだらけ」「効果もあるがADE（有害事象）もある」「有効性それ自体が別の文脈では危険性」といったことを遠慮もなく、誇張もなく、いわば数学的な明晰さで論じている。Muller’s method の David E. Muller とは別人。

2023年5月、Müller のアルゴリズムにはさらに改良の余地があることが分かった。以下のメモの内容は（大げさに言うと）最先端の成果で、まだ世界のどこでも発表されていない。あまりに限定的過ぎて、ほとんど意味がないので、ほとんど誰も研究してないってだけだけど…。

Müller のアイデアのうち、四乗非剰余の場合は Kong らによって改良されたし、以下ではそれ以外（四乗剰余の場合）についての改良を述べるので、結果的にもともとのアイデアは全部書き換えられてしまうわけだが…そのような後続の研究を刺激することそれ自体、功績といえるだろう。後から来て「こうした方が良い」と言うのは簡単だが、そのたたき台を作った先駆者が一番偉い。

まず Atkin のアイデアを振り返ると――ただし、ここでは i−1 の代わりに 1−i を使う――、p = 8k+5 のとき mod p において「−1 の平方根 i を経由して、a の平方根が分かる」。実際 (1 − i)² ≡ −2i ということ（イメージと解説）から、y² ≡ i/(2a) とすると、y²(1 − i)² ≡ i/(2a) × −2i ≡ 1/a。その両辺を a² 倍すれば:
　　[ay(1 − i)]² ≡ a　つまり　ay(1 − i) ≡ ±√a
このように、簡単に √a に到達可能。見掛け上 x² ≡ a を別の2次方程式 y² ≡ i/(2a) に「責任転嫁」したみたいだけど、Euler の基準から
　　(2a)^(p−1)/2 ≡ (2a)^4k+2 ≡ −1　ゆえに　(2a)^2k+1 ≡ i
であるから、y² ≡ i/(2a) ≡ (2a)^2k の一つの解が y ≡ (2a)^k であることは明白。y の式は簡単に解ける（逆に言うと、この y や i は与えられた a に応じて自動的に定まるものであり、この文脈では、−1 の任意の平方根を i とすることは許されない）。

§52.　p = 16k+9 の場合には (2a)^(p−1)/2 ≡ (2a)^8k+4 ≡ +1 だが（符号がプラスなのは 2 が平方剰余だから）、これは (2a)^4k+2 ≡ ±1 を含意する。複号のうちマイナスが選択されるなら、(2a)^4k+2 ≡ −1 だから、上記 Atkin の方法をそのまま流用できる（第１の場合）。では (2a)^4k+2 ≡ +1 なら、どうするか（第2の場合）――それが Müller のアイデアの核心だろう。

−1 の平方根を利用するのなら、結局のところ、何かの偶数乗 ≡ −1 の両辺を開くことになりそうだ。Euler の基準に照らせば、任意に平方非剰余 z を一つ選択して…
　　z^(p−1)/2 ≡ z^8k+4 ≡ −1
…を −1 の平方根 z^4k+2 の供給源とするのが順当だろう:
　　(2a)^4k+2 ≡ +1　なら　(2a)^2k+1 ≡ ±1
なので、複号がどっちになるにしても、(2a)^2k+1⋅z^4k+2 を i とすればいい。

s ≡ z^4k+2 を −1 の一つの平方根とした場合、もし (2a)^2k+1 ≡ +1 なら上記の i は s 自身だし、(2a)^2k+1 ≡ −1 なら上記の i は −s に当たる。

Atkin 風に進めるには y² ≡ i/(2a) ≡ (2a)^2k⋅z^4k+2 を y について解けばよく、その一つの解はもちろん y ≡ (2a)^k⋅z^2k+1。後は Atkin の方法と全く同様。

以上の方針を基に、(2a)^4k+2 ≡ −1 の場合と (2a)^4k+2 ≡ +1 の場合を一本化して説明する。まずとにかく y ≡ (2a)^k を計算し、y² ≡ i/(2a) なのだから、両辺を (2a) 倍して i ≡ (2a)y² と置く。前節で見たように y² ≡ (2a)^2k だから i ≡ (2a)^2k+1 となるが、i という名前のイメージにかかわらず、この i は −1 の平方根であるとは限らないし、「符号がどうなるか」も a によって異なる――その意味は次の通り。Euler の基準から (2a)^8k+4 ≡ +1 なので、(2a)^4k+2 は +1 の平方根、i ≡ (2a)^2k+1 はそのまた平方根――要するに +1 の四乗根。今 −1 の平方根を ±s とすると、
　　(±s)⁴ ≡ [(±s)²]² ≡ [−1]² ≡ +1
なので ±s も +1 の四乗根には違いないが、
　　(±1)⁴ ≡ [(±1)²]² ≡ [+1]² ≡ +1
も事実なので ±1 も +1 の四乗根だ。このとき (±s)² ≡ −1 だが (±1)² ≡ +1 なので s と 1 は違う種類の数。要するに x⁴ ≡ +1 (mod p) には x ≡ ±1, ±s の4種類の解があり、これは4次方程式なので、それ以外の（4種類を超える）解はない。便宜上 i と置いた数はこの4種類のどれかで、どれになるかは、やってみないと分からない。

もしこの i が +1 でも −1 でもないなら、その i は確実に −1 の平方根なので、直ちに Atkin の方法を適用でき √a ≡ ay(1 − i) となる（もちろん根号の前の符号は逆でもいい）。しかるに、もしこの i が ±1 のどちらかになってしまった場合、そのままでは Atkin の方法を適用できない。というのも、Atkin の方法は、
　　(1 − i)² ≡ −2i　　（☆）
という性質を利用している。この式に i ≡ +1 を入れれば 0 ≡ −2 になるし、i ≡ −1 を入れれば 4 ≡ 2 になるが、どちらも明らかに間違っている。（☆）は、i が −1 の平方根だからこそ成り立つ式であり、i ≡ ±1 ではこの土台が崩れてしまう。そこでこのケースでは、Tonelli のアルゴリズムの正反対のような補正を行う: もし i ≡ ±1 になったら、適当な平方非剰余 z を使って y ≡ (2a)^k⋅z^2k+1 と置き、あらためて i ≡ (2a)y² とすれば、前述の理由から、今度の i は −1 の平方根――この新しい y, i を使えば、Atkin の公式から √a ≡ ay(1 − i) とできる。実用上、この場合には最初に y ≡ (2a)^k を計算してあるので、この古い y の z^2k+1 倍をあらためて y として、その先は前半同様に進めればいい。

これで (2a)^4k+2 ≡ −1 の場合と (2a)^4k+2 ≡ +1 の場合が一本化され、その仕組みも透明になった。

〔注〕　以上はオリジナルの Müller のアルゴリズム（の Rotaru による引用）ではなく、簡単化したもの。もし内容におかしな点があったら、それはアレンジの問題で、Müller のせいではない。本質的な発見者は Atkin と Müller。

次の節では、これら2種類のケースについて、具体例を使って確かめてみたい。

§53.　16k+9型の最小の素数は p = 41 なので（k = 2 に当たる）、この節では mod 41 の例を考える。計算の流れについては前節を、途中計算については下の〔解説〕を参照。

例題1　a ≡ 10 として √10 を求める。まず y ≡ (2a)^k ≡ 20² ≡ 400 = 410 − 10 ≡ −10。次に i ≡ (2a)y² ≡ 20 × (−10)² = 20 × (82 + 18) ≡ 20 × 18 = 360 = 369 − 9 ≡ −9。この数は ≡ ±1 ではないので −1 の平方根（実際、その平方 81 は 41 の倍数 82 から見て −1）。従って、直ちに Atkin の公式を使える:
　　ay(1 − i) ≡ 10 × (−10) × (1 − (−9)) = −100 × 10 ≡ (−18) × 10 = −180 = −164 − 16 ≡ −16
符号はどちらでもいいので √10 ≡ +16 or −16 となる。

〔解説〕　途中計算では、410 = 41 × 10, 82 = 41 × 2, 369 = 41 × 9, 164 = 41 × 4 を使い 100 ≡ 18 を −100 ≡ −16 として再利用した。要するに、途中計算で（絶対値の）大きい数 B が現れたら、p（ここでは 41）の倍数 C を適宜加減して、B を（絶対値の）小さい数 B + C ないし B − C に置き換える。例えば B − C = D が小さい数なら B ≡ D なので、B を D に置き換えていい（例題2以降では、この種の説明を省く）。

〔検算〕　16² = 256 = 246 + 10 ≡ 10。ただし 246 = 41 × 6。

例題2　a ≡ 8 として √8 を求める。まず y ≡ (2a)^k ≡ 16² ≡ 256 ≡ 10。次に i ≡ (2a)y² ≡ 16 × 10² = 16 × (82 + 18) ≡ 16 × 18 = 288 = 287 + 1 ≡ 1。この数は ≡ ±1 なので −1 の平方根ではない。そこで平方非剰余 z ≡ 3 を使い（p = 41 は 12 の倍数と 1 違いでないので、3 は平方非剰余）、古い y を z^2k+1 倍つまり 3⁵ 倍することにしよう…。

3⁴ = 81 ≡ −1 なので 3⁵ ≡ −3。従って新しい y は y ≡ 10 × (−3) = −30 ≡ 11。あらためて i ≡ (2a)y² ≡ 16 × 11² = 16 × 121 ≡ 16 × (−2) = −32 ≡ 9。注意点として、例題1では i ≡ −9 だったが、ここでは i ≡ +9 となっている。i は −1 の平方根なら何でもいいのではなく、a に対応した正しい符号を選ぶ必要がある。例題1から i ≡ ±9 なので、その気になれば「符号だけ決定する手抜き」もできそうだが…。ともあれ、Atkin の公式を使える体勢が整った:
　　ay(1 − i) ≡ 8 × 11 × (1 − 9) = 88 × (−8) ≡ 6 × (−8) = −48 ≡ −7
符号はどちらでもいいので √8 ≡ +7 or −7 となる。

〔検算〕　7² = 49 ≡ 8。

例題3　a ≡ 21 ≡ −20 として √21 ≡ √−20 を求める。少しでも労力をケチるため、絶対値が小さい a ≡ −20 を使い、さらに 2a ≡ −40 ≡ +1 という関係を活用しよう。まず y ≡ (2a)^k ≡ (+1)² = 1。次に i ≡ (2a)y² ≡ (+1) × 1² = 1。この数は ≡ ±1 なので −1 の平方根ではない。例題2と同様、平方非剰余 z ≡ 3 を使い、古い y を 3⁵ 倍つまり −3 倍しよう…。

新しい y は 1 × (−3) ≡ −3。あらためて i ≡ (2a)y² ≡ (+1) × (−3)² = 9。Atkin の公式から:
　　ay(1 − i) ≡ (−20) × (−3) × (1 − 9) = 60 × (−8) ≡ 19 × (−8) = −152 = −164 + 12 ≡ 12
符号はどちらでもいいので √21 ≡ +12 or −12 となる。全く同じ平方根について、§41では Tonelli 風に求めたが、少なくともこれに関しては Müller 風の方がシンプルだ！

§54.　その気になれば、さらにショートカットができる。y ≡ (2a)^k が ±1 になってしまった場合について、その古い y（便宜上 Y とする）を使って、上では新しい y を Yz^2k+1 として i ≡ (2a)y² を求め、最終的に答えを ay(1 − i) とした。この方法だと z のべきが1回、y の平方が1回、（2a という積は無料としても）それ以外の掛け算が 4 回必要。ところが y ≡ Yz^2k+1 を i の式に代入すると:
　　i ≡ (2a)(Yz^2k+1)² ≡ (2a)Y² × (z^2k+1)²
このうち (2a)Y² の部分は古い i に他ならない。その値（ε とする）は ≡ +1 or −1 に過ぎないが、とにかく計算済みなので、符号ビットとして再利用可能。次の一手速い手順が成立する:
　　h ≡ z^2k+1　　← 平方非剰余の 2k+1 乗を h とする
　　y ≡ Yh
　　i ≡ εh²
　　答え ≡ ay(1 − i)
z のべきが1回、h の平方が1回、それ以外の掛け算が3回で済む（ε はただの符号なので無料）！

例えば、例題2では h ≡ −3, y ≡ 10 × (−3) = −30 ≡ 11, i ≡ +(−3)² ≡ 9。

既存の Müller のアルゴリズムでは（少なくとも Rotaru による引用を見る限りでは）、この最適化がないどころか、次のような遠回りをしている:
　　y ≡ Yz^2k
　　i ≡ (2a)y² × z²
　　答え ≡ ay(1 − i) × z
われわれの方法と比べ y に掛け算される z の指数が 1 小さいので――そこだけ見る限りでは 2k+1 乗より 2k 乗の方がシンプルかもしれないが――、i の計算のとき y² では足りず z² を追加で掛け算する必要が生じ、答えを求めるときにも z を追加で掛け算する必要が生じる。計算量的に無駄なばかりか、Atkin の公式を再利用できず、従って y² ≡ i/(2a) という統一的な理解もできず、内容が不透明。なぜわざわざこんな手順にしたのか不明（paywall があって Müller の原論文を自由に参照できない）。16k+9型に限定せずさらなる一般化まで考えた場合、もしかすると、この遠回りには意味があるのかもしれない。

次回は、平方根計算が「無料」になるケースについて検討したい。（続く）

2023-05-27　Müller のアルゴリズムと無料平方根　25%のラッキー

#数論　#平方根 mod p　#独自研究

1990年代、Atkin は8k+5型素数 p について、平方根 mod p の「無料」計算法を示した。「無料」とは「他の処理（強擬素数テストなど）の副産物として、ついでに求まる」という意味だが、平方根計算を主目的として Atkin のアルゴリズムを使うこともできる。

Müller はそれを16k+9型の素数 p に拡張した。mod p の平方剰余 a について、第1の場合（a が四乗非剰余）には「無料」の計算が成り立つ。この場合の「無料」は「Atkin 風」――具体的には、補助的な平方非剰余 z を必要とせず「高速」。Kong らは、第1の場合を簡単化した。第2の場合（a が四乗剰余）には、Müller のアルゴリズムは z が必要という意味で「有料」だが、2023年、われわれは第2の場合を簡単化した。

のみならず、第2の場合についても、一定確率で「無料」の計算が成り立つ。以下のアイデアの核心は自明に近いが、Müller のアルゴリズムとの組み合わせ方は自明とも言い切れず、すてきな散歩道だ！

§55.　ショートカットの仕組みは単純明快。p = 41, k = 2 として a ≡ 10 の平方根 √10 を求めたいとする（§53の例題1と同じ）。このとき…
　　a^2k+1 ≡ 10⁵ ≡ +1　　（＊）
…だが、両辺を a 倍すれば:
　　a^2k+2 ≡ 10⁶ ≡ 10
その両辺を開くと:
　　a^k+1 ≡ 10³ ≡ √10
要するに √10 ≡ 10³ ≡ 16。この例は Müller のアルゴリズムでは、もともと無料（第1の場合）。でも上の計算法は、Müller のアルゴリズムより軽快な感じがする。2a を経由せず、直接 a を扱うからだろう。

これが成り立つのは（＊）の値が +1 だから。一般に、p が16k+9型で a が平方剰余なら a^(p−1)/2 = a^8k+4 ≡ +1、従って a^4k+2 ≡ ±1。この複号がプラスになれば（確率50%）――つまり a が四乗剰余であれば――上記のようなショートカットの余地が生じる。a^4k+2 ≡ +1 の場合 a^2k+1 ≡ ±1 だが、この複号が再びプラスになれば（確率25%）――つまり a が八乗剰余であれば――平方根計算が無料になる。

補題1（16k+9型素数 p を法とする無料平方根）　a^2k+1 ≡ +1 なら √a ≡ a^k+1 (mod p)

補題1自体は、ほぼ自明。きれいな式だが、残念なことに、このショートカットが役立つのは、平均4回に1回だけ…。けれど、補題1を Müller のアルゴリズムと組み合わせる手があって、「確実に近道ができるときだけ、近道をする」「Müller のアルゴリズムが面倒になる第2の場合にだけ、近道を考える」といった選択肢もある。しかも −1 の平方根が判明している場合（例えば別の計算の副産物として）、近道ができる確率が25%から50%にアップする！

こうなると理論的に面白いばかりか、実用上、役立つ場面も考えられる。

例えば p = 89, k = 5 として a ≡ 39 の平方根 √39 (mod 89) を求めたいとしよう。2a ≡ 78 ≡ −11 (mod 89) となる。Müller のアルゴリズムに従って、まず y ≡ (2a)^k ≡ (−11)⁵ ≡ −50。次に i ≡ (2a)y² ≡ (−11)⋅(−50)² ≡ −88 ≡ +1。この数は −1 の平方根ではない。第2の場合となり、普通にやると、平方非剰余 z を探して y と i を再計算・再設定することになる。しかし…
　　i ≡ (2a)y² ≡ (2a)[(2a)^k]² ≡ (2a)^2k+1
…である。この状況において、もし 2^2k+1 ≡ 2¹¹ ≡ +1 という事実が事前に分かっていたら（あるいは、その場で計算されたとしたら）どうなるか。この例では i ≡ +1 なのだから:
　　+1 ≡ i ≡ (2a)^2k+1 ≡ 2^2k+1⋅a^2k+1 ≡ (+1)⋅a^2k+1 ≡ a^2k+1
これは補題1の前提を満たすので、
　　√39 ≡ 39⁶ ≡ 67 ≡ −22
が 39 の一つの平方根。実際 (±22)² ≡ 39。平方非剰余が不要の「無料」計算が成立！

このように i ≡ (2a)y² ≡ (2a)^2k+1 ≡ +1 の場合、もし 2^2k+1 ≡ +1 が事前に分かっているなら a^2k+1 ≡ +1 なので、100%ショートカットが成り立つ。一方、同じ状況で 2^2k+1 ≢ +1 が事前に分かっているなら、補題1のショートカットは成り立たない。

同様に i ≡ (2a)^2k+1 ≡ −1 の場合、もし 2^2k+1 ≡ −1 が事前に分かっているなら a^2k+1 ≡ +1 なので、100%ショートカットが成り立つ。一方、同じ状況で 2^2k+1 ≢ −1 が事前に分かっているなら、補題1のショートカットは成り立たない。

〔例〕　16k+9型素数として k = 17 の場合の p = 281 を考える。事前に 2^2k+1 ≡ 2³⁵ ≡ −1 が分かっているとしよう。a ≡ 35 の平方根を求めたい（2a ≡ 70）。 y ≡ 70¹⁷ ≡ 2, i ≡ 70⋅2² ≡ −1。補題1のショートカットを使える:
　　√35 ≡ 35¹⁸ ≡ 63
一方 a ≡ 29 の平方根を求めたいとしよう（2a ≡ 58）。 y ≡ 58¹⁷ ≡ −25, i ≡ 58⋅(−25)² ≡ +1。補題1のショートカットを使えない。a^2k+1 ≡ +1 ではなく a^2k+1 ≡ −1 になってしまうので。この最後の式の両辺を a 倍すると、
　　a^2k+2 ≡ −a　すなわち　a^k+1 ≡ √−a
なので、−a ≡ −29 の平方根 √−29 の計算なら、補題1のショートカットを使える（だから、もし √−1 の値が分かるなら、無料で得た √−29 に √−1 を掛け算することで、本来求めたかった √+29 も無料で計算できることになる）。

Müller のアルゴリズムの第2の場合において、ショートカットが成り立つ可能性があるのは 2 が四乗剰余の場合、つまり 2^4k+2 ≡ +1 の場合に限られる。実際、その両辺を開くと 2^2k+1 ≡ +1 or −1 なので、上記の条件で補題1を使える。

16k+9型素数を法とする無料平方根　次の場合に補題1による近道ができる:
　　2^2k+1 ≡ +1 なら Müller のアルゴリズムの i が ≡ +1 の場合
　　2^2k+1 ≡ −1 なら Müller のアルゴリズムの i が ≡ −1 の場合

他方、2 が四乗非剰余、つまり 2^4k+2 ≡ −1 の場合には、2^2k+1 は −1 の平方根であり ≡ ±1 の可能性はないので、この近道を使えない。

mod 16k+9 では 2 は平方剰余なので、50%の確率で四乗剰余（2^4k+2 ≡ ±1 の符号がプラス）。その50%のうちの半分のケース（トータルで25%の確率）で、補題1のショートカットを利用可能―― (2a)^2k+1 ≡ ±1 の符号が i ≡ ±1 の符号と一致する場合だ。

2^2k+1 が計算済みなら、第2の場合において、平均4回に1回、ショートカットを使える。ショートカットの可否も i の値から簡単に判断できる。あくまでオプションだが「条件を満たせば近道をする。満たさなければ普通に Müller のアルゴリズムを使う」といったことが考えられる。最初に1回だけ 2^2k+1 の計算コストが必要だが、その値が ≡ ±1 と分かれば、以降その mod では、本来有料のケースの半分^*が無料になる。

*　2^2k+1 が ≡ +1 or −1 になる確率は約50%。もしそうなったら（第2の場合の）半分のケースで近道ができる。50%のうちの半分のケースなので、トータルでは平均25%の確率。

補題1の変種のショートカットも考えられ、状況によっては、無料になる確率を25%より高くできる。

もし −1 の平方根（±s とする）が分かっていれば、2 が四乗剰余のケースでは、100%ショートカットが成り立つ。さらに √s と √−s も分かっているなら、2 が四乗非剰余のケースでも、100%ショートカットが成り立ち、Müller のアルゴリズムの第2の場合（本来有料）は完全に無料になる。実際には、−1 の平方根のそのまた平方根 √s などを求めるには、それなりのコストが必要。「100%無料化のためには、最初に手数料がかかる」という痛しかゆしの状況…。でも、実用性を度外視して純粋に理論的に見た場合、こんなショートカットが存在すること自体、好奇心を刺激する。

〔参考〕　2^2k+1 の値の事前計算なしに、いちかばちかで b = a^k+1 を求めてみて、b² ≡ a が成り立てば b を a の平方根として終了する手もある。理論上、一般のケースにおいて平方非剰余 z を見つけるには、成功するまで z の値を変えながら z^8k+4 を平均2回、計算する必要がある（実際には Jacobi 記号の計算でもいい）。そのコストを考えると、25%の可能性があるなら、駄目もとで b² = a^2k+2 を計算してみるという選択肢も、一応、検討に値する。ただし p が12の倍数と1違いでなければ z ≡ 3 は平方非剰余だし（ここでは p ≡ 1 mod 4 と仮定しているので p ≢ 1 mod 3 なら 3 は平方非剰余）、同様に z ≡ 5 や z ≡ 7 が使えるケースは簡単に判別できるので、非剰余探しは多くの場合、低コスト。2^2k+1 の値を知らずに、やみくもに a^2k+2 を計算するメリットがあるとすれば、非剰余を見つけにくいような限定的ケースだろう。

2023-06-01　Dirichlet の美しい証明　x⁴ ≡ 2 が解を持つ条件

#数論　#平方根 mod p　#四乗剰余

4 の倍数より 1 大きい素数 5, 13, 17, 29, … は、2個の平方数の和として（本質的には一通りに）表現される:
　　5 = 1² + 2², 13 = 3² + 2², 17 = 1² + 4², 29 = 5² + 2², …

それを p = a² + b² としよう。p は奇数の一種なので、a と b が両方奇数または両方偶数になることは、あり得ない（それじゃ p が偶数になっちゃう）。2個の数のうち、奇数の方を a、偶数の方を b としよう。もし偶数 b が 8 の倍数なら、x⁴ ≡ 2 (mod p) は解を持ち、そうでなければ、この4次式は解を持たない。

〔例1〕　p = 41 = 5² + 4²　←　b が 8 の倍数じゃないので x⁴ ≡ 2 (mod 41) は解なし。

〔例2〕　p = 73 = 3² + 8²　←　b が 8 の倍数なので x⁴ ≡ 2 (mod 73) に解あり（x ≡ ±18, ±25）。

これは本来、初等（遊び）の範囲を超えた代数的整数論の問題だが、Dirichlet による初等的な証明 [18] があるので、それを紹介したい。

われわれの文脈では、2 が四乗剰余か否かということは、ある種の平方根の計算でのショートカットの可否と関係している。別の角度から言うと、この問題は、Fermat のクリスマス定理の拡張のうち、やや高級な例（x² + ny² の形を持つ素数の特徴付けのうち、n = 64 のケース）に対応する。

§56.　p を4k+1型の（正の）素数とする。奇数 a と偶数 b を使って p = a² + b² と書ける。奇数 a を法としてこの平方和を考えると、a ≡ 0 (mod a) なので:
　　p ≡ 0² + b² ≡ b² (mod a)
つまり p は mod a の平方剰余（平方根は b）。Jacobi 記号バージョンの相互法則から、ひっくり返して a は mod p の平方剰余。Euler の基準から:
　　a^(p−1)/2 ≡ +1 (mod p)　　【ラ】

一方 mod p で考えると:
　　a² + b² ≡ 0　つまり　b² ≡ (−1)a²
mod p での a の逆数 1/a と b の積を b/a ≡ i とすると:
　　b ≡ ai (mod p)　　【リ】
両辺を平方すると b² ≡ (ai)² ≡ (i²)a²。上の式 b² ≡ (−1)a² と見比べると i² ≡ −1 だ（普通と違う方向だが、p = a² + b² と書けることを出発点に、第一補充法則が導かれた）。

さて (a + b) は、もちろん奇数だが…
　　(a + b)² = a² + b² + 2ab = p + 2ab　　【ル】
　　(a − b)² = a² + b² − 2ab = p − 2ab
上と下を足すと:
　　(a + b)² + (a − b)² = 2p
この式を mod (a + b) で考えると (a − b)² ≡ 2p となり、2p は平方剰余。Jacobi 記号で書けば:
　　(2p/(a + b)) = +1
Jacobi 記号の性質上 (2p/(a + b)) = (2/(a + b))(p/(a + b)) なので…
　　(2/(a + b)) = (p/(a + b))　　【レ】
…は、両辺とも +1 か、または両辺とも −1。【レ】の左辺は、（Jacobi 記号の）第二補充法則から (−1) の [(a+b)²−1]/8 乗に等しい。この指数は、【ル】によって、次の自然数に等しい:
　　[(a + b)² − 1]/8 = (p + 2ab − 1)/8 = (p − 1)/8 + ab/4
従って、【レ】左辺は次の値を持つ:
　　(−1)^{[(p−1)/8+ab/4]} = (i²)^{[(p−1)/8+ab/4]}
　　= (i²)^(p−1)/8⋅(i²)^ab/4 = i^(p−1)/4⋅i^ab/2
一方、【レ】の右辺に相互法則を適用すると = ((a + b)/p) だが、これは Legendre 記号とも解釈可能なので、Euler の基準から ≡ (a + b)^(p−1)/2 (mod p)。結局【レ】の右辺と左辺の関係は:
　　(a + b)^(p−1)/2 ≡ i^(p−1)/4⋅i^ab/2 (mod p)　　【ロ】

今、【ル】を mod p で考えると:
　　(a + b)² ≡ 2ab (mod p)
この両辺を (p−1)/4 乗すれば（p についての仮定から、この指数は自然数）、左辺は (a + b)^(p−1)/2 になるので、【ロ】と比較できる。まず…
　　(a + b)^(p−1)/2 ≡ (2ab)^(p−1)/4 ≡ 2^(p−1)/4⋅a^(p−1)/4⋅b^(p−1)/4
…と変形して、そこに【リ】を代入:
　　(a + b)^(p−1)/2 ≡ 2^(p−1)/4⋅a^(p−1)/4⋅a^(p−1)/4⋅i^(p−1)/4 ≡ 2^(p−1)/4⋅a^(p−1)/2⋅i^(p−1)/4
【ラ】を使って簡単にすると:
　　(a + b)^(p−1)/2 ≡ 2^(p−1)/4⋅i^(p−1)/4
これを【ロ】と比較すると:
　　2^(p−1)/4⋅i^(p−1)/4 ≡ i^(p−1)/4⋅i^ab/2 (mod p)
上の式の両辺を i^(p−1)/4 で割ると（i の累乗は ≡ ±i or ±1 で ≡ 0 じゃないんで、この割り算はOK）:
　　2^(p−1)/4 ≡ i^ab/2 (mod p)　　【ワ】

特に mod p において 2 が四乗剰余なら、もちろん 2 は平方剰余でもあるので 2^(p−1)/2 ≡ +1 で、しかも
　　2^(p−1)/4 ≡ +1
となる〔※〕。逆に、上記の関係が成り立つなら 2 は四乗剰余。ところが【ワ】右辺によれば、2^(p−1)/4 ≡ +1 が成り立つことは、i の指数 ab/2 が 4 の倍数であることと同値。 a は奇数なので、この条件は b が 8 の倍数であることと同じ。次の結論に至る。

4k+1型素数を法として 2 が四乗剰余になる条件　4k+1型の素数 p について、p を奇数 a の平方と偶数 b の平方の和として書いたとき、もし b が 8 の倍数 b = 8m なら 2 は四乗剰余、さもなければ 2 は四乗非剰余: 2 が四乗剰余になることは、p が a² + (8m)² = a² + 64m² の形を持つことと同値。

〔※〕補足　2 が四乗剰余ということは、2 の四乗根（つまり平方根の平方根）が存在。従って 2 は平方剰余で、
　　r² ≡ 2　　【ヲ】
を満たす解 r が存在し、しかもこの r 自身も平方剰余で x² ≡ r を満たす x が存在する。 r についての Euler の基準から:
　　+1 ≡ r^(p−1)/2 ≡ (r²)^(p−1)/4 ≡ 2^(p−1)/4
最後の合同記号は【ヲ】による。

クリスマス定理風に言えば: 素数 p が x² + 64y² の形を持つための必要十分条件は、p が4k+1型で、かつ 2 が法 p の四乗剰余であること。既に1750年ごろ Euler はこの事実を予想していたという。すごい直観力だ！

最初の正式な証明は Gauß の複素整数論による。Dirichlet による軽妙な簡単化（上の内容）も味わい深く、[21] では Dirichlet’s beautiful proof と呼ばれている。天下り的で特殊な奇手のように見える部分もあるが、使われているいろいろな命題の背後には、一般性のある現象が潜んでいる。

§57.　ところで p が4k+3型の場合、2 が四乗剰余か否かは自明。x⁴ ≡ 2 (mod p) の解があるとすれば、それは 2 の平方根のそのまた平方根なので、とりあえず 2 の平方根が存在することが必要。p が奇素数なら、第二補充法則から p ≡ ±1 (mod 8) が前提。従って4k+3型なら、8t+7型であることが必要（言い換えると k = 2t+1 は奇数でなければならない）。

そして、それだけで十分！　実際、その条件において Euler の基準から 2^(p−1)/2 ≡ 2^4t+3 ≡ +1 (mod p)。両辺を 2 倍すると:
　　2^4t+4 ≡ +2
今 x ≡ 2^t+1 と置けば…
　　x⁴ ≡ 2^4t+4 ≡ +2
…となって 2 は少なくとも一つの四乗根 x を持つ（四乗剰余のことを簡潔に「四乗根」と呼ぶことにする）。

例として、4k+3型の素数 p = 23 を考える（k = 5, t = 2）。 x ≡ 2^t+1 ≡ 2³ ≡ 8 は 2 の四乗根。事実 8² = 64 = 69 − 5 ≡ −5 であり、(−5)² ≡ 25 ≡ 2。法が4k+3型のときの平方根計算は簡単なので、普通に r = √2 を求めて、その r を使って ±√r と ±√−r を求めれば、4種類の解が全部得られるだろう。

4k+3型素数を法として 2 が四乗剰余になる条件　4k+3型の素数 p を法とする場合、p が8t+7型か・8t+3型かに応じて、2 は四乗剰余・四乗非剰余。証明は単純: 前者の場合 2 の四乗根を直接構成できるし、後者の場合、平方剰余ですらないので（第二補充法則）、四乗根が存在するわけない。

〔補足〕　2³ ≡ 8 が 2 の四乗根ということは (2³)⁴ = 2¹² ≡ 2 を意味するが、第二補充法則と Euler の基準から…
　　2^(p−1)/2 = 2¹¹ ≡ +1　　（p = 23 の例）
…なのだから 2¹² = 2¹¹ × 2 ≡ (+1) × 2 は必然。一般に p が 8t+7 型の素数なら t = (p−7)/8 なので:
　　t + 1 = (p−7)/8 + 8/8 = (p+1)/8　従って
　　2^t+1 = 2^(p+1)/8　の四乗は　[2^(p+1)/8]⁴ = 2^(p+1)/2 = 2^(p−1)/2 × 2 ≡ (+1) × 2 = 2
なぜなら第二補充法則と Euler の基準から 2^(p−1)/2 ≡ +1。すなわち 2^t+1 = 2^(p+1)/8 は 2 の四乗根。

ところで 4k+3型素数 p は a² + b² の形では書けない。このような素数は、Gauß 整数としても素数で、従って p = (a + bi)(a − bi) という分解が成り立たないから。

他方、4k+1型素数 p は必ず p = (a + bi)(a − bi) = a² + b² の形で書けるものの、だからといって、2 が四乗剰余になるとは限らない。まず p が8t+5型の場合、b ≡ 2 or 6 (mod 8) なので、2 は四乗剰余にならない。

mod 8 で考える。奇数の平方 a² が ≡ 1 であることに注意すると、もし p が8t+5型なら、つまり p ≡ 5 なら、a² + b² = p は 1 + b² ≡ 5 を意味する。そのとき b² ≡ 4、従って b ≡ ±2 であり、b は 8 の倍数より 2 大きいか 2 小さい。

というか、第二補充法則による制約があるので、法 p が8t+5型なら 2 は平方剰余（平方数）ですらない。ましてや四乗剰余（四乗数）のわけない。

p が8t+1型のときに限って b が 8 の倍数になるケースがあり、その場合 mod p において 2 は四乗剰余――これは b ≡ 0 (mod 8) のときだ。「p は8t+1型だけど b は 8 の倍数じゃないよ」というケースもあり、その場合 mod p において 2 は四乗非剰余―― b ≡ 4 (mod 8) のときだ。

mod 8 で考える。もし p が8t+1型つまり ≡ 1 なら、a² + b² = p は 1 + b² ≡ 1 を意味する。従って b² ≡ 0 だ。 b ≡ 0 はこの条件を満たすが b ≡ 4 も同じ条件を満たす。つまり b は 8 の倍数であるか、または 8 の倍数より 4 大きい。

法 p が8t+1型のとき、2 が四乗剰余かどうか――それが問題の核心。次の区別によって Yes か No かが決まる: p = a² + b² と書いたとき b は 8 の倍数か、それとも（4 の倍数だが）8 の倍数ではないか？　言い換えれば p は a² + 64m² 型か否か？

8t+1型の素数 p のうち、次の例では b が 8 の倍数となり、mod p において 2 は四乗剰余:
　　73 = 3² + 8²
　　89 = 5² + 8²
　　113 = 7² + 8²
　　233 = 13² + 8²
　　257 = 1² + 16²
次の例では b が 8 の倍数にならず、mod p において 2 は四乗非剰余:
　　17 = 1² + 4²
　　41 = 5² + 4²
　　97 = 9² + 4²
　　137 = 11² + 4²
　　193 = 7² + 12²
　　241 = 15² + 4²

〔文献〕
[18] Dirichlet: Ueber [Über] den biquadratischen Character der Zahl „Zwei”
 https://gdz.sub.uni-goettingen.de/id/PPN243919689_0057?tify={%22pages%22%3A[191]%2C%22view%22%3A%22info%22}
 https://gdz.sub.uni-goettingen.de/download/pdf/PPN243919689_0057/LOG_0020.pdf
[19] Lemmermeyer: Reciprocity Laws: From Euler to Eisenstein, Chap. 5, §5.1
[20] Cox: Primes of the Form x²+ny² (2nd ed.), §1D, §4B
 https://www.math.toronto.edu/~ila/Cox-Primes_of_the_form_x2+ny2.pdf
[21] Ireland & Rosen: A Classical Introduction to Modern Number Theory, Chap. 5, Exercises 26–28 (p. 64)

2023-06-24　相互法則の簡単な紹介　分かりやすく？！

#数論　#平方剰余　#心の魔法

「平方剰余（じょうよ）の相互法則」って、名前からして難しそうじゃん、チョイと「素人お断り」みたいな…？

でも、法則の「意味」を理解することは難しくないッ！

漫画やアニメや音楽と同じで、好みは人それぞれ。数論が好きな人もいれば、別のことが好きな人もいる。数論は、優雅で純粋。直接何かに（例えば学校の試験に）役立つといった、卑俗なものじゃない――「何かのため」ではなく「それ自体」の遊び。「何の役にも立たない遊びに、時間を使いたくない」という考えの人もいるだろう。「数の不思議のようなことに興味はあるんだけど、どうもよく分からない」みたいな人もいるかも…？

普通の整数の平方を「平方数」と呼ぶなら 16 や 25 や 49 などが平方数であることは、すぐ分かる。17 や 26 や 50 や 51 などが非平方数であることも、明らかだろう。これと同様のことを「普通の整数の世界」ではなく、ある種のパラレルワールド（例えば「曜日の世界」）で考えた場合――その具体的な意味については後述――、平方数・非平方数に対応する概念が、平方剰余・非剰余なんだけど、まぁ硬い話はさておき！

誰かがこう尋ねたとして、あなたなら何て答える…？！
　　 誰かの質問:　整数の平方を 3 で割ったとき、余りが 1 になることってあるの？

§58.　「整数の平方を 3 で割ったとき、余りが 1 になることがあるか？」への答えは、もちろん Yes だろう。例えば 4 っていう数の平方は 16、それを 3 で割れば（商が 5 になって）余りは 1。あるいは 5 の平方は 25 だが、それを 3 で割れば（商が 8 で）余りは 1。

実のところ、3 の倍数以外なら、どんな整数でも「平方して 3 で割ったとき 1 余る」。

〔参考〕　それを証明するには…別にしなくてもいーんだけど、参考までに。3 の倍数以外の整数は、3 の倍数より 1 大きいか・1 小さいかのどちらか。3m + 1 または 3m − 1 の形を持つ（m: 整数）。それを平方すれば、結果はこーなる…（なぜ）？
　　(3m + 1)² = (3m)² + 2⋅3m⋅1 + 1² = 3(3m² + 2m) + 1
　　(3m − 1)² = (3m)² − 2⋅3m⋅1 + 1² = 3(3m² − 2m) + 1
この右端の丸かっこ内は整数なので、右端の数は「3 × 整数、プラス1」。つまり 3 の倍数より 1 大きく、3 で割ると 1 余る。

じゃあ、次の質問への答えは…？！
　　 誰かの質問:　整数の平方を 3 で割ったとき、余りが 2 になることってあるの？

今度の答えは、まぁ No でしょうねぇ～。だって考えてる整数が 3 の倍数なら（つまり 3m の形なら）、その平方は再び 3 の倍数じゃん（ってゆーか 9 の倍数だけど、それは 3 の倍数でもある）。考えてる整数が 3 の倍数以外なら、上記のように、その平方は 3 の倍数より 1 大きい。どっちにしても 3 で割った余りは 0 か 1 のどっちか。2 には、ならない…よね？

要するに x を整数としたとき「x² を 3 で割ると余りイコール 1」という方程式（？）は解を持つけど、「x² を 3 で割ると余りイコール 2」には解がない、と。

証明できるかはともかく、実際の数値を当てはめることによって、同様に、次の事実を予想できる…
　「x² を 7 で割ると余りイコール 1」は、解あり（そういう性質を持つ x の例: 1 や 6）
　「x² を 7 で割ると余りイコール 2」は、解あり（そういう性質を持つ x の例: 3 や 4）
　「x² を 7 で割ると余りイコール 3」は、解なし
　「x² を 7 で割ると余りイコール 4」は、解あり（そういう性質を持つ x の例: 2 や 5）
　「x² を 7 で割ると余りイコール 5」は、解なし
　「x² を 7 で割ると余りイコール 6」は、解なし

解があるかないか考えている方程式（のようなもの）を例えば次のように表現する。
　　x² ≡ 2 (mod 7)
　　大ざっぱな意味:　x² を 7 で割ると余り 2
「そもそもなぜ、そんなことを考えるのか？」という点をひとまず脇に置くと、内容は単純明快。

関連して、次の例のような「分数に丸かっこを付けた変な記号」が使われる。
　　(2/7) = +1　　意味: 上記の x² ≡ 2 (mod 7) は解を持つ
　　要するに: 「平方して 7 で割ると 2 余るような整数」は存在する
　　(3/7) = −1　　意味: x² ≡ 3 (mod 7) には解がない
　　要するに: 「平方して 7 で割ると 3 余るような整数」は存在しない
別の例として…
　　(1/3) = +1　　意味: 「平方して 3 で割ると 1 余る整数」はある
　　(2/3) = −1　　意味: 「平方して 3 で割ると 2 余る整数」はない

最初に書いたのと同じ単純な内容を記号化しただけ…なのだが、この記号の意味は（慣れないうちは）微妙に分かりづらいかも…。分数のように見えるものは、本物の分数ではなく、単に便宜上の記号。

より一般的に p を 3 以上の素数として、a を p の倍数以外の数とした場合、次の記号の大ざっぱな意味は…
　　(a/p)
…平方して p で割ると余りが a になるような整数は存在するか？という質問。答えが Yes なら、この記号は +1 という値を持ち、答えは No なら、この記号は −1 という値を持つ、と。これは「なぜ？」っていう問題ではなく、「そう定義すると便利なので、そう約束しますよ」と。

この記号の名前は Legendre（ルジャンドル）記号。Legendre という数学者が発明した記号（らしい）。フランスの Legendre 大先生は、ドイツの天才青年 Gauß にばかにされてカンカンに怒った…らしいけど、高木先生の『初等整数論講義』でも「未だ平方剰余の相互法則を確定することもできず」（239ページ）、「依然として雑纂（ざっさん）的の旧態を改めていない」（371ページ）、さらに追い打ちをかけるように（Legendre と Gauß について）「ハンディキャップはあまりに大きい」と、けちょんけちょんに言われている。現代でも使われる便利な記号の発明者なのに、ちょっと気の毒…

〔注〕　HTML ページで分数みたいなものを表示するのは面倒…という実用上の理由から、このサイトの他のメモでは
　　(a/p)
の代わりに L(a, p) と表記されてることも多いです。L は Legendre の L。

§59.　「一定の数で割った余りで分類して、余りが同じ数同士を同一視する」…というと難解なようだけど、実際には日常、当たり前に行われてること。例えば「偶数」って言葉は 2 で割ると 0 余る（つまり割り切れる）ってことだし、「奇数」という言葉は 2 で割ると 1 余るという性質。2 で割った余りで分類してる。4 と 10 はどちらも偶数であるとか、3 と 7 はどちらも奇数であるという主張は、当たり前というか、だから何？って感じ。けど、あえて記号で書けば…
　　4 ≡ 0 (mod 2)　そして　10 ≡ 0 (mod 2)　だから　4 ≡ 10 (mod 2)　　大ざっぱな意味: どちらも偶数
　　3 ≡ 1 (mod 2)　そして　7 ≡ 1 (mod 2)　だから　3 ≡ 7 (mod 2)　　大ざっぱな意味: どちらも奇数

でもって「ある月の日付の 1 と 8 は曜日が同じ」とかも、極めて日常的なことでしょう…。その月の「1日」と「8日」は別の日付だけど、曜日という観点では「同じ」と同一視される。「7 で割った余りを基準とする」という考え方を mod 7 で表すなら…
　　1 ≡ 8 (mod 7)
とか、あるいは
　　2 ≡ 16 (mod 7)
のような関係式が成り立つ。この柔軟な考え方を受け入れると、前節の
　　x² ≡ 1 (mod 7) は解を持つ
という主張は、次の主張と全く同じ意味になる:
　　x² ≡ 8 (mod 7) は解を持つ
なぜなら 1 と 8 は「同じもの（例えば同じ曜日）の別名」だから。

普通の整数として「1 と 8 が等しい」という意味じゃないよ、もちろん。だけどさぁ 7 で割った余りで分類するんなら 1 と 8 は、どっちも同じ種類ってゆーか、同じ仲間じゃん？　この仲間たちについては、「1 が属する集団」という意味で 1 で表すのが最もシンプルだけど、「8 が属する集団」という意味で 8 で表しても差し支えない。その意味では 1 も 8 も、同一の集団を指す別名といえる。三本線のイコールのような記号 ≡ は「三角形の合同」なんかでも使われることがあり、図形の合同とは意味が違うけど、数論でも「mod 7 では 1 と 8 は合同」みたいな表現が使われる（この記号についてもっと詳しく）。

この一見、ちんぷんかんぷんな議論のとりあえずの重大なポイントとして、上記のように:
　　x² ≡ 8 (mod 7) は解を持つ。
前節に従って、これを「平方して 7 で割ると 8 余るような整数」が存在する…と解釈すると、訳が分からなくなっちゃう…。だって 7 で割った余りは 6 以下のはず…。8 余るなんて、おかしいじゃん。「余りが○になる」みたいな前節の話は、あくまで大ざっぱな説明。
　　質問　x² ≡ a (mod 7) は解を持つか？
の正確な意味は、「x を平方して 7 で割った余りが、a を 7 で割った余りと一致するような」そういう x は存在するか？――もう少しきっちり言うと、7 で割った余りで分類したとき、x² と a が同じ種類ってことは、あり得るか？

このコンセプト、慣れてしまうと何でもないことなんだけど、慣れないうちは、スッキリしないかも…。例えるなら、多くの人にとって 7 × 2 が 14 に等しいこと、7 × 9 が 63 に等しいことは「明らか」。でも 7 × 9 と 63 が等しいことは、なぜ明らかなんだろう…。もしも「掛け算」を知らない人から理由をしつこく尋ねられたら、あなたは「7 を 9 回足せば 63 になるでしょ」と答えるかもしれない。あるいは「7 の 10 倍は 70 で、9 倍はそれマイナス 7 だから」と答えるかも。いずれにしても「掛け算」を知らない人は混乱して、そんな計算を何で当たり前のように一瞬でできるのか、と不思議がる…。一瞬にして 9 個の 7 を足し合わせるなんて、この人は天才だろうか？　あなた自身にとっては、7 × 9 = 63 なんて、全然悩むようなことじゃないし、別に天才的な計算でもない…よねぇ？

この違いが生じる理由は単純で、あなたは掛け算に慣れてるし「掛け算は特に難しい計算でもない」と認識してる。「掛け算」を知らない人の観点では、何が何だか話が見通せない。立場を逆転させて、もしあなたが平方剰余に慣れてないなら、
　　(2/7) と (9/7) は同じ意味でどちらも +1 に等しい
…といった話について「すごく難解なこと」という印象を抱いてしまうかもしれない！

実際には、これは単に「2 ≡ 9 (mod 7) は平方数だよ。だって 3 の平方じゃん」という当たり前の主張に過ぎない…のだが、それを「当たり前」と思えるようになるためには、多少の慣れが必要。心理学的に言えば、こーゆーとき…
　　×　してはいけないこと　「難解。自分には無理」と勝手に思い込み、ネガティブな自己暗示を掛ける…　×
　　○　正しい観点　「慣れてないから、いまいちピンとこないだけかもね…」と認識する　○

「慣れてないことなので、感覚がつかめないのは当たり前。やがて慣れれば当たり前になる。わたしは今からこれを理解する。時間はかかるかもしれないけど、やがて完全に理解する」と自分に言い聞かせよう。気の持ちようの、ほんのちょっとの違いで、世界が180°変わるかも！

Everything is possible, even the impossible. — Mary Poppins

§60.　平方剰余の相互法則は、3 以上の相異なる素数 p, q が与えらたとき、
　　(p/q) と、それをひっくり返した (q/p) は一致するかしないか？
を記述したもの。p と q がどちらも 4 の倍数より 3 大きい素数なら（例えば 7 と 11 なら）、この二つは一致しない。片方が +1 なら他方は −1 になる。それ以外の場合、つまり p と q の少なくとも一方が 4 の倍数より 1 大きい素数（例えば 13）なら、この二つは一致する。片方が +1 なら他方も +1 だし、片方が −1 なら他方も −1 になる。

〔注〕　3 以上の素数はどれも「4 の倍数より 1 大きい」か、または「4 の倍数より 3 大きい」。このサイトのメモでは、非公式な愛称として、しばしば前者のタイプをバニラ素数、後者のタイプをチョコレート素数（チョコ素数）と呼んでいる。

平方剰余の相互法則については初等的な証明も可能だけど、一般向けの初等整数論入門のような本では、「読者を混乱させる」という理由で話題自体が省略されていることもあり^*、掲載されているとしても、大抵は一番最後の章に載っている。つまり「初等」の範囲では一番難しい話題の一つ（通例、次の章はもはや「初等」ではなく、ガウス整数のような「拡張された世界」の話に入る）。事実 Legendre 記号を発明した Legendre 自身、相互法則を証明できなかったし、Gauß による第一証明も、必ずしも分かりやすいものではなかった…（「吐き気を催させる」とこき下ろす人さえいる^**）。とはいえ、それは「証明がトリッキー」というだけで、法則の内容自体は、大して難しくない。Legendre 記号をひっくり返したとき、「値の ±1 の符号が変わるか変わらないか」を法則化しただけ――分数のようなものの分子・分母が「両方チョコなら、符号が変わる。少なくとも一方がバニラなら、符号が変わらない」という単純な現象となる！

* The law of quadratic reciprocity has been intentionally omitted. It is my belief that even for the serious number theory student, the law of quadratic reciprocity is better taught in a second course of number theory than a first course. After the student has worked with quadratic Diophantine equations and quadratic fields, then he can appreciate the usefulness of the quadratic reciprocity law. — Harold M. Stark, An Introduction to Number Theory

** これはガウスの第一証明についての「批判」。ガウスの第三証明は比較的分かりやすく、特にアイゼンシュタインによって簡単化されたバージョンは、現代の教科書でも広く採用されている。

相互法則がどう役立つのか先に例を見とくのも、「その証明を知りたい」と興味を抱くモチベーションになるだろう…。

【1】　相互法則の応用のほんの一例として、ここでは次の定理を導く。しげちーのスペシャルデーとの関連で、この定理は一つの鍵となる（§32）。

第三補充法則　p を 5 以上の素数とする。 p が12k±1型のとき、そしてそのときに限って、x² ≡ 3 (mod p) は解を持つ。

Legendre 記号を使えば (3/p) = +1 になる条件…ってわけ。相互法則を使って単純にひっくり返せるなら…
　　(3/p) = (p/3)
…だが、p を 3 で割った余りは 1 or 2 なので、このひっくり返した後の値（上の右辺）は、次のどちらかに等しい（§58）。
　　(1/3) = +1　さもなければ　(2/3) = −1

そんなわけで、最終的には p を 3 で割った余りによって ±1 の符号が決まる。一方、相互法則の教えによると「ひっくり返したとき何が起きるか」は p を 4 で割った余りによって決まる（p がバニラかチョコレートか）。もともとの「分子」にある 3 はチョコレートなので、「分母」もチョコレートだとストレートにひっくり返せない（分母がバニラなら単純にひっくり返せる）。従って p を 3 で割った余りと 4 で割った余りの両方を視野に入れる必要がある。すなわち p を 3 × 4 = 12 で割った余りで分類すればいい。5 以上の素数 p を 12 で割った余りは 1, 5, 7, 11 のどれかになる。

なぜ「それ以外の余りには、ならない」と言い切れるか。例えば 12 で割って余りが 2 や 4 や 6 や 8 や 10 なら、もともとの数は偶数。「5 以上の素数」という前提に反する。余りが 3 や 9 なら、もともとの数は 3 の倍数。これも前提に反する。ところで「分子・分母」という呼称は便宜上のもので、この記号は、普通の意味での分数とは何の関係もない。

p が 12 の倍数より 1 大きい場合（例えば p = 13）、その p は 3 の倍数より 1 大きい。しかも 4 の倍数より 1 大きいので、単純にひっくり返していい。
　　(3/p) = (p/3) = (1/3) = +1

〔注〕　第1の等号は相互法則による（バニラが含まれるので、分子と分母を入れ替えても、何も起きない）。第2の等号は、3 で割った余りを基準に p と 1 が同じ種類であることによる。第3の等号は x² ≡ 1 (mod 3) が解を持つこと（§58）による。ぶっちゃけ「1 = 1² は平方数」というだけだが…

p が 12 の倍数より 5 大きい場合（例えば p = 17）、その p は 3 の倍数より 2 大きい。しかも 4 の倍数より 1 大きいので、単純にひっくり返していい。
　　(3/p) = (p/3) = (2/3) = −1

p が 12 の倍数より 7 大きい場合（例えば p = 19）、その p は 3 の倍数より 1 大きい。一方、このような p は 4 の倍数より 3 大きいので、ひっくり返すと Legendre 記号全体の値は、符号が逆になる。次のように:
　　(3/p) = −(p/3) = −(1/3) = −(+1) = −1

p が 12 の倍数より 11 大きい場合（例えば p = 23）、その p は 3 の倍数より 2 大きい。一方、このような p は 4 の倍数より 3 大きいので、ひっくり返すと Legendre 記号全体の値は、符号が逆になる。次のように:
　　(3/p) = −(p/3) = −(2/3) = −(−1) = +1

「12 の倍数より 11 大きい」ということは「12 の倍数より 1 小さい」ともいえる。以上をまとめると、p が 12 の倍数 ± 1 なら x² ≡ 3 (mod p) は解を持つ。p がそれ以外なら、同じ方程式は解を持たない。

直接証明も可能な第三補充法則だが、相互法則を使うと、上記のように、機械的な単純計算。無味乾燥かもしれないが、ずいぶん楽だ。直接証明も面白いけど、効率でいえば、相互法則を経由した方が圧倒的に簡単！

【2】　最後にフィナーレにふさわしい大物（？）を…。

「第七補充法則」　p を 11 以上の素数とする。 p ≡ ±1, ±3, ±9 (mod 28) のとき、そしてそのときに限って、x² ≡ 7 (mod p) は解を持つ。

相互法則を仮定した上でこれを証明する方法については、一旦「練習問題」としたい。必要なヒントは、今回の内容に全て含まれている。勘がいい方なら、普通に【1】と同様にできるであろう。「やってみたけど、よく分からん」とか「それほど興味もないが、一応答えを知りたい」という方のために、次回に説明を掲載するが、説明を読むより自力で試してみた方が理解が深まるかと…。本当にそうなってるか、数値的に確かめてみるのも面白いだろう！

このメモの内容は「平方剰余の相互法則の説明・証明」というより、大ざっぱな「紹介」。「証明されてない定理を使うのは数学的に反則」みたいなお固い意見もあるだろうけど、Legendre も Gauß も実際に「証明」するより前に相互法則を「予想」して、その意味を考え、最終的に Gauß は証明に成功した。教科書的な「定義・定理・証明」の連鎖ばかりが、数論ではない――実際には、まず予想・直観があって、後から証明が来ることが多い。たとえ初等の範囲だって「遊び」があっていい。

2023-07-02　平方を 7 で割った余り・平方を p で割った余り　相互法則の例

#数論　#平方剰余　#平方根 mod p

整数 a = 1, 2, 3, 4, 5, … を平方した 1, 4, 9, 16, 25, … を 7 で割ると、余りは 1, 2, 4 のどれかになる。余りが 3, 5, 6 になることはない。もっとも、もともとの数 a が 7 の倍数なら a² も 7 の倍数になって 7 で割り切れる（余り 0）。

では逆に、「整数 a を平方して、7 より大きい素数 p で割ったら余りは 7 だった！」という現象が起きるのは、p がどういう素数の場合だろうか？

「整数の平方を p で割ると 7 余る」という現象の発生条件を知りたい（ここで「素数」とは、1 と自分自身でしか割り切れない 2 以上の自然数、例えば 11, 13, 17, 19, 23 を指す）。

§61.　どんな整数を平方して 11 で割っても、余りが 7 になることはない。つまり素数 p = 11 は、上記の条件を満たさない。

a = 1, 2, 3, 4, 5 の平方は、それぞれ 1, 4, 9, 16, 25 なので 11 で割った余りは 1, 4, 9, 5, 3。――ところで 11 で割った余りという観点では、a = 6 の平方は a = 5 の平方と同じ。なぜなら 6 ≡ −5 (mod 11) の両辺を平方すると 6² ≡ (−5)² ≡ 5² だから。実際 6² = 36 を 11 で割った余りは、5² = 25 を 11 で割った余りに等しい。同様に a = 7, 8, 9, 10 の平方は、それぞれ a = 4, 3, 2, 1 の平方と同じ。このことから、平方を 11 で割った余りは、上記 1, 4, 9, 5, 3 のどれかに限られる。――もっとも、例外的なケースとして、a が 11 の倍数なら a² は 11 で割り切れる（余り 0）。このケースを別にすると 12, 13, 14, … は ≡ 1, 2, 3, … なので、11 より大きい a を試しても、上記以外の新しい余りは出てこない。

同じように a の可能性を一つ一つ検討すると、p = 13 も p = 17 も、条件を満たさないことが分かる。一方、p = 19 は条件を満たす。
　　a = 8 のとき a² = 64 を 19 で割ると 7 余る！　（19 × 3 = 57）
この場合 a = −8 も平方すれば同じことなので a = 11 も「平方して 19 で割ったとき 7 余る」という性質を持つ: 11 ≡ −8 (mod 19) だからだ。
　　a = 11 のとき a² = 121 を 19 で割ると 7 余る！　（19 × 6 = 114）

整数 x = a が x² ≡ 7 (mod p) を満たすなら、0 以上 p 未満の範囲に、その a と同じ種類の数がある（等しい数はもちろん「同じ種類」だが、等しくない数も、同じ種類になることがある）。というのも mod p では、「p で割った余りが同じ」数は「同じ種類」と見なされる――言うまでもなく 0 以上 p 未満の範囲には、整数を p で割った余りが、全種類（過不足なく）存在している。

x = a が x² ≡ 7 (mod p) の解なら、x = −a も同じ式の解だが、上と同じ理由から、やはり 0 以上 p 未満の範囲に −a と同じ種類の整数がある。要するに、この範囲の中から2種類の整数 x と p − x を（2種類の解の代表例として）選択できる。

しかも、このように選んだ2種類の数のどちらか一方は p/2 より小さい: x が p の半分より小さければ、x 自身がこの性質を持つ。もし x が p の半分より大きい（けれど p よりは小さい）とすると、今度は p − x が p/2 より小さい。仮定によって p は奇数なので p/2 は整数ではなく、解の整数 x はそれより大きいか・小さいかのどっちか。結局:
　　x² ≡ 7 (mod p) が解を持てば、解の一つは 0 以上 p/2 未満の範囲にある。
　　逆に言うと、その範囲に解がないなら、解はどこにもない！

これは限定的な範囲――「有限の範囲」だ。 p が小さい場合、この範囲の整数を一つ一つ試せば、解があるかないか、簡単に確定できる。でも p がもっとでかい場合（例えば p = 8191 とか p = 131071 の場合）、p/2 未満の自然数は何千・何万種類もある。一つ一つの可能性について、「平方して p で割ると 7 余るケースがあるだろうか？」と検証するのは（たとえコンピューターを使うにしても）面倒くさい。

§62.　相互法則を利用すると、たった1回の割り算だけで、結論を出すことができる:
　　質問　x² ≡ 7 (mod p) は解を持つか？　ここで p は 7 より大きい素数。
　　答え　p を 28 で割った余りが 1, 3, 9, 19, 25, 27 のどれかなら Yes、それ以外なら No。

この答えが出る理由も、以下の通りシンプル。相互法則というのは、その理論的重要性を別にしても、この種の判定において、強力なツールとなる！

【1】　問題は、要するに (7/p) は +1 か −1 か？という判定。

ここで「分数に丸かっこを付けたような表記」は、前回紹介した Legendre 記号。

ひっくり返した (p/7) の判定なら、非常に易しい。整数の平方を 7 で割った余りは、最初に見たように 1, 2, 4 にはなり得るが 3, 5, 6 にはなり得ないのだから、p が 7 の倍数より 1 大きいか・2大きいか・4大きいなら判定結果は +1 だし、それ以外なら判定結果は −1 となる。この部分に関する限り、必要な情報は「p を 7 で割った余り」だけ。判定の内容も、簡単な算数に過ぎない（付録）。

(7/p) をひっくり返して (p/7) にしたとき、何が起きるか。この部分は、p が 4 の倍数より 1 大きいか・3 大きいか？によって変わってくる。

なぜなら、この「分数」に含まれる 7 という数は「4 の倍数よりも 3 大きい素数」だから。

もし p が「4 の倍数より 1 大きい素数」なら、単純にひっくり返すことができる。

もし p も「4 の倍数より 3 大きい素数」だと、ひっくり返したとき Legendre 記号の値の符号が反転する。

…これが平方剰余の相互法則の、表面的な意味だった。この「ひっくり返す」処理に関しては「p を 4 で割った余り」が不可欠の情報。

【2】　まとめると「p を 7 で割った余り」「p を 4 で割った余り」の両方の情報があればいい。だから p を 7 × 4 = 28 で割った余りによって、場合分けすればいい。この p という数は 7 より大きい素数（言い換えれば 11 以上の素数）なので、奇数。従って p を 28 で割った余りは、奇数に限られる。

〔補足〕　p を 28 で割った商を A、余りを B として、p = 28A + B とした場合、もしも B が偶数だとすると、p は偶数 28A と偶数 B の和なので再び偶数になってしまい、話の前提と矛盾する。つまり B は奇数でなければならない。

特に p を 28 で割った余り B が 0 になる可能性はない: 表面的には「0 も偶数だから」だけど、もっと単純に「余りが 0」とは「割り切れる」ってこと。p は素数なのだから 28 で割り切れるはずない。

同様の議論から、p を 28 で割った余りは 7 の倍数（7, 14, 21 など）にはなり得ない。

〔補足〕　もしも p = 28A + B の B が 7 の倍数だとすると、p は 7 の倍数 28A と 7 の倍数 B の和なので再び 7 の倍数になってしまい、話の前提と矛盾する: p は 11 以上の素数なので 7 の倍数ではない。つまり B は 7 の倍数以外でなければならない。

結局 p を 28 で割った余り B は 1, 3*, 5, 9, 11*, 13; 15*, 17, 19*, 23*, 25, 27* の 12 種類の数のどれか（28 より小さい正の奇数。ただし 7 の倍数 7, 21 は除外される）。この余りたちのうち * 印の 3*, 11* などは、4 の倍数より 3 大きい（言い換えると 4 の倍数より 1 小さい）。それ以外は 4 の倍数より 1 大きい。p = 28A + B のうち 28A の部分は 4 で割り切れるので、「p を 4 で割った余りが 1 になるか 3 になるか」は、「B を 4 で割った余りが 1 になるか 3 になるか」によって決まる。同様に「p を 7 で割った余り」は「B を 7 で割った余り」に等しい。

例えば B = 1 の場合（p = 29 のときなど）、p は 4 の倍数より 1 大きいので…
　　(7/p) = (p/7)
…と単純にひっくり返すことができる。この p は 7 の倍数より 1 大きいので、上記の右辺は次の値を持つ:
　　(1/7) = +1

第二の例として B = 5 の場合（p = 61 のときなど）も、p は 4 の倍数より 1 大きいので…
　　(7/p) = (p/7)
…と単純にひっくり返すことができる。この p は 7 の倍数より 5 大きいので、上記の右辺は次の値を持つ:
　　(5/7) = −1

第三の例として B = 3 の場合（p = 31 のときなど）、p は 4 の倍数より 3 大きいので…
　　(7/p) = −(p/7)
…のように、ひっくり返すと値の符号が逆転する。この p は 7 の倍数より 3 大きいので、上記の右辺は次の値を持つ:
　　−(3/7) = −(−1) = +1

他のケースも同様で、ひっくり返すとき符号が逆転するかどうかにだけ注意すると、最終的には…
　　(1/7) や (2/7) や (3/7)
…のような値の判定（それぞれ +1 か −1 か？の判定）に行き着く。

【3】　12種類のどのケースについても、上と同様の単純な議論になる。ここでは、途中経過の詳細を省く。各ケースを総合すると、次の結論に至る:
　　p を 28 で割った余り B が 1, 3, 9, 19, 25, 27 の6種類の数のどれかなら、必ず「とある整数の平方を p で割ると 7 余る」という現象が起きる。
　　B がそれ以外の数なら、どんな整数の平方を p で割っても、決して余り 7 にはならない。

B = 27 つまり「p が 28 の倍数より 27 大きい」というケース（例えば p = 83 = 28 × 2 + 27）については、あと 1 プラスすればちょうど 28 の倍数になるのだから、「p が 28 の倍数より 1 小さい」と言い換えることもできる（例えば p = 83 = 28 × 3 − 1）。

同様に「p が 28 の倍数より 25 大きい」は「p が 28 の倍数より 3 小さい」と言い換えられる。「p が 28 の倍数より 19 大きい」は「p が 28 の倍数より 9 小さい」と言い換えられる。このような言い換えを使って、上記の結論を整理し直すと:
　　素数 p が 28 の倍数 ±1, ±3, ±9 のどれかなら、「整数の平方を p で割ると 7 余る」という現象が起きる。
　　p がそれ以外の型の素数（ただし 11 以上）なら、この現象は起きない。

さらに簡潔に（少々小難しく）表現するなら、11 以上の素数 p について:
　　(7/p) = +1　⇔　p ≡ ±1, ±3, ±9 (mod 28)

これは一見、無味乾燥で、何の役に立つのか分からない「特殊な算数」のようにも思える。 Legendre は、この現象に気付いていて、何とかしてこの命題を証明しようと努力・工夫を重ねたが、「ひっくり返す」部分（相互法則）について、証明を完成させることができなかった。Gauß は最初に証明を完成させただけでなく、この法則が「整数の概念そのものを問い直すような深い観点」と結び付いていることまで見抜き、実際に「複素数の範囲の整数」という新世界を切り開いた。

§63.　（付録）　次の関係について。
　　(1/7) = +1
　　(2/7) = +1
　　(3/7) = −1
　　(4/7) = +1
　　(5/7) = −1
　　(6/7) = −1

つまり (c/7) は c = 1, 2, 4 のとき値 +1 を持ち c = 3, 5, 6 のとき値 −1 を持つ。

言い換えると x² ≡ c (mod 7) は c ≡ 1, 2, 4 なら解を持ち c ≡ 3, 5, 6 なら解を持たない。

要するに、整数（7 の倍数以外）の平方を 7 で割った余りは、必ず 1, 2, 4 のどれかになり、決して 3, 5, 6 にはならない。

証明　7 の倍数以外の整数 x は、必ず 7m ± 1 の形か、または 7m ± 2 の形か、さもなければ 7m ± 3 の形を持つ（m は整数）。

第一に、もし x = 7m ± 1 ならば:
　　x² = (7m ± 1)² = (7m)² + 2(7m)(±1) + 1²
この右端の和の第1項・第2項は 7 の倍数なので、x² を 7 で割れば、第3項の 1² = 1 が余る。

要するに 7 の倍数と 1 ずれた数を平方すると、結果は 7 の倍数より 1 大きい。 8² = 64 = 63 + 1 や 6² = 36 = 35 + 1 のように。

第二に、もし x = 7m ± 2 ならば:
　　x² = (7m ± 2)² = (7m)² + 2(7m)(±2) + 2²
この右端の和の第1項・第2項は 7 の倍数なので、x² を 7 で割れば、第3項の 2² = 4 が余る。

要するに 7 の倍数と 2 ずれた数を平方すると、結果は 7 の倍数より 4 大きい。 9² = 81 = 77 + 4 や 5² = 25 = 21 + 4 のように。

第三に、もし x = 7m ± 3 ならば:
　　x² = (7m ± 3)² = (7m)² + 2(7m)(±3) + 3²
この右端の和の第1項・第2項は 7 の倍数なので、x² を 7 で割った余りは、第3項の 3² = 9 を 7 で割った余り、つまり 2 に等しい。

要するに 7 の倍数と 3 ずれた数を平方すると、結果は 7 の倍数より 2 大きい。 10² = 100 = 98 + 2 や 4² = 16 = 14 + 2 のように。

この3種類のパターンのどれになるにしても、「7 の倍数以外の整数」の平方を 7 で割った余りは 1 or 4 or 2 であり、決して 3 or 5 or 6 ではない。具体的に:
　　x² ≡ 1 (mod 7)　は　x ≡ ±1 という解を持ち、
　　x² ≡ 4 (mod 7)　は　x ≡ ±2 という解を持ち、
　　x² ≡ 2 (mod 7)　は　x ≡ ±3 という解を持つ。
他方:
　　x² ≡ 3 (mod 7)　には解がなく、
　　x² ≡ 5 (mod 7)　にも解がなく、
　　x² ≡ 6 (mod 7)　にも解がない。□

Gauß は x² ≡ a (mod p) のような2次式が解を持つかどうか、高速に判定する方法を確立したものの、具体的な解を求める方法については、あまり高速化できなかった。最初の実用的アルゴリズムは Tonelli によって記述され、20世紀後半、Shanks によって再発見・改良された。このページでは Tonelli のアルゴリズムについて、一通り紹介した。次のページでは Shanks バージョンについて検討したい。

---