この記事は 2024/4/13 に公開した以下の動画の内容を簡易的にまとめたものです。
この動画では主に
- マチンの公式の証明
- 円周率の計算の仕方
- マチンの公式の見つけ方
について説明しています。
動画作成のためにマチンの公式の調査をしていて、いくつか気になることがありました。一つ目は wikipedia のマチンの公式の記事において、円周率の計算を
$$16 \sum_{n=0}^{3m+2} \frac{(-1)^n}{2n+1} \left(\frac{1}{5}\right)^{2n+1} -\sum_{n=0}^{m} \frac{(-1)^n}{2n+1} \left(\frac{1}{239}\right)^{2n+1}$$
としていることです。問題は和をとる項が $3m+2$ までと $m$ までで、それぞれ異なっていることです。
もう一つは、マチンの公式の見つけ方が見当たらなかったことです。正確には、ガウス整数の素因数分解を用いた方法は見当たりましたが、年代的に見てマチンがその方法を取ったとは考えにくいです。
私が調べた限りだと、マチンの公式についての情報はだいぶ偏りがあり、痒いところに手が届かない状態になっていました。
前者は誤差を評価することで解決し、後者はマチンの時代でも思いつきそうな方法を見つけたので解決し、それぞれ動画中で説明しています。動画だけだと検索に引っかかりにくいと思われるので、この記事でも簡易的に解説しようと思います。
マチンの公式について
イギリスの天文学者ジョン・マチン (1680頃 – 1751) は以下の公式
$$\frac{\pi}{4} = 4\arctan \frac{1}{5} -\arctan \frac{1}{239}$$
と、$\arctan x$ のテイラー展開
$$\arctan x = \sum_{n=0}^{\infty} \frac{(-1)^n}{2n+1} x^{2n+1} \quad (|x| \leq 1)$$
を用いて円周率を100桁目まで求めたといわれています。
マチンの公式の証明
マチンの公式の証明は $\arctan$ の加法定理
$$\arctan a \pm \arctan b = \arctan \frac{a \pm b}{1 \mp ab} \quad (|ab| < 1)$$
を用いて
$$4 \arctan \frac{1}{5} + \arctan \frac{1}{239}$$
を計算することでなされます。具体的な計算は動画をご覧ください。$\arctan$ の加法定理は $\tan$ の加法定理
$$\tan (\alpha \pm \beta) = \frac{\tan \alpha \pm \tan \beta}{1 \mp \tan \alpha \tan \beta}$$
の単純な書き換えです。
記事を書いている途中で図的に $\arctan$ の加法定理を証明する方法の解説を見つけました。
https://twitter.com/genkuroki/status/1024160463729225733
円周率の計算方法
$\arctan x$ は、
$$\arctan x = \int_0^x \frac{1}{1 + t^2} dt$$
と
$$\frac{1}{1 +t^2} = \sum_{k = 0}^n (-1)^k t^{2k} + (-1)^{n+1} \frac{t^{2n+2}}{1 +t^2}$$
を用いて
$$\arctan x = \sum_{k = 0}^n \frac{(-1)^k}{2k+1} x^{2 k + 1} + (-1)^{n+1} \int_0^x \frac{t^{2n+2}}{1 +t^2} dt $$
と表されます。これにより
$$\frac{1}{1 +x^2} \frac{x^{2n+3} }{2n + 3} \leq \left| \arctan x -\sum_{k = 0}^n \frac{(-1)^k}{2k+1} x^{2 k + 1} \right| \leq \frac{x^{2n+3} }{2n + 3}$$
がわかります。マチンの公式
$$\pi = 16\arctan \frac{1}{5} -4\arctan \frac{1}{239}$$
を用いて計算するとき、$ \frac{1}{5}$ の項と $\frac{1}{239}$ の項の両方の誤差を小さくする必要がありますが、$\frac{1}{239}$ の項の方が誤差の減りが速いので、$\frac{1}{239}$ の方が計算する項が少なくて済みます。
\begin{align}& \log_{10}\left(\frac{1}{5}\right) \fallingdotseq -0.69897 \\ &\log_{10}\left(\frac{1}{239}\right) \fallingdotseq -2.3783979 \end{align}
なので、$\frac{1}{239}$ の項の方がおおよそ $10^{3.4}$ 倍誤差の減りが速いです。
円周率を100桁まで計算する場合は、目安として
$$\frac{16}{2n_1 + 3} \frac{1}{5^{2n_1+3}} + \frac{4}{2n_2 + 3} \frac{1}{239^{2n_1+3}} < 10^{-100}$$
を満たす最小の $n_1, n_2$ を求めて実際に計算し、有効な桁数が足りなければもう少し計算します。具体的な計算方法は動画を見るか、以下の記事を見てください。
マチンの公式の見つけ方
$k$ を適当な自然数として、$\arctan$ の加法定理を用いて $2 \arctan \frac{1}{k}$, $3 \arctan \frac{1}{k}$, … と計算していき、
$$m \arctan \frac{1}{k} \leq \frac{\pi}{4} \leq (m+1) \arctan \frac{1}{k}$$
を満たす $m$ を探します。次に
\begin{align}\arctan 1 -m\arctan \frac{1}{k} &= \arctan \frac{c_m}{d_m} \\ (m+1)\arctan \frac{1}{k} -\arctan 1 &= \arctan \frac{c_{m+1}}{d_{m+1}}\end{align}
を計算して、整理すると $\arctan$ の円周率公式が得られます。$\frac{c_m}{d_m}$ と $\frac{c_{m+1}}{d_{m+1}}$ のうち小さい方が、収束が速いです。動画では図を用いて説明しています。
$k = 5$ として計算すれば、マチンの公式が得られます。マチンがこの方法でマチンの公式を発見したのかは分かりませんが、マチンは $k = 2, 3, 4$ の場合の公式も知っていたようですし、$\arctan$ の加法定理は三平方の定理さえ知っていれば導けるので、この方法を用いた可能性が高いと思われます。
この方法はプログラムで求めることができます。
参考文献
[中村] 中村 滋. 天文学者マチンが見つけた円周率公式たち