表1の出力に用いたCプログラム(文字コードはUTF-8)
#include <stdio.h>
#include <math.h>
double est(int nu, double x);
double s(int nu, double x);
int main(void)
{
double x=15.0,y,ei = -1.9186278921478669771e-8,sol;
int nu;
sol = -exp(x) * ei;
printf("%3s %18s %10s %10s\n","nu","y ","error","err est");
for(nu = 0;nu <= 19;nu++){
y = s(nu,x);
printf("%3d %-18.15g %10.3g %10.3g\n",nu,y,y-sol,est(nu,x));
}
printf("%3s %-18.15g\n","sol",sol);
return 0;
}
/* sum_{k=0}^{nu} (-1)^k k!/x^{k+1} */
double s(int nu, double x)
{
int k,sign = -1;
long long int fact=1LL;
double y = 0.0,powx = 1.0;
for(k = 0;k <= nu;k++){
powx *= x;
sign *= -1;
if (k == 0)
fact=1LL;
else
fact *= (long long int) k;
y += sign * fact /powx;
}
return y;
}
double est(int nu, double x)
{
int k;
long long int fact=1LL;
double powx = 1.0;
for(k = 0;k <= nu+1;k++){
powx *= x;
if (k==0)
fact=1LL;
else
fact *= (long long int) k;
}
return 2 * fact /powx;
}
gcc version 4.1.2 (i686-pc-linux-gnu) によるコンパイルと実行
$ gcc -Wall -lm euler.c $ ./a.out nu y error err est 0 0.0666666666666667 0.00395 0.00889 1 0.0622222222222222 -0.000498 0.00119 2 0.0628148148148148 9.45e-05 0.000237 3 0.0626962962962963 -2.4e-05 6.32e-05 4 0.0627279012345679 7.62e-06 2.11e-05 5 0.062717366255144 -2.91e-06 8.43e-06 6 0.0627215802469136 1.3e-06 3.93e-06 7 0.0627196137174211 -6.65e-07 2.1e-06 8 0.0627206625331504 3.83e-07 1.26e-06 9 0.0627200332437128 -2.46e-07 8.39e-07 10 0.0627204527700046 1.74e-07 6.15e-07 11 0.0627201451173906 -1.34e-07 4.92e-07 12 0.0627203912394818 1.12e-07 4.27e-07 13 0.0627201779336694 -1.01e-07 3.98e-07 14 0.0627203770190943 9.79e-08 3.98e-07 15 0.0627201779336694 -1.01e-07 4.25e-07 16 0.0627203902914559 1.11e-07 4.81e-07 17 0.0627201496192979 -1.29e-07 5.78e-07 18 0.0627204384258876 1.59e-07 7.32e-07 19 0.0627200726042074 -2.07e-07 9.76e-07 sol 0.0627202791074092メールアドレス 長田直樹 <osada "at" cis.twcu.ac.jp> (スパム防止のため@を"at"で置き換えています。)