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

破解求 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=d000;     }   } 少了两个变量了 ! 深入分析   好了，马上进入实质性的分析了。一定的数学知识是非常有必要的。首先 , 必须知道下面的公式可以用来求 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=d000; 　， 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));     e=d;   }   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)