[zz]破解求pi的怪异程序

来源：互联网发布：淘宝新店采集神器编辑：程序博客网时间：2024/05/17 20:31

破解求pi的怪异程序
Cong Wang
25th November,2005

Institute of Post and Telecommunication, Xi'an, PRC 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 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 wo
rld, 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)