“外星人计算圆周率的程序”原理及性能分析（上）

    在之前有人给我发了这样一段代码，让我解释一下其中的原理：

#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;}}

    这样看起来结构上就清晰很多了。

    我们每次都想输出4位数，然而圆周率只有一位整数，所以就在给数组赋值的时候前移三位，使得第一次运算的时候我们就可以得到4位整数。这也是我们把数组中的变量赋值为2000而不是2的原因。在之后的每次循环中d加上的都是f[b]*a,而a是10000，也就是说将f[b]前移四位，这也是为了让我们读取接下来的四位小数。

    或许对有些人来说从这段代码看出那个公式还是有些困难，对公式进行一下变形之后就会容易很多：

    可以看到，在程序开始运行的时候，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就可以了

    对这个公式进行欧拉级数变换则会得到下面的步骤：

   

    当然，如果利用下面这个组合恒等式再结合欧拉级数变换就可以很快地得出我们想要的结论了。

    由于篇幅有限，具体的性能分析就在下一篇博文说了。

    在本文的结尾附上一个花絮吧。

    之前有人问过我（arcsinx）^2的泰勒展开怎么做，我当时给了一个特别复杂的多项式解法。在推完上面那个和函数之后，突然想到可以把那个当成一个结论来解决这个问题。如下图所示：

 　

　转载请注明出处！

　（未完待续）