现在的位置: 首页 > 综合 > 正文

关于Pi的计算

2013年09月18日 ⁄ 综合 ⁄ 共 6001字 ⁄ 字号 评论关闭

一个看起来很牛的计算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/2 = Sum{k!/(2k+1)!!, k=0,1,2,...}
 

其实上面的程序效率并不高,下面转贴更多Pi的公式,来源于:
http://mathworld.wolfram.com/PiFormulas.html

There are many formulas of pi of many types. Among others, these include series, products, geometric constructions, limits, special values, and pi iterations.

pi is intimately related to the properties of circles and spheres. For a circle of radius r, the circumference and area are given by

C = 2pir (1)
A = pir^2. (2)

Similarly, for a sphere of radius r, the surface area and volume enclosed are

S = 4pir^2 (3)
V = 4/3pir^3. (4)

An exact formula for pi in terms of the inverse tangents of unit fractions is Machin's formula

1/4pi==4tan^(-1)(1/5)-tan^(-1)(1/(239)). (5)

There are three other Machin-like formulas, as well as thousands of other similar formulas having more terms.

GregorySeries

Gregory Eric Weisstein's World of Biography and Leibniz Eric Weisstein's World of Biography found

pi/4 = sum_(k==1)^(infty)((-1)^(k+1))/(2k-1) (6)
= 1-1/3+1/5-... (7)

(Wells 1986, p. 50), which is known as the Gregory series and may be obtained by plugging x==1 into the Leibniz series for tan^(-1)x. The error after the nth term of this series in the Gregory series is larger than (2n)^(-1) so this sum converges so slowly that 300 terms are not sufficient to calculate pi correctly to two decimal places! However, it can be transformed to

pi==sum_(k==1)^infty(3^k-1)/(4^k)zeta(k+1), (8)

where zeta(z) is the Riemann zeta function (Vardi 1991, pp. 157-158; Flajolet and Vardi 1996), so that the error after k terms is  approx (3/4)^k.

An infinite sum series to Abraham Sharp (ca. 1717) is given by

pi==sum_(k==0)^infty(2(-1)^k3^(1/2-k))/(2k+1) (9)

(Smith 1953, p. 311). Additional simple series in which pi appears are

1/4pisqrt(2) = sum_(k==1)^(infty)[((-1)^(k+1))/(4k+1)+((-1)^(k+1))/(4k-3)] (10)
= 1+1/3-1/5-1/7+1/9+1/(11)-... (11)
1/4(pi-3) = sum_(k==1)^(infty)((-1)^(k+1))/(2k(2k+1)(2k+2)) (12)
= 1/(2.3.4)-1/(4.5.6)+1/(6.7.8)-... (13)
1/6pi^2 = sum_(k==1)^(infty)1/(k^2) (14)
= 1+1/4+1/9+1/(16)+1/(25)+... (15)
1/8pi^2 = sum_(k==1)^(infty)1/((2k-1)^2) (16)
= 1+1/(3^2)+1/(5^2)+1/(7^2)+... (17)

(Wells 1986, p. 53).

In 1666, Newton Eric Weisstein's World of Biography used a geometric construction to derive the formula

pi = 3/4sqrt(3)+24int_0^(1/4)sqrt(x-x^2)dx (18)
= (3sqrt(3))/4+24(1/(12)-1/(5.2^5)-1/(28.2^7)-1/(72.2^9)-...), (19)

which he used to compute pi (Wells 1986, p. 50; Borwein et al. 1989; Borwein and Bailey 2003, pp. 105-106). The coefficients can be found from the integral

I(x) = intsqrt(x-x^2)dx (20)
= 1/4(2x-1)sqrt(x-x^2)-1/8sin^(-1)(1-2x) (21)

by taking the series expansion of I(x)-I(0) about 0, obtaining

I(x)==2/3x^(3/2)-1/5x^(5/2)-1/(28)x^(7/2)-1/(72)x^(9/2)-5/(704)x^(11/2)+... (22)

(Sloane's A054387 and A054388). Using Euler's Eric Weisstein's World of Biographyconvergence improvement transformation gives

pi/2 = 1/2sum_(n==0)^(infty)((n!)^22^(n+1))/((2n+1)!)==sum_(n==0)^(infty)(n!)/((2n+1)!!) (23)
= 1+1/3+(1.2)/(3.5)+(1.2.3)/(3.5.7)+... (24)
= 1+1/3(1+2/5(1+3/7(1+4/9(1+...)))) (25)

(Beeler et al. 1972, Item 120).

This corresponds to plugging x==1/sqrt(2) into the power series for the hypergeometric function _2F_1(a,b;c;x),

(sin^(-1)x)/(sqrt(1-x^2))==sum_(i==0)^infty((2x)^(2i+1)(i!)^2)/(2(2i+1)!)==_2F_1(1,1;3/2;x^2)x. (26)

Despite the convergence improvement, series (◇) converges at only one bit/term. At the cost of a square root, Gosper has noted that x==1/2 gives 2 bits/term,

1/9sqrt(3)pi==1/2sum_(i==0)^infty((i!)^2)/((2i+1)!), (27)

and x==sin(pi/10) gives almost 3.39 bits/term,

pi/(5sqrt(phi+2))==1/2sum_(i==0)^infty((i!)^2)/(phi^(2i+1)(2i+1)!), (28)

where phi is the golden ratio. Gosper also obtained

pi==3+1/(60)(8+(2.3)/(7.8.3)(13+(3.5)/(10.11.3)(18+(4.7)/(13.14.3)(23+...)))). (29)

A spigot algorithm for pi is given by Rabinowitz and Wagon (1995; Borwein and Bailey 2003, pp. 141-142).

More amazingly still, a closed form expression giving a digit-extraction algorithm which produces digits of pi (or pi^2) in base-16 was discovered by Bailey et al. (Bailey et al. 1997, Adamchik and Wagon 1997),

pi==sum_(n==0)^infty(4/(8n+1)-2/(8n+4)-1/(8n+5)-1/(8n+6))(1/(16))^n. (30)

This formula, known as the BBP formula, was discovered using the PSLQ algorithm (Ferguson et al. 1999) and is equivalent to

pi==int_0^1(16y-16)/(y^4-2y^3+4y-4)dy. (31)

A similar formula was subsequently discovered by Ferguson, leading to a two-dimensional lattice of such formulas which can be generated by these two formulas given by

pi==sum_(k==0)^infty((4+8r)/(8k+1)-(8r)/(8k+2)-(4r)/(8k+3)-(2+8r)/(8k+4)-(1+2r)/(8k+5)-(1+2r)/(8k+6)+r/(8k+7))(1/(16))^k (32)

for any complex value of r (Adamchik and Wagon), giving the BBP formula as the special case r==0.

There is a series of BPP-type formulas for pi in powers of (-1)^k, the first few independent formulas of which are

pi = 4sum_(k==0)^(infty)((-1)^k)/(2k+1) (33)
= 3sum_(k==0)^(infty)(-1)^k[1/(6k+1)+1/(6k+5)] (34)
= 4sum_(k==0)^(infty)(-1)^k[1/(10k+1)-1/(10k+3)+1/(10k+5)-1/(10k+7)+1/(10k+9)] (35)
= sum_(k==0)^(infty)(-1)^k[3/(14k+1)-3/(14k+3)+3/(14k+5)+4/(14k+7)+4/(14k+9)-4/(14k+11)+4/(14k+13)] (36)
= sum_(k==0)^(infty)(-1)^k[2/(18k+1)+3/(18k+3)+2/(18k+5)-2/(18k+7)-2/(18k+11)+2/(18k+13)+3/(18k+15)+2/(18k+17)] (37)
= sum_(k==0)^(infty)(-1)^k[3/(22k+1)-3/(22k+3)+3/(22k+5)-3/(22k+7)+3/(22k+9)+8/(22k+11)+3/(22k+13)-3/(22k+15)+3/(22k+17)-3/(22k+19)+1/(22k+21)]. (38)

Similarly, there are a series of BPP-type formulas for pi in powers of 2^k, the first few independent formulas of which are

pi = sum_(k==0)^(infty)1/(16^k)[4/(8k+1)-2/(8k+4)-1/(8k+5)-1/(8k+6)] (39)
= 1/2sum_(k==0)^(infty)1/(16^k)[8/(8k+2)+4/(8k+3)+4/(8k+4)-1/(8k+7)] (40)
= 1/(16)sum_(k==0)^(infty)1/(256^k)[(64)/(16k+1)-(32)/(16k+4)-(16)/(16k+5)-(16)/(16k+6)+4/(16k+9)-2/(16k+12)-1/(16k+13)-1/(16k+14)] (41)
= 1/(32)sum_(k==0)^(infty)1/(256^k)[(128)/(16k+2)+(64)/(16k+3)+(64)/(16k+4)-(16)/(16k+7)+8/(16k+10)+4/(16k+11)+4/(16k+12)-1/(16k+15)] (42)
= 1/(32)sum_(k==0)^(infty)1/(4096^k)[(256)/(24k+2)+(192)/(24k+3)-(256)/(24k+4)-(96)/(24k+6)-(96)/(24k+8)+(16)/(24k+10)-4/(24k+12)-3/(24k+15)-6/(24k+16)-2/(24k+18)-1/(24k+20)] (43)
= 1/(64)sum_(k==0)^(infty)1/(4096^k)[(256)/(24k+1)+(256)/(24k+2)-(384)/(24k+3)-(256)/(24k+4)-(64)/(24k+5)+(96)/(24k+8)+(64)/(24k+9)+(16)/(24k+10)+8/(24k+12)-4/(24k+13)+6/(24k+15)+6/(24k+16)+1/(24k+17)+1/(24k+18)-1/(24k+20)-1/(24k+21)] (44)
= 1/(96)sum_(k==0)^(infty)1/(4096^k)[(256)/(24k+2)+(64)/(24k+3)+(128)/(24k+5)+(352)/(24k+6)+(64)/(24k+7)+(288)/(24k+8)+(128)/(24k+9)+(80)/(24k+10)+(20)/(24k+12)-(16)/(24k+14)-1/(24k+15)+6/(24k+16)-2/(24k+17)-1/(24k+19)+1/(24k+20)-2/(24k+21)] (45)
= 1/(96)sum_(k==0)^(infty)1/(4096^k)[(256)/(24k+1)+(320)/(24k+3)+(256)/(24k+4)-(192)/(24k+5)-(224)/(24k+6)-(64)/(24k+7)-(192)/(24k+8)-(64)/(24k+9)-(64)/(24k+10)-(28)/(24k+12)-4/(24k+13)-5/(24k+15)+3/(24k+17)+1/(24k+18)+1/(24k+19)+1/(24k+21)-1/(24k+22)] (46)
= 1/(96)sum_(k==0)^(infty)1/(4096^k)[(512)/(24k+1)-(256)/(24k+2)+(64)/(24k+3)-(512)/(24k+4)-(32)/(24k+6)+(64)/(24k+7)+(96)/(24k+8)+(64)/(24k+9)+(48)/(24k+10)-(12)/(24k+12)-8/(24k+13)-(16)/(24k+14)-1/(24k+15)-6/(24k+16)-2/(24k+18)-1/(24k+19)-1/(24k+20)-1/(24k+21)] (47)
= 1/(4096)sum_(k==0)^(infty)1/(65536^k)[(16384)/(32k+1)-(8192)/(32k+4)-(4096)/(32k+5)-(4096)/(32k+6)+(1024)/(32k+9)-(512)/(32k+12)-(256)/(32k+13)-(256)/(32k+14)+(64)/(32k+17)-(32)/(32k+20)-(16)/(32k+21)-(16)/(32k+22)+4/(32k+25)-2/(32k+28)-1/(32k+29)-1/(32k+30)] (48)
= 1/(4096)sum_(k==0)^(infty)1/(65536^k)[(32768)/(32k+2)+(16384)/(32k+3)+(16384)/(32k+4)-(4096)/(32k+7)+(2048)/(32k+10)+(1024)/(32k+11)+(1024)/(32k+12)-(256)/(32k+15)+(128)/(32k+18)+(64)/(32k+19)+(64)/(32k+20)-(16)/(32k+23)+8/(32k+26)+4/(32k+27)+4/(32k+28)-1/(32k+31)]. (49)

F. Bellard found the rapidly converging BBP-type formula

pi==1/(2^6)sum_(n==0)^infty((-1)^n)/(2^(10n))(-(2^5)/(4n+1)-1/(4n+3)+(2^8)/(10n+1)-(2^6)/(10n+3)-(2^2)/(10n+5)-(2^2)/(10n+7)+1/(10n+9)). (50)

A related integral is

pi==(22)/7-int_0^1(x^4(1-x)^4)/(1+x^2)dx (51)

(Dalzell 1944; Dalzell 1971; Le Lionnais 1983, p. 22; Borwein et al. 2004, p. 3; Boros and Moll 2004, p. 125; Lucas 2005; Borwein et al. 2006, p. 14). This integral was known by K. Mahler in the mid-1960s and appears in an exam at the University of Sydney in November 1960 (Borwein et al. 2004, p. 3). Beukers (2000) and Boros and Moll (2004, p. 126) state this it is not clear if these exists a natural choice of rational polynomial whose integral between 0 and 1 produces pi-333/106, where 333/106 is the next convergent. However, an integral exists for the fourth convergent, namely

pi==(355)/(113)-1/(3164)int_0^1(x^8(1-x)^8(25+816x^2))/(1+x^2)dx. (52)

(Lucas 2005; Bailey et al. 2006, p. 219). In fact, Lucas (2005) gives a few other such integrals.

Backhouse (1995) used the identity

I_(m,n) = int_0^1(x^m(1-x)^n)/(1+x^2)dx (53)
= 2^(-(m+n+1))sqrt(pi)Gamma(m+1)Gamma(n+1)_3F_2(1,1/2(m+1),1/2(m+2);1/2(m+n+2),1/2(m+n+3);-1) (54)
= a+bpi+cln2 (55)

for positive integer m and n and where a, b, and c are rational constant to generate a number of formulas for pi. In particular, if 2m-n=0 (mod 4), then c==0 (Lucas 2005).

PiFormulasWagonIdentity

An even more general identity due to Wagon is given by

pi+4tan^(-1)z+2ln((1-2z-z^2)/(z^2+1))
            ==sum_(k==0)^infty1/(16^k)[(4(z+1)^(8k+1))/(8k+1)-(2(z+1)^(8k+4))/(8k+4)-((z+1)^(8k+5))/(8k+5)-((z+1)^(8k+6))/(8k+6)] (56)

(Borwein and Bailey 2003, p. 141), which holds over a region of the complex plane excluding two triangular portions symmetrically placed about the real axis, as illustrated above.

Following the discovery of the base-16 digit BBP formula and related formulas, similar formulas in other bases were investigated. Borwein and Galway have recently shown that there is no non-binary Bailey-Borwein-Plouffe-type formula, although this does not rule out a completely different scheme for digit extraction algorithms in other bases (Bailey 2002).

S. Plouffe has devised an algorithm to compute the nth digit of pi in any base in O(n^3(logn)^3) steps.

A slew of additional identities due to Ramanujan Eric Weisstein's World of Biography, Catalan, and Newton Eric Weisstein's World of Biography are given by Castellanos (1988, pp. 86-88), including several involving sums of Fibonacci numbers. Ramanujan found

sum_(k==0)^infty((-1)^k(4k+1)[(2k-1)!!]^3)/([(2k)!!]^3)==sum_(k==0)^infty((-1)^k(4k+1)[Gamma(k+1/2)]^3)/(pi^(3/2)[Gamma(k+1)]^3)==2/pi (57)

(Hardy 1923; Hardy 1924; Hardy 1999, p. 7).

Plouffe (2006) found the beautiful formula

pi==72sum_(n==1)^infty1/(n(e^(npi)-1))-96sum_(n==1)^infty1/(n(e^(2npi)-1))+24sum_(n==1)^infty1/(n(e^(4npi)-1)). (58)

An interesting infinite product formula due to Euler which relates pi and the nth prime p_n is

pi = 2/(product_(n==1)^(infty)[1+(sin(1/2pip_n))/(p_n)]) (59)
= 2/(product_(n==2)^(infty)[1+((-1)^((p_n-1)/2))/(p_n)]) (60)

(Blatner 1997, p. 119), plotted above as a function of the number of terms in the product.

A method similar to Archimedes' Eric Weisstein's World of Biography can be used to estimate pi by starting with an n-gon and then relating the area of subsequent 2n-gons. Let beta be the angle from the center of one of the polygon's segments,

beta==1/4(n-3)pi, (61)

then

pi==(2sin(2beta))/((n-3)product_(k==0)^(infty)cos(2^(-k)beta)) (62)

(Beckmann 1989, pp. 92-94).

Vieta (1593) was the first to give an exact expression for pi by taking n==4 in the above expression, giving

cosbeta==sinbeta==1/(sqrt(2))==1/2sqrt(2), (63)

which leads to an infinite product of nested radicals,

2/pi==sqrt(1/2)sqrt(1/2+1/2sqrt(1/2))sqrt(1/2+1/2sqrt(1/2+1/2sqrt(1/2)))... (64)

(Wells 1986, p. 50; Beckmann 1989, p. 95). However, this expression was not rigorously proved to converge until Rudio (1892).

A related formula is given by

pi==lim_(n->infty)2^(n+1)sqrt(2-sqrt(2+sqrt(2+sqrt(2+...+sqrt(2))))_()_(n)), (65)

which can be written

pi==lim_(n->infty)2^(n+1)pi_n, (66)

where pi_n is defined using the iteration

pi_n==sqrt((1/2pi_(n-1))^2+[1-sqrt(1-(1/2pi_(n-1))^2)]^2) (67)

with pi_0==sqrt(2) (J. Munkhammar, pers. comm., April 27, 2000). The formula

pi==2lim_(m->infty)sum_(n==1)^msqrt([sqrt(1-((n-1)/m)^2)-sqrt(1-(n/m)^2)]^2+1/(m^2)) (68)

is also closely related.

A pretty formula for pi is given by

pi==(product_(n==1)^(infty)(1+1/(4n^2-1)))/(sum_(n==1)^(infty)1/(4n^2-1)), (69)

where the numerator is a form of the Wallis formula for pi/2 and the denominator is a telescoping sum with sum 1/2 since

1/(4n^2-1)==1/2(1/(2n-1)-1/(2n+1)) (70)

(Sondow 1997).

A particular case of the Wallis formula gives

pi/2==product_(n==1)^infty[((2n)^2)/((2n-1)(2n+1))]==(2.2)/(1.3)(4.4)/(3.5)(6.6)/(5.7)... (71)

(Wells 1986, p. 50). This formula can also be written

lim_(n->infty)(2^(4n))/(n(2n; n)^2)==pilim_(n->infty)(n[Gamma(n)]^2)/([Gamma(1/2+n)]^2)==pi, (72)

where (n; k) denotes a binomial coefficient and Gamma(x) is the gamma function (Knopp 1990). Euler Eric Weisstein's World of Biography obtained

pi==sqrt(6(1+1/(2^2)+1/(3^2)+1/(4^2)+...)), (73)

which follows from the special value of the Riemann zeta function zeta(2)==pi^2/6. Similar formulas follow from zeta(2n) for all positive integers n.

An infinite sum due to Ramanujan Eric Weisstein's World of Biography is

1/pi==sum_(n==0)^infty(2n; n)^3(42n+5)/(2^(12n+4)) (74)

(Borwein et al. 1989; Borwein and Bailey 2003, p. 109; Bailey et al. 2006, p. 44). Further sums are given in Ramanujan Eric Weisstein's World of Biography (1913-14),

4/pi==sum_(n==0)^infty((-1)^n(1123+21460n)(2n-1)!!(4n-1)!!)/(882^(2n+1)32^n(n!)^3) (75)

and

1/pi = sqrt(8)sum_(n==0)^(infty)((1103+26390n)(2n-1)!!(4n-1)!!)/(99^(4n+2)32^n(n!)^3)

抱歉!评论已关闭.