円周率.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} $
Eular の公式 ()
$\dfrac{\pi}{4} = 5 \arctan \dfrac{1}{7} + 2 \arctan \dfrac{3}{79}$
Eular の公式(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 の加法定理から abn2+1 を満たす整数 ab を用いると

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

という形で 1 つの arctan 項を 2 つに分けることができるため, これを繰り返し適用することで任意の arctan 系公式を作ることができる. このとき, nanb が負になる ab を選ぶと, arctan の引数を正にするため arctan の係数を負にすることができる.

構成例
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)

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

1. m2+1 の素因数分解表を作る

m と素数 p について探索する範囲を定め, その範囲内で m2+1 を素因数分解する. このとき,

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

において,ν は 0 または 1 であり, pj mod 4 ≡ 1 が成り立つことを利用すると速い.

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

2. 使う m を決める

1. の分解に用いた素数の中から S 個 (任意に決める) の p を選び出し, その p (と 2)だけで m2+1 が分解できている mN (>S) 個抜き出す. 厳密には,各 p について,「その p の指数 e が 1 以上になっている」 m が 2 個以上ある必要がある. 条件を満たす mS 個以下の場合は p の組み合わせから作り直し.

3. 指数の行列を作る

2. で抜き出した mp について, 各 m2+1 の分解結果である指数を並べ, N×S 行列 A を作る. 各行には 1 つの m についての分解結果を, 各列には 1 つの p についての指数 e を並べているものとする. このルールさえ守っていれば順序について制限はない.

自然数の素因数分解結果なので全数とも 0 以上の値になっている, というより本当は符号が不明なので符号を決める. 符号を決めようとする e≠0 に対応する pm について, m mod pp / 2 なら e < 0, m mod pp / 2 なら e > 0 とする.

4. 係数を決める

行列の i 行目に対応する mmi を使った $\arctan\dfrac{1}{m_i}$ の係数 ci は, 行列 A から i 行目を削除した行列 Ai の行列式 (ただし i が偶数だと符号反転) になる.

5. k を決める

ここまでの手順で決まった値を使うと, $c_1\arctan\dfrac{1}{m_1} + c_2\arctan\dfrac{1}{m_2} + \cdots + c_N\arctan\dfrac{1}{m_N} = \dfrac{k\pi}{4}$ になるので,適当な精度(倍精度浮動小数で十分) で左辺を計算して k を求める.

計算方法

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

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

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

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) 式を満たすので元の条件から外れない.

Eular の展開

Eular は arctan x を以下のように展開し, Eular の公式を使って 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 \]

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