浮動小数点数
コンピュータで確率などを扱うとき、みなさん、floatやdoubleなどの浮動小数点数を使うと思います。
もしかしたら、解析的には簡単な分数で表せるような場合も、浮動小数点数を使えば、一見何がなんだか分からない数に見えます。
また、floatやdoubleなどは、2進数で表された有限桁数の浮動小数点数なので、そもそも分数を浮動小数点数では厳密には表せない場合が数多くあります。
1/3を我々が使い慣れている10進数の小数で表せば、0.333333...となります。桁数は無限にあるため、有限の桁数では、厳密に1/3を書き表すことができません。
同様に2進数の小数でもそういったことが多く起こり、0.1を厳密に表せないことは、浮動小数点数を知る人の中では、非常に有名です。
浮動小数点を分数に直したい
2進数で保存されているので、$\frac{m}{2^{k}}$のような形に直すことはできるかもしれません。
けれど、それはあまり面白くないです。できれば、より自然な形に復元したい。
勿論、情報が切り捨てられているので、完全に復元するのは不可能な場合があります。
けれど、割といい感じに復元できたら面白いと思いませんか?
浮動小数点数を連分数で近似する
連分数とは、次に示す図の右辺のような数式を指し、xを連分数に直すことを連分数展開と呼びます。なお、$a_0, a_1, a2, ...$はすべて整数で、とくに$a_1, a_2, ...$は正の整数です。
(出典: Wikipedia)
ここでは$a_3$で終わっていますが、もっと続く場合もあります。
無理数の場合は無限に続くことが知られています。
連分数展開をするアルゴリズム自体は簡単で、
- xの整数部分を求める。これが$a_i$になる。
- xから、整数部分を引く
- xが0でなければ、引いたものの逆数を求め、新たなxとする。xが0ならば終了する
- 1に戻る。ただし、十分な長さの$a_0, a_1, ...$が得られたら終了する
Pythonで書いていきましょう。
Pythonでは再帰使いたくない派なのでwhileで書いてますが、再帰で書いてもいいです。
def contfrac(x, nmax, vmin=0.000001):
'''連分数を求める。
x: 連分数を求めたい数
nmax: この個数だけaが得られたら終了する
vmin: xの絶対値がこの数を下回ったら、0になったとみなして終了する
'''
ls = []
while len(ls) < nmax:
n = int(x)
x -= n
ls.append(n)
if abs(x) < vmin:
return ls
x = 1 / x
return ls
試しに動かしてみましょう。
contfrac(0.1, 5) # => [0, 10]
contfrac(0.1388888888888889, 5) # => [0, 7, 5]
contfrac(math.sqrt(2), 5) # => [1, 2, 2, 2, 2]
連分数を分数に直す
続いて、連分数を分数にしていきます。
全部整数なので、簡単に直せそうですね。
連分数にして得られたリストを逆順に計算していけば、数値が復元できます。
この際、浮動小数点数は使わずに、分子、分母が整数の状態を保つようにしてください。
def to_frac(x, nmax, vmin=0.000001):
'数値を連分数にしてから分数に戻す。分子, 分母を返す。'
ls = contfrac(x, nmax, vmin)
m = 0
n = 1
for a in reversed(ls):
m, n = n, m
if n:
m += a * n
else:
m, n = a, 1
return m, n
やっていきましょう。
to_frac(0.1, 5) # => (1, 10) つまり0.1 = 1/10ですね
to_frac(0.1388888888888889, 5) # => (5, 36) 実はタイトルの分数は5/36でした
to_frac(math.sqrt(2), 5) # => (41, 29) 41/29 = 1.413793... 誤差はありますが、√2に近いです
to_frac(math.sqrt(2), 10) # => (3363, 2378) 3363/2378 = 1.41421362489... 連分数をより長く取ると精度は上がります
もうちょっと遊んでみましょう。
to_frac(13/17, 5) # => (13, 17) きれいに出ました
to_frac(3/15, 5) # => (1, 5) 3/15は1/5なので、約分されて出てきてしまいました
to_frac(593/193, 5) # => (212, 69) 593と193は素数なので、約分できないはず。精度が落ちて違う数字になってしまっています。連分数を長くしてみます
to_frac(593/193, 10) # => (593, 193) 無事、復元できました
常にうまくいくわけではないですが、浮動小数点数から分数を復元することができました。
面白いですね！