查人人中国名人网wow60副本:如何求圆周率？有哪些方法？
来源：百度文库 编辑：查人人中国名人网 时间：2024/04/29 01:29:56
这个比较难懂，我补充个别人的解释吧
————————————————————
一、源程序
本文分析下面这个很流行的计算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的作用就由读者自己考虑吧。
--
※ 来源:·武汉白云黄鹤站 bbs.whnet.edu.cn·[FROM: 150.12.180.101]
--
※ 来源:·计算机92 202.117.49.124·[FROM: 202.117.11.22]
发信人: Beta (.NET Beta Documentation), 信区: Algorithm
标 题: Re: 巨牛的算pi程序——绝对巨牛！
发信站: 武汉白云黄鹤站 (2001年06月25日11:28:07 星期一), 站内信件
这个程序的算法我也不明白，不过我找到了一些相关资料。
我搜索过，但是找不到对这个算法的介绍。看来想弄明白
只有找原作者Dik T. Winter了。
http://www.cs.unb.ca/~alopez-o/math-faq/mathtext/node12.html
There are essentially 3 different methods to calculate pi to many
decimals.
One of the oldest is to use the power series expansion of
atan(x) = x - x^3/3 + x^5/5 - ...
together with formulas like
pi = 16*atan(1/5) - 4*atan(1/239).
This gives about 1.4 decimals per term.
A second is to use formulas coming from Arithmetic-Geometric mean
computations. A beautiful compendium of such formulas is given in
the book pi and the AGM. They have the advantage of converging
quadratically, i.e. you double the number of decimals per iteration.
For instance, to obtain 1 000 000 decimals, around 20 iterations are
sufficient. The disadvantage is that you need FFT type multiplication
to get a reasonable speed, and this is not so easy to program.
A third one comes from the theory of complex multiplication of elliptic
curves, and was discovered by S. Ramanujan. This gives a number of
beautiful for mulas, but the most useful was missed by Ramanujan and
discovered by the ChuDnovsky's. It is the following (slightly modified
for ease of programming):
Set k_1 = 545140134; k_2 = 13591409; k_3 = 640320; k_4 = 100100025;
k_5 = 327843840; k_6 = 53360;
Then pi = (k_6 sqrt(k_3))/(S), where
S = sum_(n = 0)^oo (-1)^n ((6n)!(k_2 + nk_1))/(n!^3(3n)!(8k_4k_5)^n)
The great advantages of this formula are that
1) It converges linearly, but very fast (more than 14 decimal digits
per term).
2) The way it is written, all operations to compute S can be programmed
very simply. This is why the constant 8k_4k_5 appearing in the denominator
has been written this way instead of 262537412640768000. This is how the
Chudnovsky's have computed several billion decimals.
An interesting new method was recently proposed by David Bailey, Peter
Borwein and Simon Plouffe. It can compute the Nth hexadecimal digit of Pi
efficiently without the previous N-1 digits. The method is based on the
formula:
pi = sum_(i = 0)^oo (1 16^i) ((4 8i + 1) - (2 8i + 4) - (1 8i + 5)
- (1 8i + 6))
in O(N) time and O(log N) space.
The following 160 character C program, written by Dik T. Winter at CWI,
computes pi to 800 decimal digits.
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);}
http://numbers.computation.free.fr/Constants/TinyPrograms/tinycodes.html
Tiny programs for constants computation
(The smallest C codes to compute classical mathematical constants)
This page contains some very tiny code to compute classical mathematical
constants.
Wanted : any shorter C codes or tiny C codes for z(3), the Euler constant
g (!) or any other classical constant. Send any new tiny code to :
xavier.gourdon@free.fr or to sebah.pascal@fnac.net .
All the following programs are what is now called Spigot-Algorithm , that
that is, algorithms were the digits are calculated and printed one at the
time.
1 Computation of pi
The following tiny C code computes digits of pi. Boris Gourevitch kindly
sent me the information that this program is from Dik T. Winter (cwi
institute, Holland).
int a=10000,b,c=8400,d,e,f[8401],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);
}
The program has 158 characters long. Did you find the algorithm used ?
~~~~~~~~~~~~~faint...
Sebastian Wedeniws worked hard to find a shorter program (142 characters)
main()
{
int a=1e4,c=3e3,b=c,d=0,e=0,f[3000],g=1,h=0;
for(;b;!--b?printf("%04d",e+d/a),e=d%a,h=b=c-=15:f[b]
=(d=d/g*b+a*(h?f[b]:2e3) )%(g=b*2-1));
}
http://www.bruchez.org/erik/computer/pi/
……This code seems to have been written by Dik T.Winter, CWI, Amsterdam.
If you want to compile it, be aware that some compilers don't initialize
variables automatically, so you'll have to initialize b, d, e and g to
zero. I was told that ANSI compilers should initialize the variables to
zero.
I didn't know how this program worked.
Christoph Haenel (CHRHAENEL@aol.com) gave me this explanation:
"The program computes
1000*pi = 2000 + 1/3*(2000+2/5*(2000+3/7*(2000+ ...)))
thru a 'spigot' algorithm with spigots of 4 decimals. The main task of
the program is to convert the mixed radices in this formula to the fixed
radix 10000 that the humans like.The factor 14 is a matter of exactness:
4*k decimals need 14*k array entries. Cf. Stanley Rabinowitz and Stan
Wagon, A Spigot Algorithm for the Digits of pi, American Mathematical
Monthly, Vol 102 (1995), 3, pp 195-203."
Christoph has also written a 147 chars version of Dik's code!
A MORE READABLE FORM...
Here is a more readable version of the program:
int a = 10000, b = 0;
int c = 8400, d = 0;
int e = 0, f[8401], g = 0;
void
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)
;
}
http://forum.swarthmore.edu/dr.math/problems/piformulac.html
References :
(This is a short version for a more comprehensive list contact
Juhana Kouhia at jk87377@cc.tut.fi)
J. M. Borwein, P. B. Borwein, and D. H. Bailey, "Ramanujan,
Modular Equations, and Approximations to Pi", American Mathematical
Monthly, vol. 96, no. 3 (March 1989), p. 201 - 220.
P. Beckman
A history of pi
Golem Press, CO, 1971 (fourth edition 1977)
J.M. Borwein and P.B. Borwein
The arithmetic-geometric mean and fast computation of elementary
functions
SIAM Review, Vol. 26, 1984, pp. 351-366
J.M. Borwein and P.B. Borwein
More quadratically converging algorithms for pi
Mathematics of Computation, Vol. 46, 1986, pp. 247-253
J.M. Borwein and P.B. Borwein
Pi and the AGM - a study in analytic number theory and
computational complexity
Wiley, New York, 1987
Shlomo Breuer and Gideon Zwas
Mathematical-educational aspects of the computation of pi
Int. J. Math. Educ. Sci. Technol., Vol. 15, No. 2, 1984,
pp. 231-244
David Chudnovsky and Gregory Chudnovsky
The computation of classical constants, Columbia University,
Proc. Natl. Acad. Sci. USA, Vol. 86, 1989.
Y. Kanada and Y. Tamura
Calculation of pi to 10,013,395 decimal places based on the
Gauss-Legendre algorithm and Gauss arctangent relation
Computer Centre, University of Tokyo, 1983
Morris Newman and Daniel Shanks
On a sequence arising in series for pi
Mathematics of computation, Vol. 42, No. 165, Jan 1984,
pp. 199-217
E. Salamin
Computation of pi using arithmetic-geometric mean
Mathematics of Computation, Vol. 30, 1976, pp. 565-570
D. Shanks and J.W. Wrench, Jr.
Calculation of pi to 100,000 decimals
Mathematics of Computation, Vol. 16, 1962, pp. 76-99
Daniel Shanks
Dihedral quartic approximations and series for pi
J. Number Theory, Vol. 14, 1982, pp.397-423
David Singmaster
The legal values of pi
The Mathematical Intelligencer, Vol. 7, No. 2, 1985
Stan Wagon
Is pi normal?
The Mathematical Intelligencer, Vol. 7, No. 3, 1985
J.W. Wrench, Jr.
The evolution of extended decimal approximations to pi
The Mathematics Teacher, Vol. 53, 1960, pp. 644-650