外星人计算pi程序的数学原理

*File:      Pi800.cName:     
分析外星人计算PI的程序Author:    zyl910Blog:     
http://blog.csdn.net/zyl910/Version:   V1.0Updata:   
2006-11-5
分析外星人计算PI的程序~~~~~~~~~~~~~~~~~~~~~~原始程序：
int a=10000,b,c=2800,d,e,f[2801],g;
main()
{
or(;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公式:    pi = 2 + --- * (2 + --- * (2 + --- * (2 + ...  (2 + ---- * (2 + ... ))...)))         

                               3             5              7                      2*k+1

 

为了便于理解程序改写为：

（代码参考csdn上的某位的博客，可运行）

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

以上代码来自CSDN博客，转载处：http://blog.csdn.net/zyl910/archive/2006/11/05/1368387.aspx

 

 

用2800 long 型数组保证精确小数点后800位的精度

上面的公式可转化为以下形式：

 

              1          1    2          1    2   3               1    2    3            k    

pi = 2 + ---*2 + ---*---*2 + ---*---*---*2......+ ---*---*---*... *  ---      *2 +........；

              3          3    5          3    5   7               3    5    7         2*k+1

 

 

 1    2    3             2800               1    1   1       1           1  2799             

 ---*---*---*... *  ----------   *2 < ---*---*---*...---*2 = (---)        *2      

 3    5    7         2*2800+1          2    2   2        2           2                                         

 

 

                    

       1   10    270                1  3  270       1  810

< [(---)     ]                  <[ (---)   ]      =  [---]

       2                                10                10

 

上面那一项的值小到小数点后800为位，即在计算pi的公式中，K=2800即可保证其值精确到小数点后800位（若有人提问后面还有其他的项相加，可以用数列方法证明2800后面的所有项的和小到小数点后800位），所以仅需要用2800项

 

 

将无限项的公式转化为有限项后

 

                                          1             2              3                              2799

算Pi公式:   1000* pi = 2k+ --- * (2k+ --- * (2k+ --- * (2k+ ...  (2k+    ----   * 2k))...)))         

                                          3             5              7                          2*2799+1

 

 

只需要粗略的看看上面的程序，我们就大概知道它的确是使用的那个迭代公式来计算Pi的了，不过不知道到现在为止你是否明白了f数组的用处。如果没有明白，请继续阅读。 

余下的工作就是将这个有限项的公式计算到小数点后800位，计算方法就是上面的将数据放大10000倍再做运算，分别用f[b] 保留 小数部分用d 保留整数部分， 扩大10000倍后整数放在d中外前迭代计算，余数部分放在f[b]中，因为余数很小扩大10000倍后尚在小数点之后所以对该次需要求出的4位数字精度没有影响，只用整数部分d继续往前迭代。从上面的有限项公式逐次向前计算，每一次外层循环输出4位数字，并把原来的余数都扩大10000倍再代入运算。余数扩大10000倍后，每一次外层循环后可以省掉14次迭代（最后面14项的余数项，扩大10000倍后相当于表示的小数点后第800到803位，已不需要，所以略去），减少计算次数。这就是c = c - 14。

 

 

 

 

要是有写的不详的地方，欢迎高手讨论