Algorithms with Python

はじめに

Python は標準ライブラリ fractions を import すると、有理数 (分数, fraction) を使うことができます。今回は fractions の基本的な使い方を説明し、分数に関する簡単なプログラムを Python で作ってみましょう。

●分数の生成

Python の場合、分数はクラス Fraction のオブジェクトで、コンストラクタ Fraction() で生成します。

1. Fraction(分子, 分母)
2. Fraction(文字列)
3. Fraction(実数)

1 の場合、分母を省略すると 1 になります。引数を省略すると、値 0 の分数が生成されます。2 の場合、文字列を分数に変換します。引数の文字列は、'分子/分母' の形式や実数 (浮動小数点数) の形式になります。3 は実数を分数に変換します。ただし、実数は不正確な数なので、正確に変換できるとは限らないことに注意してください。

簡単な例を示しましょう。

>>> from fractions import Fraction
>>> Fraction(1, 2)
Fraction(1, 2)
>>> Fraction(-10, 20)
Fraction(-1, 2)
>>> Fraction(10, -15)
Fraction(-2, 3)
>>> Fraction(123)
Fraction(123, 1)
>>> Fraction()
Fraction(0, 1)

>>> Fraction('1/100')
Fraction(1, 100)
>>> Fraction('1.1')
Fraction(11, 10)
>>> Fraction('0.5')
Fraction(1, 2)
>>> Fraction('0.1')
Fraction(1, 10)
>>> Fraction('1e-2')
Fraction(1, 100)

>>> Fraction(1.1)
Fraction(2476979795053773, 2251799813685248)
>>> Fraction(0.5)
Fraction(1, 2)
>>> Fraction(0.1)
Fraction(3602879701896397, 36028797018963968)
>>> Fraction(1e-2)
Fraction(5764607523034235, 576460752303423488)

1.1 は分数で表すと 11/10 ですが、Python の実数 (規格 IEEE 754) では 1.1 を正確に表すことができません。このため、Fraction(1.1) は Fraction(11, 10) にはならず、その値に近い分数に変換されます。0.5 は実数でも正確に表すことができるので、Fraction(0.5) は Fraction(1, 2) に変換されます。Fraction() に文字列を渡すと正確に変換できるようです。

分数の分子は numerator で、分母は denominator で取得することができます。

>>> Fraction(1, 2).numerator
1
>>> Fraction(1, 2).denominator
2

関数 print() は 分子/分母 の形式で分数を表示します。

>>> print(Fraction(1, 2))
1/2
>>> print(Fraction(-1, 2))
-1/2

●分数の計算

分数の四則演算や累乗は演算子 +, -, *, /, ** で行うことができます。

>>> a = Fraction(1, 2)
>>> b = Fraction(1, 3)
>>> a + b
Fraction(5, 6)
>>> print(a + b)
5/6
>>> print(a - b)
1/6
>>> print(b - a)
-1/6
>>> print(a * b)
1/6
>>> print(a / b)
3/2
>>> print(b / a)
2/3

>>> print(a ** 2)
1/4
>>> print(b ** 3)
1/27
>>> 2 ** a
1.4142135623730951
>>> 2 ** b
1.2599210498948732

比較演算子による比較も行えます。

>>> a < b
False
>>> a > b
True
>>> a == a
True
>>> a == b
False

●分数の変換

分数は関数 int() や float() を使って整数や実数に変換することができます。

>>> int(Fraction(3/2))
1
>>> int(Fraction(4/2))
2
>>> int(Fraction(5/2))
2
>>> int(Fraction(6/2))
3
>>> float(Fraction(3/2))
1.5
>>> float(Fraction(4/2))
2.0
>>> float(Fraction(5/2))
2.5
>>> float(Fraction(6/2))
3.0

limit_denominator(max_denominator=1000000) は分母が高々 max_denominator である、もっとも self に近い分数を返します。このメソッドは実数の有理数近似を求めたり、実数で表された分数を元に戻すときに役に立ちます。

>>> import math
>>> math.pi
3.141592653589793
>>> Fraction(math.pi).limit_denominator(1000)
Fraction(355, 113)
>>> Fraction(math.pi).limit_denominator(100000)
Fraction(312689, 99532)

>>> math.sqrt(2)
1.4142135623730951
>>> Fraction(math.sqrt(2)).limit_denominator(10)
Fraction(7, 5)
>>> Fraction(math.sqrt(2)).limit_denominator(100)
Fraction(140, 99)
>>> Fraction(math.sqrt(2)).limit_denominator(1000)
Fraction(1393, 985)

>>> Fraction(1.1).limit_denominator()
Fraction(11, 10)
>>> Fraction(1.2)
Fraction(5404319552844595, 4503599627370496)
>>> Fraction(1.2).limit_denominator()
Fraction(6, 5)

●小町分数

それでは簡単な例題として、分数を使ってパズルを解いてみましょう。

このパズルの元ネタは N = 1 の場合で、参考文献 1 に掲載されています。ちなみに、3 つの分数の和が整数になる場合、その値は 1 しかありません。

プログラムは次のようになります。

リスト : 小町分数の解法

from fractions import Fraction
from itertools import permutations

# 小町分数
def komachi_solver1():
    for (a,b,c,d,e,f,g,h,i) in permutations([1,2,3,4,5,6,7,8,9], 9):
        n = Fraction(a, b * 10 + c) + Fraction(d, e * 10 + f) + Fraction(g, h * 10 + i)
        if n.numerator == 1 and a < d < g:
            print('{}/{}{} + {}/{}{} + {}/{}{} = {}'.format(a,b,c,d,e,f,g,h,i,n))

アルゴリズムは単純な生成検定法です。標準ライブラリ itertools にある関数 permutations() を使って順列を生成し、Fraction() で分数を作って足し算します。その結果 n が条件を満たしていれば、print() で解を表示します。分子の数字を a < d < g と限定することで、重複解を生成しないように工夫しています。

実行結果は次のようになります。

>>> komachi_solver1()
1/24 + 3/56 + 7/98 = 1/6
1/26 + 5/39 + 7/84 = 1/4
1/32 + 5/96 + 7/84 = 1/6
1/38 + 2/95 + 4/76 = 1/10
1/48 + 5/32 + 7/96 = 1/4
1/56 + 3/72 + 9/84 = 1/6
1/96 + 5/32 + 7/84 = 1/4
1/96 + 5/48 + 7/32 = 1/3
2/18 + 5/63 + 7/49 = 1/3
2/19 + 4/57 + 6/38 = 1/3
3/27 + 6/54 + 9/81 = 1/3
3/48 + 5/16 + 9/72 = 1/2
3/54 + 6/72 + 9/81 = 1/4
5/34 + 7/68 + 9/12 = 1

結果は全部で 14 通りになりました。

●小町分数 (2)

もうひとつ「小町分数」を出題しましょう。

このパズルの元ネタも値が 1 になる場合で、参考文献 1 に掲載されています。この問題で値が 1 / N (N は整数) になる場合は 1 と 2 の 2 通りしかないようです。

プログラムは「小町分数」を参考にすれば、簡単に作成することができます。

リスト : 小町分数 (2) の解法

def komachi_solver2():
    for (a,b,c,d,e,f,g,h,i) in permutations([1,2,3,4,5,6,7,8,9], 9):
        n = Fraction(a, b * c) + Fraction(d, e * f) + Fraction(g, h * i)
        if n.numerator == 1 and a < d < g and b < c and e < f and h < i:
            print('{}/({}*{}) + {}/({}*{}) + {}/({}*{}) = {}'.format(a,b,c,d,e,f,g,h,i,n))

>>> komachi_solver2()
1/(2*4) + 5/(3*6) + 7/(8*9) = 1/2
1/(3*6) + 5/(8*9) + 7/(2*4) = 1

●単位分数の和

パズルではありませんが、分数のお話を紹介します。分子が 1 の分数を「単位分数」といいますが、単位分数の和は古代エジプト人がよく研究していたそうです。この話は M.Kamada さん からお聞きしたのですが、参考文献 2 に「分数を単位分数の和で表す方法」がありましたので紹介します。

\(0 \lt \frac{m}{n} \lt 1\) の分数を単位分数の和で表します。まず、\(\frac{n}{m}\) の商 \(q\) を求めます。もし、割り切れれば単位分数になりますね。そうでなければ、\(\frac{m}{n}\) から \(\frac{1}{q + 1}\) を引き算して \(\frac{M}{N}\) を求めます。あとは、\(\frac{M}{N}\) に対して同じ処理を繰り返すだけです。次の式を見てください。

\(0 \lt n \bmod m \lt m\) ですから、\(M\) は必ず \(m\) より小さな値になります。つまり、この処理を繰り返していくと \(m\) は必ず 1 になるので、分数を単位分数の和で表すことができる、というわけです。なるほど納得のアルゴリズムですね。たとえば、11 / 13 を単位分数の和で表してみましょう。

11 / 13 => q = 1, 11 / 13 - 1 / 2 = 9 / 26
 9 / 26 => q = 2,  9 / 26 - 1 / 3 = 1 / 78
=> 11 / 13 = 1 / 2 + 1 / 3 + 1 / 78

このように、分子 \(m\) の値は減少していきます。このアルゴリズムを Python でプログラムすると、次のようになります。

リスト : 分数を単位分数の和で表す

def bunsu(m, n):
    if n % m == 0:
        print('1/{}'.format(n // m))
    else:
        q = n // m + 1
        print('1/{} +'.format(q), end=' ')
        bunsu(m * q - n, n * q)

関数 bunsu() は再帰定義を使ってプログラムしています。関数名は適当なので、ふさわしい名前に変更してください。あとは、アルゴリズムをそのままプログラムしただけなので、特に難しいところはないでしょう。それでは実行してみましょう。

>>> bunsu(11, 13)
1/2 + 1/3 + 1/78
>>> bunsu(12, 13)
1/2 + 1/3 + 1/12 + 1/156
>>> bunsu(19, 23)
1/2 + 1/4 + 1/14 + 1/215 + 1/138460

このほかにも、単位分数の和で表す方法は何通りもあるわけで、この方法はその中のひとつにすぎません。古代エジプトではどのような方法で求めたのでしょうか。興味深いところです。

●連分数

連分数 (continued fraction) は、分数の分母にさらに分数が含まれる形をした分数のことです。

連分数の中で、分子 \(c_1, c_2, \ldots\) がすべて 1 で、\(b_0\) が整数、\(b_1, b_2, \ldots\) がすべて正の整数であるような連分数を「正則連分数」といいます。正則連分数は \(b_i\) の情報を保持するだけで、連分数に復元することができます。そこで、以下の右辺のように記述することが多いようです。

●有理数の連分数展開

有理数 (分数) は「ユークリッドの互除法」を使って正則連分数に展開することができます。

\(a, b\) の問題を \(b, r\) の問題に帰着させるのは「ユークリッドの互除法」と同じです。この計算を繰り返し行うと、\(a, b\) はどんどん小さくなっていき、\(r = 0\) になったとき連分数展開が終了します。

簡単な例を示します。

Python で検算してみましょう。

>>> 1 + Fraction(1, 1 + Fraction(1, 1 + Fraction(1, 5)))
Fraction(17, 11)
>>> 3 + Fraction(1, 7 + Fraction(1, 16))
Fraction(355, 113)

連分数に正しく展開されていますね。

●プログラムの作成

プログラムは Fraction を使うと簡単です。次のリストを見てください。

リスト : 有理数を正則連分数に展開する

def cont_frac(n):
    if not isinstance(n, Fraction):
        raise TypeError('must be Fraction number')
    a = n.numerator
    b = n.denominator
    xs = []
    while b > 0:
        xs.append(a // b)
        a, b = b, a % b
    return xs

関数 cont_frac() は引数を正則連分数展開します。最初に引数 n が Fraction か isinstance() でチェックします。そうであれば、分子を変数 a に、分母を変数 b にセットします。展開した結果はリスト xs にセットします。あとは、ユーグリッドの互除法と同じです。b が 0 より大きければ、a // b の結果を xs に追加し、a を b に、b を a % b に更新します。最後に、return で xs を返します。

リスト : 正則連分数を有理数に戻す

def cont_frac_reduce(xs):
    i = len(xs) - 1
    a = xs[i]
    while i > 0:
        a = xs[i-1] + Fraction(1, a)
        i -= 1
    return a

# 畳み込み (reduce) を使う
def cont_frac_reduce1(xs):
    return reduce(lambda a, x: x + Fraction(1, a), reversed(xs))

(defun cont-frac-reduce1 (xs)
  (reduce (lambda (x a) (+ x (/ a))) xs :from-end t))

関数 cont_frac_reduce() は引数の正則連分数を通常の分数に直します。アルゴリズムは簡単です。変数 a を xs の最後の要素 xs[i] に初期化します。あとは、後ろから順番に a の値を xs[i - 1] + 1 / a に更新していくだけです。これで、連分数を普通の分数に戻すことができます。関数 cont_frac_reduce1() はモジュール functools の関数 reduce (畳み込み) を使ったバージョンです。

もう一つ、簡単な方法を紹介しましょう。参考 URL 4 によると、連分数 \([a_0; a_1, a_2, \ldots, a_{n-2}, a_{n-1}]\) を約分した有理数 \(\frac{p_n}{q_n}\) の分子と分母 \(p_n, q_n\) は、以下に示す漸化式で求めることができるそうです。

Python でプログラムすると、次のようになります。

リスト : 連分数を有理数に約分する (2)

# 漸化式を使う
def cont_frac_sub(xs, a0, a1):
    _, a = reduce(lambda a, x: (a[1], x * a[1] + a[0]), xs[1:], (a0, a1))
    return a

def cont_frac_reduce2(xs):
    p = cont_frac_sub(xs, 1, xs[0])
    q = cont_frac_sub(xs, 0, 1)
    return Fraction(p, q)

関数 cont_frac_sub() は漸化式を計算します。これも畳み込み reduce を使うと簡単です。関数 cont_frac_reduce2() は cont_frac_sub() を 2 回呼び出して分子 p と分母 q を求め、Fracttion(p, q) を生成して返すだけです。

簡単な実行例を示します。

>>> cont_frac(Fraction(17, 11))
[1, 1, 1, 5]
>>> cont_frac(Fraction(-17, 11))
[-2, 2, 5]
>>> cont_frac(Fraction(355, 113))
[3, 7, 16]

>>> cont_frac_reduce([1,1,1,5])
Fraction(17, 11)
>>> cont_frac_reduce([-2,2,5])
Fraction(-17, 11)
>>> cont_frac_reduce([3,7,16])
Fraction(355, 113)

>>> cont_frac_reduce1([1,1,1,5])
Fraction(17, 11)
>>> cont_frac_reduce2([1,1,1,5])
Fraction(17, 11)
>>> cont_frac_reduce1([-2,2,5])
Fraction(-17, 11)
>>> cont_frac_reduce2([-2,2,5])
Fraction(-17, 11)
>>> cont_frac_reduce1([3,7,16])
Fraction(355, 113)
>>> cont_frac_reduce2([3,7,16])
Fraction(355, 113)

>>> cont_frac_reduce([1,2,3,4,5,6,7,8])
Fraction(81201, 56660)
>>> cont_frac(Fraction(81201, 56660))
[1, 2, 3, 4, 5, 6, 7, 8]

>>> cont_frac_reduce1([1,2,3,4,5,6,7,8])
Fraction(81201, 56660)
>>> cont_frac_reduce2([1,2,3,4,5,6,7,8])
Fraction(81201, 56660)

正常に動作していますね。

●無理数の連分数展開

有理数の連分数展開は有限回で終了しますが、無理数の連分数展開は無限に続きます。たとえば、\(\sqrt 2\) の連分数展開を求めてみましょう。

係数が有理数の二次方程式の根を「二次無理数」といいます。二次無理数の正則連分数展開は循環することが知られています。

実際に試してみましょう。

>>> math.sqrt(2)
1.4142135623730951
>>> x = cont_frac_reduce([1,2,2,2,2,2,2,2,2])
>>> x
Fraction(1393, 985)
>>> float(x)
1.4142131979695431

>>> math.sqrt(3)
1.7320508075688772
>>> x = cont_frac_reduce([1,1,2,1,2,1,2,1,2])
>>> x
Fraction(265, 153)
>>> float(x)
1.7320261437908497

>>> math.sqrt(5)
2.23606797749979
>>> x = cont_frac_reduce([2,4,4,4,4,4,4,4,4])
>>> x
Fraction(219602, 98209)
>>> float(x)
2.236067977476606

>>> math.sqrt(7)
2.6457513110645907
>>> x = cont_frac_reduce([2,1,1,1,4,1,1,1,4])
>>> x
Fraction(590, 223)
>>> float(x)
2.6457399103139014

●実数の連分数展開

実数 \(\alpha\) の正則連分数展開は次のアルゴリズムで求めることができます。

プログラムは次のようになります。

リスト : 実数を正則連分数に展開する

def cont_frac_flt(x, k = 8, e = 1e-16):
    xs = []
    for _ in range(k):
        xs.append(floor(x))
        y = x - xs[-1]
        if abs(y) < e: break
        x = 1 / y
    return xs

引数 x が連分数展開を行う実数、k が段数、e が 0 と判定する基準値です。実数は不正確な数なので、演算誤差が発生します。k の値を増やしても、規則性は無限に続かずに途中で失われることがあります。

また、0.5 や 1.25 など実数でも正確に表現できる場合、途中で y の値は 0 になるので、そこで展開を終了します。ただし、演算誤差が生じる場合があるので、ぴったり 0 になるとは限りません。そこで、\(|y| \lt e\) よりも小さくなったならば、0 と判定して展開を終了します。

それでは実際に試してみましょう。

>>> import math
>>> cont_frac_flt(math.sqrt(2))
[1, 2, 2, 2, 2, 2, 2, 2]
>>> cont_frac_flt(math.sqrt(3))
[1, 1, 2, 1, 2, 1, 2, 1]
>>> cont_frac_flt(math.sqrt(5))
[2, 4, 4, 4, 4, 4, 4, 4]
>>> cont_frac_flt(math.sqrt(7))
[2, 1, 1, 1, 4, 1, 1, 1]
>>> cont_frac_flt(math.sqrt(7), k=16)
[2, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1]

>>> cont_frac_flt(0.5)
[0, 2]
>>> cont_frac_flt(0.25)
[0, 4]
>>> cont_frac_flt(0.125)
[0, 8]
>>> cont_frac_flt(1.1)
[1, 9, 1, 112589990684262, 2, 2, 6, 46912496118442]
>>> cont_frac_flt(1.1, e=1e-14)
[1, 9, 1]
>>> cont_frac_reduce([1,9,1])
Fraction(11, 10)

正常に動作しているようです。興味のある方はいろいろ試してみてください。

●平方根の連分数展開

正の整数 m の平方根は二次無理数なので連分数は循環します。ただし、循環節が長くなると浮動小数点数の演算誤差により、上記プログラムでは正確に連分数展開できない場合があります。参考 URL 6 には、平方根の連分数展開 (循環節) を正確に求めるプログラム (Python) が掲載されています。ここで、簡単にアルゴリズムを紹介しましょう。

二次無理数 \(\alpha\) を次式のように連分数展開するとします。

このとき、\(\alpha_n\) を整数 \(a_n, b_n\) を使って以下のように表すと、\(\alpha\) の連分数展開を \(a_n, b_n\) の漸化式で表すことができます。

\(a_n, b_n\) の漸化式は以下のように求めることができます。

参考 URL 6 によると、\(c_{n+1}\) を求める計算は、\( \left\lfloor \frac{c_0 + b_n}{a_{n+1}} \right\rfloor \) に置き換えることができるそうです。証明もなされているので、興味のある方はお読みくださいませ。

アルゴリズムの説明に合わせて、参考 URL 6 のプログラムを改変すると、次のようになります。

リスト : 平方根の連分数展開

def cont_frac_root(m):
    c0 = floor(sqrt(m))
    cf = [c0]
    if c0 * c0 != m:
        a, b, c = 1, 0, c0
        while True:
            b = c * a - b
            a = (m - b * b) // a
            c = (c0 + b) // a
            cf.append(c)
            if a == 1: break
    return cf

連分数展開の結果はリスト (1 次元配列) cf に格納します。漸化式の変数 \(a_n, \ b_n, \ c_n\) を変数 a, b, c で表しています。最初に floor(sqrt(m)) で m の平方根の整数部を求めて変数 c0 にセットし、cf を [c0] に初期化します。c0 の二乗が m と等しい場合、m は平方数なので cf をそのまま返します。

m が平方数でなければ、\(\sqrt m\) を漸化式で連分数に展開します。変数 a, b, c を初期化して、while ループで漸化式を計算します。これはアルゴリズムをそのままプログラムするだけです。変数 a の値が 1 に戻ったら循環したので、break で while ループを脱出し、return で cf を返します。

簡単な実行例を示します。

>>> for m in range(1, 17): print(m, cont_frac_root(m))
...
1 [1]
2 [1, 2]
3 [1, 1, 2]
4 [2]
5 [2, 4]
6 [2, 2, 4]
7 [2, 1, 1, 1, 4]
8 [2, 1, 4]
9 [3]
10 [3, 6]
11 [3, 3, 6]
12 [3, 2, 6]
13 [3, 1, 1, 1, 1, 6]
14 [3, 1, 2, 1, 6]
15 [3, 1, 6]
16 [4]

>>> cont_frac_root(97)
[9, 1, 5, 1, 1, 1, 1, 1, 1, 5, 1, 18]
>>> cont_frac_flt(math.sqrt(97), k = 24)
[9, 1, 5, 1, 1, 1, 1, 1, 1, 5, 1, 18, 1, 5, 1, 1, 1, 1, 1, 1, 5, 1, 11, 2]

>>> cont_frac_root(1234)
[35, 7, 1, 3, 1, 4, 4, 2, 9, 1, 1, 2, 3, 1, 1, 34, 1, 1, 3, 2, 1, 1, 9, 2, 4, 4, 1, 3, 1, 7, 70]
>>> cont_frac_flt(math.sqrt(1234), k = 32)
[35, 7, 1, 3, 1, 4, 4, 2, 9, 1, 1, 2, 3, 1, 1, 33, 1, 1, 4, 1, 1, 1, 2, 1, 11, 10, 4, 1, 4, 18, 1, 1]

cont_frac_flt() では循環する前に演算誤差により規則性が失われることがありますが、cont_frac_root() は循環節を正確に求めることができます。

●無理数の有理数近似

無理数を有理数で近似することを考えます。たとえば、\(\sqrt 2 = 1.4142135623730951\) の近似値として、以下のような有理数を考えることができます。

これらの近似値の誤差を Python を使って求めると、次のようになります。

>>> from fractions import Fraction
>>> import math
>>> a = math.sqrt(2)
>>> a
1.4142135623730951
>>> for q in [10, 100, 1000, 10000]:
...     r = Fraction(math.floor(a * q), q)
...     print(abs(a - r), 1 / q, 1 / (q*q))
...
0.014213562373095234 0.1 0.01
0.004213562373095225 0.01 0.0001
0.00021356237309522186 0.001 1e-06
1.3562373095243885e-05 0.0001 1e-08

近似値を \(r = \frac{p}{q}\) とすると、誤差は \(\frac{1}{q}\) の範囲に収まりますが、\(\frac{1}{q^2}\) には収まっていません。実数 \(\alpha\) の近似値 \(\frac{p}{q}\) が次の不等式を満たすとき、有理数 \(\frac{p}{q}\) を「\(\alpha\) の良い近似」といいます。

誤差を \(\frac{1}{q^2}\) の範囲内に収めるのは難しいように思いますが、「\(\alpha\) が無理数なら \(\alpha\) の良い近似が無限個ある」ことが証明されています。これを「ディリクレのディオファントス近似定理」といいます。無理数 \(\alpha\) の正則連分数展開を途中で打ち切ったもの \(\frac{p_n}{q_n} = [b_0, b_1, \ldots, b_n]\) は、無理数 \(\alpha\) の良い近似になることが知られていて、その誤差は以下の式で表されます。

\(\sqrt 5\) で試してみると、次のようになります。

>>> cont_frac_reduce([2,4])
Fraction(9, 4)
>>> abs(math.sqrt(5) - Fraction(9, 4))
0.013932022500210195
>>> 1/16
0.0625

>>> cont_frac_reduce([2,4,4])
Fraction(38, 17)
>>> abs(math.sqrt(5) - Fraction(38, 17))
0.0007738598527309293
>>> 1/(17**2)
0.0034602076124567475

>>> cont_frac_reduce([2,4,4,4])
Fraction(161, 72)
>>> abs(math.sqrt(5) - Fraction(161, 72))
4.31336113213554e-05
>>> 1/(72**2)
0.00019290123456790122

誤差は \(\frac{1}{{q_n}^2}\) よりも小さくなることがわかります。

●ペル方程式

変数の個数が式の本数よりも多い方程式を「不定方程式」といいます。一次不定方程式 \(ax + by = c \ (c \ne 0)\) は、\(a, b, c\) が整数で、\(c\) が \(a\) と \(b\) の最大公約数の倍数であるとき、整数解を持つことが知られています。これを「ベズーの等式」といいます。二次の不定方程式では、以下に示す「ペル方程式 (Pell's Equation)」が有名です。

\(m\) が平方数でない正の整数とき、ペル方程式は自明の解 (\(x = 1, y = 0\)) 以外の整数解を無限個持つことが知られています。その中で \(x + y\sqrt m\) を最小にする解を \(x_1, y_1\) とすると、ペル方程式の解は次式で表すことができます。

ペル方程式の解はこれで全部です。最小解 \(x_1, y_1\) は \(\sqrt m\) の連分数展開で求めることができます。

簡単な例を示しましょう。

循環節の長さが偶数の場合、\(x^2 - m y^2\) は \(+1\) になり、奇数の場合は \(-1\) になります。

これを Python でプログラムすると、次のようになります。

リスト : ペル方程式

def pells_equation(m):
    cf = cont_frac_root(m)[:-1]
    if len(cf) == 0:
        raise ValueError("pells-equation: domain error")
    p = cont_frac_sub(cf, 1, cf[0])
    q = cont_frac_sub(cf, 0, 1)
    if p * p - m * q * q == 1:
        return p, q
    return p * p + m * q * q, 2 * p * q

最初に cont_frac_root() で \(\sqrt m\) の連分数展開を求め、最後の要素をスライスで削除します。もし、cf が空リストであれば、m は平方数なのでエラーを送出します。あとは、cont_frac_sub() で近似値の分子と分母を求め、変数 p, q にセットします。式の値が 1 ならば、return で p と q を返します。そうでなければ、\(p + q\sqrt m\) の二乗を計算して return で返します。

それでは実行してみましょう。

>>> pells_equation(2)
(3, 2)
>>> pells_equation(3)
(2, 1)
>>> pells_equation(7)
(8, 3)
>>> pells_equation(61)
(1766319049, 226153980)
>>> pells_equation(109)
(158070671986249, 15140424455100)

正常に動作しているようです。興味のある方はいろいろ試してみてください、

●循環小数

次は、分数を小数に直すことを考えてみましょう。1/2, 1/5, 1/8, 1/10 などは 0.5, 0.2, 0.125, 0.1 のように途中で割り切れて、有限な桁で表すことができます。これを「有限小数」といいます。ところが、1/3, 1/6, 1/7 は、次のように有限な桁では表すことができません。

1/3 は 3 を無限に繰り返し、1/6 は 0.1 のあと 6 を無限に繰り返し、1/7 は数列 142857 を無限に繰り返します。このような少数を「循環小数 (repeating decimal)」といい、繰り返される数字の列を「循環節」といいます。分数を小数に直すと、有限小数か循環小数のどちらかになります。

循環小数は次のように循環節の始めと終わりを点で示します。

このほかに、循環節に下線を引いたり、カッコで囲む表記法もあります。

●分数を循環小数に変換

分数を循環小数に直すのは簡単です。筆算と同じように計算していくだけです。次の図を見てください。

     0 1 4 2 8 5 7
   ----------------
 7 ) 1
     0
    ----- 
     1 0  <-- 余り 1
       7
    ------- 
       3 0
       2 8
      -------
         2 0
         1 4
        -------
           6 0
           5 6
          -------
             4 0
             3 5
            -------
               5 0
               4 9
              -----
                 1  <-- 余り 1 

7 で割った余り 1 が 2 回現れていますね。これから先は同じ数列を繰り返すことがわかります。7 の剰余は 0 から 6 まであり、剰余が 0 の場合は割り切れるので循環小数にはなりません。現れる余りの数が限られているので、割り切れなければ循環することになるわけです。また、このことから循環節の長さは分母の値よりも小さいことがわかります。

それでは Python3 でプログラムを作ってみましょう。

リスト : 循環小数を求める

# 探索
def position(n, xs):
    for i, x in enumerate(xs):
        if x == n: return i
    return -1

# 分数を循環小数に変換する
def repdec(x):
    if not isinstance(x, Fraction):
        raise TypeError('must be Fraction number')
    m, n = x.numerator, x.denominator
    xs = []
    ys = []
    while True:
        p = m // n
        q = m % n
        if q == 0:
            ys.append(p)
            return ys, [0]
        else:
            x = position(q, xs)
            if x >= 0:
                ys.append(p)
                return ys[:x+1], ys[x+1:]
            xs.append(q)
            ys.append(p)
            m = q * 10

関数 repdec(x) は分数 x を循環小数に変換します。返り値は 2 つのリストを格納したタプルで、先頭のリストが冒頭の小数部分を、次のリストが循環節の部分を表します。途中で割り切れる場合は循環節を [0] とします。これ以降、冒頭の小数部分を有限小数部分と記述することにします。

変数 xs が余りを格納するリスト、変数 ys が商を格納するリストです。最初に変数 m と n に x の分子と分母をセットし、その商と余りを計算して、変数 p と q にセットします。q が 0 ならば割り切れたので有限小数です。p を append() で ys に追加して、それと [0] をタプルに格納して返します。

割り切れない場合、余り q が xs にあるか関数 position() でチェックして、その位置を変数 x にセットします。同じ値を見つけた場合、ys の先頭から x 番目までの要素が有限小数部分で、それ以降が循環部になります。p を append() で ys に追加してから、ys を 2 つに分けて返します。見付からない場合は p, q を append() で xs, ys に追加して、m の値を q * 10 に更新してから割り算を続行します。

それでは実行してみましょう。

>>> for x in range(2, 12): print(repdec(Fraction(1, x)))
...
([0, 5], [0])
([0], [3])
([0, 2, 5], [0])
([0, 2], [0])
([0, 1], [6])
([0], [1, 4, 2, 8, 5, 7])
([0, 1, 2, 5], [0])
([0], [1])
([0, 1], [0])
([0], [0, 9])

>>> repdec(Fraction(355, 113))
([3], [1, 4, 1, 5, 9, 2, 9, 2, 0, 3, 5, 3, 9, 8, 2, 3, 0, 0, 8, 8, 4, 9, 5, 5, 7, 
5, 2, 2, 1, 2, 3, 8, 9, 3, 8, 0, 5, 3, 0, 9, 7, 3, 4, 5, 1, 3, 2, 7, 4, 3, 3, 6, 
2, 8, 3, 1, 8, 5, 8, 4, 0, 7, 0, 7, 9, 6, 4, 6, 0, 1, 7, 6, 9, 9, 1, 1, 5, 0, 4, 
4, 2, 4, 7, 7, 8, 7, 6, 1, 0, 6, 1, 9, 4, 6, 9, 0, 2, 6, 5, 4, 8, 6, 7, 2, 5, 6, 
6, 3, 7, 1, 6, 8])

正常に動作していますね。

●循環小数を分数に変換

循環小数を分数に直すことも簡単にできます。循環小数 - Wikipedia によると、有限小数部分を a、循環節の小数表記を b、節の長さを n とすると、循環小数 x は次の式で表すことができる、とのことです。

ここで、カッコの中は初項 1, 公比 \(\frac{1}{10^n}\) の無限等比級数になるので、その和は次の公式より求めることができます。

簡単な例を示しましょう。

プログラムを作る場合、次のように考えると簡単です。

有限小数部分を表すリストを xs とすると

有限小数部分 = p0 / q0
ただし p0 = xs を整数に変換
      q0 = 10 ^ (len(xs) - 1)

循環節を表すリストを ys とすると

循環節 = p1 / q1
ただし p1 = ys を整数に変換
      q1 = 10 ^ len(ys) - 1

            p0        p1       p0 * q1 + p1
循環小数 = ---- + --------- = --------------
            q0     q0 * q1        q0 * q1

冒頭の有限小数部分を分数に変換するのは簡単ですね。先頭が整数部分になるので、小数部分の桁はリスト xs の長さ - 1 になります。循環節だけを分数に変換する場合、たとえば 1/7 の循環節は [1, 4, 2, 8, 5, 7] になりますが、分子 p' は \(0.142857 \times 10^6 = 142857\) となるので、循環節を表すリストを整数に変換するだけでよいことがわかります。有限小数部分がある場合は、その桁数だけ循環節の部分を右シフトすればよいので、p1/q1 に 1/q0 を掛けるだけです。

これを Python でプログラムすると、次のようになります。

リスト : 循環小数を分数に直す

from functools import reduce

# 循環小数を分数に直す
def from_repdec(xs, ys):
    # 有限小数の部分を分数に直す
    p0 = reduce(lambda a, x: a * 10 + x, xs, 0)
    q0 = 10 ** (len(xs) - 1)
    # 循環節を分数に直す
    p1 = reduce(lambda a, x: a * 10 + x, ys, 0)
    q1 = (10 ** len(ys) - 1)
    # 有限小数 + 循環節
    return Fraction(q1 * p0 + p1, q0 * q1)

アルゴリズムをそのままプログラムしただけなので、とくに難しいところは無いと思います。

それでは実行してみましょう。

>>> for x in range(2, 12): print(from_repdec(*repdec(Fraction(1, x))))
...
1/2
1/3
1/4
1/5
1/6
1/7
1/8
1/9
1/10
1/11
>>> from_repdec(*repdec(Fraction(355, 113)))
Fraction(355, 113)

正常に動作していますね。

●参考文献

1. 芦ヶ原伸之,『超々難問数理パズル 解けるものなら解いてごらん』, 講談社, 2002
2. 奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
3. 連分数展開とその計算方法, (マスオさん)
4. 連分数 - Wikipedia
5. 循環小数 - Wikipedia
6. 平方根の連分数とペル方程式 (PDF), (有澤健治さん)

●プログラムリスト

#
# frac.py : 分数
#
#           Copyright (C) 2024 Makoto Hiroi
#
from fractions import Fraction
from itertools import permutations
from functools import reduce
from math import floor, sqrt

#
# 小町分数
#
def komachi_solver1():
    for (a,b,c,d,e,f,g,h,i) in permutations([1,2,3,4,5,6,7,8,9], 9):
        n = Fraction(a, b * 10 + c) + Fraction(d, e * 10 + f) + Fraction(g, h * 10 + i)
        if n.numerator == 1 and a < d < g:
            print('{}/{}{} + {}/{}{} + {}/{}{} = {}'.format(a,b,c,d,e,f,g,h,i,n))

def komachi_solver2():
    for (a,b,c,d,e,f,g,h,i) in permutations([1,2,3,4,5,6,7,8,9], 9):
        n = Fraction(a, b * c) + Fraction(d, e * f) + Fraction(g, h * i)
        if n.numerator == 1 and a < d < g and b < c and e < f and h < i:
            print('{}/({}*{}) + {}/({}*{}) + {}/({}*{}) = {}'.format(a,b,c,d,e,f,g,h,i,n))

# 単位分数の和
def bunsu(m, n):
    if n % m == 0:
        print('1/{}'.format(n // m))
    else:
        q = n // m + 1
        print('1/{} +'.format(q), end=' ')
        bunsu(m * q - n, n * q)

#
# 正則連分数
#
def cont_frac(n):
    if not isinstance(n, Fraction):
        raise TypeError('must be Fraction number')
    a = n.numerator
    b = n.denominator
    xs = []
    while b > 0:
        xs.append(a // b)
        a, b = b, a % b
    return xs

# 連分数を約分する
def cont_frac_reduce(xs):
    i = len(xs) - 1
    a = xs[i]
    while i > 0:
        a = xs[i-1] + Fraction(1, a)
        i -= 1
    return a

def cont_frac_reduce1(xs):
    return reduce(lambda a, x: x + Fraction(1, a), reversed(xs))

def cont_frac_sub(xs, a0, a1):
    _, a = reduce(lambda a, x: (a[1], x * a[1] + a[0]), xs[1:], (a0, a1))
    return a

def cont_frac_reduce2(xs):
    p = cont_frac_sub(xs, 1, xs[0])
    q = cont_frac_sub(xs, 0, 1)
    return Fraction(p, q)

# 実数を正則連分数に展開する
def cont_frac_flt(x, k = 8, e = 1e-16):
    xs = []
    for _ in range(k):
        xs.append(floor(x))
        y = x - xs[-1]
        if abs(y) < e: break
        x = 1 / y
    return xs

# 平方根の連部数展開
def cont_frac_root(m):
    c0 = floor(sqrt(m))
    cf = [c0]
    if c0 * c0 != m:
        a, b, c = 1, 0, c0
        while True:
            b = c * a - b
            a = (m - b * b) // a
            c = (c0 + b) // a
            cf.append(c)
            if a == 1: break
    return cf

# ペル方程式
def pells_equation(m):
    cf = cont_frac_root(m)[:-1]
    if len(cf) == 0:
        raise ValueError("pells-equation: domain error")
    p = cont_frac_sub(cf, 1, cf[0])
    q = cont_frac_sub(cf, 0, 1)
    if p * p - m * q * q == 1:
        return p, q
    return p * p + m * q * q, 2 * p * q

#
# 循環小数
#

# 探索
def position(n, xs):
    for i, x in enumerate(xs):
        if x == n: return i
    return -1

# 分数を循環小数に変換する
def repdec(x):
    if not isinstance(x, Fraction):
        raise TypeError('must be Fraction number')
    m, n = x.numerator, x.denominator
    xs = []
    ys = []
    while True:
        p = m // n
        q = m % n
        if q == 0:
            ys.append(p)
            return ys, [0]
        else:
            x = position(q, xs)
            if x >= 0:
                ys.append(p)
                return ys[:x+1], ys[x+1:]
            xs.append(q)
            ys.append(p)
            m = q * 10

# 循環小数を分数に直す
def from_repdec(xs, ys):
    # 有限小数の部分を分数に直す
    p0 = reduce(lambda a, x: a * 10 + x, xs, 0)
    q0 = 10 ** (len(xs) - 1)
    # 循環節を分数に直す
    p1 = reduce(lambda a, x: a * 10 + x, ys, 0)
    q1 = (10 ** len(ys) - 1)
    # 有限小数 + 循環節
    return Fraction(q1 * p0 + p1, q0 * q1)