解读传说中计算π的超牛的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)