表1の出力に用いたCプログラム(文字コードはUTF-8)
/* The E-algorithm for Leibniz */ #include <stdio.h> #include <math.h> #define NMAX 10 int main(void){ int n,k,m,j,sign=-1; double e[NMAX+1][NMAX+1]; // e[n][k]=E_k^{(n)}(s) double g[NMAX+1][NMAX+1][NMAX+1]; // g[n][k][j]=g^{(n)}_{k,j} double sum=0.0,tv; double c[NMAX+1][NMAX]; // c[n][k]=g[n][k-1][k]/(g[n+1][k-1][k]-g[n][k-1][k]) tv=4.0*atan(1.0); /* pi */ for(n=1;n<=NMAX;n++){ sign *= -1; sum += sign*4.0 /((double) 2*n-1); e[n][0]=sum; if (n==1) continue; j=1; if (n==2){ g[1][0][j]=g0(1,j); g[2][0][j]=g0(2,j); } else{ g[n][0][j]=g0(n,j); } c[n-1][1]=g[n-1][0][1]/(g[n][0][1]-g[n-1][0][1]); e[n-1][1]=e[n-1][0]-c[n-1][1]*(e[n][0]-e[n-1][0]); for(j=2;j<=n-1;j++){ if (j==n-1){ for(m=1;m<=n;m++){ g[m][0][j]=g0(m,j); } for(k=1;k<=n-2;k++){ for(m=1;m<=n-k;m++){ g[m][k][n-1]=g[m][k-1][n-1]-c[m][k]*(g[m+1][k-1][n-1]-g[m][k-1][n-1] ); } // endfor m } // endfor k c[1][n-1]=g[1][n-2][n-1]/(g[2][n-2][n-1]-g[1][n-2][n-1]); e[1][n-1]=e[1][n-2]-c[1][n-1]*(e[2][n-2]-e[1][n-2]); break; } // endif g[n][0][j]=g0(n,j); for(k=1;k<=j-1;k++){ g[n-k][k][j]=g[n-k][k-1][j]-c[n-k][k]*(g[n-k+1][k-1][j]-g[n-k][k-1][j]); } c[n-j][j]=g[n-j][j-1][j]/(g[n-j+1][j-1][j]-g[n-j][j-1][j]); e[n-j][j]=e[n-j][j-1]-c[n-j][j]*(e[n-j+1][j-1]-e[n-j][j-1]); } // endfor j } // endfor n printf("e_algor.c\n"); printf("\n"); for(n=1;n<=NMAX;n++){ printf("%3d",n); for(k=0;k<=n-1;k++){ printf("%19.15f",e[n-k][k]); if ((k+1)%4==0 && k<n-1){ printf("\n"); printf(" "); } } printf("%19.15f",e[n-k][k]); if ((k+1)%4==0 && k<n-1){ printf("\n"); printf(" "); } } printf("\n"); } for(n=1;n<=NMAX;n++){ for(k=0;k<=n-1;k++){ printf("%10.2e",e[n-k][k]-tv); } printf("\n"); } return 0; } double g0(int n,int j) { return pow(-1,n-1)*pow((double) n,1-2*j); }
お気づきのことがありましたら、下記のメールアドレスにご一報ください。
メールアドレス 長田直樹 <osada "at" cis.twcu.ac.jp> (スパム防止のため@を"at"で置き換えています。)