円周率の求め方の公式は多数存在するが、ここでは以下の綺麗な連分数の公式を利用する。
$$
\frac{4}{\pi} = 1 + \dfrac{1^2}
{3 + \dfrac{2^2}
{5 + \dfrac{3^2}
{7 + \dfrac{4^2}
{9 + \dfrac{5^2}
{11 + \dfrac{6^2} {\ddots} }}}}}
$$
この書き方は美しくて良いのだが、どんどん縦にも横にも伸びてきて、書くのが困難である。
それで、以下のように書くことがある。
$$
\frac{4}{\pi} = 1 +
\dfrac{1^2}{3+}\:
\dfrac{2^2}{5+}\:
\dfrac{3^2}{7+}\:
\dfrac{4^2}{9+}\:
\dfrac{5^2}{11+}\:
$$
分子が$n^{2}$になっているが、これが全部１の場合の連分数を、正則連分数という。
ネット上の連分数の解説のほとんどは正則連分数であるが、ここではより一般的な、そして美しい連分数の話をする。
畳み込み関数
連分数は再帰的になっているから、再帰呼び出しすればできるが、それはエレガントでない。
せっかくJuliaを使うのだから、畳み込み関数を使おう。
1, 1/2, 1/3, 1/4, 1/5 の5項の間に – を入れて計算してみる。
julia> 1 - 1/2 - 1/3 - 1/4 - 1/5
-0.2833333333333333これは、次のように、先頭から順番に () で括って計算するのと同じである。
julia> (((1 - 1/2) - 1/3) - 1/4) - 1/5
-0.2833333333333333このように、先頭（左）から順に演算する関数がfoldlである。
julia> foldl(-,map(x->1/x,1:5))
-0.2833333333333333今度は結合の順番を、右側から順に行ってみよう。
julia> 1 - (1/2 - (1/3 - (1/4 - 1/5)))
0.7833333333333333これを、右結合というが、演算子が – という非対称な演算子の場合、右結合するか、左結合するかで結果が異なる。
これを、Juliaの右畳み込み関数 foldr を使うと、次のようになる。
julia> foldr(-,map(x->1/x,1:5))
0.7833333333333333連分数
ということで、最初の円周率の連分数を、右畳み込み関数で、次のようにして円周率を求めることができる。
foldrの第1引数は、2項演算子または引数が２の関数を与える。
すると、次のようになる。
julia> 4/foldr((n, x) -> 2n-1 + n*n/x, 1:20)
3.1415926535897793でも、これで終わってはつまらない。
連分数は分数なんだから、分数のままで求めてみよう。
Juliaは、浮動小数点数だけではなく、有理数（分数）にも対応している。
julia> 1//2 + 1//3
5//6除算演算子(/)の場合には、結果が浮動小数点数になってしまうが、有理数除算演算子(//)の場合には、分数表示のままである。もちろん、約分できる場合には、約分してくれる。
では、/を//に変えて、円周率を求めてみよう。
julia> 4/foldr((n, x) -> 2n-1 + n*n//x, 1:20)
6759794961465344//2151709564809805次に、項数を2から20まで増やしながら、円周率の近似値がどうなるか、見てみよう。
julia> for k in 2:20
my_pi = 4/foldr((n, x) -> 2n-1 + n*n//x, 1:k)
@printf( "%-20.15g\t%s\n", my_pi,my_pi)
end
2.66666666666667 8//3
3.25 13//4
3.12121212121212 103//33
3.14529914529915 368//117
3.14093376764387 14464//4605
3.14170854271357 3126//995
3.14157240608854 112484//35805
3.14159617607035 2297088//731185
3.14159204257068 52481536//16705395
3.1415927593578 345371//109935
3.14159263530947 1836284096//584507385
3.14159265674555 18431540224//5866941465
3.14159265304551 22207688704//7068926865
3.1415926536836 34020159616//10828953135
3.14159265357363 1864534365952//593499721815
3.14159265359258 24727311941632//7870947849765
3.14159265358931 240822550134784//76656198523905
3.14159265358988 2672379400192//850644782715
3.14159265358978 6759794961465344//2151709564809805超高精度演算
もっと精度を上げるために、110項までの計算をしてみよう。
julia> 4/foldr((n, x) -> 2n-1 + n^2//x, 1:110)
ERROR: OverflowError: 10201 * 21317595526195387 overflowed for type Int64
Stacktrace:
・・・以下省略分子や分母が64ビット整数の範囲を超えてしまうとオーバーフローが発生しダメなようだ。
しかし、次のようにbig()を使うと計算できる。
julia> 4/foldr((n, x) -> 2n-1 + n^2//big(x), 1:110)
127429759805932607441365489711200209651957424850172829134567731871622154823776662448269808802257401020416
//40562152340254191526972571181283917950850886703397826665197980233035381516064301728895113563704773029205でも、これでは分かりづらいので浮動小数点表示にしてみた。
julia> my_pi = 4.0/foldr((n, x) -> 2n-1 + n^2//big(x), 1:110)
3.141592653589793238462643383279502884197169399375105820974944592307816406286198先頭の4(整数)を4.0(浮動小数点数)に変えただけである。
円周率は、もともと定数πとして用意されているので、それと比べてみよう。
julia> my_pi - π
0.0完全に一致している、と喜んでよいのだろうか。
これは、現在の BigFloat型の精度の範囲で一致しているというだけである。
ちなみに、現在の精度は次で求まる。
julia> precision(BigFloat)
256256は、256ビット浮動小数点数ということを意味していることに注意しよう。
超高精度演算をしたい場合には、BigFloatで精度指定をすることで、何万桁もの浮動小数点同士の演算もできるようだが、これ以上立ち入らないことにする。
＃人気の記事
＃タグ
