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