Maximaで学ぶ連分数

数学の具体的な計算にMaximaを使って、数学もMaximaも同時に学んでしまいましょう。今回はMaximaを使って連分数の計算をしてみたいと思います。計算自体はシンプルなものですが、桁数が大きくなるにつれ、また累乗根などを含むにつれ、手計算では大変になっていきます。実装自体は簡単ながらも、計算機の便利さを実感しやすい題材ではないかと思います。

連分数

連分数

連分数（れんぶんすう、continued fraction）とは、分数の分母のなかに更に分数が含まれているような数のことを指します。より詳しい説明についてはWikipediaなどを参照してください：

ja.wikipedia.org

ここでは具体的な有理数、無理数の連分数を求める方法をMaximaを使って実演してみます。最終的には、インプットの数を適当に変えることで、任意の有理数、無理数に対する連分数展開を求めるMaximaのプログラムを与えます。

具体例（有理数 8/5 の連分数展開）

有理数 $8/5$ を整数部分とそれ以外に分解してみます。

$\qquad\displaystyle{\frac{8}{5}=1+\frac{3}{5}}$

さて、右辺第２項の有理数 $3/5$ ですが、これを書き換えると

$\qquad\displaystyle{\frac{3}{5}=\frac{1}{\frac{5}{3}}}$

です。この右辺の分母に現れた有理数 $5/3$ を、上でやったように整数部分とそれ以外に分解してみると

$\qquad\displaystyle{\frac{5}{3}=1+\frac{2}{3}}$

ですね。つまり

$\qquad\displaystyle{\frac{3}{5}=\frac{1}{\frac{5}{3}}=\frac{1}{1+\frac{2}{3}}}$

と書けることがわかります。同じ手順を踏んで

$\qquad\displaystyle{\frac{2}{3}=\frac{1}{\frac{3}{2}}=\frac{1}{1+\frac{1}{2}}}$

ですから、結局

$\qquad\displaystyle{\frac{8}{5}=1+\frac{1}{1+\frac{1}{1+\frac{1}{2}}}}$

と表せることがわかります。

連分数を求める手順

以上の操作をMaximaで実装するために、少し抽象的に手順を追っていきます。まず出発点となる数 $f_0$ を整数部分 $I_0$ とそれ以外の部分 $p_0$ に分解します。

$\qquad\displaystyle{f_0=I_0+p_0}$

そして右辺に現れた数 $p_0$ の逆数を

$\qquad\displaystyle{f_1=\frac{1}{p_0}}$

として、再び $f_1$ を整数部分 $I_1$ とそれ以外の部分 $p_1$ に分解します。

$\qquad\displaystyle{f_1=I_1+p_1}$

この手順を同様に繰り返して数列 $\{f_0,f_1,f_2,\dots\},\ \{I_0,I_1,I_2,\dots\},\ \{p_0,p_1,p_2,\dots\}$ が得られます。この手順は最終的に（有理数の場合には） $p_{n+1}=0$ となるまで続きます。このとき分数 $p_n$ は

$\qquad\displaystyle{p_n=\frac{1}{I_{n+1}}}$

で与えられます。なぜなら $p_{n+1}=0$ から

$\qquad\displaystyle{0=p_{n+1}=f_{n+1}-I_{n+1}=\frac{1}{p_n}-I_{n+1}}$

となるためです。すると数 $f_0$ の連分数が以下のように与えられます：

$\qquad\displaystyle{\begin{aligned}f_0&=I_0+\frac{1}{f_1}\\ &=I_0+\frac{1}{I_1+\frac{1}{f_2}}\\&=\cdots\\&=I_0+\frac{1}{I_1+\frac{1}{I_2+\frac{1}{I_3+\cdots}}}\end{aligned}}$

ただし最後の部分は

$\qquad\displaystyle{I_n+\frac{1}{I_{n+1}}}$

という分数で終わる形をしていて、数列 $\{I_0,I_1,I_2,\dots,I_{n+1}\}$ の情報から求めることができます。

Maximaによる実装

まず手始めに、与えらた数に対して、整数部分とそれ以外とに分解する方法をみてみます。上の具体例に沿って $f_0=8/5$ としてみましょう。まず、入力された数に対して整数部分とそれ以外を求める関数を定義します。

kill(all)$
Int(x):=fix(x)$     /*  fix は Maxima の組み込み関数で、ガウス記号です */
P(x):=x-Int(x)$

を使って

input_num: 8/5$
f[0]: input_num$
I[0]: Int(f[0]);
p[0]: P(f[0]);

より

$\qquad\displaystyle{I_0=1 \qquad\quad p_0=\frac{3}{5}}$

です。ここから

$\qquad\displaystyle{f_1=\frac{1}{p_0}=\frac{5}{3}}$

と求まり、この操作を再び $f_1=5/3$ に施して、それを逐次的に繰り返せばよいことがわかります。すなわち

N: 10$         /*  繰り返し操作の上限回数  */
n:1$ while n < N do(
f[n]: 1/P(p[n-1]),
I[n]: Int(f[n]),
p[n]: P(f[n]),
n: if p[n] = 0    /*  繰り返し操作の終了  */
  then block(n_max: n, n+2*N) 
  else n+1        /*  繰り返し操作の継続  */
);

とし、得られた数列 $I=\{I_0,I_1,I_2,\dots,I_{n_{max}}\}$ を表示すると

makelist(I[n],n,0, n_max);

により

$\qquad\displaystyle{I=\{1,1,1,2\}}$

であることがわかります。よってこの数列の情報から $f_0=8/5$ の連分数が

$\qquad\displaystyle{\frac{8}{5}=1+\frac{1}{1+\frac{1}{1+\frac{1}{2}}}}$

と求まります。当然ですが、これは上の具体的な説明で求めた結果に一致します。

Maximaの組み込み関数

さて、実はわざわざ自分で上記のようなプログラムを組まなくても、Maximaには組み込み関数として、上記の数列 $I$ を求めるコマンドが用意されています：

cf(8/5);

から、リスト

$\qquad\displaystyle{[1,1,1,2]}$

が返ってきます。ただ、この計算出力を正しく利用する上で、上記のように計算手順を自分で実装し、理解しておくことも大切なことではないかなと思います。

具体例（無理数の連分数展開）

上の例では有理数を考えましたが、連分数は無理数に対しても求めることができます。ここでは具体的な無理数として $\sqrt{2}$ を考えて、やはり整数部分とそれ以外に分解してみます。

$\qquad\displaystyle{\sqrt{2}=1+(\sqrt{2}-1)}$

ここで右辺第２項を書き換えると

$\qquad\displaystyle{\sqrt{2}-1=\frac{1}{\frac{1}{\sqrt{2}-1}}=\frac{1}{\frac{\sqrt{2}+1}{(\sqrt{2}-1)(\sqrt{2}+1)}}=\frac{1}{1+\sqrt{2}}}$

です。少し面倒なのが有理化の操作が入っているという点です。さて、右辺の分母内にある $\sqrt{2}$ を再び整数部分とそれ以外に分解して

$\qquad\displaystyle{\sqrt{2}-1=\frac{1}{2+(\sqrt{2}-1)}}$

がわかります。ここで $t= \sqrt{2}-1$ とおくと、この関係式は

$\qquad\displaystyle{t=\frac{1}{2+t}}$

という形をしていることがわかります。分母を払って移項すると、関係式は結局

$\qquad\displaystyle{t^2+2t-1=0}$

という $t$ の２次方程式になるわけですが、解の公式から

$\qquad\displaystyle{t=-1\pm\sqrt{2}}$

その解の一つが $t=\sqrt{2}-1$ であることがわかります。さて、連分数に戻って、上記の関係式を逐次的に用いると

$\qquad\displaystyle{\sqrt{2}-1=\frac{1}{2+(\sqrt{2}-1)}=\frac{1}{2+\frac{1}{2+(\sqrt{2}-1)}}=\frac{1}{2+\frac{1}{2+\frac{1}{2+(\sqrt{2}-1)}}}=\cdots}$

という連分数が求まることがわかります。ただし、有理数の場合とは異なり、この操作は無限に続くことになります。

これをMaximaで実装してみます。といっても、有理化の処理と、操作が打ち切りにならないことに起因する処理を加えるだけです：

kill(all)$
input_num: sqrt(2)$
N: 10$         /*  繰り返し操作の上限回数  */
algebraic: ture$  /*  有理化を処理するためのフラグ  */
Int(x):=fix(x)$     /*  fix は Maxima の組み込み関数で、ガウス記号です */
P(x):=x-Int(x)$
f[0]: input_num$
I[0]: Int(f[0])$
p[0]: P(f[0])$
n:1$ while n < N do(
f[n]: rat(1/P(p[n-1])),  /*  有理化の処理  */
I[n]: Int(f[n]),
p[n]: P(f[n]),
n: if p[n] = 0    /*  繰り返し操作の終了  */
  then block(n_max: n, n+2*N) 
  else 
     if n = N-1 then n+1                    /*  繰り返し操作の継続  */
     else block(n_max: N-1, n +1)   /*  n_max の取り扱い */
)$
makelist(I[n],n,0, n_max);

から数列が

$\qquad\displaystyle{I(\sqrt{2})=[1,2,2,2,2,2,2,2,2,2,2]}$

と求まります。繰り返し操作の上限 N をもっと大きくしても、2が続いていくことがわかります（実際に試してみてください）。ちなみに、Maximaの実装コマンド

cf(sqrt(2));

を使うと、リストとして

$\qquad\displaystyle{[1,2]}$

が返ってきます。他の場合として、例えば input_num = $7^{1/2},\ 3^{1/3}$ などで試してみると、Maximaの組み込み関数コマンドは２乗根までしかうまく対応していないことがわかります。実際

cf(7^(1/2));
cf(3^(1/3));

ですが、今回実装したプログラムでは（ここでは x = input_num , N = Nmax を引数として数列のリストを返す関数として Cf を定義しました）

kill(all)$
algebraic: ture$ 
Int(x):=fix(x)$
P(x):=x-Int(x)$
Cf(x,N):= block(
   f[0]: x,
   I[0]: Int(f[0]),
   p[0]: P(f[0]),
   n:1,
      while n < N do(
      f[n]: rat(1/P(p[n-1])),
      I[n]: Int(f[n]),
      p[n]: P(f[n]),
      n: if p[n] = 0 then block(n_max: n, n+2*N) 
          else 
             if n = N-1 then n+1 
             else block(n_max: N-1, n +1)
),
makelist(I[n],n,0, n_max)
)$

として

Cf(7^(1/2), 20);
Cf(3^(1/3), 100);

となり、今回の記事で作成したプログラムのほうがより汎用性の高い実装になっていることがわかります。この実装にさまざまな値を入れてみると、２乗根については同じ配列の数字が周期性をもって出現する一方、３乗根、４乗根などにはそのような性質がみられないことが観察できます。

キーワード：Maxima、連分数、代数学、方程式

pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts