圆周率π的精确计算

一、圆周率π计算 

#include <bits/stdc++.h>
using namespace std;
long a=10000,b,c=56000,d,e,f[56001],g;
int main(){
    for(;b-c ; )f[b++]=a/5;
    for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
    for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
    //cout << "Hello world!" << endl;
    return 0;
}

二、数学公式

π = 2 + 1/3 * (2 + 2/5 * (2 + 3/7 * (2 + ...  (2 + k/2k+1 * (2 + ... ))...)))

三、分析

#include <stdio.h>
int a=10000,b,c=2800,d,e,f[2801],g;
main() {
int i;
for(i=0;i<c;i++)
     f[i]=a/5;
while(c!=0)
     {
         d=0;
         g=c*2;
         b=c;
         while(1)
            {
                d=d+f[b]*a;
                g--;
                f[b]=d%g;
                d=d/g;
                g--;
                b--;
                if(b==0) break;
                d=d*b;
            }
         c=c-14;
         printf("%.4d",e+d/a);
         e=d%a;
    }
}

要想计算出无限精度的PI，我们需要上述的迭代公式运行无数次，并且其中每个分数也是完全精确的，这在计算机中自然是无法实现的。那么基本实现思想就是迭代足够多次，并且每个分数也足够精确，这样就能够计算出PI的前n位来。上面这个程序计算800位，迭代公式一共迭代2800次。

int a=10000,b,c=2800,d,e,f[2801],g;

这句话中的2800就是迭代次数。

由于float或者double的精度远远不够，因此程序中使用整数类型（实际是长整型），分段运算（每次计算4位）。我们可以看到输出语句printf("%.4d",e+d/a); 其中%.4就是把计算出来的4位输出，我们看到c每次减少14（ c=c-14;），而c的初始大小为2800，因此一共就分了200段运算，并且每次输出4位，所以一共输出了800位。

由于使用整型数运算，因此有必要乘上一个系数，在这个程序中系数为1000，也就是说，公式如下：

1000*pi = 2K+ 1/3 * (2K+ 2/5 * (2K+ 3/7 * (2K+ ... (2K+ k/2k+1 * (2K+ ... ))...)))

这里的2K表示2000，也就是f[2801]数组初始化以后的数据,a=10000,a/5=2000,所以下面的程序把f中的每个元素都赋值为2000：
for(i=0;i<c;i++)
f[i]=a/5;

你可能会觉得奇怪，为什么这里要把一个常数储存到数组中去，请继续往下看。我们先来跟踪一下程序的运行：

while(c!=0)  //假设这是第一次运行，c=2800，为迭代次数
     {
         d=0;
         g=c*2;  //这里的g是用来做k/(2k+1)中的分子
         b=c;  //这里的b是用来做k/(2k+1)中的分子
         while(1)
         while(1)
            {
                d=d+f[b]*a;  //f中的所有的值都为2000，这里在计算时又把系数扩大了a=10000倍。这样做的目的稍候介绍，你可以看到输出的时候是d/a,所以这不影响计算
                g--;
                f[b]=d%g;  //先不管这一行
                d=d/g;  //第一次运行的g为2*2799+1，你可以看到g做了分母
                g--;
                b--;
                if(b==0) break;
                d=d*b;  //这里的b为2799，可以看到b做了分子。
            }
         c=c-14;
         printf("%.4d",e+d/a);
         e=d%a;
    }
 

只需要粗略的看看上面的程序，我们就大概知道它的确是使用的那个迭代公式来计算Pi的了，不过不知道到现在为止你是否明白了f数组的用处。如果没有明白，请继续阅读。
d=d/g,这一行的目的是除以2k+1，我们知道之所以程序无法精确计算的原因就是这个除法。即使用浮点数，答案也是不够精确的，因此直接用来计算800位的Pi是不可能的。那么不精确的成分在哪里？很明显：就是那个余数d%g。程序用f数组把这个误差储存起来，在下次计算的时候使用。现在你也应该知道为什么d=d+f[b]*a;中间需要乘上a了吧。把分子扩大之后，才好把误差精确的算出来。
d如果不乘10000这个系数，则其值为2000，那么运行d=d/g；则是2000/(2*2799+1)，这种整数的除法答案为0，根本无法迭代下去了。
现在我们知道程序就是把余数储存起来，作为下次迭代的时候的参数，那么为什么这么做就可以使得下次迭代出来的结果为接下来的数字呢？
这实际上和我们在纸上作除法很类似：

　0142
　　/———
　7/　1
　　　10
　　　 7
——————
　　　 30
　　　 28
——————
　　　　20
　　　　14
——————
　　　　 6
.....

我们可以发现，在做除法的时候，我们通常把余数扩大之后再来计算，f中既然储存的是余数，而f[b]*a;则正好把这个余数扩大了a倍，然后如此循环下去，可以计算到任意精度。
这里要说明的是，事实上每次计算出来的d并不一定只有4位数，例如第一次计算的时候，d的值为31415926，输出4位时候，把低四位的值储存在e中间，e=d%a，也就是5926。

最后，这个c=c-14不太好理解。事实上没有这条语句，程序计算出来的仍然正确。只是因为如果迭代2800次，无论分数如何精确，最后Pi的精度只能够达到800。
你可以把程序改为如下形式尝试一下：

for(i=0;i<800;i++)
     {
     {
         d=0;
         g=c*2;
         b=c;
         while(1)
            {
                d=d+f[b]*a;
                g--;
                f[b]=d%g;
                d=d/g;
                g--;
                b--;
                if(b==0) break;
                d=d*b;
            }
       //c=c-14;  //不要这句话。
         printf("%.4d",e+d/a);
         e=d%a;
    }

最后的答案仍然正确。
不过我们可以看到内循环的次数是c次，也就是说每次迭代计算c次。而每次计算后续位数的时候，迭代次数减少14，而不影响精度。