円周率.jp > 公式 > arctan系公式

arctan系公式

arctan は三角関数 tan の逆関数であり, 有理数からπを求める計算の代表例である. 詳しくはarctan とは? のページを見てほしい.

代表例

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

Machin の公式
\frac{\pi}{4} = 4 arctan \frac{1}{5} − arctan \frac{1}{239}
Klingenstierna の公式
\frac{\pi}{4} = 8 arctan \frac{1}{10} − arctan \frac{1}{239} − 4 arctan \frac{1}{515}
Eular の公式 ()
\frac{\pi}{4} = 5 arctan \frac{1}{7} + 2 arctan \frac{3}{79}
Eular の公式(2)
\frac{\pi}{4} = 4 arctan \frac{1}{5} − arctan \frac{1}{70} + arctan \frac{1}{99}
Gauß の公式
\frac{\pi}{4} = 12 arctan \frac{1}{18} + 8 arctan \frac{1}{57} − 5 arctan \frac{1}{239}
Störmer の公式
\frac{\pi}{4} = 6 arctan \frac{1}{8} + 2 arctan \frac{1}{57} + arctan \frac{1}{239}
Störmer の公式(2)
\frac{\pi}{4} = 44 arctan \frac{1}{57} + 7 arctan \frac{1}{239} − 12 arctan \frac{1}{682} + 24 arctan \frac{1}{12943}
高野喜久雄の公式
\frac{\pi}{4} = 12 arctan \frac{1}{49} + 32 arctan \frac{1}{57} − 5 arctan \frac{1}{239} + 12 arctan \frac{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 を素因数分解する. このとき,

m2+1 = 2ν p1e1 p2e1ps es

において,ν は 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\frac{1}{m_i} の係数 ci は, 行列 A から i 行目を削除した行列 Ai の行列式 (ただし i が偶数だと符号反転) になる.

5. k を決める

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

計算方法

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

arctan x\sum_{k=0}%5e{\infty} (-1)k\frac{1}{2k+1} x2k+1x\frac{1}{3} x3\frac{1}{5} x5\frac{1}{7} x7 + …

となるが,x を整数の逆数にとった

arctan \frac{1}{n}\sum_{k=0}%5e{\infty} (-1)k\frac{1}{2k+1} \frac{1}{n%5e{2k+1}}\frac{1}{n}\frac{1}{3} \frac{1}{n%5e3}\frac{1}{5} \frac{1}{n%5e5}\frac{1}{7} \frac{1}{n%5e7} + …

の形だと計算が容易に(且つ,n が大きければ収束が速く) なるため,大半のarctan系公式は

\frac{\pi}{4}\sum_{k=0}%5en ak・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 & 1 == 1 ) {
       A = A - C;
     } else {
       A = A + C;
     }
  }
  return A;
}

一応動くサンプルプログラムを置いておく.コンパイルの仕方や実行手順,桁数の設定についての質問は受け付けないのでご自身でコードを読んで理解してほしい.:long.c long.h arctan.c

計算項数

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

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

10-d\frac{1}{2n+1} x2n+1

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

d < log (2n + 1) − (2n + 1) log x     …(1)

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

\frac{1}{2} \Bigl( \frac{d}{\log x} + 1\Bigr)n
となる. この条件を満たす n は符号の関係から必ず (1) 式を満たすので元の条件から外れない.

Eular の展開

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

arctan x \frac{x}{1+x%5e2} \Bigl( 1 + \sum_{k=1}%5e{\infty} \prod_{j=1}%5e{k} \frac{2j}{2j+1} \Bigl(\frac{x%5e2}{1+x%5e2}\Bigr) k \Bigr)
\frac{y}{x} \Bigl( 1+ \frac{2}{3} y\frac{2\cdot4}{3\cdot5} y2\frac{2\cdot4\cdot6}{3\cdot5\cdot7} y3\Bigr)     ただし y\frac{x%5e2}{1+x%5e2}
\frac{y}{x} \Bigl( 1+ \frac{2}{3} y \Bigl( 1+ \frac{4}{5} y \Bigl( 1+ \frac{6}{7} y \Bigl( 1+ … \Bigr) \Bigr)\Bigr)

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

x\frac{1}{7}y\frac{2}{100} ,   x\frac{3}{79}y\frac{144}{100000}

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

\frac{1}{x} \prod_{j=1}%5e{n} \frac{2j}{2j+1} \Bigl(\frac{x%5e2}{1+x%5e2}\Bigr) n+1 < 10-d
log x\sum_{j=1}%5e{n} log \Bigl(1+\frac{1}{2j}\Bigr) + (n + 1) log \Bigl(1+\frac{1}{x%5e2}\Bigr)d

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

d−log x\Bigl/ log \Bigl(1+\frac{1}{x%5e2}\Bigr) − 1 < n

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