８章　連分数とPell方程式

2．連分数展開を利用した解

intermezzo

ここまでの計算を見てほとんどの人は驚いたと思う。
そして、その驚き方には2通りあると思う。

まず第1は、解が確かに必ずあるということ、
及び n=94 の場合のように、ひじょうに大きな数になる場合があるということ。
100以下だったから求められたようなものの、1000以下だったらたぶんこの方法ではだめだろう。

しかし、pell方程式について元々知っていた人は、もっと驚いたのではないか。
つまり、こんな単純な方法で、本当に解を求めてしまった、という驚きである。
pell方程式は初期解が n に比べてひじょうに大きくなることで有名であり、
こんな方法で解が得られるのであればなんの苦労もない。
だから普通は、しらみつぶしで解をもとめようとは思わない。
特に手計算では絶対こんなことはしないだろう。

しかし、計算機は n=94 の場合の解を 221064回の繰り返しという愚直な方法で、
しかも数秒で求めてしまった。それもスーパーコンピュータではなく、
そのへんにころがっている普通のパソコンによってである。
実際、何の苦労もいらない。

まったく我々は何という力を手にいれてしまったのか。

連分数

pell方程式は、一般に連分数展開を利用して解くが、ここではまず連分数について説明する。
例として、sqrt(2)を、以下のように変形する。

sqrt(2) = 1 + sqrt(2) -1

                     1
        = 1 + ---------------
                     1
                -----------
                 sqrt(2)-1

                     1
        = 1 + ---------------
                 sqrt(2)+1
                -----------
                    2-1

                   1
        = 1 + -----------
               1+sqrt(2)

以下 sqrt(2) のところに、自分自身を代入していくと、

                      1
sqrt(2) = 1 + ----------------
                        1
             2 + -------------
                         1
                 2 + ---------

                     2 + ....

となる。このように、分子が 1の分数が重なる形の表示を連分数という。
このままでは書くのが大変なので、次のように簡略化して書くことにする。

sqrt(2) = 1 [2]
( [ ] は、かっこ内が繰り返されることを表す)

sqrt(n)の連分数展開は、必ずこのように循環的になることが知られている（例えば参考文献）。
もうひとつサンプルとして、sqrt(3) の連分数展開を求めてみる。

sqrt(3) = 1 + sqrt(3) -1

                     1
        = 1 + ---------------
                     1
                -----------
                 sqrt(3)-1

                     1
        = 1 + ---------------
                 sqrt(3)+1
                -----------
                    3-1

                       1
        = 1 + ------------------
                    sqrt(3)-1
               1 + -----------
                        2

                        1
        = 1 + ---------------------
                         3-1
               1 + --------------
                    2(sqrt(3)+1)

                       1
        = 1 + ------------------
                        1
               1 + -----------
                    1+sqrt(3)

ここでsqrt(3)が現れたので自分自身を代入する。
連分数展開は、

sqrt(3) = 1 [1,2]

となる。

Pell方程式及び連分数展開を求めるためのアルゴリズム

ここでは "Recreations in the Theory of Numbers" で解説されている解法を紹介する。
a_nを連分数展開の各桁の数、 p_n/q_nを n番目で打ち切ったときの部分和とし、
補助変数として、P_n, Q_nを導入する（大文字／小文字を区別する）。
これらの数について、以下の漸化式により、第 n 項目を計算する。

a_n a₁ = int(sqrt(n))
……
a_n = int((a₁+P_n)/Q₁)

P_n P₁ = 0
P₂ = a₁
……
P_n = a_n-1 * Q_n-1-P_n-1

Q_n Q₁ = 1
Q₂ = n-a₁²
……
Q_n = (n-P_n²) / Q_n-1

p_n p₁ = a₁
p₂ = a₁ * a₂+1
……
p_n = a_n * p_n-1+p_n-2

q_n q₁ = 1
q₂ = a₂ ……
q_n = a_n * q_n-1+q_n-2

a_n	a₁ = int(sqrt(n)) …… a_n = int((a₁+P_n)/Q₁)
P_n	P₁ = 0 P₂ = a₁ …… P_n = a_n-1 * Q_n-1-P_n-1
Q_n	Q₁ = 1 Q₂ = n-a₁² …… Q_n = (n-P_n²) / Q_n-1
p_n	p₁ = a₁ p₂ = a₁ * a₂+1 …… p_n = a_n * p_n-1+p_n-2
q_n	q₁ = 1 q₂ = a₂ …… q_n = a_n * q_n-1+q_n-2

Q_n = 1となる n について、 p_n-1, q_n-1が Pell方程式の解となる。

プログラム

 10   ' pell equation's solution by continue fraction method
 20   input N
 30   ' n=1
 40   Lp1=0:Lq1=1:A1=isqrt(N):Sp1=A1:Sq1=1
 50   ' n=2
 60   Lp2=A1:Lq2=N-A1^2:A2=int((A1+Lp2)/Lq2):Sp2=A1*A2+1:Sq2=A2
 70   if Lq2=1 then F=1:goto 150
 80   ' n
 90   Lp3=A2*Lq2-Lp2:Lq3=(N-Lp3^2)\Lq2:A3=int((isqrt(N)+Lp3)/Lq3)
100   Sp3=A3*Sp2+Sp1:Sq3=A3*Sq2+Sq1
110   if Lq3=1 then F=2:goto 150
120   Lp1=Lp2:Lp2=Lp3:Lq1=Lq2:Lq2=Lq3:Sp1=Sp2:Sp2=Sp3:Sq1=Sq2:Sq2=Sq3:A2=A3
130   goto 90
140   '
150   if F=1 then print Sp1,Sq1,Sp1^2-N*Sq1^2:end
160   print Sp2,Sq2,Sp2^2-N*Sq2^2:end

変数はそれぞれ、

Lp := P
Lq := Q
Sp := p
Sq := q

としている。

このプログラムを実行して、例えば n = 94 の場合を求めてみると、

2143295 221064 1

という正しい結果が返って来る。

前	この章の目次	次

E-mail : kc2h-msm@asahi-net.or.jp

三島　久典