分析外星人计算Pi的程序

有一个只用4行代码就实现的计算Pi的程序，被称为外星人计算Pi的程序。

有许多人讨论分析了该程序的实现原理，如：http://blog.csdn.net/panqiaomu/archive/2006/05/07/711776.aspx

但我总感觉它分析得不够透彻，于是自己分析了一下。

1.将原程序修改成更易看懂的形式；

2.采用同样的算法，用Excel表格将Pi算了出来。

下载（注意修改下载后的扩展名）
 

/*
File:      Pi800.c
Name:      分析外星人计算PI的程序
Author:    zyl910
Blog:      http://blog.csdn.net/zyl910/
Version:   V1.0
Updata:    2006-11-5

分析外星人计算PI的程序
~~~~~~~~~~~~~~~~~~~~~~

原始程序：
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);
}

算Pi公式:
          1          2          3                    k
pi = 2 + --- * (2 + --- * (2 + --- * (2 + ...  (2 + ---- * (2 + ... ))...)))
          3          5          7                   2k+1

*/

int main(void)
{
 long a=10000;// 缩放系数
 long c=2800; // 迭代次数
 long f[2801];// 中间计算结果
 long d;      // 内循环误差累计项
 long e=0;    // 外循环误差累计项
 long b, g;   // 分子 与 分母. k/(2k+1)
 int i, j;

 for(i=0; i<c; i++) // 若完全依照公式，循环条件应为“for(i=1; i<=c; i++)”。幸好f[2800]的贡献非常小，取任意值都不会对精度造成太大的影响
  f[i] = 2 * a/10; // 2就是公式中的系数2。 a/10的意思是保留一个十进制位，因为输出是从个位3开始的
 while(c > 0)
 {
  d = 0;
  g = (c*2 + 1) - 2;      // 分母。因为每次循环都输出了4位，所以在后面运算时乘以了a，所以这里得 -2
  b = c;                  // 分子
  while(b > 0)
  {
   /* 根据公式，乘以分子 */
   d *= b;
   d += f[b]*a;  // 因为每次外循环都输出了4位
   /* 根据公式，除以分母 */
   f[b] = d % g; // 带分数的 分子部分
   d /= g;       // 带分数的 整数部分
   /* Next */
   g -= 2;
   b--;
  }
  printf("%.4d", e+d/a);
  e = d % a;
  c -= 14; // 因为精度固定为800位，每输出4位后，相当于精度需求降低了4位，所以每次可以少算14项
 }
 printf("/n");
 return 0;
}