円周率 π を精度よく計算する関数。

<関数名>
  pi —- 円周率 π の値を多桁数で算出する

<形式>
  int pi(int *data, int keta);

<引数>
  data (出力) 円周率の値が10進4桁ずつ入った整数の配列
  keta (入力) 小数点以下の桁数の指定

<関数値>
  正常に計算できたときは 0、エラーのときは -1。

用例
  int data[1000];
  pi(data, 1000);

  π の値 (小数点以下1万桁)

<関数本体>
  pi.c

<説明>
  πの値を求めるには、いろいろな式があるが、arctan(逆正接関数)を無限級数で表した式がよく使われる。式が単純な形で扱いやすく、ある程度の桁数の値を得るには適している。しかし欠点は、非常に長い桁数を得るには時間がかかりすぎることである。
  ここで使う式は、最もよく知られている、J. Machinの式
     π/4 = 4 arctan (1/5) – arctan (1/239)
  である。彼は1706年にこの公式でπを小数点以下100桁を求めた。同じ公式で、1949年にはENIACにより2037桁が約70時間で、1958年にはIBM704により1万桁が約100分で求められている。
  arctanをによる古典的な公式にほかに以下のものがある。
  ガウスの公式
     π/4 = 12 arctan (1/18) + 8 arctan (1/57) – 5 arctan (1/239)
  クリンジェンシェルナの公式
     π/4 = 8 arctan (1/10) – arctan (1/239) – 4 arctan (1/515)
  シュテルマーの公式
     π/4 = 6 arctan (1/8) + 2 arctan (1/57) + arctan (1/239)
  高野喜久雄の公式
     π/4 = 12 arctan (1/49) + 32 arctan (1/57) – 5 arctan (1/239) + 12 arctan (1/110443)
  柴田昭彦の公式
     π/4 = 17 arctan (1/22) + 3 arctan (1/172) – 2 arctan (1/682) – 7 arctan (1/5357)

  シュテルマーの公式に8の割り算を工夫すれば高速化期待できるし、高野の公式を利用して並列処理をすれば良い成績が期待できよう。

  arctanを使った公式でπを計算すると、桁数を2倍にするのに計算時間は4倍の増加になるが、以下の方法では2.1〜2.2倍で済む。

  ガウス・ルジャンドルの公式

     A = 1、B = (1/2)1/2、T = 1/4、X = 1

  として、AとBの差が必要とする精度より大きい間、次の計算を繰り返し実行する。

     Y = A、A = (A + B)/2、B = (B * Y)1/2
     T = T – X * (Y-A)2、X = 2 * X

  すると、πの値は(A + B)2/(4 * T) となる。ただし、A、B、T、Yの各値は、求めようとする精度以上で計算しておく必要がある。

Comments are closed.

Post Navigation