圆周率π的计算方法

2024年4月6日发(作者：2021年江西高考数学试卷分析)

圆周率π的计算方法

圆周率的计算方法

古人计算圆周率，一般是用割圆法。即用圆的内接或外切正多边形来逼近圆的周长。

Archimedes用正96边形得到圆周率小数点后3位的精度；刘徽用正3072边形得到5位

精度；Ludolph Van Ceulen用正262边形得到了35位精度。这种基于几何的算法计算

量大，速度慢，吃力不讨好。随着数学的发展，数学家们在进行数学研究时有意无意地发

现了许多计算圆周率的公式。

1、 Machin公式

这个公式由英国天文学教授John Machin于1706年发现。他利用这个公式计算到了

100位的圆周率。Machin公式每计算一项可以得到1.4位的十进制精度。因为它的计算

过程中被乘数和被除数都不大于长整数，所以可以很容易地在计算机上编程实现。

用马青公式计算Pi至小数点后100位程序

program Pi_Value;

{$APPTYPE CONSOLE}

1

//将Pi计算精确小数点后100位

//Machin公式

//Pi=16arctan(1/5)-4arctan(1/239)

uses

SysUtils;

const

N=100;

S=2*N+50;

aNum=5;

bNum=239;

type

Num=array [1..S] of byte;

//初始化数组

2

procedure AZero(var arr:Num);

var

i:smallint;

begin

for i:=1 to S do

arr:=0;

end;

//除法

procedure Division(var arr:Num;const b:smallint);

var

c,y,i:smallint;

begin

c:=0;

3

for i:=1 to S do

begin

y:=arr+c*10;

c:=y mod b;

arr:=y div b;

end;

end;

//加法

procedure Addition(var arr:Num;const b:Num);

var

i,y,c:smallint;

begin

c:=0;

4

for i:=S downto 1 do

begin

y:=arr+b+c;

if y>=10 then

begin

c:=1;

arr:=y-10;

end

else

begin

c:=0;

arr:=y;

end;

5

end;

end;

//减法

procedure Minus(var arr:Num;const b:Num);

var

i,y,c:smallint;

begin

c:=0;

for i:=S downto 1 do

begin

y:=arr-b-c;

if y<0 then

begin

6

c:=1;

arr:=10+y;

end

else

begin

c:=0;

arr:=y;

end;

end;

end;

var

tag:boolean;

a,b,Ra,Rb,t:Num;

7

i,j:smallint;

begin

AZero(t);

Ra:=t;Rb:=t;

tag:=true;

writeln(\'计算中，请等待......\');

for i:=1 to N do

begin

a:=t;b:=t;

a[1]:=16;b[1]:=4;

for j:=1 to i*2-1 do

begin

Division(a,aNum);

8

DiVision(b,bNum);

end;

Division(a,i*2-1);

Division(b,i*2-1);

if tag then

begin

tag:=false;

Addition(Ra,a);

Addition(Rb,b);

end

else

begin

tag:=true;

9

Minus(Ra,a);

Minus(Rb,b);

end;

end;

Minus(Ra,Rb);

writeln(\'计算结果如下：\');

writeln(Ra[1],\'.\');

for i:=2 to N+1 do

write(Ra);

readln;

End.

还有很多类似于Machin公式的反正切公式。在所有这些公式中，Machin公式似乎

是最快的了。虽然如此，如果要计算更多的位数，比如几千万位，Machin公式就力不从

心了。下面介绍的算法，在PC机上计算大约一天时间，就可以得到圆周率的过亿位的精

度。这些算法用程序实现起来比较复杂。因为计算过程中涉及两个大数的乘除运算，要用

10

FFT(Fast Fourier Transform)算法。FFT可以将两个大数的乘除运算时间由O(n2)缩短为

O(nlog(n))。

2、 Ramanujan公式

1914年，印度数学家Srinivasa Ramanujan在他的论文里发表了一系列共14条圆周

率的计算公式，这是其中之一。这个公式每计算一项可以得到8位的十进制精度。1985

年Gosper用这个公式计算到了圆周率的17,500,000位。

1989年，David & Gregory Chudnovsky兄弟将Ramanujan公式改良成为：

这个公式被称为Chudnovsky公式，每计算一项可以得到15位的十进制精度。1994

年Chudnovsky兄弟利用这个公式计算到了4,044,000,000位。Chudnovsky公式的另一

个更方便于计算机编程的形式是：

3、AGM(Arithmetic-Geometric Mean)算法

11

Gauss-Legendre公式：

初值：

重复计算：

最后计算：

这个公式每迭代一次将得到双倍的十进制精度，比如要计算100万位，迭代20次就

够了。1999年9月Takahashi和Kanada用这个算法计算到了圆周率的206,158,430,000

位，创出新的世界纪录。

4、Borwein四次迭代式：

初值：

重复计算：

最后计算：

这个公式由Jonathan Borwein和Peter Borwein于1985年发表，它四次收敛于圆

周率。

12

5、 Bailey-Borwein-Plouffe算法

这个公式简称BBP公式，由David Bailey, Peter Borwein和Simon Plouffe于1995

年共同发表。它打破了传统的圆周率的算法，可以计算圆周率的任意第n位，而不用计算

前面的n-1位。这为圆周率的分布式计算提供了可行性。1997年，Fabrice Bellard找到

了一个比BBP快40％的公式：

参考百度文库、人民网-数学故事欣赏

13