円周率.jp (http://xn--w6q13e505b.jp/formula/arctan.html)

$\arctan$ は三角関数 $\tan$ の逆関数であり、 有理数から $\pi$ を求める計算の代表例である。 詳しくはarctan とは？ のページを見てほしい。

[JB07][JW04][JW05] などでたくさん紹介されているので、 ここでは有名なモノ、計算記録に使われたことがあるモノを示す。 なお Machin の公式のように $\arctan$ の引数が整数の逆数であり 項数が 2 つしかない公式は 4 種類しかないことが証明されている。

Machin の公式

$\dfrac{\pi}{4} = 4 \arctan \dfrac{1}{5} - \arctan \dfrac{1}{239}$

Klingenstierna の公式

$\dfrac{\pi}{4} = 8 \arctan \dfrac{1}{10} - \arctan \dfrac{1}{239} - 4 \arctan \dfrac{1}{515}$

Euler の公式 (※)

$\dfrac{\pi}{4} = 5 \arctan \dfrac{1}{7} + 2 \arctan \dfrac{3}{79}$

Euler の公式(2)

$\dfrac{\pi}{4} = 4 \arctan \dfrac{1}{5} - \arctan \dfrac{1}{70} + \arctan \dfrac{1}{99}$

Gauß の公式

$\dfrac{\pi}{4} = 12 \arctan \dfrac{1}{18} + 8 \arctan \dfrac{1}{57} - 5 \arctan\dfrac{1}{239}$

Störmer の公式

$\dfrac{\pi}{4} = 6 \arctan \dfrac{1}{8} + 2 \arctan \dfrac{1}{57} + \arctan \dfrac{1}{239}$

Störmer の公式(2)

$\dfrac{\pi}{4} = 44 \arctan \dfrac{1}{57} + 7 \arctan \dfrac{1}{239} - 12 \arctan \dfrac{1}{682} + 24 \arctan \dfrac{1}{12943}$

高野喜久雄の公式

$\dfrac{\pi}{4} = 12 \arctan \dfrac{1}{49} + 32 \arctan \dfrac{1}{57} - 5 \arctan\dfrac{1}{239} + 12 \arctan \dfrac{1}{110443}$

連鎖探索法 [JB07][JW05][JM05]

$\arctan$ は $\tan$ の逆関数なので

\[ \frac{\pi}{4} = \arctan 1 \]

が成り立つ。また $\tan$ の加法定理から $ab=n^2+1$ を満たす整数 $a$、$b$ を用いると

\[ \arctan \frac{1}{n} = \arctan \frac{1}{n+a} + \arctan \frac{1}{n+b} \]

という形で 1 つの $\arctan$ を 2 つに分けることができるため これを繰り返し適用することで多くの $\arctan$ 系公式を作ることができる。 このとき、$n+a$ や $n+b$ が負になる $a$、$b$ を選ぶと $\arctan$ の引数を正にするため $\arctan$ の係数を負にすることができる。

\[ \begin{eqnarray} \arctan 1 &=& \arctan(1/2) + \arctan(1/3)\\ \arctan(1/2) &=& \arctan(1/3) + \arctan(1/7)\\ \arctan(1/3) &=& \arctan(1/4) + \arctan(1/13) &=& \arctan(1/5) + \arctan(1/8)\\ \arctan(1/7) &=& \arctan(1/8) + \arctan(1/57) &=& \arctan(1/12) + \arctan(1/17)\\ \arctan(1/57) &=& \arctan(1/70) + \arctan(1/307) &=& \arctan(1/32) - \arctan(1/73)\\ \arctan(1/73) &=& \arctan(1/47) - \arctan(1/132) \end{eqnarray} \]

素数探索法 (Störmer の方法[JM05][JT05])

こちらは素因数分解を利用して、複数の公式を纏めて求める事ができる手法である。 背景となっている数学、計算原理などについては特に解説しない。

1. $m^2+1$ の素因数分解表を作る

パラメータとして適当に決めた $M$、$P$ に対して 自然数 $m\lt M$ と素数 $p\lt P$ の範囲内で $m^2+1$ を素因数分解する。 このとき

\[ m^2 + 1 = 2^{\nu} p_1^{e_1} p_2^{e_2} \cdots p_s^{e_s} \]

において $\nu$ は 0 または 1 であり、 $p_j \bmod 4 \equiv 1$ が成り立つ。

この素因数分解において適度に広範囲に渡った分解をしていれば、 同じ分解表を以後の手順で繰り返し利用することができ、 何度も因数分解しなくてすむ。

2. 使う $p$、$m$ の組み合わせを決める

この $n$ 項の $\arctan$ が含まれる公式を作るとしよう。 1. の分解に用いた素数の中から $n-1$ 個 の $p$ を選び出し、 その $p$ (と 2) だけで $m^2+1$ が分解できている $m$ を $n$ 個抜き出す。

具体例として $n=3$ として $\{p\}=\{5, 13\}$、 $\{m\} = \{18, 57, 239\}$ を選び出す。 ちなみにそれぞれの $m^2+1$ の分解は

\[ \begin{eqnarray} 18^2+1 &=& 325 &=& 5^2 \cdot 13\\ 57^2+1 &=& 3250 &=& 2 \cdot 5^3 \cdot 13\\ 239^2+1 &=& 57122 &=& 2 \cdot 13^4 \end{eqnarray} \]

となっている。

3. 指数の行列を作る

2. で抜き出した $\{m\}$ について 各 $m^2+1$ の分解結果である指数を並べ、 $n \times (n-1)$ 行列 $C$ を作る。 各行には 1 つの $m$ についての分解結果が、 各列には 1 つの $p$ についての指数が対応する形になる。 ただし、$m \bmod p \gt p/2$ なら符号を負にする。

上記の例では $(m,p)=(18,5), (239,5)$ で $m \bmod p \gt p/2$ となるが、$(239,5)$ の指数は 0 なので $(18,5)$ の指数だけ符号を負にして

\[ C=\begin{pmatrix} -2 & 1 \\ 3 & 1 \\ 0 & 4 \end{pmatrix} \]

となる。

4. 公式の係数を決める

行列の $i$ 行目に対応する $m=m_i$ を使った $\arctan\dfrac{1}{m_i}$ の係数 $c_i$ は、 行列 $C$ から $i$ 行目を削除した行列 $C_i$ の行列式 (ただし $i$ が偶数だと符号反転) になる。 先の計算例でいえば

\[ c_1 = \det(C_1) = \begin{vmatrix} 3 & 1 \\ 0 & 4 \end{vmatrix} = 12 \] \[ c_2 = -\det(C_2) = -\begin{vmatrix} -2 & 1 \\ 0 & 4 \end{vmatrix} = 8 \] \[ c_3= \det(C_3) = \begin{vmatrix} -2 & 1 \\ 3 & 1 \end{vmatrix} = -5 \]

から

\[ 12\arctan\dfrac{1}{m_1} + 8\arctan\dfrac{1}{m_2} - 5\arctan\dfrac{1}{m_3} \]

という基本形ができる。

5. $k$ を決める

ここまでの手順で決まった式の値 ${\displaystyle \sum_{i=1}^{n} c_i\arctan\dfrac{1}{m_i}}$ は $\dfrac{\pi}{4}$ の整数倍になっている。 そこで、適当な精度 (例えば倍精度浮動小数) で計算して その倍数 $k$ を求める。 先ほどの計算例では

\[ 12\arctan\dfrac{1}{18} + 8\arctan\dfrac{1}{57} - 5\arctan\dfrac{1}{239} = k\dfrac{\pi}{4} \]

とすると、左辺の値が 0.7853981633974$\cdots$ なので $k=1$ となり、 Gauß の公式

\[ \dfrac{\pi}{4} = 12\arctan\dfrac{1}{18} + 8\arctan\dfrac{1}{57} - 5\arctan\dfrac{1}{239} \]

ができたことがわかる。

$\arctan x$ の値は Taylor 展開(Maclaurin 展開) を利用して求めることができる。 具体的な式展開は

\[ \arctan x = \sum_{k=0}^{\infty} \frac{(-1)^k}{2k+1} x^{2k+1} = x - \frac{1}{3} x^3 + \frac{1}{5} x^5 - \frac{1}{7} x^7 + \cdots \]

となる。$x$ を整数の逆数にとった

\[ \arctan \frac{1}{n} = \sum_{k=0}^{\infty} \frac{(-1)^k}{2k+1}\frac{1}{n^{2k+1}} = \frac{1}{n} - \frac{1}{3}\frac{1}{n^3} + \frac{1}{5} \frac{1}{n^5} - \frac{1}{7} \frac{1}{n^7} + \cdots \]

の形であれば計算が容易に (そして $n$ が大きければ収束が速く) なるため、大半の $\arctan$ 系公式は

\[ \frac{\pi}{4} = \sum_{k=0}^n a_k \arctan \frac{1}{b_k} \]

という形になっている。 (ただし Binary Splitting や DRM を使う場合、 整数の逆数であるメリットは特に無い。) なお、Euler の計算では別の形の展開を求め計算したため、 その方法は別に書いている。

コンピュータを用いて計算する定番の方法として、$x$ → $x^3$ → $x^5$ →… を計算していき、 そのそれぞれの段階で $x^3$ → $x^3/3$、$x^5$ → $x^5/5$ を求め、 加算・減算を交互に繰り返すことで求める方法がある。 下に擬似コードでもってその具体的な手順を示す。 この計算手順では $k$ と $n$ 以外の変数に関わる計算を多倍長固定小数点数で行う必要がある。

1. 多倍長の固定小数点数 $A = 0$, $B = x$, $C$ を用意する。
2. 計算したい桁数と $x$ から計算項数 $n$ を求める。
3. $k=0, 1, \cdots, n-1$ で以下を繰り返す。
   1. $B = B / x^2$
   2. $C = B / (2k + 1)$
   3. もし $k$ が偶数なら
      1. $A = A - C$
   4. $k$ が奇数なら
      1. $A = A + C$
4. $A$ を返す

計算項数

上記のアルゴリズムで曖昧にしている計算項数 $n$ の方法について考えよう。 本当なら $C$ に入る値が表現する範囲で 0 になったらそれ以降は計算する必要はない。 が、毎回 0 チェックをするのは計算コストが高いので計算を行う前に必要な計算項数を見積もる必要がある。

計算の例外を無くす前提条件として、必要な項数 $n$ は 4 以上、$x$ は 2 以上の整数の逆数であるとする。 この前提条件自体は $\pi$ 計算においてはキツい制限とはならない。 求める桁数を小数点以下 $d$ 桁としたとき、$n$ の満たす条件は

\[ \frac{1}{(2n+1)x^{2n+1}} \lt 10^{-d} \]

となる。対数をとって計算すると

\[ d \lt \log (2n+1) + (2n+1) \log x \tag{1} \]

ここで $O(\log n) < O(n)$ から $\log(2n+1)$ の項を無視して考えると

\[ \frac{1}{2} \left(\frac{d}{\log x} - 1\right) \lt n \]

Euler は $\arctan x$ を以下のように展開し、 Euler の公式を使って 1 時間程度で 20 桁を求めた。[JB02][FB02]

\[ \begin{eqnarray} \arctan x &=& \dfrac{x}{1+x^2} \left(1+ \sum_{k=1}^{\infty} \prod_{j=1}^{k} \dfrac{2j}{2j+1} \left( \frac{x^2}{1+x^2}\right)^k \right)\\ &=& \frac{y}{x} \left(1+\dfrac{2}{3}y + \dfrac{2\cdot4}{3\cdot5} y^2 + \dfrac{2\cdot4\cdot6}{3\cdot5\cdot7}y^3 + \cdots \right) & \quad ({\rm ただし\ }y = \dfrac{x^2}{1+x^2})\\ &=& \dfrac{y}{x} \left(1 + \dfrac{2}{3}y \left(1 + \dfrac{4}{5}y \left(1 + \dfrac{6}{7}y \left(1 + \cdots \right) \right) \right) \cdots \right) \end{eqnarray} \]

この展開において、x に 1/7 や 3/79 を代入すると

\[ x = \frac{1}{7} \Rightarrow y = \frac{2}{100} \] \[ x = \frac{3}{79} \Rightarrow y = \frac{144}{100000} \]

となるので 10 進数での計算が行いやすい。 なお、この展開を用いた場合に $d$ 桁まで求めるのに必要な計算項数 $n$ は

\[ \begin{eqnarray} \dfrac{1}{x} \prod_{j=1}^{n} \dfrac{2j}{2j+1} \left(\dfrac{x^2}{1+x^2}\right)^{n+1} \lt 10^{-d}\\ \log x + \sum_{j=1}^{n} \log \left(1+ \dfrac{1}{2j}\right) + (n+1) \log \left(1+\dfrac{1}{x^2}\right) \gt d \end{eqnarray} \]

ここで Σ の項を無視して計算すると

\[ \frac{d-\log x}{\log \left(1+\frac{1}{x^2}\right)} - 1 \lt n \]

となる。Σ の項をある程度計算しておけばより正確な項数を求められる。