AOYAMA Koji's プログラミングブログ - プログラミングを楽しく体験

円周率算出の高速化

2025/04/20
円周率算出の高速化

　本記事では、この記事で円周率を算出するために用いた、 高速に円周率を求める式について、できるだけ簡潔に解説します。 図を用いることで、公式の意味を視覚的に理解しやすいようにしています。
　この記事で解説したことが前提になっていますので、ご注意ください。



[PR]

マチンの公式


　結論としては、以下のマチンの公式を用いました。
\( \frac{\pi}{4} = 4\arctan\frac{1}{5} - \arctan\frac{1}{239} \)

マクローリン展開


　\(\arctan x\)のマクローリン展開は以下になります。
\( \arctan x = x - \frac{1}{3}x^{3} + \frac{1}{5}x^{5} - \frac{1}{7}x^{7} \cdots \)

当てはめる


　元の式に当てはめて4倍すると以下が得られます。
\( \pi = 4(\frac{4}{5} - \frac{1}{239}) - \frac{4}{3}(\frac{4}{5^{3}} - \frac{1}{239^{3}}) + \frac{4}{5}(\frac{4}{5^{5}} - \frac{1}{239^{5}}) - \frac{4}{7}(\frac{4}{5^{7}} - \frac{1}{239^{7}}) \cdots \)

　この式は、この記事で解説したライプニッツの公式に比べて、項が小さくなっていく速度が圧倒的に早いです。 そのため収束が早く、高速に、指定の桁数までの計算結果が求まります。
　あとはこれを、できるだけ再計算が少なくなるよう、プログラミングします。
[PR]

式の意味


　マチンの公式の意味するところを図で解説します。

\(\arctan\frac{1}{5}\)の三角形


円周率算出の高速化 fig.2
　\(\arctan\frac{1}{5}\) は、直角三角形で表すとこの図になります。 左下の角度が \(\arctan\frac{1}{5}\) です。

4つ並べる


円周率算出の高速化 fig.3
　それを4つ重ねていくとこうなります。 角度は合計して \( 4\arctan\frac{1}{5} \) になります。
　最後に並べた三角形の頂点は、右斜45度に近いところにあるのがわかります。

時計回りに45度回転させる


円周率算出の高速化
　わかりやすくなるように、図の直角三角形4つを、時計回りに45度回転させて、その頂点の座標を計算してみます。 すると…なんということでしょう！ その頂点のx座標とy座標の比率が 239:1 という、整数比になりました。

\(\arctan\frac{1}{239}\)の三角形が登場


円周率算出の高速化 fig.4
　つまり、\(\arctan\frac{1}{5}\)の三角形を4つ積み重ねると、45度に対して、この図のような直角三角形ぶん行き過ぎます。 これは、左下の角度が \(\arctan\frac{1}{239}\) の直角三角形です。 この数字、ついに登場しました。
　なお、かなり横に細長い三角形になるため、ここでは、比率を変更して図示しています。

マチンの公式(再掲)


　以上がマチンの公式の意味です。 公式を再掲します。 角度 \(\arctan\frac{1}{5}\) ぶん4回回転させて、\(\arctan\frac{1}{239} \) ぶん戻すと45度、すなわち \(\frac{\pi}{4}\) になる、というのがこの公式です。
\( \frac{\pi}{4} = 4\arctan\frac{1}{5} - \arctan\frac{1}{239} \)

通常はなかなか整数比にならない


円周率算出の高速化 fig.5
　直角三角形の横線と縦線の長さの比率は、通常、整数比にはなりません。
　この図の直角三角形の左下の角度は30度ですが、\(\sqrt{3}:1\) です。 これを用いて \(\pi\) を求める場合は、先に \(\sqrt{3}\) を、高精度で算出する必要があり、計算量が増えます。
　マチンさんは、整数比率になる組み合わせをよく発見したなあ、、、と驚きます。

回転行列について


　各座標を回転させるには、回転行列を用いることで、その座標が計算できます。
　以下に、\((x_0,y_0)\)を、原点まわりに角度\(\theta\)だけ回転させて\((x_1,y_1)\)を求める、回転行列の図と式を記します。 ここでは式と例を提示するだけに留めますが、本記事の内容を、ご自身で計算したい場合には、ぜひ、ご確認ください。
円周率算出の高速化 fig.6
\(\begin{pmatrix} x_1 \\ y_1 \end{pmatrix} = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} x_0 \\ y_0 \end{pmatrix}\)


マルチスレッド


　本記事の元になったベンチマークアプリケーションでは、マルチスレッドによる実装も行っています。
　本記事で解説している式の各項は、他の項の影響を受けないので、別々に計算することができます。 それを利用して、複数の項をグループ化し、分散して計算します。
　マルチスレッドについては、題材は別ですがこの記事で詳しく解説していますので、併せてご覧ください。
[PR]

まとめ


　本記事では、マチンの公式を用いた、高速に円周率を求める式について、簡潔に解説しました。 図を用いることで、式に登場する\(\arctan\frac{1}{239}\)などの意味が、直感的に理解できるようにしています。

補足


・\(\arctan x\)のマクローリン展開は、厳密には、\(|x|\le{1}\)のときのみですが、今回は条件に合致します。
・画像内のラスタライズ文字フォントにOpen Font LicenseのZen Antiqueを使用しております。
・数式表現にMathJaxを使用しております。助かります！

カテゴリー：プログラミング解説,円周率計算
[PR]