単純連分数

Haskell

http://www.sampou.org/cgi-bin/haskell.cgi?Everyday%3a2004-12-19&l=jpより。

TODO

  • 実数ωとはなんだろう？
追記(2004-12-20)

このエントリを書いた後にhttp://www.sampou.org/cgi-bin/haskell.cgi?Everyday%3a2004-12-19&l=jpが更新された。特別な数ではなくてただ任意の実数に名前をつけただけだった。

  • [1,1..]のときが一番誤差が入りやすいと思って、toRational (SCF [1,1..]) :: Doubleの値が変化しなくなった50桁に制限したのは正しいか？
追記(2004-12-20)

だいたい正しい。IEEEのdoubleは仮数部53bitだったと思うので56桁にしておく。

56桁の理由*1

単純連分数 a = [c_0, c_1, \ldots ]に対し f_n(x) = \frac{1}{c_n + f_{n+1}(x)}, f_{N + 1}(x) = x, f_0(x) = c_0 + f_1(x)により[0,1]区間上の関数 \{ f_j \}_{0 \le j \le N + 1}を定める。aをN桁までで切ったものはf_0(0)である。以後cは1以上の整数をあらわす。
0 \lt \frac{1}{c + x} \lt 1より0 \lt \exists b \lt 1, a = f_0(b)である。
\min{f_1} \le f_1(b) \le \max{f_1}したがって、r := (\max{f_1}/\min{f_1} - 1) \lt 1 / 2^{53}となるようなNをとれば、仮数53bitの浮動小数で表現するには十分である。
\frac{\max{f_n}}{\min{f_n}} - 1 = \frac{ c_n + \max{f_{n + 1}}}{ c_n + \min{f_{n+1}}} - 1 \lt \frac{1}{\frac{c_n}{\min{f_{n+1}}} + 1} ( \frac{\max{f_{n+1}}}{\min{f_{n+1}}} - 1) \lt (1/2)( \frac{\max{f_{n+1}}}{\min{f_{n+1}}} - 1)(c_n = 1でおさえる)。(c + x)/(c + y) \lt 2, \text{if 0 \lt y \le x \lt 1}より r \lt (1/2)^{N - 1} \frac{\max{f_N}}{\min{f_N}} \lt (1/2)^{N - 2}。 N - 2 = 53。

ついでに \sqrt{2}の単純連分数について

[n,n...]と同じ数字だけの単純連分数を考えてみる。
x = n + 1/xと書ける。二次方程式だから解けてx = \frac{n + \sqrt{n^2 + 4}}{2}。
[m,n,n...]という連分数を考える。 m + \frac{2}{n + \sqrt{n^2 + 4}} = (1/2)(\sqrt{n^2 + 4} + (2 m - n))。ここでn = 2 mと置くと、 = \sqrt{ 1 + m^2}。
\sqrt{2} = \sqrt{1 + 1^2}だから、 \sqrt{2}の連分数は[1,2,2,2,..]( m = 1, n = 2 m = 2)とわかる。

Code

module CFrac where

newtype SCFrac = SCF [Int] deriving (Show, Eq)

instance Ord SCFrac where
instance Num SCFrac where
instance Fractional SCFrac where
  fromRational a = SCF cs
    where 
      as = a:map (1 /) (takeWhile (/= 0) 
              $ zipWith (-) as $ map fromIntegral cs)
      cs = map floor as
instance Real SCFrac where
  toRational (SCF xs) 
    = foldr1 (\x y -> x + 1 / y) $ map fromIntegral $ take 56 xs

Example

CFrac> realToFrac (SCF (1:repeat 2))
1.4142135623731
CFrac> realToFrac (sqrt 2) :: SCFrac
SCF [1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,1,1,2,7,1,2,33,2,7,5,2,1,1,16,2]

*1:何回も書き直した。この程度のものでもいきなりタイプというのはよくないのか