小数を近似する有理数の計算法

質問

0から1までの範囲に入る、小数点以下K桁までの数xが与えられたとき、つまり、（A=10のK乗　と書くことにすると）
x = q/A, qは整数, 0＜q＜A
が与えられたときに、「xが有理数 m/n (0＜m＜n) の値の、小数点以下K桁までの四捨五入による近似値になっている」という条件、つまり
q/A - 1/(2A) ≦ m/n ＜ q/A + 1/(2A)
を満たすような、最小のn を求めるための、速い計算方法が分かりません。

aをbで割った余りをrem(a,b)と書くとき、上記の条件は
2A ＜ rem( (2q-1)n, 2A) + 2n
を満たすことと同値であることまでは考えたんですが、この先工夫しても、結局Aに比例する程度の手間を掛ける計算方法しか思いつかないんです。
もっと速く（ひょっとしてKに比例するぐらいの手間でも）計算できそうな気がするのですが…

なお、この問題の動機は、統計データなどを見ていると母数が少ないはずなのに、頻度がやたら桁の多い小数で表示してあることがあります。そういう場合に、元の有理数（特に分母である母数）が知りたいと考えたことです。（もちろん、「nの倍数」としか分からないですが。）なので、実用上はKは3～8程度です。

回答No.1

UKY · ベストアンサー

例えばこんな風に……

0.444444
= 1 / 2.25
= 1 / (2 + 0.25)
= 1 / (2 + 1 / 4)
= 4 / 9

0.538462
= 1 / 1.85714
= 1 / (1 + 0.85714)
= 1 / (1 + 1 / 1.1667)
= 1 / (1 + 1 / (1 + 0.1667))
= 1 / (1 + 1 / (1 + 1 / 6))
= 7 / 13

誤差を無視した計算なので絶対にうまくいくという保証はありません。
計算の手間は n の桁数に大体比例します。

補足

ありがとうございます。
なっ、なるほど。連分数ってこう使うんですかー

さっそくExcelをつかって、やってみました。連分数で展開して、j回目の展開の端数の部分を四捨五入する、というのをj=1,2,...で試して行って、誤差が問題の条件を満たす範囲内になったら終わるようにしましたら、大変うまく行くようです。（この場合の手間はKの2乗のオーダです。）

実用上はこれでバッチリOKなのですが、でも、本当にこれが最小のnなのかは、難しくて証明できませんでした。うーむむ…

回答No.6

UKY

さっき
> 「分母の桁数が K の半分以下なら必ず分数に直せる」ということが言えそうです。
と書きましたが、桁落ち誤差を考えるとやっぱりこれは必ずしも正しくなさそうです。

やはりろくに考えもせずにものを言うものじゃありませんね。

補足

No.4で教えて戴いた最適近似の意味がやっと飲み込めて、解決したように思います。

a を連分数展開のj回目で近似したものを、n[j]/m[j]とするとき、
∀m（m＜m[j] ならば ∀n (|n/m - a | ＞|n[j]/m[j] - a|))
だと思い込んでいたので分からなくなっていました。
正しくは
∀m（m＜m[j] ならば ∀n (|n/m - a | ＞|n[j-1]/m[j-1] - a|))
ということが、最適という意味ですね。

これなら、1＜q＜A-1のとき、
(2q-1)/(2A) を連分数展開の2j-1回目で近似したものをu[j]/v[j]とすれば、
(i) u[j]/v[j]は既約
(ii) u[j]/v[j]≧(2q-1)/(2A)
(iii) u[j]/v[j]＜u[j-1]/v[j-1]
(iv) ∃jmax (u[jmax]/v[jmax]=(2q-1)/(2A)）
だから、
(2q-1)/(2A)≦u[J]/v[J]＜(2q+1)/A≦u[J-1]/v[J-1]
となるJが存在し、
∀n∀m ((m＜v[J] かつ n/m≧(2q-1)/A) ならば n/m＞u[J-1]/v[J-1])
である。
一方、
(2q+1)/(2A) を連分数展開の2h回目で近似したものをs[h]/t[h]とすれば、同様に
s[H-1]/t[H-1]＜(2q-1)/(2A)≦s[H]/t[H]≦(2q+1)/(2A)
となるHが存在し、
∀n∀m((m＜t[H] かつ n/m＜(2q+1)/(2A)) ならば n/m＜s[H-1]/t[H-1])
である。
だから、v[J]＜t[H] ならば u[J]/v[J]、さもなければ s[H]/t[H]が、範囲[(2q-1)/(2A), (2q+1)/(2A))に入る既約分数の中で、分母が最小である。

さらに、No.5の補足に書いた「別の問題」のほうも、
分母が(2A)以下であって、範囲[(2q-1)/(2A), (2q+1)/(2A))に入る既約分数は、
{u[J]/v[J], ....., u[jmax]/v[jmax]}∪{s[H]/t[H], ...., s[hmax]/t[hmax]}
で尽くされている。
という形で解決できたように思います。

皆様、ご指導ありがとうございました。ご回答下さった全員に届くよう、補足欄を使ってお礼申し上げます。

もうしばらくだけこの質問を開けておきますので、もしわたしが勘違いしているようであれば、ご指摘ください。

回答No.5

UKY

> 切り捨てより四捨五入の方が誤差が小さくなる気がする

あ、やっぱり見落としてた……。
でも、最終項を四捨五入した場合でも特に問題なさそうです。

例えば、
0.615384615
を
1 / (1 + 1 / 1.6)
まで展開したところで四捨五入して、
1 / (1 + 1 / 2)
にしたとします。でもこれは結局
1 / (1 + 1 / (1 + 1 / 1))
に等しいので、更にもう一項多く展開した
1 / (1 + 1 / (1 + 1 / 1.666666))
を四捨五入した
1 / (1 + 1 / (1 + 1 / 2))
よりは誤差の絶対値は大きい(か等しい)はずです。

> 3. この列に現れる整数niを分母とする有理数の誤差の絶対値は、全てのn(n＜ni)について、nを分母とする有理数の誤差の絶対値よりも小さい
という命題は
> 0. j 回展開するより j+1 回展開したほうが、誤差の絶対値が小さくなる
と
> 2. 分母の列 n0, n1, ...,nj, ..が単調増加である
より自明と考えていいと思います。
任意の仮分数は連分数に展開できるので、アルゴリズムが誤差をできるだけ小さくする連分数展開を求め続ける限りは解を「見落とす」ことはないと思います。

> xが実は分母がK桁以上ある既約分数を小数点以下K桁目まで近似した値だった、という場合でもうまく行くんでしょうか。

さすがにそれは無理です。
例えば既約分数 333333 / 1000000 を小数点以下3桁で丸めてしまうと 0.333 になり、どっから見ても 1 / 3 にしか見えません。せいぜい 333 / 1000 が限界です。
私もこのアルゴリズムがうまく働く桁数の限界についてはよく分かっていないのですが、分母が K 桁以内でも近似小数値から分数に直せるとは限らないことは確かです。
例えば
9998 / 9999 = 0.99989998999899989998999899989999……
9997 / 9998 = 0.99989997999599919983996799359872……
の二つを近似小数値から判別するには、少なくとも小数点以下 8 桁必要です。(実際には連分数展開をするときに「桁落ち」という誤差が出るので、もっとたくさんの桁がいると思います)

……と書いているうちに、専門家の方が出てきてくださったようです。
「誤差は1/(m(j)^2)以下」ということなので、とりあえず「分母の桁数が K の半分以下なら必ず分数に直せる」ということが言えそうです。

補足

ありがとうございます。No.3に補足をした後すぐに、「四捨五入をやった場合、j回目の展開時に五入が生じたとすると、切り捨てで展開した場合のj+1回目と同じ結果になるだけ」と気がつきました。ご指摘の通りであり、切り捨てだけ考えれば十分でした。

>> 3. この列に現れる整数niを分母とする有理数の誤差の絶対値は、全てのn(n＜ni)について、nを分母とする有理数の誤差の絶対値よりも小さい
に関しては、
> アルゴリズムが誤差をできるだけ小さくする連分数展開を求め続ける限りは解を「見落とす」ことはないと思います。
と仰る通り、その証明については、No.4で概要を示して頂いたのでなんとかなりそうですが、No.4の補足に書いたように、まだ難点が残っていました。

> 例えば既約分数 333333 / 1000000 を小数点以下3桁で丸めてしまうと 0.333 になり、どっから見ても 1 / 3 にしか見えません。

に関しても、わざわざありがとうございます。ご指摘の通り、最初の質問に「（もちろん、「nの倍数」としか分からないですが。）」と書いたのは間違いでした。ですからこの質問は、元の動機から見るとズレた問題になってしまっています。動機の方を満たすには別の問題：「K桁目以下を四捨五入したときq/Aとなるような、分母が2A(?)未満の全ての既約分数を求む」を考えるべきでした。もちろんここでは、最初の質問についてのみお尋ねしますが、関連が非常に強いのは明らかだと思います。

回答No.4

nakaizu

連分数を使えば目的の近似値を得られるのは既に解答されている通りです。
これが最良の近似であることをこの場で示すのはスペースの関係上困難ですが、次のことが成り立ちますのでそれを利用すると良いでしょう。
a を j回までの連分数で近似したものを n(j)/m(j)と書くことにします。

(1) n(j)/m(j)<a<n(j+1)/m(j+1) またはn(j+1)/m(j+1)<a<n(j)/m(j) が成り立つ(aが有理数のときは等号成立も有り得る)
(2) n(j)/m(j)-n(j+1)/m(j+1)=1/(m(j)m(j+1))または
n(j)/m(j)-n(j+1)/m(j+1)=-/(m(j)m(j+1))
が成り立つ

このことから誤差は1/(m(j)^2)以下であることがわかります。分母がm(j)以下の分数 x/yとn(j)/m(j)との差は最小でも1/(ym(j))ですから誤差の方が小さいことがわかります。
実際には |x/y-n(j)/m(j)|=1/ym(j) となるのはx/y=n(j-1)/m(j-1)の時しかないので 他の分数 x/yとn(j)/m(j)との差は2/(ym(j))以上になります。
よってn(j)/m(j)が最良の近似となります。

(1)と(2)の証明はここに書ききるのはちょっと無理なので、本(たとえば高木貞治の初等整数論講義)で勉強してみてください。

補足

ありがとうございます。なるほど。No.3の補足に書いた命題3. はこれで証明できて、最良の近似になっていることが分かります。そうすると、残った問題は、以下のように言えるでしょうか。

「j回目までの連分数展開で近似したものn(j)/m(j)の誤差が初めて所定の範囲内に入った(a-1/(2A)≦n(j)/m(j)＜a+1/(2A))とします。このとき、分母がm(j)以下の分数 x/yであって、なおかつ誤差が所定の範囲内に入るものはあり得えないのかどうか」

> (1) n(j)/m(j)<a<n(j+1)/m(j+1) またはn(j+1)/m(j+1)<a<n(j)/m(j)

および

> 実際には |x/y-n(j)/m(j)|=1/ym(j) となるのはx/y=n(j-1)/m(j-1)の時しかないので 他の分数 x/yとn(j)/m(j)との差は2/(ym(j))以上になります。

というヒントから、j回目の展開で初めて誤差が所定の範囲内に入ったとするとき、もしy<m(j)であってx/y も同じ範囲内に入るようなyが存在したとすると、
(i) y=m(j-1), |n(j)/m(j)-x/y|= 1/(ym(j)) または
(ii) y≠m(j-1), |n(j)/m(j)-x/y|≧ 2/(ym(j))
である。
(i)の場合は、j-1回目の展開では誤差が所定の範囲内に入らなかったのだからあり得ない。
(ii)の場合は
(iia) n(j-1)/m(j-1) ＞n(j)/m(j) ≧ x/y+2/(ym(j)), x/y＞n(j)/m(j)-1/A または
(iib) x/y-2/(ym(j)) ≧ n(j)/m(j) ＞ n(j-1)/m(j-1), n(j)/m(j)+1/A＞x/y
のどちらかであり、いずれにせよ、
ym(j)＞2A
であるから y＜m(j)より
m(j)＞√(2A)
が(ii)が生じるための必要条件であり、従って、m(j)≦√(2A) が示せれば解決すると思うんですが…どうも難しいです。

回答No.3

UKY

最小の n になることの証明はちょっと手間が掛かりますが可能です。

連分数を仮分数に直したとき、
1. 約分できないこと (つまり、m と n が「互いに素」)
2. j 回展開するより j+1 回展開したほうが n が大きいこと
の二つを示します。

一つ目の条件は、数学的帰納法を使います。
1 / (p + s / t) = t / (pt + s)
において、s と t が「互いに素」であるならば、t と (pt + s) も「互いに素」になりますので、これに帰納法を適用すると連分数展開がどんなに複雑になっても、仮分数が既約であることが示せます。

二つ目のほうはけっこう複雑です。

1 / a

1 / (a + 1 / b)
= b / (ab + 1)

1 / (a + 1 / (b + 1 / c))
= (bc + 1) / (a(bc + 1) + c)

1 / (a + 1 / (b + 1 / (c + 1 / d)))
= (b(cd + 1) + d) / (a(b(cd + 1) + d) + (cd + 1))

とりあえず、3 回目までの展開の分母を次のように置きます。
n0(a) = a
n1(a, b) = ab + 1
n2(a, b, c) = a(bc + 1) + c
n3(a, b, c, d) = a(b(cd + 1) + d) + (cd + 1)
ついでに分子も
m0() = 1
m1(b) = b
m2(b, c) = bc + 1
m3(b, c, d) = b(cd + 1) + d
で、証明すべきことは
n0(a) < n1(a, b) < n2(a, b, c) < n3(a, b, c, d) < ……
です。

まず n0(a) < n1(a, b) を証明してみます。
n0(a) の a を (a + 1 / b) に置き換えると

n0(a + 1 / b)
= a + 1 / b
= (ab + 1) / b
= n1(a, b) / b

よって、n0(a + 1 / b) * b = n1(a, b)
したがって、
n0(a)
< n0(a + 1 / b)
<= n0(a + 1 / b) * b
= n1(a, b)
となり、n0(a) < n1(a, b) が示せました。

これを一般化して、任意の [j] について
n[j](a, b, c, ……, x) < n[j+1](a, b, c, ……, x, y)
が証明できれば OK ということになります。
そのために示すべきことは、
1. n[j](……, x) < n[j](……, x + 1 / y)
2. n[j](……, x + 1 / y) <= n[j](……, x + 1 / y) * y
3. n[j](……, x + 1 / y) * y = n[j+1](……, x, y)
です。
一つ目と二つ目は直感的に明らかだと思います。

三つ目は、n[j](……, x) が x に関する一次式であることがポイントです。
つまり、ax や bcx などの項があっても x^2 や x^3 の項はないということです。
これにより、n[j](……, x + 1 / y) * y は y に関しても一次式になります。
つまり、引数の (1 / y) は y を掛けることで余すところなく整数になり、
(n[j](……, x + 1 / y) * y) 全体も整数になります。

同じく分子についても、(m[j](……, x + 1 / y) * y) が整数になることを示すと、
(n[j](……, x + 1 / y) * y) が (n[j+1](……, x, y)) になることが示せます。

……これで三つの条件が揃ったので、
n0(a) < n1(a, b) < n2(a, b, c) < n3(a, b, c, d) < ……
が証明できました。
……ふぅ。確かにこれは難しいですね。どこか見落としてないかちょっと心配です (^^;

補足

まだ全部理解できたわけではないんですが、凄い証明ですね。ありがとうございます。
証明なさっていないけれど自明と思われるのが
0. j 回展開するより j+1 回展開したほうが、誤差の絶対値が小さくなる
であって、その上に
1. 連分数展開をj回目(j=0,1,....)で打ち切って得られるのは既約分数である
2. 分母の列 n0, n1, ...,nj, ..が単調増加である
従って、j回目の展開で誤差の絶対値が初めて1/(2A)以下になったとき、njは n0, n1, ...,nj, ...の中では、「誤差の絶対値が1/(2A)以下」を満たす最小の整数である。この方針は納得しました。
あとは、
3. この列に現れる整数niを分母とする有理数の誤差の絶対値は、全てのn(n＜ni)について、nを分母とする有理数の誤差の絶対値よりも小さい
が示せたら解決するのかな、と思いますが、これは自明なんでしょうか。
証明の中ではj回目の展開において端数を(No.1の補足で書いたように四捨五入するんじゃなくて、）切り捨てていらっしゃるようです。 でも直感的に、切り捨てより四捨五入の方が誤差が小さくなる気がするんで、ちょっとパラドックスの気分です。

回答No.2

venus_j

ちなみに先ほどの小数0.538462を正確に連分数展開すると
1/(1+1/(5+1/(1+1/(12819+1/(1+1/(1+1/2)))))
となります。
ここで12819のところを無限大で近似してみてはどうでしょうか。計算には自信がありませんが。

補足

ありがとうございます。
確かにこの例ではここで切るっきゃない感じです。ご指摘は、No.1の補足で「端数」と言った部分が0に近いか1に近いときを目安にして、そこで切ると判断すれば、Kに比例する手間になる、という意味であるように思います。「近い」と判断する簡単な基準が何か具体的に出せればアルゴリズムとして完成するんでしょうね。xが実は分母がK桁以上ある既約分数を小数点以下K桁目まで近似した値だった、という場合でもうまく行くんでしょうか。