最近兩天在看圓周率的算法,網(wǎng)上流傳這么一個(gè)經(jīng)典的算法,

?1 #include < stdio.h >
?2 long ?a = 10000 ,b = 0 ,c = 2800 ,d,e,f[ 2801 ],g;?
?3 void ?main()?
?4 {??
?5 ????printf( " %d " ,b);
?6 ???? for (;b != c;)?
?7 ???? {?
?8 ????????f[b] = a / 5 ;?
?9 ????????b ++ ;?
10 ????} ?
11 ???? for (;?d = 0 ,g = c * 2 ;?c -= 14 ,printf( " %.4d " ,e + d / a),e = d % a)?
12 ???? for (b = c;d += f[b] * a,f[b] = d %-- g,d /= g -- , -- b;d *= b);??
13 } ?
14

??????? 實(shí)際上,我看f[2080]和e都是沒有初始化的,f[2080]應(yīng)初始化為a/5,e初始化為0.
上機(jī)驗(yàn)證了一下,算出來打印256位,對(duì)比圓周率表(附后),可知答案是正確的,但是打印256位之后,程序卡在那兒不動(dòng),不知道是哪陷入了死循環(huán).
?????? 至于到底有什么問題,還望大蝦幫忙看一看.
?????? 先不管它,這個(gè)程序雖然風(fēng)格惡劣(“短小精悍”讓人看不懂),但這個(gè)算法還是值得研究的.畢竟,它是確定性算法,比各種隨機(jī)版本(比如蒙特卡羅法模擬)的算法要可靠.

一、源程序
本文分析下面這個(gè)很流行的計(jì)算PI的小程序。下面這個(gè)程序初看起來似乎摸不到頭腦，
不過不用擔(dān)心，當(dāng)你讀完本文的時(shí)候就能夠基本讀懂它了。
程序一：很牛的計(jì)算Pi的程序

int ?a = 10000 ,b,c = 2800 ,d,e,f[ 2801 ],g;
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);
}

二、數(shù)學(xué)公式
數(shù)學(xué)家們研究了數(shù)不清的方法來計(jì)算PI,這個(gè)程序所用的公式如下：
??

至于這個(gè)公式為什么能夠計(jì)算出PI，已經(jīng)超出了本文的能力范圍。
下面要做的事情就是要分析清楚程序是如何實(shí)現(xiàn)這個(gè)公式的。
我們先來驗(yàn)證一下這個(gè)公式：


程序二：Pi公式驗(yàn)證程序
#include "stdio.h"
void main()
{
?? float pi=2;
?? int? i;
?? for(i=100;i>=1;i--)
????? pi=pi*(float)i/(2*i+1)+2;
?? printf("%f\n",pi);
?? getchar();
}
上面這個(gè)程序的結(jié)果是3.141593。

三、程序展開
在正式分析程序之前，我們需要對(duì)程序一進(jìn)行一下展開。我們可以看出程序一都是使用for循環(huán)來完成計(jì)算的，這樣做雖然可以使得程序短小，但是卻很難讀懂。根據(jù)for循環(huán)的運(yùn)行順序，我們可以把它展開為如下while循環(huán)的程序：
程序三：for轉(zhuǎn)換為while之后的程序

?1 int ?a = 10000 ,b,c = 2800 ,d,e,f[ 2801 ],g;
?2 main()? {
?3 int ?i;
?4 for (i = 0 ;i < c;i ++ )
?5 ?????f[i] = a / 5 ;
?6 while (c != 0 )
?7 ????? {
?8 ?????????d = 0 ;
?9 ?????????g = c * 2 ;
10 ?????????b = c;
11 ????????? while ( 1 )
12 ???????????? {
13 ????????????????d = d + f[b] * a;
14 ????????????????g -- ;
15 ????????????????f[b] = d % g;
16 ????????????????d = d / g;
17 ????????????????g -- ;
18 ????????????????b -- ;
19 ???????????????? if (b == 0 )? break ;
20 ????????????????d = d * b;
21 ????????????}
22 ?????????c = c - 14 ;
23 ?????????printf( " %.4d " ,e + d / a);
24 ?????????e = d % a;
25 ????}
26 }
27

下面我們就針對(duì)展開后的程序來分析。

四、程序分析
要想計(jì)算出無限精度的PI，我們需要上述的迭代公式運(yùn)行無數(shù)次，并且其中每個(gè)分?jǐn)?shù)也是完全精確的，這在計(jì)算機(jī)中自然是無法實(shí)現(xiàn)的。那么基本實(shí)現(xiàn)思想就是迭代足夠多次，并且每個(gè)分?jǐn)?shù)也足夠精確，這樣就能夠計(jì)算出PI的前n位來。上面這個(gè)程序計(jì)算800位，迭代公式一共迭代2800次。
int a=10000,b,c=2800,d,e,f[2801],g;
這句話中的2800就是迭代次數(shù)。由于float或者double的精度遠(yuǎn)遠(yuǎn)不夠，因此程序中使用整數(shù)類型（實(shí)際是長(zhǎng)整型），分段運(yùn)算（每次計(jì)算4位）。我們可以看到輸出語句 printf("%.4d",e+d/a); 其中%.4就是把計(jì)算出來的4位輸出，我們看到c每次減少14（ c=c-14;），而c的初始大小為2800，因此一共就分了200段運(yùn)算，并且每次輸出4位，所以一共輸出了800位。
要是想每次計(jì)算3位或5位,怎么辦?首先得unsigned long才可以一次計(jì)算5位.
c的大小隨便改,舉個(gè)例子: 比如,
a=100000,c=c-14,迭代2800次,共循環(huán)200次,輸出1000位
a=1000,c=c-7,迭代2800次,共循環(huán)400次,輸出1200位
由于使用整型數(shù)運(yùn)算，因此有必要乘上一個(gè)系數(shù)，在這個(gè)程序中系數(shù)為1000，也就是說，公式如下：

?
這里的2k表示2000，也就是f[2801]數(shù)組初始化以后的數(shù)據(jù),a=10000,a/5=2000,所以下面的程序把f中的每個(gè)元素都賦值為2000：

?1 for (i = 0 ;i < c;i ++ )
?2 ?????f[i] = a / 5 ;你可能會(huì)覺得奇怪，為什么這里要把一個(gè)常數(shù)儲(chǔ)存到數(shù)組中去，請(qǐng)繼續(xù)往下看。
?3 我們先來跟蹤一下程序的運(yùn)行：
?4 while (c != 0 )??假設(shè)這是第一次運(yùn)行，c = 2800 ，為迭代次數(shù)
?5 ????? {
?6 ?????????d = 0 ;
?7 ?????????g = c * 2 ;??這里的g是用來做k / (2k + 1 )中的分母
?8 ?????????b = c;????這里的b是用來做k / (2k + 1 )中的分子
?9 ????????? while ( 1 )
10 ????????? while ( 1 )
11 ???????????? {
12 ????????????????d = d + f[b] * a;?f中的所有的值都為2000，這里在計(jì)算時(shí)又把系數(shù)擴(kuò)大了a = 10000倍。這樣做的目的稍候介紹，你可以看到輸出的時(shí)候是d / a,所以這不影響計(jì)算.
13 ????????????????g -- ;
14 ????????????????f[b] = d % g;?先不管這一行
15 ????????????????d = d / g;???第一次運(yùn)行的g為2 * 2799 + 1 ，你可以看到g做了分母
16 ????????????????g -- ;
17 ????????????????b -- ;
18 ???????????????? if (b == 0 )? break ;
19 ????????????????d = d * b;?這里的b為2799，可以看到d做了分子。
20 ????????????}
21 ?????????c = c - 14 ;
22 ?????????printf( " %.4d " ,e + d / a);
23 ?????????e = d % a;
24 ????}
25

只需要粗略的看看上面的程序，我們就大概知道它的確是使用的那個(gè)迭代公式來計(jì)算Pi
的了，不過不知道到現(xiàn)在為止你是否明白了f數(shù)組的用處。如果沒有明白，請(qǐng)繼續(xù)閱讀。
?d=d/g,這一行的目的是除以2k+1，我們知道之所以程序無法精確計(jì)算的原因就是這個(gè)除法。即使用浮點(diǎn)數(shù)，答案也是不夠精確的，因此直接用來計(jì)算800位的Pi是不可能的。那么不精確的成分在哪里？很明顯：就是那個(gè)余數(shù)d%g。程序用f數(shù)組把這個(gè)誤差儲(chǔ)存起來，再下次計(jì)算的時(shí)候使用。現(xiàn)在你也應(yīng)該知道為什么d=d+f[b]*a;中間需要乘上a了吧。
把分子擴(kuò)大之后，才好把誤差精確的算出來。
d如果不乘10000這個(gè)系數(shù)，則其值為2000，那么運(yùn)行d=d/g；則是2000/(2*2799+1)，這種整數(shù)的除法答案為0，根本無法迭代下去了。
現(xiàn)在我們知道程序就是把余數(shù)儲(chǔ)存起來，作為下次迭代的時(shí)候的參數(shù)，那么為什么這么做就可以使得下次迭代出來的結(jié)果為接下來的4位數(shù)呢？
這實(shí)際上和我們?cè)诩埳献鞒ê茴愃疲?br />?


我們可以發(fā)現(xiàn)，在做除法的時(shí)候，我們通常把余數(shù)擴(kuò)大之后再來計(jì)算，f中既然儲(chǔ)存的是余數(shù)，而f[b]*a;則正好把這個(gè)余數(shù)擴(kuò)大了a倍，然后如此循環(huán)下去，可以計(jì)算到任意精度。
這里要說明的是，事實(shí)上每次計(jì)算出來的d并不一定只有4位數(shù)，例如第一次計(jì)算的時(shí)候，d的值為31415926，輸出4位時(shí)候，把低四位的值儲(chǔ)存在e中間，e=d%a，也就是5926。最后，這個(gè)c=c-14不太好理解。事實(shí)上沒有這條語句，程序計(jì)算出來的仍然正確。只是因?yàn)槿绻?800次，無論分?jǐn)?shù)如何精確，最后Pi的精度只能夠達(dá)到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;
????}

最后的答案仍然正確。
不過我們可以看到內(nèi)循環(huán)的次數(shù)是c次，也就是說每次迭代計(jì)算c次。而每次計(jì)算后續(xù)位
數(shù)的時(shí)候，迭代次數(shù)減少14，而不影響精度。
附:圓周率小數(shù)點(diǎn)后1000位：

1415926535 8979323846 2643383279 5028841971 6939937510
5820974944 5923078164 0628620899 8628034825 3421170679
8214808651 3282306647 0938446095 5058223172 5359408128
4811174502 8410270193 8521105559 6446229489 5493038196
4428810975 6659334461 2847564823 3786783165 2712019091
4564856692 3460348610 4543266482 1339360726 0249141273
7245870066 0631558817 4881520920 9628292540 9171536436
7892590360 0113305305 4882046652 1384146951 9415116094
3305727036 5759591953 0921861173 8193261179 3105118548
0744623799 6274956735 1885752724 8912279381 8301194912
9833673362 4406566430 8602139494 6395224737 1907021798
6094370277 0539217176 2931767523 8467481846 7669405132
0005681271 4526356082 7785771342 7577896091 7363717872
1468440901 2249534301 4654958537 1050792279 6892589235
4201995611 2129021960 8640344181 5981362977 4771309960
5187072113 4999999837 2978049951 0597317328 1609631859
5024459455 3469083026 4252230825 3344685035 2619311881
7101000313 7838752886 5875332083 8142061717 7669147303
5982534904 2875546873 1159562863 8823537875 9375195778
1857780532 1712268066 1300192787 6611195909 2164201989

想要1000萬位的發(fā)email給我就是.