GMPを使用して円周率πを計算するCプログラムの例

Machinの公式により円周率πを計算するC言語プログラムの例。多倍長計算を行うために、GMP(GNU Multiple Precision Arithmetic Library)を使用しています。

Machinの公式とは

Machin(マチン)の公式とは、以下のような円周率とarctan(x)との関係式です:

π/4=4*arctan(1/5)-arctan(1/239)

この式と、arctanの冪級数展開

arctan(x)=x-x^3/3+x^5/5-x^7/7……

を利用して、πの値を計算します。

プログラム例

以下のプログラムをコンパイルして実行すると、円周率を小数点以下10,000桁程度まで計算します。

/* machin.c */

#include <stdio.h>
#include <gmp.h>

#define BASE 10     /* 基数(10進数) */
#define LAST 7200   /* arctan の冪級数展開の項数 */
#define PREC 33220  /* 実数計算の精度(bit) */

/* arctan(x) の x^{2*last-1} の項までの冪級数展開 */
void arcTan(mpq_t sum, mpq_t x, unsigned int last);

int main(void)
{
  mpq_t res_q, x, atan, term, coeff;
  mpf_t res_f;

  /* 初期化 */
  mpq_inits(res_q, x, atan, term, coeff, NULL);

  /* Machinの公式を利用して円周率を計算 */
  mpq_set_ui(x, 1, 5);          /* x = 1/5 */
  mpq_set_si(coeff, 4 * 4, 1);  /* coeff = 4 * 4 */
  arcTan(atan, x, LAST);
  mpq_mul(term, coeff, atan);
  mpq_add(res_q, res_q, term);

  mpq_set_ui(x, 1, 239);           /* x = 1/239 */
  mpq_set_si(coeff, 4 * (-1), 1);  /* coeff = 4 * (-1) */
  arcTan(atan, x, LAST);
  mpq_mul(term, coeff, atan);
  mpq_add(res_q, res_q, term);  

  /* 分数を小数に変換 */
  mpf_set_default_prec(PREC);  /* 実数計算の精度を指定 */
  mpf_init(res_f);
  mpf_set_q(res_f, res_q);

  /* 結果出力 */
  mpf_out_str(stdout, BASE, 0, res_f);

  /* 後始末 */
  mpq_clears(res_q, x, atan, term, coeff, NULL);
  mpf_clear(res_f);

  return 0;
}

void arcTan(mpq_t sum, mpq_t x, unsigned int last)
{
  mpq_t term, sq_x, coeff;
  unsigned int i;

  mpq_inits(term, sq_x, coeff, NULL);

  /* Horner法で級数の和を後ろから計算 */
  mpq_mul(sq_x, x, x);
  mpq_set_ui(sum, 0, 1);
  for(i = 2*last-1; i > 1; i -= 2) {
    mpq_set_ui(coeff, i, 1);
    mpq_div(term, x, coeff);
    if (i % 4 == 3) {
      mpq_neg(term, term);  /* term = -term */
    }
    mpq_add(sum, sum, term);
    mpq_mul(sum, sum, sq_x);
  }
  mpq_add(sum, x, sum);  /* 最初の項を加える */

  mpq_clears(term, sq_x, coeff, NULL);
}

出力結果は以下のようになります:

0.31415926535897932384626433832795028841971693……
(中略)
……678566722797e1

GMPは有理数を扱えるので、上のプログラムでは有理数の範囲で計算しています。そして、最終的な計算結果を有理数から小数に変換して表示しています。

GMPのmpf_out_str関数で小数を出力するとき、"0."から始まります。最後の"e1"は、ここでは10倍という意味です。

最後の桁は誤差が発生することがあります。上の出力結果でも発生しています(「…7966…」→「…797」)。なお、小数点以下10,000桁の最後は「…678」であり、そこまでは合っています。

関連記事

Windows XPにおけるGMPのセットアップ作業のメモ
GMPを使用して自然対数の底を計算するCプログラムの例

【theme : プログラミング
【genre : コンピュータ

プロフィール

よしいず

Author:よしいず
MATHEMATICS.PDFというウェブサイトを運営しています。

管理の都合上、トラックバックとコメントはオフにしてあります。ブログ経験者なら分かっていただけると思いますが、スパム(アダルトやその他の宣伝)ばかりなのが現実です。

リンクは自由です。当サイトの記事に対する間違いの指摘・意見・感想などを述べた記事からのリンクは歓迎です。ただし、ブログ記事アップ直後はミスが多く、頻繁に修正します。場合によっては削除する可能性もあります。その際、何も断りもなく修正・削除しますがご了承ください。内容を参考にする場合には投稿後一週間ほど様子を見てからにしてください(笑)。

記事の間違いを指摘するときは、その具体的箇所、理由(仕様に反するなど)・根拠(参考にした文献など)、代替案(同じ結果を得るための正しいやり方)も教えてください。そうしないと、(指摘される側および第三者はその時点では無知の状態なので、)どこが間違いなのか分かりませんし、本当に間違っているのかどうかが判断・検証できません。実際、間違いだと指摘されたことが結局は正しかったというケースもありますので。

このブログのタイトル一覧

リンク
月別アーカイブ
カテゴリ
最新記事
検索フォーム
RSSリンクの表示