*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。
要是有写的不详的地方,欢迎高手讨论