誤差について

[いろは]

最近C++で数値解析をしてみようと思い、一先ず数学関数を自作してみようとしました。
三角関数等はかなりいい精度が取れるのですが、指数関数では誤差が発生し精度があまりよくありません。
一応べき級数と極限の2つの場合で計算しているのですが、極限での実装は流石に避けたいです。

これ以上どう精度に対してどう対策をすればいいか、また、他の計算方法があれば教えてください。

コード:

#include <cmath>
#include <stdio.h>

//打切り誤差
#define MATHLIB_EPSILON			1e-15

namespace Math {
	double exp1(double x) {
		double x1 = 1.0, x2 = 1.0, buf = 0;
		for (int i = 1;; ++i) {
			buf = x2;
			x1 = x1*(x / i);
			x2 += x1;
			if (std::abs(x2 - buf) < MATHLIB_EPSILON)
				break;
		}
		return x2;
	}
	double exp2(double x) {
		return std::pow(1.0 + x / 10000000, 10000000);
	}
}

void main() {

	for (int i = -20; i < 21; ++i) {
		printf("%4d:%25.15lf\n     %25.15lf\n", i, Math::exp1(i), Math::exp2(i));
		printf("     %25.15lf\n", std::exp(i));
	}

	getchar();
	return;
}

[みけCAT]

オフトピック

とりあえず実行してみようと思いましたが、
C言語としてはcmathという非標準のヘッダが使用されているためコンパイルできず、
C++としてはグローバルで戻り値がintでないmain関数が定義されているためコンパイル出来ませんでした。
C++で、"-Dmain=mein();int main(){mein();}void mein"というオプションを使用すればコンパイルできました。

[みえ]

いい試みですね。円周率については簡単なものを自分で書いてみたことがあるのですが、ネイピア数を試したことはなかったです。

これ以上どう精度に対してどう対策をすればいいか、また、他の計算方法があれば教えてください。

以下のサイトによると、去年 e を 5 兆桁計算した人は e^1 のテイラー展開を使ったようなので、e そのものを計算するアルゴリズムとしては Math::exp1 も悪くなさそうです。

e
http://www.numberworld.org/digits/E/

収束が遅いアルゴリズムや、そもそも得たい結果が double の精度を超えていると double は使えません。このような場合は実数を保持するデータ構造を自分で書く必要があります。多倍長演算とか任意精度演算とか呼ばれるテクニックです。以下は wikipedia ですが、検索するとほかのサイトもたくさんヒットします。基本的には整数の計算をうまく積み重ねることで実現しているはずです。イメージ言えば、人間が筆算で手計算をするような感じでしょうか。

任意精度演算 - Wikipedia
https://ja.wikipedia.org/wiki/%E4%BB%BB ... 4%E7%AE%97

e の値はそれで極めるとして、より現実的な実装であるはずの std::exp(double) の計算方法も気になったので調べてみました。

Ubuntu の g++ 5.4.0 で試してみると、最終的には glibc の中の __ieee754_exp という関数にたどり着きました (確かめていませんが、Visual C++ では Microsoft 実装の似たようなものがあるはずです)。パッと見ただけでは何をやっているのか理解できないほど複雑ですが、任意精度演算ではないはずです。これをコピペして自分で実行してみても面白いかもしれません。

sourceware.org Git - glibc.git/blob - sysdeps/ieee754/dbl-64/e_exp.c
https://sourceware.org/git/?p=glibc.git ... ads/master

[lbfuvab]

ベキ級数を使う実装の場合、
exp(a) = exp(a/n)^n
の類の等式を使えば収束は速くなりますね。