Haskellの畳み込み関数foldrで連分数を計算する

こんにちは｜こんばんは。カエルのアイコンで活動しております @kyamaz です。

はじめに

連分数とは

「連分数」とは、ある数を分数表記したときに、分母の中に更に分数が含まれる次のような形の分数のことを指します。

$$
a_0 + \cfrac{b_1}{a_1 + \cfrac{b_2}{a_2 + \cfrac{b_3}{a_3 + \cfrac{b_4}{\ddots}}}}
$$

特に、分子が全て$1$である場合$(b_1,b_2,b_3,\dots=1)$を「単純連分数」と言います。また単に、単純連分数を「連分数」と呼ぶ場合が多いです。

Haskellではcontinued-fractionsというHackageがあり、連分数を扱うことができますが、本稿では、ワンライナーに近い実装で計算します。参考にした記事は脚注^{1 2}に紹介します。

Haskellの畳み込み関数 foldr

Haskellの畳み込み関数に関しては、以下の記事で詳しく紹介されています。

この記事に解説があるとおり、Haskellではリストは遅延評価されるため、foldrであれば無限リストを扱うことができます。
連分数の定義から、無限リストを扱う方法が自然ですので、連分数を扱う計算では foldr関数を使って計算してみましょう。

連分数の定義を実装する

連分数展開したときのパタメータを、$(a_0,b_1),(a_1,b_2),\dots$というタプルで準備します。$n$回まで繰り返すことを前提にして実装します。

ghci> import Data.Ratio
ghci> continuedFraction n = foldr (\(a, b) c -> fromIntegral a + fromIntegral b / c) 1%1 . take n

この関数の型を見てみましょう。

ghci> :t continuedFraction
continuedFraction
  :: (Integral a1, Integral a2, Integral a3) =>
     Int -> [(a2, a3)] -> Ratio a1

Haskellは型推論がしっかりしており、素直に実装してみると型がこのように決まり、とても便利です。

無理数を連分数で計算する

無理数を計算する

Wikipediaの連分数に取り上げられている既知の無理数の計算を幾つか試してみましょう。まずは、タプルのリストを整えます。

sqrt2 = zip (1 : a) [1,1 ..] where a = repeat 2
sqrt3 = zip (1 : a) [1,1 ..] where a = cycle [1,2]
sqrt5 = zip (2 : a) [1,1 ..] where a = repeat 4
golden = zip (1: a) a where a = repeat 1
napier = zip (2:1: a) [1,1 ..] where a = xs 2 where xs n = [n,1,1] ++ ( xs (n+2))

$\sqrt{2}, \sqrt{3}, \sqrt{5},\varphi (=\frac{\sqrt{5}+1}{2}), e $ について、関数を使って計算した値と連分数で計算した値をそれぞれ出力した例を以下に示します。

ghci> sqrt 2
1.4142135623730951
ghci> fromRational $ continuedFraction 30 sqrt2
1.4142135623730951

ghci> sqrt 3
1.7320508075688772
ghci> fromRational $ continuedFraction 30 sqrt3
1.7320508075688772

ghci> sqrt 5
2.23606797749979
ghci> fromRational $ continuedFraction 15 sqrt5
2.23606797749979

ghci> (sqrt 5 + 1) / 2
1.618033988749895
ghci> fromRational $ continuedFraction 30 golden -- まだ誤差が大きい
1.6180339887496482
ghci> fromRational $ continuedFraction 40 golden
1.618033988749895

ghci> exp 1
2.718281828459045
ghci> fromRational $ continuedFraction 20 napier -- 少し誤差あり
2.7182818284590398
ghci> fromRational $ continuedFraction 21 napier
2.718281828459045

πを計算する

$\pi$を計算する連分数の式として、次のエントリを参考にします。

この２つの連分数を試してみましょう。係数のタプルは次のようになります。

pi'' = zip (3 : [6,6 ..]) ((^ 2) <$> [1,3 ..])
rpi'4 = zip [1,3 ..] ((^ 2) <$> [1 ..])

$$
\pi = 3 + \cfrac{1^2}{6 + \cfrac{3^2}{6 + \cfrac{5^2}{6 + \cfrac{7^2}{6 + \cfrac{9^2}{6 + \cfrac{11^2}{\ddots}}}}}}
$$

ただ、この連分数展開は収束性が悪いです。以下のように$1,000$回目でも小数点以下4桁までしか合わないほどの精度です。

ghci> pi
3.141592653589793
ghci> fromRational $ continuedFraction 40 pi''
3.1415885465407083
ghci> fromRational $ continuedFraction 1000 pi''
3.1415926533392926

そこで逆数で定義される次式は収束が速いです。
$$
\frac{4}{\pi} = 1 + \cfrac{1^2}{3 + \cfrac{2^2}{5 + \cfrac{3^2}{7 + \cfrac{4^2}{9 + \cfrac{5^2}{11 + \cfrac{6^2}{\ddots}}}}}}
$$

ghci> pi
3.141592653589793
ghci> fromRational $ 4 / (continuedFraction 22 rpi'4 )
3.1415926535897936
ghci> fromRational $ 4 / (continuedFraction 23 rpi'4 )
3.141592653589793

試してみると、23回目で良い精度の計算ができております。

小数から連分数を求める

小数から連分数のリストを求める関数toCFracを次のように定義します。

ghci> toCFrac x | x < 0 = [] | otherwise = (i : toCFrac (1/(x - fromIntegral i))) where i = floor x

型推論は以下のようになりました。

ghci> :t toCFrac
toCFrac :: (Integral a, RealFrac t) => t -> [a]

この関数toCFracを使って、前述の値を試してみましょう。

ghci> take 20 $ toCFrac $ sqrt 2
[1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2]
ghci> take 20 $ toCFrac $ sqrt 3
[1,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1]
ghci> take 20 $ toCFrac $ sqrt 5  -- 最初の13個までが正しく計算されています
[2,4,4,4,4,4,4,4,4,4,4,4,4,2,1,10,2,1,7,8]
ghci> take 20 $ toCFrac $ exp 1
[2,1,2,1,1,4,1,1,6,1,1,8,1,1,10,1,1,12,1,1]
ghci> take 20 $ toCFrac $ (sqrt 5 + 1)/2
[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]

つぎに$\pi$も試してみました。その答えは円周率.jp の連分数のページで確認しましたが、最初の13個までが正しく計算されました。

ghci> take 20 $ toCFrac pi  -- 最初の13個までが正しく計算されています
[3,7,15,1,292,1,1,1,2,1,3,1,14,3,3,23,1,1,7,4]

おわりに

最後に、$\pi$を計算できたところまで連分数で表記して、本稿を締めたいと思います。
$$
\pi = 3+\dfrac{1}{7+\dfrac{1}{15+\dfrac{1}{1+\dfrac{1}{292+\dfrac{1}{1+\dfrac{1}{1+\dfrac{1}{1+\dfrac{1}{2+\dfrac{1}{1+\dfrac{1}{3+\dfrac{1}{1+\dfrac{1}{14+\dfrac{1}{\ddots}}}}}}}}}}}}}
$$

なお、お手元の環境で検証する際に、動作が異なるときには参考になるかもしれませんので、毎度ではございますが最後に、本稿を記載するために検証したHaskellの環境を記しておきます。

本稿の環境

本稿のために使用した環境は以下となります。
macOS: Sonoma 14.4 (chip: Apple M1)
GHCup: 0.1.22.0
GHC: 9.6.4

ご一読いただきまして誠に有り難うございます。
（●）（●） Happy Hacking!
/""　＿＿""＼

1. @tsujimotterさんのはてなブログの記事『連分数展開とその計算方法【連分数計算アプリの紹介付き】』 ↩

2. Chris Smith氏のエントリ『Continued Fractions: Haskell, Equational Reasoning, Property Testing, and Rewrite Rules in Action』 ↩