一、圓周率π計算 

/*
*@Author:   STZG
*@Language: C++
*/
#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，而不影響精度。