円周率の公式は、ウィキペディアの記事「円周率の近似」で確認できます。コンパクトで効率的な計算が約束され、10 進法に特化しているため、この式に惹かれました。式は
pi = -3 + SUM(n=0,oo): n*(2^n)*(n!)^2/(2*n)!
私のCコードは以下です。それは非常に簡単に思えます。すべての中間ステップを計算しますが、完全に収束しません。バグは何ですか？
#include <stdio.h>
#include <math.h>
#define sq(x) ((x)*(x))
int nfac(int);
int main()
{
double term, denom, sum;
double w, x, y, z, pi;
int n;
/* Plouffe's 1996 algorithm
(see http://en.wikipedia.org/wiki/Approximations_of_π: */
sum = -3.;
for(n=1; n <= 11; n++)
{
printf("n= %d\n", n);
printf("n! = %.0f\n", w = nfac(n));
printf("(n!)^2 = %.0f\n", x = sq(w));
printf("2^n = %.0f\n", y = pow(2,n));
printf("(2*n)! = %.0f\n", z = nfac(2*n));
printf("n*2^n)*((n!)^2)/(2*n)! = %f\n", term = n*y*x/z);
printf("sum = %f\n\n", sum += term);
}
printf("pi = %.10f\n\n", pi = sum);
}
int nfac(int n)
{
int i, nn;
if(n==0) return 1;
nn = 1;
for(i=1; i<=n; i++)
nn= i*nn;
return nn;
}