分享
 
 
 

外星人计算Pi的程序

王朝other·作者佚名  2006-01-09
窄屏简体版  字體: |||超大  

外星人计算Pi的程序

一、源程序

本文分析下面这个很流行的计算PI的小程序。下面这个程序初看起来似乎摸不到头脑,

不过不用担心,当你读完本文的时候就能够基本读懂它了。

程序一:很牛的计算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

至于这个公式为什么能够计算出PI,已经超出了本文的能力范围。

下面要做的事情就是要分析清楚程序是如何实现这个公式的。

我们先来验证一下这个公式:

程序二:Pi公式验证程序

#include "stdio.h"

void main()

{

float pi=2;

int i;

for(i=100;i>=1;i--)

pi=pi*(float)i/(2*i+1)+2;

printf("%f\n",pi);

getchar();

}

上面这个程序的结果是3.141593。

三、程序展开

在正式分析程序之前,我们需要对程序一进行一下展开。我们可以看出程序一都是使用

for循环来完成计算的,这样做虽然可以使得程序短小,但是却很难读懂。根据for循环

的运行顺序,我们可以把它展开为如下while循环的程序:

程序三:for转换为while之后的程序

int a=10000,b,c=2800,d,e,f[2801],g;

main() {

int i;

for(i=0;i<c;i++)

f[i]=a/5;

while(c!=0)

{

d=0;

g=c*2;

b=c;

while(1)

{

d=d+f[b]*a;

g--;

f[b]=d%g;

d=d/g;

g--;

b--;

if(b==0) break;

d=d*b;

}

c=c-14;

printf("%.4d",e+d/a);

e=d%a;

}

}

注:

for([1];[2];[3]) {[4];}

的运行顺序是[1],[2],[4],[3]。如果有逗号操作符,例如:d=0,g=c*2,则先运行d=0,

然后运行g=c*2,并且最终的结果是最后一个表达式的值,也就是这里的c*2。

下面我们就针对展开后的程序来分析。

四、程序分析

要想计算出无限精度的PI,我们需要上述的迭代公式运行无数次,并且其中每个分数也

是完全精确的,这在计算机中自然是无法实现的。那么基本实现思想就是迭代足够多次

,并且每个分数也足够精确,这样就能够计算出PI的前n位来。上面这个程序计算800位

,迭代公式一共迭代2800次。

int a=10000,b,c=2800,d,e,f[2801],g;

这句话中的2800就是迭代次数。

由于float或者double的精度远远不够,因此程序中使用整数类型(实际是长整型),分

段运算(每次计算4位)。我们可以看到输出语句 printf("%.4d",e+d/a); 其中%.4就是

把计算出来的4位输出,我们看到c每次减少14( c=c-14;),而c的初始大小为2800,因

此一共就分了200段运算,并且每次输出4位,所以一共输出了800位。

由于使用整型数运算,因此有必要乘上一个系数,在这个程序中系数为1000,也就是说

,公式如下:

1 2 3 k

1000*pi = 2k+ --- * (2k+ --- * (2k+ --- * (2k+ ... (2k+ ---- * (2k+ ... )).

..)))

3 5 7 2k+1

这里的2k表示2000,也就是f[2801]数组初始化以后的数据,a=10000,a/5=2000,所以下面

的程序把f中的每个元素都赋值为2000:

for(i=0;i<c;i++)

f[i]=a/5;

你可能会觉得奇怪,为什么这里要把一个常数储存到数组中去,请继续往下看。

我们先来跟踪一下程序的运行:

while(c!=0) 假设这是第一次运行,c=2800,为迭代次数

{

d=0;

g=c*2; 这里的g是用来做k/(2k+1)中的分子

b=c; 这里的b是用来做k/(2k+1)中的分子

while(1)

{

d=d+f[b]*a; f中的所有的值都为2000,这里在计算时又把系数扩大了

a=10000倍。

这样做的目的稍候介绍,你可以看到

输出的时候是d/a,所以这不影

计算

g--;

f[b]=d%g; 先不管这一行

d=d/g; 第一次运行的g为2*2799+1,你可以看到g做了分母

g--;

b--;

if(b==0) break;

d=d*b; 这里的b为2799,可以看到d做了分子。

}

c=c-14;

printf("%.4d",e+d/a);

e=d%a;

}

只需要粗略的看看上面的程序,我们就大概知道它的确是使用的那个迭代公式来计算Pi

的了,不过不知道到现在为止你是否明白了f数组的用处。如果没有明白,请继续阅读。

d=d/g,这一行的目的是除以2k+1,我们知道之所以程序无法精确计算的原因就是这个除

法。即使用浮点数,答案也是不够精确的,因此直接用来计算800位的Pi是不可能的。那

么不精确的成分在哪里?很明显:就是那个余数d%g。程序用f数组把这个误差储存起来

,再下次计算的时候使用。现在你也应该知道为什么d=d+f[b]*a;中间需要乘上a了吧。

把分子扩大之后,才好把误差精确的算出来。

d如果不乘10000这个系数,则其值为2000,那么运行d=d/g;则是2000/(2*2799+1),这

种整数的除法答案为0,根本无法迭代下去了。

现在我们知道程序就是把余数储存起来,作为下次迭代的时候的参数,那么为什么这么

做就可以使得下次迭代出来的结果为

接下来的4位数呢?

这实际上和我们在纸上作除法很类似:

0142

/——------

7 / 1

10

7

---------------

30

28

---------------

20

14

---------------

60

.....

我们可以发现,在做除法的时候,我们通常把余数扩大之后再来计算,f中既然储存的是

余数,而f[b]*a;则正好把这个余数扩大了a倍,然后如此循环下去,可以计算到任意精

度。

这里要说明的是,事实上每次计算出来的d并不一定只有4位数,例如第一次计算的时候

,d的值为31415926,输出4位时候,把低四位的值储存在e中间,e=d%a,也就是5926。

最后,这个c=c-14不太好理解。事实上没有这条语句,程序计算出来的仍然正确。只是

因为如果迭代2800次,无论分数如何精确,最后Pi的精度只能够达到800。

你可以把程序改为如下形式尝试一下:

for(i=0;i<800;i++)

{

d=0;

g=c*2;

b=c;

while(1)

{

d=d+f[b]*a;

g--;

f[b]=d%g;

d=d/g;

g--;

b--;

if(b==0) break;

d=d*b;

}

// c=c-14; 不要这句话。

printf("%.4d",e+d/a);

e=d%a;

}

最后的答案仍然正确。

不过我们可以看到内循环的次数是c次,也就是说每次迭代计算c次。而每次计算后续位

数的时候,迭代次数减少14,而不影响精度。为什么会这样,我没有研究。另外最后的

e+d/a,和e=d/a的作用就由读者自己考虑吧。

--

 
 
 
免责声明:本文为网络用户发布,其观点仅代表作者个人观点,与本站无关,本站仅提供信息存储服务。文中陈述内容未经本站证实,其真实性、完整性、及时性本站不作任何保证或承诺,请读者仅作参考,并请自行核实相关内容。
2023年上半年GDP全球前十五强
 百态   2023-10-24
美众议院议长启动对拜登的弹劾调查
 百态   2023-09-13
上海、济南、武汉等多地出现不明坠落物
 探索   2023-09-06
印度或要将国名改为“巴拉特”
 百态   2023-09-06
男子为女友送行,买票不登机被捕
 百态   2023-08-20
手机地震预警功能怎么开?
 干货   2023-08-06
女子4年卖2套房花700多万做美容:不但没变美脸,面部还出现变形
 百态   2023-08-04
住户一楼被水淹 还冲来8头猪
 百态   2023-07-31
女子体内爬出大量瓜子状活虫
 百态   2023-07-25
地球连续35年收到神秘规律性信号,网友:不要回答!
 探索   2023-07-21
全球镓价格本周大涨27%
 探索   2023-07-09
钱都流向了那些不缺钱的人,苦都留给了能吃苦的人
 探索   2023-07-02
倩女手游刀客魅者强控制(强混乱强眩晕强睡眠)和对应控制抗性的关系
 百态   2020-08-20
美国5月9日最新疫情:美国确诊人数突破131万
 百态   2020-05-09
荷兰政府宣布将集体辞职
 干货   2020-04-30
倩女幽魂手游师徒任务情义春秋猜成语答案逍遥观:鹏程万里
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案神机营:射石饮羽
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案昆仑山:拔刀相助
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案天工阁:鬼斧神工
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案丝路古道:单枪匹马
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:与虎谋皮
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:李代桃僵
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案镇郊荒野:指鹿为马
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案金陵:小鸟依人
 干货   2019-11-12
倩女幽魂手游师徒任务情义春秋猜成语答案金陵:千金买邻
 干货   2019-11-12
 
推荐阅读
 
 
 
>>返回首頁<<
 
靜靜地坐在廢墟上,四周的荒凉一望無際,忽然覺得,淒涼也很美
© 2005- 王朝網路 版權所有