圆周率π的计算方法

2012-12-16 00:07阅读：

http://blog.sina.cn/dpool/blog/u/1459730044

圆周率π的计算方法
圆周率π的计算方法，是一个饶有趣味，值得探讨的问题。最直观的计算方法自然是从几何上着手，历史上也正是如此，这便是割圆法。设一半径为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 <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',j);
printf('%f',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公式在形势上仍显复杂，且全部乘除算法也难以提高计算机计算效率，最好是有乘除项之和，如：
圆周率π的计算方法