元祖ワシ的日記

調子こいてもう一発。単なるメモという話もありますが。

4n+1型の素数（すなわち4で割ると1あまる素数）は平方数の和に出来ます。それも只一通りです。不思議なことですが、これは紛れもない事実なのです。

5=1+4
13=4+9
17=1+16

ちなみに4n+3型の素数は平方の和には出来ません。

小さい素数は手当たりしだいに計算していけばいいのですが、巨大な素数となると結構厄介です。（例えば120710677とか）

そこでこれを計算するアルゴリズムを紹介します。

まずB=1として

A^2+B^2=kp

となるような数Aを探します（ようするに二乗して1を足したらpの倍数になる数を探す）。このようなAはpが4n+1型なら必ず存在します。（平方剰余の相互法則）

もしk=1ならA^2+1^2 = pだからこれが平方和です。よって、k>1とします。 次のような u, vを探します。

u = A mod k
v = B mod k

ここでu, vは -k/2からk/2の間に入るようにとります。すると u, vはkを法としてA, Bに合同ですから、A^2+B^2 ≡ u^2+v^2 ≡ 0 (mod k)。すなわち u^2+v^2はkで割り切れることに注意します。

ここで次の恒等式を考えます。

(uA + vB)^2 + (vA - uB)^2 = (A^2 + B^2)(u^2 + v^2)

右辺はk^2で割り切れます。

また、(vA - uB) ≡ 0 (mod k)ですから、(vA - uB)^2も k^2で割り切れます。

故に恒等式が成立するためには (uA + vB)^2も k^2で割り切れなければなりません。

したがって、 (vA - uB)、(uA + vB)はともにkで割り切れます。

A' = (uA + vB)/k
B' = (vA - uB)/k
k' = (u^2 + v^2)/k
とすれば

A'2 + B'2 = k'p

u, vの条件より k' < kだから A, Bからはじめてより小さいk'を得ることができました。以下これをkが1になるまで繰り返せば平方和が求まります。

参考までに120710677の平方和をこのアルゴリズムで求めると以下の様になります。

46505545^2 + 1^2 == 17916938 * 120710677
18805958^2 + 3^2 == 2929849 * 120710677
7874929^2 + 18^2 == 513745 * 120710677
2586742^2 + 270^2 == 55432 * 120710677
866197^2 + 12690^2 == 6217 * 120710677
283914^2 + 31516^2 == 676 * 120710677
14455^2 + 107238^2 == 97 * 120710677
48346^2 + 8768^2 == 20 * 120710677
18011^2 + 16708^2 == 5 * 120710677
3081^2 + 10546^2 == 1 * 120710677

このアルゴリズムを使用せずに、これを手当たり次第で計算していくのはすごい時間がかかります。このアルゴリズムが如何に便利かがよくわかるかと。

ところで、最初の 46505545^2 + 1^2 をどうやって求めるのでしょう？これを手当たり次第に計算していくのは、偉い時間がかかります。これは原始根をうまく使うと簡単に求められます。

120710677(=p)の原始根は5ですから 5^(p-1) ≡ 1 (mod p) 120710676乗すれば1なのだから、その半分すなわち60355338乗したときの剰余は-1か1です。しかし5は原始根だから (p-1)/2乗したときの剰余は-1。故にその半分 (p-1)/4乗した数が求めるAです。 なぜなら A^2 ≡ (5^((p-1)/4))^2 ≡ 5^((p-1)/2) ≡ -1 (mod p) 故に A^2 + 1 ≡ 0 (mod p) ここでpが4n+1型というのが効いてきます。なぜなら4n+1型でなければ、(p-1)/4は整数ではないからです。 最後に原始根をどうやって求めるのか？という問題もありますが、 これは手当たり次第に計算していくしかありません。ただそれも、全ての冪(2～p-1)をチェックする必要はなく、 p-1の約数乗だけチェックすればいいのです。この場合ですと p-1 = 2^2 3^1 17^2 34807^1 だから36通りのべき乗で1に合同になるかどうかを確認すればよいことになります。 (もっともp-1は絶対1に合同になるし、1乗もあまり意味がないですから、34通りについてどれも1に合同にならないことを示せばOKです） これでチェックが120710676通りから35通りに減らすことができます。

参考
二個の平方数の和 - Wikipedia

追記
A^2 + 1 ≡ 0 (mod p)

の求め方ですが原始根をワザワザ求める必要はありません。適当に数を選んで(p-1)/4乗してみます。すると1/2の確率でAになります。

まず半数の数は(p-1)/2乗すると1、もう半数は-1になります。
-1になる数の平方根は(p-1)/4乗した数です。これがAになります。
確率的なアルゴリズムですが、p-1が大きな約数を持つときはこちらのほうが圧倒的に簡単です。