怒了，求高人解释程序算法，很简短的一个程序_CMS教程

外星人计算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=c2;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 "stdioh"

void main()

{

float pi=2;

int i;

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

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

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

getchar();

}

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

三、程序展开

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

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=c2;

b=c;

while(1)

{

d=d+f[b]a;

g--;

f[b]=d%g;

d=d/g;

g--;

b--;

if(b==0) break;

d=db;

}

c=c-14;

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

e=d%a;

}

注：

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

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

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

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

四、程序分析

要想计算出无限精度的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

1000pi = 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=c2; 这里的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为22799+1，你可以看到g做了分母

g--;

b--;

if(b==0) break;

d=db; 这里的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/(22799+1)，这

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

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

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

接下来的4位数呢？

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

0142

/——------

7 / 1

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

我们可以发现，在做除法的时候，我们通常把余数扩大之后再来计算，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=c2;

b=c;

while(1)

{

d=d+f[b]a;

g--;

f[b]=d%g;

d=d/g;

g--;

b--;

if(b==0) break;

d=db;

}

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

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

e=d%a;

}

最后的答案仍然正确。

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

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

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

!高斯消去法

subroutine agaus(a,b,n,x,l,js)

dimension a(n,n),x(n),b(n),js(n)

double precision a,b,x,t

l=1 !逻辑变量

do k=1,n-1

d=00

do i=k,n

do j=k,n

if (abs(a(i,j))>d) then

d=abs(a(i,j))

js(k)=j

is=i

end if

end do

end do !把行绝对值最大的元素换到主元位置

if (d+10==10) then

l=0

else !最大元素为0无解

if(js(k)/=k) then

do i=1,n

t=a(i,k)

a(i,k)=a(i,js(k))

a(i,js(k))=t

end do !最大元素不在K行，K行

end if

if(is/=k) then

do j=k,n

t=a(k,j)

a(k,j)=a(is,j)

a(is,j)=t !交换到K列

end do

t=b(k)

b(k)=b(is)

b(is)=t

end if !最大元素在主对角线上

end if !消去

if (l==0) then

write(,100)

return

end if

do j=k+1,n

a(k,j)=a(k,j)/a(k,k)

end do

b(k)=b(k)/a(k,k) !求三角矩阵

do i=k+1,n

do j=k+1,n

a(i,j)=a(i,j)-a(i,k)a(k,j)

end do

b(i)=b(i)-a(i,k)b(k)

end do

if (abs(a(n,n))+10==10) then

l=0

write(,100)

return

end if

x(n)=b(n)/a(n,n)

do i=n-1,1,-1

t=00

do j=i+1,n

t=t+a(i,j)x(j)

end do

x(i)=b(i)-t

end do

100 format(1x,'fail')

js(n)=n

do k=n,1,-1

if (js(k)/=k) then

t=x(k)

x(k)=x(js(k))

x(js(k))=t

end if

end do

return

end

program main

dimension a(4,4),b(4),x(4),js(4)

double precision a,b,x

real m1,m2,j

open(1,file="laiyitxt")

read(1,)m1,m2,j

close(1)

n=4

print,m1,m2,j

a(1,1)=m1cos(314159j/180)

a(1,2)=-m1

a(1,3)=-sin(314159j/180)

a(1,4)=0

a(2,1)=m1sin(314159j/180)

a(2,2)=0

a(2,3)=cos(314159j/180)

a(2,4)=0

a(3,1)=0

a(3,2)=m2

a(3,3)=-sin(314159j/180)

a(3,4)=0

a(4,1)=0

a(4,2)=0

a(4,3)=-cos(314159j/180)

a(4,4)=1

b(1)=0

b(2)=m198

b(3)=0

b(4)=m298

call agaus(a,b,n,x,l,js)

if (l/=0) then

write(,)"a1=",x(1),"a2=",x(2) ,"n1=",x(3),"n2=",x(4)

end if

end

!逆矩阵求解

SUBROUTINE qiuni(A,N,L,IS,JS)

DIMENSION A(N,N),IS(N),JS(N)

DOUBLE PRECISION A,T,D

L=1

DO K=1,N

D=00

DO I=K,N

DO J=K,N

IF(ABS(A(I,J))GTD) THEN !把最大的元素给D

D=ABS(A(I,J))

IS(K)=I

JS(K)=J

END IF

END DO

IF (D+10EQ10)THEN

L=0

WRITE(,200)

RETURN

END IF

200 FORMAT(1X,'ERRNOT INV')

DO J=1,N

T=A(K,J)

A(K,J)=A(IS(K),J)

A(IS(K),J)=T

END DO

DO I=1,N

T=A(I,K)

A(I,K)=A(I,JS(K))

A(I,JS(K))=T

END DO

A(K,K)=1/A(K,K)

DO J=1,N

IF(JNEK)THEN

A(K,J)=A(K,J)A(K,K)

END IF

END DO

DO I=1,N

IF(INEK)THEN

DO J=1,N

IF(JNEK)THEN

A(I,J)=A(I,J)-A(I,K)A(K,J)

END IF

END DO

END IF

END DO

DO I=1,N

IF(INEK)THEN

A(I,K)=-A(I,K)A(K,K)

END IF

END DO

DO K=N,1,-1

DO J=1,N

T=A(K,J)

A(K,J)=A(JS(K),J)

A(JS(K),J)=T

END DO

DO I=1,N

T=A(I,K)

A(I,K)=A(I,IS(K))

A(I,IS(K))=T

END DO

RETURN

END

SUBROUTINE BRMUL(A,B,N,C)

DIMENSION A(N,N),B(N),C(N)

DOUBLE PRECISION A,B,C

DO I=1,N

DO J=1,N

C(I)=00

DO L=1,N

C(I)=C(I)+A(I,L)B(L)

END DO

RETURN

END

program main

DIMENSION A(4,4),B(4,1),C(4,1),IS(4),JS(4)

DOUBLE PRECISION A,B,C

REAL M1,M2,JD

OPEN(1,FILE='LAIYITXT')

READ(1,) M1,M2,JD

PRINT,M1,M2,JD

CLOSE(1)

A(1,1)=M1COS(314JD/180)

A(1,2)=-M1

A(1,3)=-SIN(314JD/180)

A(1,4)=0

A(2,1)=M1SIN(314JD/180)

A(2,2)=0

A(2,3)=COS(314JD/180)

A(2,4)=0

A(3,1)=0

A(3,2)=M2

A(3,3)=-SIN(314JD/180)

A(3,4)=0

A(4,1)=0

A(4,2)=0

A(4,3)=-COS(314JD/180)

A(4,4)=1

B(1,1)=0

B(2,1)=M198

B(3,1)=0

B(4,1)=M298

CALL QIUNI(A,4,L,IS,JS)

CALL BRMUL(A,B,4,C)

WRITE(,) (C(I,1),I=1,4)

END

画图

USE MSFLIB

INTEGER status

TYPE(xycoord) xy

status=SETCOLORRGB(#FFFFFF)

status1=SETCOLORRGB(#0000FF)

OPEN(1,FILE="GTXT")

READ(1,) G1,G2,G3,G4

OPEN(2,FILE="NTXT")

READ(2,) N1,N2,N3,N4

CALL MOVETO(INT(20),INT(20),xy)

status=LINETO(INT(40),INT(G1))

status=LINETO(INT(80),INT(G2))

status=LINETO(INT(120),INT(G3))

status=LINETO(INT(160),INT(G4))

CALL SETLINESTYLE(#FF00)

CALL MOVETO(INT(20),INT(20),xy)

status1=LINETO(INT(40),INT(N1))

status1=LINETO(INT(80),INT(N2))

status1=LINETO(INT(120),INT(N3))

status1=LINETO(INT(160),INT(N4))

READ(,)

END

如果对您有帮助，请记得采纳为满意答案，谢谢！祝您生活愉快！

虽然算法与计算机程序密切相关，但二者也存在区别：计算机程序是算法的一个实例，是将算法通过某种计算机语言表达出来的具体形式；同一个算法可以用任何一种计算机语言来表达。

算法列表

图论

路径问题

0/1边权最短路径

BFS

非负边权最短路径（Dijkstra）

可以用Dijkstra解决问题的特征

负边权最短路径

Bellman-Ford

Bellman-Ford的Yen-氏优化

差分约束系统

Floyd

广义路径问题

传递闭包

极小极大距离 / 极大极小距离

Euler Path / Tour

圈套圈算法

混合图的 Euler Path / Tour

Hamilton Path / Tour

特殊图的Hamilton Path / Tour 构造

生成树问题

最小生成树

第k小生成树

最优比率生成树

0/1分数规划

度限制生成树

连通性问题

强大的DFS算法

无向图连通性

割点

割边

二连通分支

有向图连通性

强连通分支

2-SAT

最小点基

有向无环图

拓扑排序

有向无环图与动态规划的关系

二分图匹配问题

一般图问题与二分图问题的转换思路

最大匹配

有向图的最小路径覆盖

0 / 1矩阵的最小覆盖

完备匹配

最优匹配

稳定婚姻

网络流问题

网络流模型的简单特征和与线性规划的关系

最大流最小割定理

最大流问题

有上下界的最大流问题

循环流

最小费用最大流 / 最大费用最大流

弦图的性质和判定

组合数学

解决组合数学问题时常用的思想

逼近

递推/动态规划

概率问题

Polya定理

计算几何 / 解析几何

计算几何的核心：叉积 / 面积

解析几何的主力：复数

基本形

点

直线，线段

多边形

凸多边形 / 凸包

凸包算法的引进，卷包裹法

Graham扫描法

水平序的引进，共线凸包的补丁

完美凸包算法

怒了，求高人解释程序算法，很简短的一个程序

发表评论

评论列表（0条）