円周率を求める式

円周率の公式

プログラム中で円周率が必要となった時に使われる式

プログラム中に円周率が欲しくなった場合、以下の式を使うことが多い。これらの式は「凡そ 3」のような近似ではなく、関数の返す精度で正しい値を得る。

\begin{align}
\tan45^∘ &= \tan\frac\pi4 = 1 より\\
\pi &= 4\arctan(1)\\_.\\
\sin90^∘ &= \sin\frac\pi2 = 1 より\\
\pi &= 2\arcsin(1)\\_.\\
\cos180^∘ &= \cos\pi = -1 より\\
\pi &= \arccos(-1)
\end{align}

マクローリン展開で求めているであろうことを考えると引数の値が一番引数が小さい $\arctan$ が有利ですね。ただ、$45^∘$ で折り返してる可能性を考えると $\arcsin$ や $\arccos$ でも悪くなさそうですが。

四則演算を用いて計算する方法

上記のプログラム中で円周率を求めるのも内部的では四則演算で計算を行い値を求めている。以下では四則演算で求めることのできる方法を挙げる。

以下、$\arctan$は主値 $\left( -\frac\pi2 \le \arctan x \le \frac\pi2 \right)$ とする。

ライプニッツ (Leibniz) の式

遅いけども超有名な式。

\frac \pi4 = \arctan1 = \sum_{n=0}^\infty \frac {\left(-1\right)^n}{2n+1}

$\tan\dfrac\pi4 = 1$ なのでそれの逆 $\arctan 1 = \dfrac\pi4$ から求められる。
右辺の総和 ($\sum$) の式は$\arctan1$ をマクローリン展開したもの。
ちなみに $\arctan x$ のマクローリン展開は

\begin{align}
\arctan x &= \sum_{n=1}^\infty \frac{\left(-1\right)^{n-1}}{2n-1} x^{2n-1}\\
&= \frac {x^1}1 - \frac {x^3}3 + \frac{x^5}5 - \frac {x^7}7 + \frac {x^9}9 \dotsb\\
\end{align}

これは 2項を 1つに纏めて以下のようにも書ける。

\arctan x = \sum_{n=1}^\infty \left( \frac1{4n-3} - \frac1{4n-1}x^2\right)x^{4n-3} 

これは符号の計算をしなくて良くなるなど少しの手間で計算量を減らすことができる。
これでライプニッツの式を書くと$x=1$なので、

\begin{align}
\arctan 1 = \frac \pi4
&= \sum_{n=1}^\infty \frac{\left(-1\right)^{n-1}}{2n-1} &\rightarrow ① \\
&= \sum_{n=1}^\infty \left( \frac1{4n-3} - \frac1{4n-1}\right) &\rightarrow ②\\
&= \sum_{n=1}^\infty \frac 2{16n\left(n-1\right)+3} &\rightarrow ③\\
\end{align}

割り算はコストが高いので、プログラミングするなら後の式のほうが少し早いかな？

GNU Octave で ①～③ の計算時間を比較してみた。

多少早くなるけど、冪乗を消した方が遥かに効いてますね。

arctan系 多項式 (マチンの式など)

マクローリン展開は 0 付近の数値の収束は早いが 0 から離れるにしたがって収束が遅くなる。そのためライプニッツの式の $\arctan$ の引数を小さくする式の変形を繰り返していったもの。

\begin{array}{rll}
\frac \pi4 &= 4 \arctan \frac 15 - \arctan \frac 1{239} &\mbox{(by John Machin:1706)}\\
&= \arctan \frac 12 - \arctan \frac 13 &\mbox{(by John Machin)}\\
&= 2 \arctan \frac 12 - \arctan \frac 17 &\mbox{(by John Machin)}\\
&= 2 \arctan \frac 13 + \arctan \frac 17 &\mbox{(by John Machin)}\\
&= 3 \arctan \frac 14 + \arctan \frac 5{99} &\mbox{(by John Machin)}\\
&= 3 \arctan \frac 14 + \arctan \frac 1{20} + \arctan \frac 1{11985} &\mbox{(by John Machin)}\\
&= 8 \arctan \frac 1{10} - \arctan \frac 1{100} - \arctan \frac {11}{5637} + \arctan \frac {10893}{41480222056636} &\mbox{(by John Machin)}\\

&= 5 \arctan \frac17 + 2 \arctan \frac 3{79} &\mbox{(by Leonhard Euler:17xx)}\\

&= 12 \arctan \frac 1{18} + 8 \arctan \frac 1{57} - 5 \arctan \frac 1{239} &\mbox{(by Johann Carl Friedrich Gauß:1863)}\\

&= 8 \arctan \frac 1{10} - 4 \arctan \frac 1{515} - \arctan \frac 1{239} &\mbox{(by Robert Simson)}\\

&= 6 \arctan \frac 1{8} + 2 \arctan \frac 1{57} + \arctan \frac 1{239} &\mbox{(by Fredrik Carl Mülertz Størmer:1896)}\\

&= 44 \arctan \frac 1{57} - 5 \arctan \frac 1{239} + 12 \arctan \frac {369}{128467} &\mbox{(—)}\\

&= 17 \arctan \frac 1{23} + 8 \arctan \frac 1{182} + 5 \arctan \frac {127}{228636} &\mbox{(—)}\\

&= 44 \arctan \frac 1{57} + 7 \arctan \frac 1{239} - 12 \arctan \frac 1{682} + 24 \arctan \frac 1{12943} &\mbox{(by Fredrik Carl Mülertz Størmer:1896)}\\

&= 22 \arctan \frac 1{28} + 2 \arctan \frac 1{443} - 5 \arctan \frac 1{1393} - 10 \arctan \frac 1{11018}  &\mbox{(by E.B.Escott:1896)}\\

&= 17 \arctan \frac 1{23} + 8 \arctan \frac 1{182} + 10 \arctan \frac 1{5118} + 5 \arctan \frac 1{6072}   &\mbox{(by Jörg Arndt:1993)}\\

&= 12 \arctan \frac 1{49} + 32 \arctan \frac 1{57} - 5 \arctan \frac 1{239} + 12 \arctan \frac 1{110443} &\mbox{(by 高野喜久雄:1982)}\\
\end{array}

この高野喜久雄の式は 2002年に 1兆2441億桁まで計算したときに使われ、検算は $24 \arctan \frac 1{12943}$ で終わる Størmer の式で行われた。

\begin{array}{rll}
\frac \pi4
&= 83 \arctan \frac 1{107} + 17 \arctan \frac 1{1710} - 44 \arctan \frac 1{225443} - 68 \arctan \frac 1{2513489} + 22 \arctan \frac 1{42483057} + 34 \arctan \frac 1{7939642926390344818}  &\mbox{(by 清宮 宏:？)}\\

&= 83\arctan\frac{1}{107} + 17\arctan\frac{1}{1710} - 22\arctan\frac{1}{103697}
- 24\arctan\frac{1}{2513489} - 44\arctan\frac{1}{18280007883}
+ 12\arctan\frac{1}{7939642926390344818}
+ 22\arctan\frac{1}{3054211727257704725384731479018} &\mbox{(by M. Wetherfield:2004)}\\
\end{array}

これらは加法定理を変形した関係式などを使い導出された。(他の関係を使うこともある)

ライプニッツの式からマチンの式の導出

\tan\alpha = \frac15 ― ① の両辺を \arctan をとって\\
\alpha = \arctan\frac15 ― ①'\\
\begin{array}{rl}
ここで \tan2\alpha を求める。加法定理より\\
\tan\left(\alpha + \alpha\right) &= \frac{\tan\alpha + \tan\alpha}{1-\tan\alpha \tan\alpha}\\
&= \frac{2\tan\alpha}{1-\tan^2\alpha}\\
①を代入して　　　　　　　\\
\tan2\alpha &=  \frac{2\times\frac15}{1-\left(\frac15\right)^2}
= \frac{\frac25}{1-\frac1{25}}
= \frac{\frac25}{\frac{24}{25}}
= \frac{10}{24}\\
\therefore \tan2\alpha &= \frac5{12} ― ②\\
次に \tan4\alpha を求める。同じく加法定理より&\\
\tan\left(2\alpha + 2\alpha\right) &= \frac{\tan2\alpha + \tan2\alpha}{1-\tan2\alpha \tan2\alpha}\\
&= \frac{2\tan2\alpha}{1-\tan^22\alpha}\\
②を代入して　　　　　　　\\
\tan4\alpha &= \frac{2\times\frac5{12}}{1-\left(\frac5{12}\right)^2}\
= \frac{\frac56}{1-\frac{25}{12^2}}
= \frac{\frac56}{\frac{119}{12^2}}
= \frac{\frac56\times12^2}{119}
= \frac{120}{119}\\
\therefore \tan4\alpha &= 1+\frac1{119} ― ③\\
更に \tan\left(4\alpha-\frac\pi4\right) を加法定理により\\
\tan\left(4\alpha-\frac\pi4\right)
&= \frac{\tan4\alpha + \tan\left(-\frac\pi4\right)}{1 - \tan4\alpha \times \tan\left(-\frac\pi4\right)}
= \frac{\tan4\alpha + \left(-1\right)}{1-\tan4\alpha\times\left(-1\right)}
= \frac{\tan4\alpha -1}{1+\tan4\alpha}\\
③を代入　　　　　　　　\\
&= \frac{1 + \frac1{119} -1}{1 + 1 + \frac1{119}}
= \frac{119+1-119}{119+119+1}
= \frac{1}{239}\\
\tan\left(4\alpha-\frac\pi4\right) &= \frac1{239}\\
両辺の \arctan をとって　　\\
4\alpha-\frac\pi4 &= \arctan\frac1{239}\\
\frac\pi4 &= 4\alpha - \arctan\frac1{239}\\
①'を代入　　　　　　　　\\
&= 4\arctan\left(\frac15\right) - \arctan\left(\frac1{239}\right)\\
\end{array}

図で書くと結構分かり易いんですけど、面倒ですねぇ。

こんな関係式から導出したりもするようです。

\tan(\alpha + \beta) = \frac {\tan\alpha+\tan\beta}{1-\tan\alpha\tan\beta}\quad (加法定理)\\

\left\{
\begin{array}{rr}
p = \tan \alpha &\left( \leftrightarrow \arctan p = \alpha \right)\\
q = \tan \beta  &\left( \leftrightarrow \arctan q = \beta \right)
\end{array}
\right.\\

\begin{align}
とおき、これを加法定理に代入する\\
まず右辺の \tan を書き換える\\
\tan(\alpha + \beta) &= \frac {p+q}{1-pq}\\
両辺の \arctan をとる\\
\alpha + \beta &= \arctan \frac {p+q}{1-pq}\\
左辺の \alpha と \beta を書き換える\\
\arctan p + \arctan q &= \arctan \frac {p+q}{1-pq}\\
\arctan p &= \arctan \frac {p+q}{1-pq} - \arctan q
\end{align}

これを使って、マチンの公式を変形してみると

\begin{align}
\frac \pi4 &= 4 \arctan \frac 15 - \arctan \frac 1{239}\\
適当に \\
p = \frac 15, & q = \frac {1}{25} とおくと \frac {p+q}{1-pq} = \frac {15}{62} になるので \\
\frac \pi4 &= 4 \left( \arctan \frac {15}{62} -  \arctan \frac {1}{25}\right) - \arctan \frac 1{239}\\
&= 4 \arctan \frac {15}{62} -  4\arctan \frac {1}{25} - \arctan \frac 1{239}
\end{align}

という関係式を作ることができる。ただし、分子を １にした方が計算量が減るのであまり良い式ではない。それを考慮すると、

\left\{
\begin{align}
p &> \frac {p+q}{1-pq} \\
p &> q > 0 \\
\left( 1-pq \right) &= n\left( p+q \right) : nは正の整数\\
\end{align}
\right. 

の条件を満たす必要がある。もう少し変形して $\arctan$ の式

\begin{align}
\arctan p &= \arctan \frac {p+q}{1-pq} - \arctan q \\
& を少し書き換えて \\
\arctan p &= \arctan \frac 1s - \arctan \frac 1r \\
& とすると \\
r &= \frac 1q : r は正の整数 \\
s &=  \frac {1-pq}{p+q}\\
 &=  \frac {1-\frac pr}{p+\frac 1r}\\
 &=  \frac {r-p}{rp+1} : s は正の整数 \\
\end{align}

よって条件は

p を与えたとき \\
\left\{
\begin{align}
s &= \frac {r-p}{rp+1} \\
p &< s \\
p &> \frac 1r > 0 \\
r, &sは正の整数\\
\end{align}
\right. 

となるので、wolframalpha に突っ込んでみる。

これの後ろに p = 1/10

他にも以下のような関係式がある

\begin{align}
\arctan \frac1n &= 2 \arctan \frac 1{2n} - \arctan \frac1{4n^3+3n}\\
\arctan \frac2n &= 2 \arctan \frac 1n - \arctan \frac 2{n^3+3n}
\end{align}

等々
他の $\arctan$ 系式としては $\tan \frac\pi6 = \frac 1{\sqrt3} \Rightarrow \arctan \frac 1{\sqrt3} = \frac\pi6 より$

\frac \pi6 = \frac 1{\sqrt3} \left(1 - \frac 1{3\times3} + \frac 1{5 \times 3^2} - \frac 1 {7\times 3^3} + \cdots \right)\qquad\mbox{(by Sharp:1699)}

などがある。

参考 計算効率の良いarctan関係式の探索の試み.pdf,
arctan 関係式一覧 松元隆二
逆正接関数による円周率の公式 清宮 宏
清宮さんの評価によると、上記 $\arctan$ 系の式で最良であると考えられるのは 6項の清宮さん発見の式らしい。

無限乗積

\left\{
\begin{align}
x_1 &= \sqrt{\frac12}\\
x_{n+1} &= \sqrt{\frac{1+x_n}2}
\end{align}\right. \quad として\\
\frac{\pi}{2} = \prod_{n=1}^\infty x_n \qquad\mbox{(by François Viète:1579)}\\
\begin{array}{l}
これを総積の記号を使わずに書くと\\
\frac 2\pi =
 \sqrt{\frac12} \cdot
 \sqrt{\frac12 + \frac12\sqrt{\frac12}} \cdot
 \sqrt{\frac12 + \frac12\sqrt{\frac12 + \frac12\sqrt{\frac12}}} \cdot
 \sqrt{\frac12 + \frac12\sqrt{\frac12 + \frac12\sqrt{\frac12 + \frac12\sqrt{\frac12}}}} \cdots\\
\end{array}

項が増える度に $+ \frac12\sqrt{\frac12}$ が中に入っていく。→円周率.jp 多角形の利用
ウォリス積(by John Wallis:1655) 収束は遅いらしい。

\begin{align}
\frac{\pi z}{\sin \pi z} & = \prod_{n=1}^\infty \frac{n^2}{n^2-z^2}\\
この式に& z=\frac12 を代入して \\
\frac{\pi}{2} & = \prod_{n=1}^\infty \frac{4n^2}{4n^2-1}\\
& = \prod_{n=1}^\infty \frac{2n \cdot 2n}{(2n-1)(2n+1)}\\
& = \frac{2 \cdot 2}{1 \cdot 3} \times \frac{4 \cdot 4}{3 \cdot 5} \times \frac{6 \cdot 6}{5 \cdot 7} \times \frac{8 \cdot 8}{7 \cdot 9} \cdots
\end{align}

数字の並びは凄く分かりやすいけど、こんなので円周率が求まるんですねぇ。(考えたウォリスさんの証明は甘く、その後に三角関数の積分を使って証明された。)

別の形式で書くと

\begin{align}
正弦関数に関する&乗積公式 \\
\frac{\sin(\pi x)}{\pi x} &= \prod_{n=1}^\infty (1-\frac{x^2}{n^2}) \\
から、ガンマ関数&の反転 (相反) 公式 \\
\Gamma(x)\Gamma(1-x) &= \frac{\pi}{\sin(\pi x)} \\
を経由して　　　 \\
\sqrt{\pi} &= \Gamma(1/2) \\
を得る。これを展開&して \\
&= \lim_{n \to \infty} \frac{2^{2n}(n!)^2}{(2n)!\sqrt{n}}
\end{align}

となるそうな。
複素関数論における無限積の公式 青山学院大学 江添 優里

Newton / Euler の式

\frac\pi2 =
 \sum^{\infty}_{k=0} \frac{k!}{(2k+1)!!}  =
 \sum^{\infty}_{k=0} \frac{2^kk!^2}{(2k+1)!} =
 1 + \frac13 \left( 1+ \frac25 \left( 1 + \frac37 \left( 1 + \cdots\right) \right) \right)

ブラウンカーの公式 (William Brouncker, 2nd (1620-1684))

計算内容はライプニッツの式とほぼ同じで、収束がとても遅く実用的ではありませんが、無理数を規則性のある連分数で書き表した式。

\begin{align}
\frac4\pi &= 1 + \frac{1^2}{
  \displaystyle 2 + \frac{3^2}{
  \displaystyle 2 + \frac{5^2}{
  \displaystyle  2 + \frac{7^2}{
  \displaystyle   2 + \frac{9^2}{
  \displaystyle    2 + \frac{11^2}{
       2 + \cdots
     }
    }
   }
  }
 }
} \\
この式は&場所を取るため以下のように書かれることもある\\
 &= 1 + \frac{1^2}{
 2 +} \frac{3^2}{
 2 +} \frac{5^2}{
 2 +} \frac{7^2}{
 2 +} \frac{9^2}{
 2 +} \frac{11^2}{
 2 +} \cdots \\

\pi &= \frac4{
 \displaystyle 1 + \frac{1^2}{
 \displaystyle  3 + \frac{2^2}{
 \displaystyle   5 + \frac{3^2}{
 \displaystyle    7 + \frac{4^2}{
 \displaystyle     9 + \frac{5^2}{
 \displaystyle     11 + \frac{6^2}{
 \displaystyle      13 + \frac{7^2}{
        15 + \cdots
       }
      }
     }
    }
   }
  }
 }
} \\
同様に\\
 &= \frac4{
 1 +} \frac{1^2}{
 3 +} \frac{2^2}{
 5 +} \frac{3^2}{
 7 +} \frac{4^2}{
 9 +} \frac{5^2}{
 11+} \frac{6^2}{
 13+} \frac{7^2}{
 15+} \cdots \\
\end{align}

ラマヌジャン系 (モジュラー関数)

\begin{align}
\frac{1}{\pi} &= \frac{2\sqrt{2}}{99^2} \sum_{k=0}^\infty \frac{(4k)!(1103+26390k)}{(4^k99^kk!)^4}\\
&=\frac{2\sqrt{2}}{9801} \sum^\infty_{k=0} \frac{(4k)!(1103+26390k)}{396^{4k}k!^4}
&\mbox{(by Srinivasa Aiyangar Ramanujan:1914)}\\
\frac{4}{\pi} &= \sum_{n=0}^\infty \frac{(-1)^n(4n)!(1123+21460n)}{882^{2n+1}(4^nn!)^4}
&\mbox{(by Srinivasa Aiyangar Ramanujan)}\\
\\
\frac{1}{\pi} &= 12\sum^\infty_{k=0} \frac{(-1)^k (6k)!(13591409+ 545140134k)}{(3k)!k!^3 640320^{3k+\tfrac{3}{2}}}
&\mbox{(by Chudnovsky:1994)}\\

\frac1\pi &= \frac{12}{\sqrt C} \sum_{n=0}^\infty \frac{(-1)^n(6n)!(A+Bn)}{(3n)!(n!)^3C^n}\\

&\hspace{3em} \left\{
\begin{array}{l}
A=13591409\\
B=545140134\\
C=640320^3
\end{array}
\right.
&\mbox{(by Chudnovsky)}\\

&\mbox{Chudnovsky と同じ関係式を用いて}\\
&\hspace{3em} \left\{
\begin{array}{l}
A=1657145277365 + 212175710912\sqrt{61}\\
B=107578229802750 + 13773980892672\sqrt{61}\\
C=(5280(236674+30303\sqrt{61}))^3
\end{array}
\right.
&\mbox{(by Borwein)}\\
\end{align}

ラマヌジャンの式を使い William Gosper が 1985年に 1752万6200桁の計算をしたが、当然ラマヌジャンの式が証明されている訳はなく、当時は本当に円周率の値なのかどうか分からなかった。(その後証明された)

Chudnovsky の式の変形

\begin{align}
\frac1\pi &= \frac{12}{\sqrt C} \sum_{n=0}^\infty \frac{(-1)^n(6n)!(A+Bn)}{(3n)!(n!)^3C^n}\\
K_n の項 &= \frac{(-1)^n(6n)!}{(3n)!(n!)^3C^n} と置くと\\

K_n の次の項は &= \frac{(-1)^n(-1)(6n)!(6n+1)(6n+2)(6n+3)(6n+4)(6n+5)(6n+6)}{(3n)!(3n+1)(3n+2)(3n+3)((n!)(n+1))^3C^{n+1}}\\
\frac{K_{n+1}}{K_n} &= \frac{(-1)^n(-1)(6n)!(6n+1)(6n+2)(6n+3)(6n+4)(6n+5)(6n+6)}{(3n)!(3n+1)(3n+2)(3n+3)(n!)^3(n+1)^3C^{n+1}} \times \frac{(3n)!(n!)^3C^n}{(-1)^n(6n)!}\\
&= \frac{-(6n+1)(6n+2)(6n+3)(6n+4)(6n+5)(6n+6)}{(3n+1)(3n+2)(3n+3)(n+1)^3C}\\
&= \frac{-72(6n+1)(3n+1)(2n+1)(3n+2)(6n+5)(n+1)}{3(3n+1)(3n+2)(n+1)(n+1)^3C}\\
&= \frac{-24(6n+1)(2n+1)(6n+5)}{(n+1)^3C}\\
& 元の式に戻すと\\
\dfrac1\pi &=\frac{12}{\sqrt C} \sum_{n=0}^\infty Kn(A+Bn)\\
&= \frac{12}{\sqrt C} \sum_{n=0}^\infty (A+Bn) \prod_{k=0}^{n-1}\frac{K_{n+1}}{K_n}\\
&= \frac{12}{\sqrt C} \sum_{n=0}^\infty (A+Bn) \prod_{k=0}^{n-1}\frac{-24(6n+1)(2n+1)(6n+5)}{(n+1)^3C}\\
&= \frac{12}{\sqrt C}
\left(
  \underset{n=0 の項}
  {A + \frac{K_1}{K_0}}
\left(
 \underset{n=1 の項}
 {A + B + \frac{K_2}{K_1}} \cdots
\left(
 \underset{n の一般項}
 {A + Bn + \frac{K_{n+1}}{K_n}}
\left(\cdots\right) \right) \cdots \right.\right)\\
\end{align}

この変形を使うことで計算が楽になります。
Chudnovsky の式を使った計算方法については Chudnovsky の公式を用いた円周率の計算用メモをどうぞ。

C での実装例 Computing billions of π digits using GMP から gmp-chudnovsky.c
C++ での実装例 C++ - 円周率計算(Chudnovsky の公式使用)!
※コンパイルは apt install libgmp-dev してから

四次収束 (Eugene Salamin & Richard Brent)

算術幾何平均 (AGM:Arithmetic Geometric Mean) を使った式。いや、良く分からないんですが。

n=0\,のとき\left\{
\begin{array}{l}
\alpha_0 = \dfrac12 \left( 1 + \dfrac1{\sqrt[4]2} \right)\\
\beta_0 = \dfrac12 \left( 1-\dfrac1{\sqrt[4]2} \right)\\
\gamma_0 = \dfrac12
\end{array}\right.\quad
それ以降
\left\{
\begin{align}
\alpha_{n+1} &= \dfrac12 \left( \alpha_n + \sqrt[4]{\left(\alpha_n^2 + \beta_n^2\right) \left(\alpha_n^2 - \beta_n^2\right)} \right)\\

\beta_{n+1} &= \dfrac12 \left( \alpha_n - \sqrt[4]{\left(\alpha_n^2 + \beta_n^2\right) \left(\alpha_n^2 - \beta_n^2\right)} \right)\\

\gamma_{n+1} &= \left(2\alpha_n^2 + \beta_n^2\right)\beta_n^2
\end{align}
\\
漸化式\quad
\pi = \lim_{n \to \infty} \dfrac{2\alpha_n^４ - \beta_n^4}{\displaystyle 1 - \sum_{j=0}^{n+1}4^j\gamma_j}\\
\right.\\

らしい。よくわからない。
このへん？この辺も理論としては繋がってそう。

Gauß-Legendre

これも算術幾何平均 (AGM) 系 の一つ。
Super PI 採用のアルゴリズムで、これも先ほどのと同じような感じで以下のように表される。

n=0\,のとき\left\{
\begin{align}
a_0 &= 1\\
b_0 &= \dfrac1{\sqrt2}\\
t_0 &= \dfrac14\\
p_0 &= 1
\end{align}\right.\qquad それ以降
\left\{
\begin{array}{rcl}
a_{n+1} &=& \dfrac{a_n + b_n}{2}\\

b_{n+1} &=& \sqrt{a_nb_n}\\

t_{n+1} &=& t_n - p_n\left(a_n - a_{n+1}\right)^2\\

p_{n+1} &=& 2p_n
\end{array}\right.\\
\\
漸化式\quad
\pi = \lim_{n \to \infty} \dfrac{\left(a+b\right)^2}{4t}\\

C 言語の実装 super pi を作ろう

ここらへんから円周率を計算するプログラムを作ってみようかな。

Bailey–Borwein–Plouffe (BBP) の式

16進数や 2進数で任意桁数の値を求める方法。Computation of the n'th digit of pi in any base in O(n^2)
普通は小数点以下を上の桁から順番に求めていくが、この式は途中の位をすっ飛ばして欲しい桁の値を直接求めることができる。円周率.jp - BBP系公式, Wikipedia

\pi = \sum_{k = 0}^\infty \left[ \frac{1}{16^k} \left( \frac{4}{8k + 1} - \frac{2}{8k + 4} - \frac{1}{8k + 5} - \frac{1}{8k + 6} \right) \right]\quad \mbox{(David H. Bailey, Peter Borwein, Simon Plouffe:1995)}

Bellard の式

Bailey-Borwein-Plouffe の式を約 43% 高速化したもの。Wikipedia

\begin{align}
\pi = \frac1{2^6} \sum_{n=0}^\infty \frac{(-1)^n}{2^{10n}} \, \left(-\frac{2^5}{4n+1} \right. & {} - \frac1{4n+3} + \frac{2^8}{10n+1} - \frac{2^6}{10n+3} \left. {} - \frac{2^2}{10n+5} - \frac{2^2}{10n+7} + \frac1{10n+9} \right)
\end{align}\quad\mbox{(by Fabrice Bellard:1997)}

最近の円周率の記録は y-cruncher とつかってるのかな＿

一つの ‘‘‘math‘‘‘ 環境の中には align を二つ入れられないんですねぇ。

その他

バーゼル問題

平方数の逆数の和を求める問題が提起されて100年近く経った 1735年オイラーが解くことに成功。→ バーゼル問題の初等的な証明

\frac{\pi^2}6 = \sum_{n=1}^\infty\frac1{n^2} = \frac1{1^2} + \frac1{2^2} + \frac1{3^2} + \cdots \quad\mbox{(by Leonhard Euler:1735)}

arcsin の展開 (建部賢弘)

図形の問題を解くため建部は $(\arcsin^2x)$ の無限級数展開を求めた。円周率を求めたわけではないが、この級数を用いれば以下のような式もできる。

\begin{align}
\arcsin^2x &= x^2 + \frac{2^2}{3\cdot4} x^4 + \frac{2^2\cdot4^2}{3\cdot4\cdot5\cdot6} x^6 + \frac {2^2\cdot4^2\cdot6^2}{3\cdot4\cdot5\cdot6\cdot7\cdot8} x^8 + \cdots \quad\mbox{(by 建部賢弘:1722)}\\
&= \sum_{n=1}^{\infty} \frac{2^{2n-1}((n-1)!)^2}{(2n)!}x^{2n}\\
この式を使うと　\\
\sin \frac\pi2 &=1\\
\frac\pi2&=\arcsin1\\
両辺を2乗して　\\
\frac{\pi^2}4 &=\arcsin^21\\
これを最初の式に代入\\
\frac{\pi^2}4 &= 1 + \frac{2^2}{3\cdot4} + \frac{2^2\cdot4^2}{3\cdot4\cdot5\cdot6} + \frac {2^2\cdot4^2\cdot6^2}{3\cdot4\cdot5\cdot6\cdot7\cdot8} +\cdots\\
&= \sum_{n=1}^{\infty} \frac{2^{2n-1}((n-1)!)^2}{(2n)!}
\end{align}

森本 光生, 小川 束, 建部賢弘の数学, 数学, 2004, 56 巻, 3 号, p. 308-319, 公開日 2008/12/25, Online ISSN 1883-6127, Print ISSN 0039-470X の (9) 式

参考　[円周率の公式集]
(https://web.archive.org/web/20111005043044/http://www.pluto.ai.kyutech.ac.jp/plt/matumoto/pi_small/_pi_small.html)

世界記録を破る“パイ”の作り方 : Google Japan Blog
Google Cloud の良い宣伝ですね(^^;

円周率とarctan型公式 - 小林健太 一橋大学商学研究科

Chrome で編集しているが、数式部分のプレビューがほぼ機能しなくて結構つらい...