古城

@ 2007.11.21

　

　

　

首先看一段程序：

#include <stdlib.h> 
#include <stdio.h> 

long 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); 
} 

运行一下，输出结果为：

31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185 



一、源程序 
本文分析下面这个很流行的计算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的作用就自己考虑吧。 

　

　

　

上一篇：ADO[转]

下一篇：TCP三次握手及原理