解讀傳說中計算π的超牛的C程式

在我上大學的時候就流傳著這樣一個超牛的C程式，只用三行代碼就能計算π到小數點後800位，還有的地方開玩笑說是外星人寫的，的確是牛的不得了。那個時候大家一起研究都搞不懂，昨天看了一篇文章解釋這段代碼，今天自己實驗了很久，終於弄明白了，所以記下來和大家一起交流。
這段C代碼是這樣的：
#include "stdio.h"
long a=10000, b, c=2800, d, e, f[2801], g;
void 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) ;
}
我把解釋這段代碼的文章附在最後面了，作者叫王聰，看他的部落格是很高手的樣子。我的主要分析都是從他的文章裡看到的。謝謝這位高手！這麼多年的疑惑終於解開了。
從我的角度來理解，我不習慣把for語句拆開變成while來分析，我覺著現在這個樣子就挺好的，當然個別語句還是要調整一下。我覺得說明以下幾個問題就能完全搞明白了：演算法（兩點）、錯誤（兩點）、其他（四點）。
一、演算法
1、π的計算公式
π/2=1+1!/3!!+2!/5!!+3!/7!!+...+k!/(2*k+1)!!+...
這個公式不記得了，好像是高數裡學過。註：5!=1*2*3*4*5，5!!=1*3*5
2、公式的程式實現（數組儲存餘數）
把上面的公式做一下展開和調整
π/2 = 1 + 1!/3!! + 2!/5!! + 3!/7!! + ... + k!/(2*k+1)!!
= 1 + 1/3 + (1*2)/(3*5) + (1*2*3)/(3*5*7) + ... + (1*2*...*k)/(3*5*...*(2k+1))
= 1 + 1/3 + (1/3)*(2/5) + (1/3)*(2/5)*(3/7) + ... + (1/3)*(2/5)*...*(k/(2k+1))
= (1/3)*(2/5)*...*(k/(2k+1)) + ... + (1/3)*(2/5)*(3/7) + (1/3)*(2/5) + 1/3 + 1
= (k/(2k+1))*...*(2/5)*(1/3) + ... + (3/7)*(2/5)*(1/3) + (2/5)*(1/3) + 1/3 + 1
= (((((k/(2k+1)+1)*((k-1)/(2(k-1)+1)+1)*...)*3/7+1)*2/5+1)*1/3+1)/1
= (((((1/(2k+1)*k+1)/(2(k-1)+1)*(k-1)+1)/...)/7*3+1)/5*2+1)/3*1+1)/1
可以了，我們要做的就是做除法、做乘法、加1，做除法、做乘法、加1，...，這樣直到做除法除以1。那麼我們就可以在一個迴圈中完成它，用一個計數器i，1除以2i+1、乘以i、加1，迴圈。但是我們要輸出800位小數，就只能用大數除法的辦法了：分為n趟執行，這一趟把公式中每一次除法的餘數儲存下來，把商累計後輸出；下一趟把儲存的餘數提高精度再做除法，把餘數儲存下來，把商累計後輸出，以此類推。根據數學運算，當k=2800時，可以精確到小數點後800位（實際上是精確到小數點後844位），那麼我們就設定一個2800的數組，全部初始化為1（也就是設定餘數為1，就是公式中加的1），然後在迴圈中不斷的做除法、乘法、加餘數，過程中還要儲存餘數，和手工除法一樣。
要說明的一點就是提高精度的辦法，比如f[2800]，初始值為1，除以5601、乘以2800，把商的整數部分加1放到下一項的運算中，餘數儲存下來，在下一趟乘以10，再除以5601、乘以2800，得到的商的整數部分就是上一趟運算結果的後一位元字。我們如果要每次輸出4位元字，那麼就要乘以10000，得到的商的整數部分不會超過4位元，就是上一趟運算結果的後四位元字。注意到上述公式的結果是π/2，那麼我們就把數組都初始化為2，由於要一次輸出4位，所以擴大10000倍。
可以由如下的一段代碼實現：
#define M 2800
#include <stdio.h>
long b, d, f[M+1], i;
void main() {
for( ; b<=M; ) f[b++]=1000*2;
for( ; i<800/4; i++, printf("%.4d",f[0]+d/10000), f[0]=d%10000)
for(d=0, b=M; d+=f[b]*10000, f[b]=d%(2*b+1), d/=(2*b+1), b; b--)
d*=b;
}
如果寫成這樣就再明白不過了：
#define M 2800
#include <stdio.h>
void main()
{
long yushu[M+1], k = 0, i = 0, j = 0, sum = 0, PI = 0;
for (k = 0; k <= M; k++)
yushu[k] = 1000 * 2;
for (i = 1; i <= 800 / 4; i++)
{
for (j = M, sum = 0; j > 0; j--)
{
sum += yushu[j] * 10000;
yushu[j] = sum % (2 * j + 1);
sum /= (2 * j + 1);
sum *= j;
}
sum += yushu[j] * 10000;
yushu[j] = sum % (2 * j + 1);
printf("%.4d", PI + sum / 10000);
PI = sum % 10000;
}
}
把第一段代碼稍加變形，就是把b替換為b-1，就可以得到
#define M 2800
#include <stdio.h>
long b, d, e, f[M+1], i;
void main() {
for( ; b-M; ) f[++b]=1000*2;
for( ; i<800/4; i++, printf("%.4d",e+d/10000), e=d%10000)
for(d=0, b=M; d+=f[b]*10000, f[b]=d%(2*b-1), d/=(2*b-1), b>1; b--)
d*=(b-1);
}
再對照著繼續變形就能得到那段傳說的原始代碼了。這個過程中用了一些障眼法做混淆作用，尤其是故意設下兩個錯誤，但卻沒有影響計算結果的精度，下面就要提到。
二、錯誤
1、數組初始化
for( ;b-c; ) f[b++] =a/5;
這一行代碼把數組f[0]至f[2799]設定為2000，f[2800]沒有初始化，還是0，那麼在實際運算中這一項是沒有起到作用的，由於這一項的結果實在太小了。但是從計算邏輯上，卻有很大的混淆作用，實際上f[0]是從來沒有用到的，所以初始化應該是把f[1]至f[2800]初始化為2000，這樣精度在小數點後848位才出現變化，由於結果是精確到第844位，所以這點精度實際上是沒有意義的，但是演算法卻正確了。也就是
for( ;b-c; ) f[++b] =a/5;
2、除法的運算過程（除法的計算起點）
如果細心可以發現在第三個for迴圈的初始化中有一點不同，原始代碼用的是b=c，我們用的是b=M，也就是b=2800，不同點在於我們按照公式把每一次除法都做了，而原始代碼是每次漏掉一項，由於每趟計算的結果儲存在餘數數組裡，所以這種遺漏對結果的影響非常小以致於可以忽略不計。但這樣的混淆作用太大了，原來的演算法幾乎面目全非。所以這一點錯誤對於理解代碼至關重要。
三、其他
1、輸出結果（printf）
代碼中的d就是每趟運算中商的累計，就是要輸出的結果，由公式及我們的演算法處理可知，d不會超過8位，而且其左4位與前次餘數e的和也不會超過8位，所以我們用printf("%.4d",e+d/10000)把左4位輸出，把d的右4位儲存在e中e=d%10000。printf("%.4d",…)輸出4位整數，且顯示最左面的0，就是那個.的作用了。
2、兩個乘數（1000和10000）
在公式的程式實現中已經解釋過這兩個數位來曆了。
3、c的作用
在第二個for迴圈中用c-=14來控制迴圈次數，實際上的影響在於第三個for迴圈的初始化部分b=c，上面分析過這樣的處理（錯誤）對精度幾乎沒有影響。實際上c-=14沒有特殊的含義，只是因為2800/(800/4)=14，我們完全可以用計數器來替代，就像我們的代碼中一樣。
4、g的作用
g就是公式中的2k+1，在我們的程式中完全可以用2*b+1來替代。
在我們的程式中，2800就是公式中的k，通過它來調整計算結果的精度，如果為了偵錯工具當然多少都可以了，我調試的時候用的是6。還有用到計數器i，i<800/4就是顯示800位的意思，實際上我們顯示的是小數點後799位，加上小數點前面的3總共是800位。
下面的文章中對f[2800]的分析實際上是不對的，對c-=14的說明如上“c的作用”，第二個for迴圈中b=c的錯誤作者沒有提到，而且公式的展開還不到位，和程式中的演算法實現對不上。
解讀求pi的怪異代碼
Cong Wang
25th November,2005
Institute of Post and Telecommunications,Xi'an, P.R.China
Network Engineering Dep.
引言
網上流傳著一個怪異的求pi程式，雖然只有三行卻能求出pi值連小數點前共800位。這個程式如下：
/*某年Obfuscated C Contest佳作選錄:*/
#include <stdio.h>
long 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);
}
/* (本程式可算出pi值連小數點前共800位)
(本程式錄自sci.math FAQ,原作者未詳)*/
咋一看，這程式還挺嚇人的。別慌，下面就告訴你它是如何做到的，並且告訴你寫怪異C程式的一些技巧。^_^
展開化簡
我們知道，在C語言中，for迴圈和while迴圈可以互相代替。
for(statement1;statement2;statement3){
statements;
}
上面的for語句可以用下面的while語句來代替：
statement1;
while(statement2){
statements;
statement3;
}
而且要寫怪異的C程式，逗號運算子無疑是一個好的助手，它的作用是：
從左至右依次計算各個運算式的值，並且返回最右邊運算式的值。
把它嵌入for迴圈中是寫怪異代碼的常用技巧之一。所以，上面的程式可以展開為：
#include <stdio.h> /*1*/
/*2*/
long a=10000, b, c=2800, d, e, f[2801], g; /*3*/
main(){ /*4*/
while(b-c!=0){ /*5*/
f[b]=a/5; /*6*/
b++; /*7*/
} /*8*/
d=0; /*9*/
g=c*2; /*10*/
while(g!=0){ /*11*/
b=c; /*12*/
d+=f[b]*a; /*13*/
f[b]=d%--g; /*14*/
d=d/g--; /*15*/
--b; /*16*/
while(b!=0){ /*17*/
d=d*b+f[b]*a; /*18*/
f[b]=d%--g; /*19*/
d=d/g--; /*20*/
--b; /*21*/
} /*22*/
c-=14; /*23*/
printf("%.4d",e+d/a); /*24*/
e=d%a; /*25*/
d=0; /*26*/
g=c*2; /*27*/
} /*28*/
} /*29*/
現在是不是好看一點了？
進一步化簡
你應該能注意到a的值始終是10000，所以我們可以把a都換成10000。再就是，仔細觀察g，在外層迴圈中，每次迴圈用它做除法或取餘時，它總是等於2*c-1，而b總是初始化為c。在內層迴圈中，b每次減少1，g每次減少2。你這時會想到了吧？用2*b-1代替g！代進去試試，沒錯！另外，我們還能做一點化簡，第26行的d=0是多餘的，我們把它合并到第13行中去，第13行可改寫為：d=f[b]*a;　。所以程式可以改為：
#include <stdio.h>
long b, c=2800, d, e, f[2801];
main(){
while(b-c!=0){
f[b]=2000;
b++;
}
while(c!=0){
b=c;
d=f[b]*10000;
f[b]=d%(b*2-1);
d=d/(b*2-1);
--b;
while(b!=0){
d=d*b+f[b]*10000;
f[b]=d%(b*2-1);
d=d/(b*2-1);
--b;
}
c-=14;
printf("%.4d",e+d/10000);
e=d%10000;
}
}
少了兩個變數了!
深入分析
好了，馬上進入實質性的分析了。一定的數學知識是非常有必要的。首先,必須知道下面的公式可以用來求pi:
pi/2=1+1!/3!!+2!/5!!+3!/7!!+...+k!/(2*k+1)!!+...
只要項數足夠多，pi就有足夠的精度。至於為什麼，我們留給數學家們來解決。
寫過高精度除法程式的人都知道如何用整數數組來進行除法用法，而避免小數。其實很簡單，回想一下你是如何徒手做除法的。用除數除以被除數，把得數放入結果中，餘數乘以10後繼續做下一步除法，直到餘數是零或達到了要求的位元。
原程式使用的數學知識就那麼多,之所以複雜難懂是因為它把演算法非常巧妙地放到迴圈中去了。我們開始具體來剖析器。首先，我們從數學公式開始下手。我們求的是pi，而公式給出的是pi/2。所以，我們把公式兩邊同時乘以2得:
pi=2*1+2*1!/3!!+2*2!/5!!+2*3!/7!!+...+2*k!/(2*k+1)!!+...
接著，我們把它改寫成另一種形式,並展開：
pi=2*1+2*1!/3!!+2*2!/5!!+2*3!/7!!+...+2*n!/(2*n+1)!!
=2*(n-1)/(2*n-1)*(n-2)/(2*n-3)*(n-3)/(2*n-5)*...*3/7*2/5*1/3
+2*(n-2)/(2*n-3)*(n-3)/(2*n-5)*...*3/7*2/5*1/3
+2*(n-3)/(2*n-5)*...*3/7*2/5*1/3
+2*3/7*2/5*1/3
+2*2/5*1/3
+2*1/3
+2*1
對著公式看看程式，可以看出，b對應公式中的n，2*b-1對應2*n-1。b是從2800開始的，也就是說n=2800。（至於為什麼n=2800時,能保證pi的前800位準確不在此討論範圍。）看程式中的相應部分：
d=d*b+f[b]*10000;
f[b]=d%(b*2-1);
d=d/(b*2-1);
d用來存放除法結果的整數部分，它是累加的，所以最後的d將是我們要的整數部分。而f[b]用來存放計算到b為止的餘數部分。
到這裡你可能還不明白。一是，為什麼數組有2801個元素？很簡單，因為作者想利用f[1]~f[2800]，而C語言的數組下標又是從0開始的，f[0]是用不到的。二是，為什麼要把數組元素除了f[2800]都初始化為2000？10000有什麼作用？這倒很有意思。因為從printf("%.4d",e+d/10000);看出d/10000是取d的第4位以前的數字，而e=d%10000;　，e是d的後4位元字。而且，e和d差著一次迴圈。所以列印的結果恰好就是我們想要的pi的相應的某4位！開始時之所以把數組元素初始化為2000，是因為把pi放大1000倍才能保證整數部分有4位，而那個2就是我們公式中兩邊乘的2！所以是2000！注意，餘數也要相應乘以10000而不是10！f[2800]之所以要為0是因為第一次乘的是n-1也就是2799而不是2800!每計算出4位，c都要相應減去14，也就保證了一共能夠列印出4*2800/14=800位。但是，這竟然不會影響結果的準確度！本人數學功底不高，無法給出確切答案。（要是哪位達人知道，請發email到xiyou.wangcong@gmail.com告訴我哦。）
偶然在網上見到一個根據上面的程式改寫的“準確”（這個準確是指沒有漏掉f[]數組中的的任何一個元素。）列印2800位的程式，如下：
long b,c=2800,d,e,f[2801],g;
int main(int argc,char* argv[])
{
for(b=0;b<c;b++)
f[b] = 2;
e=0;
while(c > 0)
{
d=0;
for(b=c;b>0;b--)
{
d*=b;
d+=f[b]*10;
f[b]=d%(b*2-1);
d/=(b*2-1);
}
c-=1;
printf("%d",(e+d/10)%10);
e=d%10;
}
return 0;
}
不妨試試把上面的程式壓縮成3行。
結論
以Knuth圖靈演講中的一句話結束全文：
We have seen that computer programming is an art, because it applies accumulated knowledge to the world, because it requires skill and ingenuity, and especially because it produces objects of beauty. A programmer who subconsciously views himself as an artist will enjoy what he does and will do it better.
與大家共勉！^_^
(EOF)