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

arctan系公式

arctan は三角関数 tan の逆関数であり, 有理数からπを求める計算の代表例である. 詳しくは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 の逆関数なので

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

が成り立つ.また,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$ の組み合わせを決める

1. の分解に用いた素数の中から $n-1$ 個 の $p$ を選び出し, その $p$ (と 2) だけで $m^2+1$ が分解できている $m$ を $n$ 個抜き出す。この $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$ の値は,arctan の 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} \]

という形になっている. なお,Euler の計算では別の形の展開を求め計算したため, その方法は別に書いている.

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

function arctan(int x) {
  long float A = 0, B = x, C;
  for (int k = 0; k < n; ++k) {
    B = B / x; B = B / x;   // もしくは B = B / (x * x);
    C = B / (2 * k + 1);
    if (k % 2 == 1) {
      A = A - C;
    } else {
      A = A + C;
    }
  }
  return A;
}

計算項数

上記のプログラム例で表している n は計算項数を表している. 本来なら c に入る値が表現する範囲で 0 になったらそれ以降は計算する必要はないが,毎回 0 チェックをするのは計算コストが高いため,計算を行う前に必要な計算項数を見積もる必要がある.

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

\[ 10^{-d} \gt \frac{1}{2n+1} x^{2n+1} \]

を満たす n を求める. まずは両方の常用対数を取る.(ここで log は log10 を示す)

\[ 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 \]
となる. この条件を満たす n は符号の関係から必ず (1) 式を満たすので元の条件から外れない.

Euler の展開

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 \]

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