連分数 ― 塩田研一覚書帳 ―
連分数
整数 $a_0$ と、自然数 $a_1$, $a_2$, ... から次のように表される式を連分数と呼びます: $$a_0 +\frac{1}{a_1 +\displaystyle{\frac{1}{a_2 +\displaystyle{\frac{1}{a_3 +\cdots }}}}} $$ ただしこれでは幅を取りすぎるので $$ a_0 \,+\, \frac{1}{a_1} \,\raise -7pt {+}\, \frac{1}{a_2} \,\raise -7pt {+}\, \frac{1}{a_3} \,\raise -7pt {+}\, \raise -7pt {\cdots} $$ と書いたり $$ [\,a_0,\, a_1,\, a_2,\, a_3,\,\cdots\,] $$ と書いたりもします。実数の連分数展開
実数を連分数の形に書き表すことを連分数展開と言います。
連分数展開
・入力: 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 の連分数に展開されます:
$$
\frac{1+\sqrt{5}}{2}
=
1
+\frac{1}{1
+\displaystyle{\frac{1}{1
+\displaystyle{\frac{1}{1
+\cdots
}}}}}
$$
近似分数
連分数展開を適当な有限項で打ち切ると近似分数が得られます。 これは「本当は有理数なのだけど分母が分からない」という近似値が与えられたときに、 分母を推定するための強力な武器となります。 例えば近似値： 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$ に対してユークリッドの互除法を実行したときの各ステップの商を並べたものと等しく、有限項で終了します。 また、近似分数の分母 ･ 分子は剰余列の数字が並びます。循環連分数
有限項で終わる連分数を有限連分数、無限に続く連分数を無限連分数と呼びます。 無限連分数のうち、十分番号の大きい項以降で数列が循環するものを循環連分数と呼びます。 連分数の顕著な性質は連分数展開とペル方程式
自然数 $n$ と整数の未知数 $x$, $y$ についての $$ x^2 － n y^2 = ±1 $$ という方程式をペル方程式と呼びます。 $n$ が平方数でなければペル方程式は必ず解を持ち、しかも $\sqrt{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$ を、$\sqrt{n}$ の連分数展開を用いて素因数分解する方法を連分数法と言います。 $p / q$ が $\sqrt{n}$ の近似分数ならば $$ p^2 － n q^2 = \mbox{比較的小さい数} $$ ですので、この式を $\bmod{n}$ すると $$ p^2 \equiv \mbox{比較的小さい数} \pmod{n} $$ という合同式が得られます。 右辺に現れる「比較的小さい数」を掛け合わせて平方数が作れれば $$ u^2 \equiv v^2 \pmod{n} \quad \cdots (\ast) $$ という形の式が作れ、$\gcd(u － v, n)$ が $n$ の真の約数になる可能性が大きい、というアイデアです。 ただし、$\sqrt{n}$ の連分数展開しか材料が無い、というのは欠点で、 $(\ast)$ という関係式を導くアイデアは、2次ふるい法や、複数多項式2次ふるい法、数体ふるい法へと発展していきます。