在之前有人给我发了这样一段代码,让我解释一下其中的原理:
#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);
}
乍一看根本不知所云,于是自己运行了一下,发现这个主体只有三行的程序居然能算出圆周率的前800位!这也激起了我的兴趣。
当然,懒惰的天性还是促使着我先去网上搜索了一下相关的东西,于是在贴吧里面发现了这样一篇帖子:
http://tieba.baidu.com/p/173404246
网上还有许多与这个程序相关的博文,不过仔细读完之后发现很多都是由这篇帖子改编的,其他的也都大同小异。这些帖子和博文中都选择性地回避了几个重要的问题,这几个问题也将是本文重点讨论的部分。
#include <stdio.h>
long a=10000,b,c=2800,d,e,f[2801],g;
main()
{
for(int i=0;i<c;i++)//将数组中的变量都赋值为2000
{
f[i]=a/5;
}
while(c)
{
b=c;//做分子
d=0;
g=2*c;//做分母
while(b)
{
d=d+f[b]*a;
g--;
f[b]=d%g;//保留余数
d=d/g;
g--;
b--;
if(b)
d*=b;
}
c-=14;
printf("%.4d",e+d/a);
e=d%a;
}
}
可以看到,在程序开始运行的时候,b被赋值为2800,g被赋值为2*2800=5600,所以在g第一次参与运算时,g的值为5600-1=2*2799+1,b第一次参与运算时的值为2799,所以第一次循环恰如上面公式中最小的括号所示,是2799/(2*2799+1),当进行下一次循环时,就用之前得出来的值加上2000再乘上2798/(2*2798+1)。即while(b)中的每一循环都是之前的结果加上2000再向外拓展一层括号。
由于只需要精确到前800位,所以2800次迭代就可以满足所需要的精确度。
当然,在C语言中是不存在能够精确到小数点后800位的数据类型,所以在计算的时候需要用高精度计算的方法。while(c)中的每次循环都将数组中剩余的数据(即之前除以分母时的余数)前移四位,这样可以得到接下来的四位小数。当然,在运算的时候可能会出现进位的现象,即d/a是大于0的情况,所以要将这个部分加到之前保留的四位数据中,即每次输出的是e+d/a.
在输出那里%.4d的意思是输出4位整数,即若整数少于四位时,从高位开始将缺少的位补全为0。这个的作用在那篇帖子和很多博文中解释的都是输出的整数不止4位的时候取前四位输出,并把后四位保留到e中。这是错误的。%.4d没有截取最高四位的功能,如果整数的位数大于等于4,输出的会是整数本身,而当整数不足4位时它才会起作用。这里之所以要用%.4d,是因为e是d对a的取余,但这个余数可能小于1000,如果直接输出%d就会导致小数中一些0的缺失。 所以要用0进行占位,即%.4d是每次都输出4位整数的保障。
至于c的作用,帖子里面和其他的博文基本上都忽略掉了。在这里c的步长为14,即每次运算都舍弃掉最后的14项,这样可以保证整个运算过程的精确。若c的步长过小,则会因为出现冗余导致误差;若c的步长过大,则会因为忽略掉不该忽略的项而导致误差。在原帖中说这个对精确度没有影响,这也是不对的。不是对精确度没有影响,而是它的影响不足以影响到前800位小数。如果步长为1,最多只能精确到九百多位;步长为14的时候,则可以保证前16275位小数的准确程度,这也是仅通过c的改变能达到的最高的精确度。一般来讲,步长过大比步长过小的精确度会高一些。具体会在后面的性能分析中给出详细的分析。
程序实现的过程大体上就是这样,下面说一下帖子和博文普遍回避的两个问题:
1.那个计算圆周率的公式是怎样得来的。
2.这个程序的性能如何?究竟能达到多高的精确度。
QQ952866274 http://blog.csdn.net/qq_31474315
一、公式推导
先把这个公式做一下变形:
可以发现右边变成了级数的形式,由于本人数学水平有限,想不到特别酷炫的方法。只能用两种比较常规的方法去解决它
1.和函数法
我做级数求和的题目时,一般会习惯性地想到和函数。这个我一开始想着直接乘上一个x的幂,结果怎么也凑不出来。联想到arcsinx的麦克劳林展开后x的次幂都是奇数,所以我把次幂改成n+1/2试了一下,于是就有了下面的步骤:
2.欧拉级数变换
在求完和函数之后,我又想着利用欧拉级数变换解决这个问题。
这个需要利用Leibniz公式:
这个公式的推导很容易,对arctanx泰勒展开然后令x=1就可以了