橋本先生の本：
には明示的には載っていないのですが、math.stackexchange.comの関連質問記事：
やそのほかの記事で出てくる方法があります。それはガウス整数での互除法を用いた方法です。ちなみにガウス整数とは\(a, b \in Z\), \(i\)は虚数単位として\(a+b\,i\)の形の複素数を指します。ガウス整数同士の割り算が出来て、商と余りを求めることが出来ます。
この方法は、本質的には以前の記事
で書いた普通の互除法を用いる方法と同じです。ガウス整数の互除法が理論的には簡明で、その最適化が整数の互除法を使うやり方なのだと思います。
普通の互除法の時と同じように、
\(p \equiv 1 \left( mod 4\right)\)となる素数\(p\)に対して整数\(a\,\left(0\lt a \lt \frac{p}{2}\right) \)で\(a^2 \equiv -1 \left( mod p \right)\)を満たすものがただ１つ定まります。
の\(a\)を求めます。するとガウス整数としての\(p\)と\(a+i\)の最大公約数の実部の２乗と虚部の２乗を足すと\(p\)になるのです。
早速確認して見ましょう。Maximaにはガウス整数の互除法のコマンドはないのですが、ガウス整数の素因数分解はあります。一旦そちらで確認して見ましょう。
まず上記\(a\)を求める関数find_sqrtmone(p)を定義します。また使う素数は888777677としましょう。
(%i1) find_sqrtmone(p):=block([b],for i:3 thru p-1 do if jacobi(i,p)=-1 and (b:power_mod(i,(p-1)/4,p))<(p-1)/2 then return(b))$
(%i2) p:888777677;
$$ \tag{%o2} 888777677 $$
これが本当に素数でかつ4で割ると1余ることを確認して見ます。
(%i3) [primep(p),mod(p,4)];
$$ \tag{%o3} \left[ \mathbf{true} , 1 \right] $$
では\(p\)を素因数分解して見ましょう。gcfactor(z)というコマンドを使えば一発です。
(%i4) gcfactor(p);
$$ \tag{%o4} -i\,\left(12989+26834\,i\right)\,\left(26834+12989\,i\right) $$
次に\(a\)を求めてgcfactorで素因数分解します。
(%i5) a:find_sqrtmone(p);
$$ \tag{%o5} 304219384 $$
(%i6) gcfactor(a+%i);
$$ \tag{%o6} -\left(5+2\,i\right)\,\left(10+19\,i\right)\,\left(30+83\,i\right)\,\left(26834+12989\,i\right) $$
明らかに共通素因数が\(26834+12989\,i\)です。そして確かに\(26834^2+12989^2=888777677\)が成り立ちます。
\(a\)の定義から\(a^2=-1\,mod\,p\)が成り立ちます。つまり\(a^2=-1+k\,p, \,k \in Z\)ということです。\(-1\)を移項してガウス整数の範囲で素因数分解すると、\(\left(a+i\right)\,\left(a-i\right)=k\,p\)となります。\(p\)が左辺のどちらの因数も割り切らないのは自明なので、\(p\)のガウス整数としての２つの素因数が左辺の２つの素因数を割り切ることになります。従って例えば\(\left(a+i\right)\)と\(p\)の公約数が\(p\)の素因数の一つになるわけです。そして\(p\)の２つの素因数は互いに複素共役ですからその実部の２乗と虚部の２乗を足すと\(p\)に戻るのも自明です。
ここでは自前でガウス整数の互除法を実装してみます。小道具として、ガウス整数のノームを求めるnorm_gir(z), 与えられた複素数に最も近いガウス整数を求めるround_gir(z), ガウス整数同士の割り算divide_gir(z,w)（商と余りを求める）を定義します。
(%i7) norm_gir(z):=realpart(z)^2+imagpart(z)^2;
$$ \tag{%o7} \mathrm{norm\underline{\quad}gir}\left(z\right):=\mathrm{realpart}\left(z\right)^2+\mathrm{imagpart}\left(z\right)^2 $$
(%i8) round_gir(z):=round(realpart(z))+round(imagpart(z))*%i;
$$ \tag{%o8} \mathrm{round\underline{\quad}gir}\left(z\right):=\mathrm{\%round}\left(\mathrm{realpart}\left(z\right)\right)+\mathrm{\%round}\left(\mathrm{imagpart}\left(z\right)\right)\,i $$
(%i9) divide_gir(z,w):=block([quo,rem],quo:round_gir(z*conjugate(w)/norm_gir(w)),rem:ratsimp(z-quo*w),[quo,rem]);
$$ \tag{%o9} \mathrm{divide\underline{\quad}gir}\left(z , w\right):=\mathbf{block}\;\left(\left[ \mathrm{quo} , \mathrm{rem} \right] , \mathrm{quo}:\mathrm{round\underline{\quad}gir}\left(\frac{z\,w^{\star}}{\mathrm{norm\underline{\quad}gir}\left(w\right)}\right) , \mathrm{rem}:\mathrm{ratsimp}\left(z-\mathrm{quo}\,w\right) , \left[ \mathrm{quo} , \mathrm{rem} \right] \right) $$
これだけの小道具を定義すると、２つのガウス整数の最大公約数を求める関数gcd_gir(z,w)を定義するのは簡単です。
(%i10) gcd_gir(z,w):=block([quo,rem],
if norm_gir(z)<norm_gir(w) then [z,w]:[w,z],
[quo,rem]:divide_gir(z,w),
if rem=0 then w else gcd_gir(w,rem));
$$ \tag{%o10} \mathrm{gcd\underline{\quad}gir}\left(z , w\right):=\mathbf{block}\;\left(\left[ \mathrm{quo} , \mathrm{rem} \right] , \mathbf{if}\;\mathrm{norm\underline{\quad}gir}\left(z\right)<\mathrm{norm\underline{\quad}gir}\left(w\right)\;\mathbf{then}\;\left[ z , w \right] :\left[ w , z \right] , \left[ \mathrm{quo} , \mathrm{rem} \right] :\mathrm{divide\underline{\quad}gir}\left(z , w\right) , \mathbf{if}\;\mathrm{rem}=0\;\mathbf{then}\;w\;\mathbf{else}\;\mathrm{gcd\underline{\quad}gir}\left(w , \mathrm{rem}\right)\right) $$
早速、\(p\)と\(a+i\)の最大公約数を求めます。
(%i11) gcd_gir(p,a+%i);
$$ \tag{%o11} 12989\,i+26834 $$
検算として、実部と虚部をとってそれぞれ２乗して足します。
(%i12) realpart(%)^2+imagpart(%)^2;
$$ \tag{%o12} 888777677 $$
