連分数   ― 塩田研一覚書帳 ―


連分数

 整数 a0 と、自然数 a1, a2, ... から次のように表される式を連分数と呼びます:
ただしこれでは幅を取りすぎるので
と書いたり
と書いたりもします。

実数の連分数展開

 実数を連分数の形に書き表すことを連分数展開と言います。
連分数展開
・入力: x
・出力: x の連分数展開
    
a[0] = ( x の切捨て )
n = 0
for ( 適当な回数 ) {
    x = 1 / x
    n += 1
    a[n] = ( x の切捨て )
    x = x - a[n]
}
[a[0], a[1], a[2], ...] を出力
例えば、黄金比は全ての項が 1 の連分数に展開されます:

近似分数

 連分数展開を適当な有限項で打ち切ると近似分数が得られます。 これは「本当は有理数なのだけど分母が分からない」という近似値が与えられたときに、 分母を推定するための強力な武器となります。 例えば
近似値:
x = 0.665513264

連分数展開:
[0, 1, 1, 1, 95, 1, 1, 1, 10297, 1, 1, 2, 1, 21, 3, 1, ...

近似分数:
0 / 1 = 0.000000000000000
1 / 1 = 1.000000000000000
1 / 2 = 0.500000000000000
2 / 3 = 0.666666666666667
191 / 287 = 0.665505226480836
193 / 290 = 0.665517241379310
384 / 577 = 0.665511265164645
577 / 867 = 0.665513264129181
5941753 / 8928076 = 0.665513263999993
5942330 / 8928943 = 0.665513264000005 ...
という具合で、この場合分母は 867 かな、と思う訳です。 Python で書けばこの位のプログラムで求められます:
# ApproximateRatio.py

def ContFrac(x, M = 32):
    a = []
    b = int(x)
    for i in range(M):
        a.append(b)
        x -= b
        if x == 0:
            break
        x = 1 / x
        b = int(x)
    return a
    
def ApproximateRatio(x, M = 32):
    a = ContFrac(x, M)
    p0, q0 = a[0], 1
    p1, q1 = a[1] * p0 + 1, a[1]
    t = [[p0, q0], [p1, q1]]
    del a[0: 2]
    while len(a) > 0:
        p2, q2 = a[0] * p1 + p0, a[0] * q1 + q0
        t.append([p2, q2])
        p0, p1 = p1, p2
        q0, q1 = q1, q2
        del a[0]
    return t

print "Input a deciaml number"
x = input()
a = ContFrac(x)
t = ApproximateRatio(x)
print "x = %15.15lf" % x
print "Continued Fraction Expansion:"
print a
print "Approximate Ratio:"
for p, q in t:
    print "%d / %d = %15.15lf" % (p, q, p * 1.0 / q)
print
整数論も実験科学ですので、近似値を計算して分母を推定する、という作業が必要な研究もあるのです。

連分数展開とユークリッドの互除法

 有理数 a / b の連分数展開は、 入力 a, b に対してユークリッドの互除法を実行したときの各ステップの商を並べたものと等しく、有限項で終了します。 また、近似分数の分母 ・ 分子は剰余列の数字が並びます。

循環連分数

 有限項で終わる連分数を有限連分数、無限に続く連分数を無限連分数と呼びます。 無限連分数のうち、十分番号の大きい項以降で数列が循環するものを循環連分数と呼びます。 連分数の顕著な性質は
  • 有限連分数に展開されること ⇔ 有理数であること
  • 循環連分数に展開されること ⇔ 2次の実数であること
ということです。ここで、2次の実数とは、有理数係数の2次方程式の根になっている実数のことで、 ( a + b √c ) / d の形の数のことです。
 有限小数または循環小数であることと、有理数であることが同値でしたが、 それに対応する定理です。

連分数展開とペル方程式

 自然数 n と整数の未知数 x, y についての
x2 - n y2 = ±1
という方程式をペル方程式と呼びます。 n が平方数でなければペル方程式は必ず解を持ち、しかも √n の連分数展開を用いて最小解を求めることができます。 連分数展開から求めた近似分数 p / q が解を与えるのです。

 解が小さければ、上のプログラムを改良して解けます。
# pellsmall.py

import math

def ContFrac(x, M = 32):
    同上
    
def ApproximateRatio(x, M = 32):
    同上

print "Input a non-square positive integer"
while 1:
    n = input("n = ")
    s = int(math.sqrt(n))
    if s * s != n:
        break
t = ApproximateRatio(math.sqrt(n), 256)
for p, q in t:
    w = p * p - n * q * q  
    if abs(w) == 1:
        print "%d^2 - %d * %d^2 = %d" % (p, n, q, w)
        break
実行例:
170^2 - 19 * 39^2 = 1
1520^2 - 31 * 273^2 = 1
2143295^2 - 94 * 221064^2 = 1
1764132^2 - 193 * 126985^2 = -1
 しかし、ペル方程式の解は思いがけず大きな数になることもあり、 上のプログラムでは精度不足になることも少なくありません。 実数計算を回避して次のように書けば、かなり使い物になります:
# pell.py

# floor of sqrt
def longsqrt(a):
    if a < 0:
        print 'sqrt for negative number'
        exit()
    else:
        if a <= 1:
            return a
        s = a
        while s * s > a:
            s = (s + a / s) / 2
    return s

# bit length
def bitlen(n):
    return len(bin(n)) - 2

# The smallest positive solution for  abs( x^2 - n * y^2 ) = 1'
def pell(n):
    s = longsqrt(n)

    a1, b0, c0, p0, q0 = s, 0, 1, 1, 0
    a2, b1, c1, p1, q1 = (s + b0) / c0, s, n - s * s, s, 1

    while c1 != 1:
        a2 = ( s + b1 ) / c1
        b2 = a2 * c1 - b1
        c2 = c0 + a2 * ( b1 - b2 )
        p2 = p0 + a2 * p1
        q2 = q0 + a2 * q1
        a0, b0, c0, p0, q0 = a1, b1, c1, p1, q1
        a1, b1, c1, p1, q1 = a2, b2, c2, p2, q2
    return [p1, q1]

while 1:
    n = input("n = ")
    s = longsqrt(n)
    if s * s != n:
        break
x, y = pell(n)
print 'n = %d ( %d-bit )' % (n, bitlen(n))
print 'x = %d ( %d-bit )' % (x, bitlen(x))
print 'y = %d ( %d-bit )' % (y, bitlen(y))
print 'x^2 - n * y^2 =', x ** 2 - n * (y ** 2)
実行例:
n = 123 ( 7-bit )
x = 122 ( 7-bit )
y = 11 ( 4-bit )
x^2 - n * y^2 = 1

n = 1234 ( 11-bit )
x = 586327869067265 ( 50-bit )
y = 16691023073856 ( 44-bit )
x^2 - n * y^2 = 1

n = 12345 ( 14-bit )
x = 1196823028442576899590849641 ( 90-bit )
y = 10771703481902106796084652 ( 84-bit )
x^2 - n * y^2 = 1

n = 123456 ( 17-bit )
x = 32153667637494049 ( 55-bit )
y = 91511235212695 ( 47-bit )
x^2 - n * y^2 = 1

n = 1234567 ( 21-bit )
x = 20371567825887579087969922203933358793493846332810110
    6974127231916998110712447355624 ( 277-bit )
y = 18334417735362515888338401277540899079617602694996527
    9326746283914164614149841725 ( 267-bit )
x^2 - n * y^2 = 1

n = 12345678 ( 24-bit )
x = 84798178834254206501069091776063801828858638442491294
    17328362269333115794116905372168214650084077804079568
    14267997463737746006791945771905787959250226099632516
    97702932038103575417349220734335472506700750376403558
    08406717611730435396100932042631926106061387158808852
    82843935289983300226959591157572728353245598926332234
    67646161202348129574120811322088113694783 ( 1193-bit )
y = 24133985779040703487042001168546778233886602226457226
    20173249543953984275822758639086925826348272191944555
    90043334608533759158774990106847117663387640744380349
    49056494509471323897949013056855661657290867499070288
    51068559649123641801036739929261482249920820622857447
    19149238640257078627722120620236367595257737252253551
    47035172813431107542245055578655655064 ( 1181-bit )
x^2 - n * y^2 = 1
n = 123456789, n = 12345678901, n = 111111111 ( 11桁 ) などで実行すると解は驚くほど大きい数字になります。 お試しあれ。

連分数法

 合成数 n を、√n の連分数展開を用いて素因数分解する方法を連分数法と言います。 p / q が √n の近似分数ならば
p2 - n q2 = 比較的小さい数
ですので、この式を mod n すると
p2 ≡ 比較的小さい数 ( mod n )
という合同式が得られます。 右辺に現れる「比較的小さい数」を掛け合わせて平方数が作れれば
u2 ≡ v2 ( mod n )    ..... (*)
という形の式が作れ、gcd(u - v, n) が n の真の約数になる可能性が大きい、というアイデアです。 ただし、√n の連分数展開しか材料が無い、というのは欠点で、 (*) という関係式を導くアイデアは、2次ふるい法や、複数多項式2次ふるい法、数体ふるい法へと発展していきます。

戻る