C言語でf(x) = e^x-1の高精度な数値を求めたいです

Question

C言語でf(x) = e^x-1の高精度な数値を求めたいですのですが、10^-15の位からどうしても値が合いません。
以下のfor文の書き方で積み残しは対処できたと思っています。
自分は桁落ちが対処できていないのかなと思います。四倍精度小数点型を使用すれば対処できるのかなとは思うのですが、使用方法がわかりません。
以下に実行結果を掲載させていただきます
【実行結果（項数100の時）】
項数 = 100
f(x) = 　0.0000100000500001666682125655433166500075
解析値 =0.0000100000500000696490587870357558131218
誤差 = 　0.0000000000000000970191537785075608368857

#include<stdio.h>
#include<math.h>

double kaizyou(int);

int main(void){
  int N;
  int n; // 項数
  int k; // 項数カウント
double x; // e^xのx
double xk; // x^n/n!
double e; // テイラー展開の部分和  
double fx;//解析値
double er; // exp(x)との差

  printf("N = ");
  scanf("%d", &N); 
  x = pow(10, -5);
  n = pow(10, N);

  printf("項数 = %d\n", n);

  e = 0.0;
  xk = 1.0;

  for(k = n ; k >= 1 ; --k){
    e += kaizyou(k);
  }

  fx = exp(x)-1;
  er = fabs(e-fx);
  printf("f(x) =\t%.40f\n", e);
  printf("解析値 =%.40f\n", fx);
  printf("誤差 =\t%.40f\n", er);

  return 0;
}

double kaizyou(int n){
  int i;
  long double x = 1, d = pow(10, -5);

  for(i = n ; i  >= 1 ; --i){
    x *= d/(double)i;
  }

  return x;
}

Answer 1

他の方もおっしゃっている通り、本来 15 桁ほどあるはずの double の有効桁数が exp(x) - 1 の部分で 10 桁ほどに減ってしまうのが原因です。
実は exp(x) - 1 を (桁落ちさせずに精度よく) 求める専用の関数があります: https://en.cppreference.com/w/c/numeric/math/expm1
fx = exp(x)-1; の行を fx = expm1(x); に書き換えれば e == fx になることが確かめられます。

とはいえ、有効数字 15 桁 (小数点以下 20 桁) を超えて完全に一致しているのは偶然でしかないので、もし有効数字を増やしたいのであれば

• コンパイラの独自拡張の 4 倍精度を使う
• 多倍長実数のライブラリ (mpfr またはそれの boost ラッパ) を使う

のどちらかですね。

Answer 2

どれだけの精度が必要なのかわかりませんが、多倍長浮動小数点数ライブラリを使う手もあります。

ちょっと調べてみたところ、Boostにもあるみたいですね。URLを貼っておきます。

多倍長浮動小数点数

Answer 3

まず前提としてdoubleの精度は十進数表現で約15桁というのはご認識されていると思います。
exp(1e-5)は約15桁の精度を持っていますが、それは先頭の整数部の1から15桁です。
ということで、

解析値 =0.0000100000500000696490587870357558131218

とお書きの数値のうち、有効なのは、0.00001000005000だけです。
ほぼ同じ値同士の引き算で、桁落ち（正しい意味）が発生し、有効精度が減ってしまいました。

任意精度計算ツールのbcを使うと、

$ echo 'scale=50;e(0.00001)-1' | bc -l
.00001000005000016666708333416666805555753968501984

なので、お書きのプログラムで計算された値は有効桁15桁で正しいようです。

追記

gccだとデフォルトではlong doubleは80bit拡張倍精度（19桁）のようなので（16バイトなのに）、オプションで4倍精度（34桁）にします。
printfは4倍精度に対応していないようです。また、浮動小数点関数は倍精度なので使えません。

gcc -mlong-double-128 -lm -lquadmath -o foo foo.c

#include <stdio.h>
#include <math.h>
#include <quadmath.h>

typedef long double FLOAT;
FLOAT kaizyou(int);

int main(void){
  int N;
  int n; // 項数
  int k; // 項数カウント
  FLOAT e; // テイラー展開の部分和
  FLOAT er; // exp(x)との差
  FLOAT fx = .00001000005000016666708333416666805555753968Q;
  char buf[100];

  printf("N = ");
  scanf("%d", &N);
  n = pow(10, N);
  printf("項数 = %d\n", n);
  e = 0.0;
  for(k = n ; k >= 1 ; --k){
    e += kaizyou(k);
  }
  er = e-fx;
  if(er<0) er = -er;
  quadmath_snprintf (buf, sizeof buf, "%.40Qf", e);
  printf("f(x) =\t%s\n", buf);
  quadmath_snprintf (buf, sizeof buf, "%.40Qf", fx);
  printf("解析値=\t%s\n", buf);
  quadmath_snprintf (buf, sizeof buf, "%.40Qf", er);
  printf("誤差=\t%s\n", buf);
  return 0;
}

FLOAT kaizyou(int n){
  int i;
  FLOAT x = 1, d = 0.00001Q;

  for(i = n ; i  >= 1 ; --i){
    x *= d/(FLOAT)i;
  }
  return x;
}

Comments

[otn]

scale=50と書いてますが、e()を計算する際の内部の各計算が50桁で行われると言うことだと思いますので、e()の結果が50桁の精度を持っているわけではないでしょうね。

[raiboz1115]

回答ありがとうございます。とても分かりやすくて助かりました！
有効桁を15桁以上にするにはどうすればよいでしょうか？

[otn]

doubleをlong doubleにするのでしょうか。
4倍精度やったことなかったので、やってみました。