解密那一段神奇的计算圆周率的程序(一）

网络上广为流传的一段程序，计算圆周率的程序，看起来是这样的：

#include<stdio.h>
//2015.8.12
int 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);
}

我把它重新写过：

#include   "stdio.h"

int a=10000, bit, cur=2800, sum, rem, fill[2801], q;

void cal_to_next(){
  sum += fill[bit]*a;
  --q;
  fill[bit] = sum % q;
  sum /=q;
  --q;
  --bit;
}

main()   
{   
  while(bit-cur){    
    fill[bit++]=a/5;  
  }

  do{
    sum=0;
    q=cur*2;
    bit = cur;
    cal_to_next();    
     
    while(bit){  
      sum *= bit;
      cal_to_next();
    };
    
    printf("%.4d",rem + sum/a);
    
    rem  = sum % a;
    cur  = cur-14;    
  } while(cur);
return 0;   
}  

这样应该好懂了。如果还是不明白，那么，lisp版本就是这样的：

(defun frac-pi(n)
  (do ((s 0 (+ s v))
       (a 1 (+ a 1))
       (v 2 (* v (/ a b)))
       (b 3 (+ b 2)))
      ((= a n) s)))

如果还是不明白，那么，excel版本应该是这样的：

ppi.png

其中，唯一的公式是 D3=D2*C3，向下复制。D列所有数字的和就是要求的圆周率的近似值。

学数学的人是这样写的：

pi.png

（虽然数学符号很混乱的说，但果然看起来紧凑，比lisp看起来还要紧凑。）
这么多解释了，你一定能看懂这个公式了。
那么，这个公式是如何推导出来的呢？且听下回分解。