スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

GMPを使用して自然対数の底を計算するCプログラムの例

expの冪級数展開により自然対数の底を計算するC言語プログラムの例。多倍長計算を行うために、GMP(GNU Multiple Precision Arithmetic Library)を使用しています。

概要

expの冪級数展開

exp(x)=1+x+x^2/2!+x^3/3!+……

を利用して、自然対数の底 e=exp(1) の値を計算します。

プログラム例

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

/* e.c */

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

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

/* exp(x) の x^{last} の項までの冪級数展開 */
void expSeries(mpq_t sum, mpq_t x, unsigned int last);

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

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

  /* 冪級数展開を利用して自然対数の底を計算 */
  mpq_set_ui(x, 1, 1);
  expSeries(res_q, x, LAST);

  /* 分数を小数に変換 */
  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, NULL);
  mpf_clear(res_f);

  return 0;
}

void expSeries(mpq_t sum, mpq_t x, unsigned int last)
{
  mpq_t tmp;
  unsigned int i;

  mpq_init(tmp);

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

  mpq_clear(tmp);
}

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

0.27182818284590452353602874713526624977……
(中略)
……788567430286e1

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

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

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

関連記事

Windows XPにおけるGMPのセットアップ作業のメモ
GMPを使用して円周率πを計算するCプログラムの例

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

プロフィール

よしいず

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

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

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

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

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

リンク
月別アーカイブ
カテゴリ
最新記事
検索フォーム
RSSリンクの表示
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。