计算 π

日期: 05 月 20日, 2014
标签:

计算公式

取近似值

即

该公式可以方便地用迭代法实现，代码如下

#include <stdio.h>

int main()
{
	int i;
	double sum = 2.0;

	for (i = 100; i >= 1; i--)
	{
		sum = sum * i / (i * 2 + 1);
		sum += 2.0;
	}
	printf("%1.15lf\n", sum);

	return 0;
}

程序运行结果为

gcc -o pi pi.c
./pi
3.141592653589793

高精度实现

高精度实现的代码如下

#include <stdio.h>

#define ARRAY_SIZE 2502
#define ITERATIONS 33223

int pi[ARRAY_SIZE+1];

int main()
{
    int i, j;

    // pi = 2.0;
    pi[0] = 2000;
    for(i = ITERATIONS; i > 0 ; i--)
    {
        // pi *= i;
        for(j = ARRAY_SIZE; j >= 0; j--)
        {
            pi[j] *= i;
        }
        for(j = ARRAY_SIZE-1; j > 0; j--)
        {
            pi[j-1] += pi[j] / 10000;
            pi[j] %= 10000;
        }

        // pi /= (i * 2 + 1);
        for(j = 0; j < ARRAY_SIZE; j++)
        {
            pi[j + 1] += pi[j] % (i * 2 + 1) * 10000;
            pi[j] /= (i * 2 + 1);
        }
        // pi += 2.0;
        pi[0] += 2000;
    }
    for(i = 0; i < 2500; i++)
    {
        printf("%04d", pi[i]);
    }
    return 0;
}

误差分析

截断误差

易得

所以截断误差

舍入误差

设数组pi的大小为k+1，则每次运算的舍入误差不超过$\frac{1}{10^{4k}}$，则多次运算积累的舍入误差

参数选取

若要计算$\pi$的前N位有效数字，须满足

又

故取n和k满足下式即可

不妨取

取