有一个只用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;
}
|
