表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"で置き換えています。)