圆周率π的计算方法

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

I.

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

http://s5/mw690/5701b67cgd0e65b603a74&690

图1 割圆法

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

                                 http://s11/mw690/5701b67cgd0e65c43351a&690                                               (1)

进而有

http://s7/bmiddle/5701b67cgd0e65cd39656&690                            (2)

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

#include <stdio.h>

#include <math.h>

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公式”，有几种表达方式。如下：

http://s4/bmiddle/5701b67cgd0e65db15643&690                        (3)

或

http://s4/bmiddle/5701b67cgd0e65edc82c3&690                        (4)

或

http://s5/bmiddle/5701b67cgd0e65f653204&690                            (5)

下面证明这个公式：

令

http://s16/bmiddle/5701b67cgd0e6601cf9bf&690                           (6)

利用分部积分法

http://s6/bmiddle/5701b67cgd0e660897cc5&690             

于是有关系式

http://s8/bmiddle/5701b67cgd0e6611fbbc7&690                             (7)

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

http://s3/bmiddle/5701b67cgd0e661a77702&690                         (8)

或者

http://s14/bmiddle/5701b67cgd0e66267f43d&690                       (9)

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

http://s12/bmiddle/5701b67cgd0e6639881fb&690                           (10)

将(10)式代入(9)式

http://s14/bmiddle/5701b67cg7b4a3d42674d&690                     

http://s14/bmiddle/5701b67cgd0e66541e60d&690         

即

http://s14/bmiddle/5701b67cgd0e666388b4d&690                         (11)

其中

http://s11/bmiddle/5701b67cgd0e666d5016a&690                   

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

http://s6/bmiddle/5701b67cg7b4a3d8d15a5&690              

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

http://s3/bmiddle/5701b67cgd0e668183d82&690                                  

Wallis公式得证。

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

III.

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

http://s2/bmiddle/5701b67cgd0e668ae1b31&690                                  

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

http://s8/bmiddle/5701b67cgd0e669bbe107&690            

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

http://s5/bmiddle/5701b67cgd0e66af4a814&690                             (12)

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

对(12)式做一个小变换

http://s8/bmiddle/5701b67cgd0e66bd97797&690           

于是有关系式

http://s14/bmiddle/5701b67cgd0e66d21dced&690                            (13)

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

http://s14/bmiddle/5701b67cgd0e66e82519d&690

 

 

...

总结规律得

http://s1/bmiddle/5701b67cgd0e66fb9be20&690           (14)

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

http://s7/bmiddle/5701b67cg7b4a3e83a7d6&690                               

结合(14)式，得到

http://s2/bmiddle/5701b67cgd0e671f67791&690              (15)

或者

http://s13/bmiddle/5701b67cgd0e672df22fc&690                          (16)

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

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

http://s14/bmiddle/5701b67cgd0e6739934ad&690         (17)

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

 

 参考资料：

http://episte.math.ntu.edu.tw/articles/sm/sm_27_05_1/index.html