をお持ちの方はp281に記載されている定理13.10をご覧になってください。
定理13.10 \(p \equiv 1 \left(mod \,4\right)\)となる素数を\(p=4\,n+1\)と表す。このとき整数\(a,b \left(-\frac{p}{2}<a, \,b<\frac{p}{2}\right)\)を
\(a\equiv\frac{{{2\,n}\choose{n}}}{2}, \,b\equiv a\,\left(2\,n\right)! \left(mod\,p\right)\)
を満たすようにとるとき\(p=a^2+b^2\)が成立する。
この定理、短いし、二項係数と階乗しか出てこないし、、、もっと広く世間に知られても良いと思うのですが、大抵の整数論の本には載っていません。その理由は後回しにして、この定理が述べている\(a,\,b\)の求め方をMaximaで実際にいくつかの\(4\,n+1\)型素数に対して確認して見ましょう。
まず準備として\(gnorm\left(k,p\right)\)という関数を定義します。\(k\,mod\,p\)という数を、\(-\frac{p}{2}\)から\(\frac{p}{2}\)の範囲に調整するためのものです。
(%i1) gnorm(k,p):=block([m:mod(k,p)], if m>(p-1)/2 then m-p else m)$
では早速実際にやってみます。\(10037\)は4で割ると1余る素数です。
(%i2) p:10037;
$$ \tag{%o2} 10037 $$
次に\(p=4\,n+1\)の\(n\)を求めます。
(%i3) n:(p-1)/4;
$$ \tag{%o3} 2509 $$
後は定理の通りにやれば良いのです。aを二項係数を使って求め、gnormを使って範囲を調整します。
(%i4) a:gnorm(binomial(2*n,n)/2,p);
$$ \tag{%o4} 89 $$
bを階乗を使って求め、gnormを使って範囲を調整します。
(%i5) b:gnorm((2*n)!*a,p);
$$ \tag{%o5} -46 $$
最後に\(a,\,b\)を平方して足すと\(p\)に戻ることを確認します。
(%i6) a^2+b^2;
$$ \tag{%o6} 10037 $$
もう一つくらいやってみましょう。
(%i7) p:28813;
$$ \tag{%o7} 28813 $$
(%i8) n:(p-1)/4;
$$ \tag{%o8} 7203 $$
(%i9) a:gnorm(binomial(2*n,n)/2,p);
$$ \tag{%o9} 93 $$
(%i10) b:gnorm((2*n)!*a,p);
$$ \tag{%o10} 142 $$
(%i11) a^2+b^2;
$$ \tag{%o11} 28813 $$
ではなぜこの定理は紹介されないのか、、、おそらく証明も定理の成り立つ背景の説明も面倒だからなのだと思います。橋本先生の本では「この定理のガウスによる証明は、初等的であっても決して単純ではなく、その本質を理解することはむしろ困難です」と述べ、ガウスによる証明自身は紹介されていません。
そして、その代わりに素数\(p=4\,n+1\)に対して\(F_p\)上の楕円曲線\(E : y^2=x^3-x\)のHasse不変量\(\hat{H}_p\left(E\right)\)を考えると、
\(\hat{H}_p\left(E\right)=\left(-1\right)^n\,\binom{2\,n}{n}\)
が成立すること、\(F_p\)上の楕円曲線\(E\)の解の個数を\(N_p\left(E\right)\)として、
\(N_p\left(E\right)-1-p\equiv\hat{H}_p\left(E\right) \left(mod\,p\right)\)
が成立すること、この左辺が実はヤコブスタール和と一致することなどからこの二項係数が\(a\)を計算していることを説明しています。
問題もアルゴリズムも超初等的なのに、その背景や証明に楕円曲線やハッセ不変量が出てくるようでは初等整数論の本では説明できないでしょう。
ここで紹介したMaximaのコードは非常に簡単ですが、遅いです。階乗を求める際の乗算を\(mod\,p\)で計算するように変更し、二項係数もその階乗\(mod\,p\)で定義しなおせばずっと早くなると思います。