Maxima で綴る数学の旅

紙と鉛筆の代わりに、数式処理システムMaxima / Macsyma を使って、数学を楽しみましょう

-数学- 4で割ると1余る素数を2つの平方数の和に分解する方法 - 互除法-

探検! 数の密林・数論の迷宮

探検! 数の密林・数論の迷宮

 

今度は、4n+1型素数を二つの平方数の和に分割するのに、互除法を使ってみます。この方法は「数論の迷宮」p283を参照してください。

\(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)\)を満たすものがただ1つ定まります。このとき明らかに\(p\)と\(a\)は互いに素なので、\(r_0=p, r_1=a\)の2数間で互除法を行うと、

\( r_i=k_i\,r_{i+1}+r_{i+2}, 0\leq r_{i+2} \lt r_{i+1} \, \left(0 \leq i \leq n-1 \right)\)

\( r_n=1, r_{n+1}=0 \)

を満たす正整数の単調減少列\(\{r_i\}_{i=0}^n\)が得られます。したがって、

\( r_{k+1} \lt r_k \lt \sqrt{p} \lt r_{k-1} \)

を満たす正整数の組\(\left(r_k,r_{k+1}\right)\)が確定します。

 

定理13.12 上記の正整数の組\(\left(r_k,r_{k+1}\right)\)について次の等式が成立する。
\(p=r_{k}^2+r_{k+1}^2\)

 

(%i1) debug:false$

まず、1ステップ互除法を行う関数eucを定義します。euc([a,b])とすると[b, mod(a,b)]が戻ってきます。大域変数debugにtrueを代入してから呼び出すと\(\frac{a}{b}\)の商が表示されます。

(%i2) euc(pair):=block([qo,rem],[qo,rem]:divide(pair[1],pair[2]),if debug then print(qo),[pair[2],rem])$

次にsearchmn()を定義します。この互除法を繰り返し呼び出します。余りの単調減少列\(\{r_i\}_{i=0}^n\)が\(\sqrt{p}\)を跨いだところで、再帰呼び出しは終了します。

(%i3) searchmn(pair,p):=block([sq],sq:ev(sqrt(p),numer), if(pair[1]>sq) and (sq>pair[2]) then euc(pair) else searchmn(euc(pair),p))$

ところでこの互除法を起動するためには-1の平方根で範囲の条件を満たす\(a\)を見つける必要があります。その関数をfind_sqrtmone(p)として定義します。

(%i4) 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))$

この定義はtsujimotter氏のブログで紹介されていた方法をmaximaに書き換えたものでした。

 

ここまでの準備が整えば、後は試してみるだけです。ここでは\(4\,n+1\)型の素数として28813を使います。

(%i5) p:28813;
$$ \tag{%o5} 28813 $$

次にmod pでの-1の平方根\(a\)を求めます。

(%i6)a:find_sqrtmone(p);
$$ \tag{%o6} 5885 $$

互除法を再帰呼び出しします。

(%i7) searchmn([p,a],p);
$$ \tag{%o7} \left[ 142 , 93 \right] $$

結果が得られました。確認のために答えのリストの最初の数と2番目の数を取り出して2乗して足してみます。

(%i8) %[1]^2+%[2]^2;

$$ \tag{%o8} 28813 $$

ちゃんと元に戻りました!!

 

この方法は実は連分数と深い関係があります。それをみるために、互除法の計算の際に現れる商を次々に表示してみます。

(%i9) searchmn([p,a],p),debug=true;
\( \tag{*} 4\)
\(  \tag{*} 1\)
\(  \tag{*} 8\)
\(  \tag{*} 1\)
\(  \tag{*} 1\)
\(  \tag{*} 1\)

$$ \tag{%o9} \left[ 142 , 93 \right] $$

一方、有理数\(\frac{p}{a}\)を連分数表示してみます。

(%i10) p/a=cfdisrep(cf(p/a));
$$ \tag{%o10} \frac{28813}{5885}=4+\frac{1}{1+\frac{1}{8+\frac{1}{1+\frac{1}{1+\frac{1}{1+\frac{1}{1+\frac{1}{1+\frac{1}{1+\frac{1}{8+\frac{1}{1+\frac{1}{4}}}}}}}}}}} $$

なんだかみたような数字が並んでいます。これをリスト表示に戻すと、

(%i11) cf(p/a);
$$ \tag{%o11} \left[ 4 , 1 , 8 , 1 , 1 , 1 , 1 , 1 , 1 , 8 , 1 , 4 \right] $$
となります、、、。最初の6つの数字は互除法の商の表示と全く一致し、そのあとの6つの数字はそれを逆から読んだものになっています!!

 

残念ながら橋本先生の本にはこの辺の事情は載っていなかったので、ネット検索してみました。その結果次の文献に行き当たりました。

Brillhart, John. "Note on Representing a Prime as a Sum of Two Squares." Mathematics of Computation 26, no. 120 (1972): 1011-013. doi:10.2307/2005889.
この互除法のアルゴリズムは元々HermiteとSerretという二人の人が独立に見つけた連分数展開ベースのアルゴリズムをBrillhartが互除法に書き直したものだったのです。上記文献の中に次のような記述があります。

  1. \(\frac{p}{a}\)の連分数展開はリスト表示するとその長さは偶数で、全体が回文になっている。
    \(\frac{p}{a}=[q_0,q_1,\cdots,q_k,q_k,\cdots,q_1,q_0]=\frac{A_{2\,k+1}}{B_{2\,k+1}}\)
  2. \(A_{2\,k+1}=p\) かつ \(A_{2\,k}=a\)
  3. \(p=A_k^2+A_{k-1}^2\)
  4. 途中打ち切りの分子\(A_k\)に関して次の再帰的な関係式が成り立つ
    \(p=q_0\,a+A_{2\,k-1},\, a=q_1\,A_{2\,k-1}+A_{2\,k-2},\cdots\)

上記1,3より、まず回文の前半を分数で表示してみます。

(%i12) cfdisrep([4,1,8,1,1,1]);
$$ \tag{%o12} 4+\frac{1}{1+\frac{1}{8+\frac{1}{1+\frac{1}{1+\frac{1}{1}}}}} $$

(%i13) rat(%);
$$ \tag{%o13} \frac{142}{29} $$

確かに、分子\(A_5\)は142で互除法で求めた数字の一つと一致しています。

さらに回文の前半を一つ短くして分数表示してみます。

(%i14) cfdisrep([4,1,8,1,1]);
$$ \tag{%o14} 4+\frac{1}{1+\frac{1}{8+\frac{1}{1+\frac{1}{1}}}} $$

(%i15) rat(%);
$$ \tag{%o15} \frac{93}{19} $$

確かに、分子\(A_4\)は93で互除法で求めたもう一つの数字と一致しています。

 

ところで今回紹介したfind_sqrtmone(p)及びsearchmn([p,a])は非常に高速です。Maximaでちょっと頑張れば300桁程度の\(4\,n+1\)型素数を作るのも簡単ですが、aをfind_sqrtmone(p)で求め、searchmn([p,a])で2つの整数(平方して和をとるとpになる)を求めるのはほとんど一瞬です。