π的计算方法有哪些

圆周率π的计算方法

圆周率π的计算方法，是一个饶有趣味，值得探讨的问题。最直观的计算方法自然是从几何上着手，历史上也正是如此，这便是割圆法。设一半径为1的圆，作这个圆的内接正n边形，用此正n边形的周长去近似圆的周长。显然当n→∞时，正n边形的周长就无限趋近于圆周长，求得正n边形周长后除以直径便求出了圆周率。

I.

从几何上观察，可知：正n边形周长随n递增而递增，但始终是个有限值。割法如图1：

向左转|向右转

图1 割圆法

设圆半径为1，令半弦长AB=2a，AC=2c，OG和OD分别是等腰△OAB和△OAC 的中线。则我们要做的只是求出c关于a的表达式c=c（a）.令GC=b，根据勾股定理有:

向左转|向右转

(1)

进而有

向左转|向右转

(2)

得到此式后，编写计算机程序就很容易了，C语言程序如下：

#include

#include

main()

{double a,b,c,d,pi;double sqrt(double);int i,j,n;a=0.5;b=0;c=0;d=0.5;scanf("%d",&n);

for(i=1;i<=n;i++)

{b=sqrt(1-a*a);c=(1-b)*0.5;d=sqrt(c);a=d;}

j=pow(2,n)*3;pi=2*d*j;printf("%d\n",j);printf("%f\n",pi)}

这里有一个问题就是a的初值如何选择？显然越简单直观越好，而已知对于圆内接正六边形的每一条边长等于圆的半径。所以取a=0.5，程序中参数n是对正六边形分割的次数，d的作用是当输入n=0（正六边形）的时候，得到π=3，此所谓的“径圆一三”。将这个文件保存为文本，在linux下用“gcc -lm”命令编译后，打开编译后得到的文件就能执行。

在古代可没有电子计算机，而祖冲之利用割圆法算得圆周率介于3.1415926和

3.1415927之间，可见古人之伟大！

II.

上面的方法简单直观，但是缺点也很明显。计算机在底层只能做“加减乘除四则整数运算”，显然开根号运算还是要通过转化为整数运算（级数展开等）才最后到硬件级计算。那么我们能否直接用整数的四则运算得到π的值？有！而且方法是多样的，其中一种叫作“Wallis公式”，有几种表达方式。如下：

向左转|向右转

(3)

或

向左转|向右转

(4)

或

(5)

下面证明这个公式：

令

向左转|向右转

(6)

利用分部积分法

于是有关系式

向左转|向右转

(7)

从上式可知I0=1，I1=π/4.根据这两个初值条件有向左转|向右转

(8)

或者

向左转|向右转

(9)

其中m=0,1,2,...而由(7)式也可知

向左转|向右转

(10)

将(10)式代入(9)式

向左转|向右转

即

向左转|向右转

(11)

其中

由式(11)可知Wm>0且有上限，而

说明Wm随着m的增大递增，所以如下极限存在，且由夹逼准则得其值

Wallis公式得证。

实际上Wallis公式的发现在微积分建立之前，其探寻过程限于篇幅不在这里给出，这也反映出同一个问题可以有不同的论证方法，也令我们不得不佩服古人的智慧。

III.

虽然Wallis公式比割圆法要易于计算得多，但是Wallis公式在形势上仍显复杂，且全部乘除算法也难以提高计算机计算效率，最好是有乘除项之和，如：

反观(6)式，实际上令x=cosθ，则有dx=-sinθdθ.式(6)变为

如果令x=sinθ，则只变换形式不影响结果。那我们设想利用其它的三角函数能否得到同样的结果？令

向左转|向右转

(12)

注意这里的积分上限改成了π/4，因为π/2>θ>π/4的时候tanθ>1，将导致积分发散。

对(12)式做一个小变换

于是有关系式

向左转|向右转

(13)

而初值T0=π/4，观察规律有

...

总结规律得

向左转|向右转

(14)

其中m=1,2,3,...而从式(12)中可知

结合(14)式，得到

向左转|向右转

(15)

或者

向左转|向右转

(16)

显然这种方法形式上比前两种方法要简单得多，计算机执行的时候也能更高效。

而在我前面的文章中讲过幂级数的应用，arctanθ展开为幂级数（泰勒级数）后表达式为

向左转|向右转

(17)

该级数的收敛域为[-1,1]，将x=1代入，则得到式(15)，这又是一个殊途同归的例子！