文档章节

分析-外星人计算Pi的程序

RoyceInWh
 RoyceInWh
发布于 2016/06/20 16:04
字数 1967
阅读 2
收藏 0

FROM http://hyry.dip.jp/tech/blog/index.html?month=200201



分析-外星人计算Pi的程序

RY  程序人生 2002/01/22

* 本文作于2002年,最早发表于网易论坛。现在用Google搜索一下,真还有不少转载的。

源程序

本文分析下面这个很流行的计算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的作用就由读者自己考虑吧。

本文转载自:http://blog.csdn.net/jingxia2008/article/details/49493365

RoyceInWh

RoyceInWh

粉丝 5
博文 240
码字总数 1282
作品 0
武汉
程序员
私信 提问
Python学习之路12-外星人

本系列是对入门书籍《Python编程:从入门到实践》的笔记整理,属于初级内容。标题顺序采用书中标题。 本章主要是对上一篇的继续,添加“外星人”,“外星人”与飞船的交互。 1. 回顾项目 开发...

vpointer701
2018/04/16
0
0
VUE+WebPack前端游戏设计:实现外星人攻击建筑物时的冒烟效果

玩过红警或是星际的玩家都知道,当子弹或对手攻击建筑物时,建筑物会产生冒烟效果,并且逐步变形,当攻击足够大后,建筑物会爆炸毁灭,这种动态特效极大的增强了游戏的视觉观赏性和娱乐性,本...

望月从良
2018/03/23
0
0
计算器 abacus 技术文档之二----初步设计

======================================= 计算器 abacus 的下载地址:http://www.oschina.net/code/snippet73693213725 如果你有关于 abacus 的问题或者建议,请发邮件至 zhoucosin@163.co......

zhcosin
2012/10/30
599
0
例2.2 圆柱体的表面积

/* 例2.2 圆柱体的表面积    输入底面半径r和高h,输出圆柱体的表面积,保留3位小数,格式见样例。     样例输入:3.5 9 样例输出:274.889 【分析】    圆柱体的表面积由3部分组成:...

技术小美
2017/11/16
0
0
外星人大战----------------------游戏开发(四)

前面已经实现外星人的移动,现在开始完成射杀外星人。我们将要使用sprite.groupcollide()检测两个编组的碰撞。我们要在碰撞的时候立马就让外星人消失,所以在更新位置的时候就检测有无碰撞。...

shinhwa96
2018/11/24
0
0

没有更多内容

加载失败,请刷新页面

加载更多

只需一步,在Spring Boot中统一Restful API返回值格式与统一处理异常

统一返回值 在前后端分离大行其道的今天,有一个统一的返回值格式不仅能使我们的接口看起来更漂亮,而且还可以使前端可以统一处理很多东西,避免很多问题的产生。 比较通用的返回值格式如下:...

晓月寒丶
今天
59
0
区块链应用到供应链上的好处和实际案例

区块链可以解决供应链中的很多问题,例如记录以及追踪产品。那么使用区块链应用到各产品供应链上到底有什么好处?猎头悬赏平台解优人才网小编给大家做个简单的分享: 使用区块链的最突出的优...

猎头悬赏平台
今天
27
0
全世界到底有多少软件开发人员?

埃文斯数据公司(Evans Data Corporation) 2019 最新的统计数据(原文)显示,2018 年全球共有 2300 万软件开发人员,预计到 2019 年底这个数字将达到 2640万,到 2023 年达到 2770万。 而来自...

红薯
今天
61
0
Go 语言基础—— 通道(channel)

通过通信来共享内存(Java是通过共享内存来通信的) 定义 func service() string {time.Sleep(time.Millisecond * 50)return "Done"}func AsyncService() chan string {retCh := mak......

刘一草
今天
57
0
Apache Flink 零基础入门(一):基础概念解析

Apache Flink 的定义、架构及原理 Apache Flink 是一个分布式大数据处理引擎,可对有限数据流和无限数据流进行有状态或无状态的计算,能够部署在各种集群环境,对各种规模大小的数据进行快速...

Vincent-Duan
今天
59
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部