解读传说中计算解读传说中计算π的超牛的的超牛的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;
评论0