円周率公式のまとめ
今日の話は円周率公式の話です。
円周率公式と一口に言っても、円周率に関わる公式はたくさんあります。
その中でも今回は、コンピューターで円周率を計算する際に実際に使われている、あるいは使われていたものについて取り扱っていきたいと思います。
ただ紹介するだけでは味気ないので、詳細な性質・証明も添えていきます。
意外と円周率公式の紹介みたいな記事はあっても、深い内容まで書いてあるものってあまり無いんですよね。
そんなわけで、できる限りself-containedを目指したいと思います。
1. Machinの公式
Machinの公式の発見は1706年のことだそうですが、まずはそれ以前に用いられていた関係式を紹介します。それがライプニッツの公式です。
\begin{align}
\frac{\pi}{4}=\sum_{n=0}^{\infty} \frac{(-1)^{n}}{2 n+1}=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9} - \cdots
\end{align}
式の形としてはとても綺麗で、また証明も簡潔です。
この式の正体は のマクローリン展開です。
\begin{align}
\sum_{n=0}^{\infty} \frac{(-1)^n}{2 n + 1} x^{2 n+1} = x - \frac{1}{3} x^{3}+\frac{1}{5} x^{5} - \frac{1}{7} x^{7} + \frac{1}{9} x^{9} - \cdots \quad (|x|<1)
\end{align}
上に書いてあるとおり収束半径が1なので を代入できるか議論が必要ですが、この場合は代入することができ*1ライプニッツの公式を得ます。
ただこの公式の欠点は形がキレイなものの収束が非常に遅いことにあります。
もっと速く正確に円周率を計算することはできないか、と人々は考えるわけです。
その一つの解決策として現れたのが、Machinの公式でした。
この式は、ライプニッツの公式のある意味での応用によって作られます。
証明はガウス整数（2つの有理整数a,bを用いてa+biのように書ける複素数）を用いることによって（ライプニッツの公式より少し難しくなりますが）簡単に得ることができます。また、証明を追うことによってなぜ5だったり239といった一見意味のなさそうな整数が出てくるのかも理解できるので、一度読んでみてください。
Machinの公式の証明
図形的に考えるために、複素数平面を用いて証明を行います。
以下 は虚数単位です。
したがって、偏角を考えることによって、Machinの公式は以下と同値になります。
\begin{align}
\frac{(5+i)^4 (1-i)}{239+i} \in \mathbb{R}
\end{align}
は
回転させるためだけのもので、本質的には
の偏角が
であればよい、ということです。
ここで と
をガウス整数環の範囲で素元に分解します。
\begin{align}
5 + i &= (3-2i)(1+i) \\
239 + i &= (3-2i)^4 (1+i) i
\end{align}
なんと、素元分解すると同じ素元が含まれていることが分かります。
これが と
が現れる理由でした！
を4乗すると
が4つになってちょうど割ると打ち消し合ってくれるわけです。
実際に計算してみると、以下のようになります。
\begin{align}
\frac{(5+i)^4 (1-i)}{239+i} = \frac{(3-2i)^4 (1+i)^4 (1-i)}{(3-2i)^4 (1+i) i} = 4
\end{align}
初めの仮説通り、実数になったのでMachinの公式の証明ができました。
この議論を応用すれば自分でMachin-like formulaを生成することもできるのですが、それはまた別の機会にでもしましょう。
（ヒントはノルムと素因数分解です）
さて、Machinの公式はライプニッツの公式の収束の遅さを解決できると書きました。
実際に計算するとライプニッツの公式は本当に遅く、マクローリン級数のうち600項以上を足さないと3.14まですらたどり着きません。
（実際にPythonで計算してみたら620項ほど必要でした）
一方で、Machinの公式の方は3.14までならたったの2項で事足りてしまいます。
この差は何から生まれるのかというと、arctanの引数の大きさです。
ライプニッツの公式では、arctanの引数は1でした。1は何乗しようが1です。
しかし、Machinの公式ではそれぞれ1/5と1/239で1より小さいため、べき乗を繰り返していけば小さくなっていきます。
また、分母が大きいほど0に近づくスピートも速くなるので、できるだけ分母は大きいほうがいいことが分かります。
Machin自身はこの公式を用いて円周率を100桁まで求めることに成功しています。
この公式の類似は非常に多くの形が見つかっていて、中でも以下の高野喜久雄による公式は2002年の円周率1兆桁計算に用いられました。
「おや？」と思った人もいるかもしれませんが、そうです、このarctan型の公式は20年ほど前まで使われていたんです。実際にこの型が主流だったので1970年ごろまでだそうですが、意外と最近の数学なんです。
学校などでは数百年前くらいまでの数学しかほとんどやりませんから、特別数学に馴染みのない人には意外に聞こえるかもしれません。
2. ラマヌジャン型公式
さて、Machin型の次に主流となり今でも用いられているのが、みんな大好きラマヌジャン型の円周率公式です。
\begin{align}
\frac{1}{\pi} = \frac{2 \sqrt{2}}{99^{2}} \sum_{n=0}^{\infty} \frac{(4n)!}{(n!)^4} \frac{1103+26390n}{396^{4n}}
\end{align}
上の式は1914年にラマヌジャンによって発見されたものです。
残念ながら証明は前の2つよりも格段に難しいですが、解説記事を書いていますのでよかったら読んでいただけると嬉しいです。（証明はありませんが……）
ラマヌジャン自身が発見した公式もとても収束が速く1項ごとに8桁ずつ正確に計算が進むのですが、実際に計算に用いられているのはこれをより高速にしたChudnovsky兄弟の公式です。1994年に発表されました。
\begin{align}
\frac{1}{\pi} = 12 \sum_{n=0}^{\infty} (-1)^{n} \frac{(6 n) ! (13591409 + 545140134 n)}{(3 n) !(n !)^{3} (640320^3)^{n+1/2}}
\end{align}
これは、上の解説記事で言う の場合で、
とすると得られます。
こちらの方は、1項計算するたびに14桁ずつ進んでいきます。最初のライプニッツの公式と比べたらかなり高速です。
Chudnovsky兄弟の公式の発表の少し前、Borwein兄弟という人たちもラマヌジャン型円周率公式を発表しました。
Borwein兄弟はラマヌジャンの円周率公式がモジュラー関数と関係することを証明し、Chudnovsky兄弟はその功績のおかげで公式を改良できたのです。
Borwein兄弟が発見した公式を参考までに載せておきます。（1989）
\begin{align}
\frac{1}{\pi} = 12 \sum_{n=0}^{\infty} (-1)^{n} \frac{(6 n) !(A+B n)}{(3 n) !(n !)^{3} C^{n+1/2}}
\end{align}
\begin{align}
\begin{array}{l}{A=1657145277365+212175710912 \sqrt{61}} \\ {B=107578229802750+13773980892672 \sqrt{61}} \\ {C=(5280 \times (236674+30303 \sqrt{61}))^{3}}\end{array}
\end{align}
こちらは一項ごとに25桁進むそうですが、あまり主流ではないようです。
無理数があまり機械的には嬉しくないのでしょうか。 の場合なことは明らかですが
に何を代入したかはパッと見ではわかりませんでした……。
3. BBP型公式
略さずに言うと「ベイリー（Bailey）・ボールウェイン（Borwein）・プラウフ（Plouffe）の公式」です。
Borwein兄弟はこちらでも登場です。
どういう公式かと端的に言えば、「円周率を2(あるいは16)進数表記にしたときの小数点以下第n桁目をそれまでの数字を計算することなく求める」式です。
計算の現場では、検算に用いられているようです。
\begin{align}
\pi = \sum_{k=0}^{\infty} \frac{1}{16^{k}}\left(\frac{4}{8 k+1}-\frac{2}{8 k+4}-\frac{1}{8 k+5}-\frac{1}{8 k+6}\right)
\end{align}
まだ勉強中なため詳しくは話せませんが、証明自体は高校範囲の積分計算で済むそうです。ただ、そもそもの発想の部分は未だ謎で、誰かご存じの方いらっしゃれば教えていただきたいです。
\begin{align}
\pi = \frac{1}{2^{6}} \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2^{10 n}}\left(-\frac{2^{5}}{4 n+1}-\frac{1}{4 n+3}+\frac{2^{8}}{10 n+1}-\frac{2^{6}}{10 n+3}-\frac{2^{2}}{10 n+5}-\frac{2^{2}}{10 n+7}+\frac{1}{10 n+9}\right)
\end{align}
参考
*1:アーベルの連続性定理により示される
