piの算出

はじめに

今回もアルゴリズムの体操をしていきたいと思います。
今日のお題は「 」
円周率ってことは知ってますが、算出方法までは気にしたことなかったのでコレを機に。
ここでは何十万桁もの値を求めるためのアルゴリズムではなく、
伝統的なアプローチについてのお話となります。

C言語による最新アルゴリズム事典 (ソフトウェアテクノロジー)

円周率

はい。では円周率を言ってみましょう。皆さん義務教育なので小数点10桁までは暗記していますね。
せーのっ、「さんてんいちよんいちごきゅうにろくごさんご」

それ以降のモノについては 円周率百万桁 こちらでご確認ください♫

$ \pi $ を求めて

この円周率($ \pi $)さんですが、どうやって求めていこう。。。となるわけですが、 $ \pi $を使って計算するものがありましたよね。そうです三角関数です。

そのなかでも$ tan(x) $を用います。(理由は後付で。。。) $ tan(x) $と$ \pi $には$ \eqref{eq:tan-pi} $ , $ \eqref{eq:arc-tan-pi} $のような関係が成り立ちます。

$$ \begin{eqnarray} tan(1) = \frac{\pi}{4} \label{eq:tan-pi} \end{eqnarray} $$

$$ \begin{eqnarray} \pi = 4 \times a rctan(1) \label{eq:arc-tan-pi} \end{eqnarray} $$

この$ a rctan(x) $は、級数展開することができ、$ \eqref{eq:arc-tan-ex} $のように表すことができます。

$$ \begin{eqnarray} a rctan(x) = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + \cdots \label{eq:arc-tan-ex} \end{eqnarray} $$

このように、$ a rctan(x) $を使うことで、$ \pi $を「シンプルな関数」として表現することができるようになるのです。
これが$ tan(x) $を使う理由だと考えています。

これを実装したものが以下になります。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
from matplotlib import pyplot as plt

def arctan(x, n):
    val = float(0)
    for i in range(n):
        r = 2 * i + 1
        val += ((-1) ** i) * (x ** r) / r
    return val

def calcPi(n):
    return 4 * arctan(1, n)

# loop count
n = 1000

x = np.arange(0, n, 1)
y = np.array(list(map(calcPi, x)))

# draw result
plt.plot(x, y)
plt.show()

この結果は次の図のようになります。 少々結果がわかりづらいですが、第1000項まで計算した結果になります。
1000項目の値は、「3.14259365434」となり、小数点2桁までは得られていますね。

もっと欲しがる

「3.14」で納得できませんよ。もっとほしいよ。と人は考える生き物です。 昔の賢い人達もモチロン考えました。

そして出てきたのがJ.Machin(マチン, 1680-1752)の公式$ \eqref{eq:machin} $です。 $$ \begin{eqnarray} \frac{\pi}{4} = 4 \times a rctan\bigl(\frac{1}{5}\bigr) - a rctan\bigl(\frac{1}{239}\bigr) \label{eq:machin} \end{eqnarray} $$

先ほど同様、第1000項までの値の変遷を表示すると、以下のようになります。 なにやら直線のようなものが書かれていますが、拡大してみるとちゃんと値が出ていますね。 ただ、先程のものと比べると圧倒的に早く収束していることが見て取れるかと思います。 その証拠に、第10項の時点ですでに「3.1415926535898362」とかなりの制度を出すことができています。

もっともっと

R.P.Brent, E.Salaminがそれぞれ独立に考えだした、相加相乗平均を用いる方法もあり、コレはSalamin-Brentアルゴリズムと呼ばれています。
このアルゴリズムの特徴として、その収束スピードがあります。
先ほどのMachinの公式では、計算する項が1つ増えるごとに精度が1桁増えるものでした。(少なくとも序盤はそう見える)

しかし、Salamin-Brentアルゴリズムは、各繰り返しによって、精度が3→8→19→41→84桁というスピードで収束していくことが知られています。 アルゴリズムは $$ \begin{align} \begin{split} a_0 &= 1 \\ b_0 &= \frac{1}{\sqrt{2}} \\ a_{n + 1} &= \frac{a_n + b_n}{2} \\ b_{n + 1} &= \sqrt{a_n b_n} \\ \pi_n &= \frac{(2 a_{n + 1})^2}{1 - 4\sum^{n}_{k = 0}{2^{k} ( a^2_k - b^2_k )}} = \frac{(2 a_{n + 1})^2}{1 - 4\sum^{n}_{k = 0}{2^{k} ( a_k - a_{k - 1} )^2}} \end{split} \end{align} $$ となります。 (※C言語による最新アルゴリズム事典では別の記述がありますが、おそらく数式の作成時にミスがあったのだと思います。)

C言語による最新アルゴリズム事典 (ソフトウェアテクノロジー)

def salamin(n):
    a = 1.0
    b = 1.0 / math.sqrt(2.0)
    s = 1.0
    for i in range(n + 1):
        next_a = (a + b) / 2
        b = math.sqrt(a * b)
        s -= 4 * (2 ** i) * ((next_a - a) ** 2)
        a = next_a
    return ((a + b) ** 2) / s

ここまで収束が早くなると、もはや図も必要なく、 結果は、

3.1405792505221686
↓
3.1415926462135428
↓
3.141592653589794

と収束します。 繰り返し3回にして、ここまでの精度を求めることができます。(最後の4は丸め誤差です)

おわりに

$ \pi $をさらにもっと桁の大きいところまで高精度で計算したい場合には、更に効率のよい方法があるようです。
(多倍長演算etc...)
ですがそれはまた別のお話ということで。

三角関数の級数展開からの近似、相加相乗平均を用いた極限値での収束など、
ほんとに頭のいい人はよく思いつきますね(遠い目)