解读传说中计算π的超牛的C程序

来源:本站
导读:目前正在解读《解读传说中计算π的超牛的C程序》的相关信息,《解读传说中计算π的超牛的C程序》是由用户自行发布的知识型内容!下面请观看由(电工技术网 - www.9ddd.net)用户发布《解读传说中计算π的超牛的C程序》的详细说明。
简介:在我上大学的时候就流传着这样一个超牛的C程序,只用三行代码就能计算π到小数点后800位,还有的地方开玩笑说是外星人写的,的确是牛的不得了。那个时候大家一起研究都搞不懂,昨天看了一篇文章解释这段代码,今天自己试验了很久,终于弄明白了,所以记下来和大家一起交流。

这段C代码是这样的:

#include"stdio.h"

longa=10000,b,c=2800,d,e,f[2801],g;

voidmain(){

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倍。

可以由如下的一段代码实现:

#defineM2800

#include<stdio.h>

longb,d,f[M+1],i;

voidmain(){

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;

}

如果写成这样就再明白不过了:

#defineM2800

#include<stdio.h>

voidmain()

{

longyushu[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,就可以得到

#defineM2800

#include<stdio.h>

longb,d,e,f[M+1],i;

voidmain(){

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的怪异代码

CongWang

25thNovember,2005

InstituteofPostandTelecommunications,Xi'an,P.R.China

NetworkEngineeringDep.

引言

网上流传着一个怪异的求pi程序,虽然只有三行却能求出pi值连小数点前共800位。这个程序如下:

/*某年ObfuscatedCContest佳作选录:*/

#include<stdio.h>

longa=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.mathFAQ,原作者未详)*/

咋一看,这程序还挺吓人的。别慌,下面就告诉你它是如何做到的,并且告诉你写怪异C程序的一些技巧。^_^

展开化简

我们知道,在C语言中,for循环和while循环可以互相代替。

for(statement1;statement2;statement3){

statements;

}

上面的for语句可以用下面的while语句来代替:

statement1;

while(statement2){

statements;

statement3;

}

而且要写怪异的C程序,逗号运算符无疑是一个好的助手,它的作用是:

从左到右依次计算各个表达式的值,并且返回最右边表达式的值。

把它嵌入for循环中是写怪异代码的常用技巧之一。所以,上面的程序可以展开为:

#include<stdio.h>/*1*/

/*2*/

longa=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>

longb,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位的程序,如下:

longb,c=2800,d,e,f[2801],g;

intmain(intargc,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;

}

return0;

}

不妨试试把上面的程序压缩成3行。

结论

以Knuth图灵演讲中的一句话结束全文:

Wehaveseenthatcomputerprogrammingisanart,becauseitappliesaccumulatedknowledgetotheworld,becauseitrequiresskillandingenuity,andespeciallybecauseitproducesobjectsofbeauty.Aprogrammerwhosubconsciouslyviewshimselfasanartistwillenjoywhathedoesandwilldoitbetter.

与大家共勉!^_^

提醒:《解读传说中计算π的超牛的C程序》最后刷新时间 2024-03-14 00:59:47,本站为公益型个人网站,仅供个人学习和记录信息,不进行任何商业性质的盈利。如果内容、图片资源失效或内容涉及侵权,请反馈至,我们会及时处理。本站只保证内容的可读性,无法保证真实性,《解读传说中计算π的超牛的C程序》该内容的真实性请自行鉴别。