π等于多少

圆周率π的⼏种计算⽅法与分析

圆周率（）的计算⽅法

慈悲⼼待⼈，时时有平安；智慧⼼安⼰，处处有⾃在！

From©

圆周率介绍

从代数上来讲，就是⼀个⽆理数（irrationalnumber，即⽆限不循环⼩数），也是⼀个超越数(transcendentalnumbers,即不能作

为有理系数多项式⽅程的根的数)，与三⾓函数有密切关系，并且与⽆穷级数的和有很⼤关系。

从⼏何上来说，就是⼀个⽐值，特殊图形圆的周长与其直径的⽐值，也就是单位圆的半周长，还可以是单位圆的⾯积或圆⾯积与其半

径的平⽅之⽐；也有⼈提出：如果⽤圆的周长与其半径的⽐值作为圆周率的话，也是可以的。

是⼀个数学常数，互联⽹上有很多关于圆周率的介绍，⽐较权威的参见

2021年8⽉19⽇，由瑞⼠格⼒孙应⽤技术⼤学(UniversityofAppliedSciencesoftheGrisons—UASGrisons)打破记录，历时仅

114天(4⽉28⽇~8⽉19⽇)，计算到了万亿位(Trilliondigits),精确位数为。最后10位数字为

。

打破了2020年1⽉29⽇由美国阿拉巴马州的TimothyMullican耗时8个⽉（303天）创造的万亿位!

⾕歌⼯程师EmmaHarukaIwao于2019年1⽉宣布已经计算到了万亿位，为了纪念圆周率⽇-

PAIDay（每年的3⽉14⽇）。

下⾯主要阐述现代分析⽅法对计算圆周率的贡献。

泰勒级数理论（Taylor’stheoremandTaylorries）

多项式逼近任意函数

多项式(polynomial)逼近(approximate)未知函数在给定点附近的值，实际上就是找到多项式的系数(coefficients):

靠近给定点逼近函数,此处假定(assuming)函数在有微分(differentiable).从⽽有,下

⾯试着找出系数.

>>>f=Function('f')

>>>ries(f(x))

f(0)+x*Subs(Derivative(f(_x),_x),_x,0)+x**2*Subs(Derivative(f(_x),(_x,2)),_x,0)/2+x**3*Subs(Derivative(f(_x),(_x,3)),_x,0)/6+x**4*Subs(Derivativ

e(f(_x),(_x,4)),_x,0)/24+x**5*Subs(Derivative(f(_x),(_x,5)),_x,0)/120+O(x**6)

实函数或复函数(arealorcomplex-valuedfunction)在有⽆限阶微分，则可以表⽰为幂级数(powerries):

此处表⽰的阶乘(denotesthefactorial),表⽰在处的微分(derivative).可以写成更加紧凑的求和式⼦：

泰勒级数

的0阶导数就是它⾃⼰，且.当,此时的泰勒级数⼜叫做麦克劳林级数.

举例(Forinstance)

π

π

62.862831853071750′′′′

7817924264

50

31.4=31.4×10=123.14×1013

fa,a,⋯,a01n

P(x)=na+0a(x−1a)+a(x−2a)+2⋯+a(x−na)+nR((x−a))n+1

x=af(x)f(x)x=anthf(x)≈P(x)n

a,n=n0,1,2,⋯,n

f(x)x=a

f(x)≈f(a)+(x−

1!

f(a)′

a)+(x−

2!

f(a)′′

a)+

2(x−

3!

f(a)′′′

a)+

3⋯

n!nf(a)(n)fx=anth

f(x)=(x−a).n=0

∑∞

n!

f(a)(n)

n

ff(x−a)=00!=1a=0

指数函数(exponentialfunction)

指数函数是重要的基本初等函数之⼀。⼀般地，函数为常数且叫做指数函数。

指数函数的定义域是;

值域;

单调递增：时,单调递减：时。

注意，在指数函数的定义表达式中，在前的系数必须是数1,

指数函数应⽤到值上的这个函数写为。还可以等价的写为，这⾥的是数学常数，就是⾃然对数的底数，近似等于

，还称为欧拉数.它的严格定义是。

作为实数变量的函数，的图像总是正的(在轴之上)并递增(从左向右看)。它永不触及轴，尽管它可以⽆限程度地靠近轴(所

以，轴是这个图像的⽔平渐近线。它的反函数是⾃然对数，其定义域是所有正数。

指数函数在处的泰勒级数展开式为：

对求导，始终不变，即，

且.

剩下的分⼦(numerator)为,分母(denominator)为.

欧拉公式(Eulerformula):

>>>ries(sin(x),n=10)

x-x**3/6+x**5/120-x**7/5040+x**9/362880+O(x**10)

>>>ries(cos(x),n=10)

1-x**2/2+x**4/24-x**6/720+x**8/40320+O(x**10)

>>>ries(exp(x),n=10)

1+x+x**2/2+x**3/6+x**4/24+x**5/120+x**6/720+x**7/5040+x**8/40320+x**9/362880+O(x**10)

指数函数在的泰勒级数展开式为

令(纯虚数),此处虚数单位(imaginaryunit),满⾜.

右边为

所以与正弦函数和余弦函数的泰勒展开式对⽐，得到

y=ax(aa>0，a=1)

x∈R

(Domain)(0,+∞)

a>10

ax

a=eexp(x)exe

2.718281828e=1+n→+∞

lim(

n1)n

xexxxx

xln(x)x>0

exx=00

n=0

∑

n!

xn

=++++++⋯

0!

x0

1!

x1

2!

x2

3!

x3

4!

x4

5!

x5

=1+x+++++⋯

2

x2

6

x3

24

x4

120

x5

exxe=

dxk

dk

xex

e=01

(x−0)nn!

e=ixcos(x)+isin(x)

ezz=00

ez=n=0

∑

n!

zn

=++++++⋯

0!

z0

1!

z1

2!

z2

3!

z3

4!

z4

5!

z5

=1+z+++++⋯

2

z2

6

z3

24

z4

120

z5

z=ixii=4k1,i=4k+1i,i=4k+2−1,i=4k+3−i,k∈N+

1+ix−−i++i−−i+⋯

2!

x

3!

x

4!

x

5!

x

6!

x

7!

x

上⾯也是函数和的泰勒级数展开式.(请参照,)

右边可以简写为,简写为,三个字母的缩写，产⽣⼀个新的函数。

故有

写成紧凑的求和式为：

可见，是偶函数，是奇函数。

显然，我们也可以通过泰勒展开得到上述式⼦。多阶导数有下列规律：

在处的偶次导数为0,在奇次导数为

欧拉恒等式(Euleridentity)的由来

只要在欧拉公式中令,即可得到上帝恒等式——欧拉恒等式

神奇的是：这个欧拉恒等式包含了这五个神奇的数,还有这两个基本运算符，故称为上帝公式。

多项式逼近理论应⽤(Polynomialapproximationtheorem)

下⾯应⽤多项式逼近理论，得到⼏个常⽤的三⾓函数（）的乘积展开式。主要还是代数学中的求根⽅法。

问题:求⽅程的解。

cos(x)=1−+

2!

x2

−

4!

x4

+

6!

x6

⋯（偶函数）

sin(x)=x−+

3!

x3

−

5!

x5

+

7!

x7

⋯（奇函数）

cos(x)sin(x)∀x

cos(x)+isin(x)cis(x)

e=cos(x)+isin(x)ix

cos(x)=n=0

∑∞

(2n)!

(−1)xn2n

=1−+−+⋯.

2!

x2

4!

x4

6!

x6

sin(x)=n=0

∑

(2n+1)!

(−1)x

=x−+−+⋯.

3!

x3

5!

x5

7!

x7

cos(x)sin(x)

sin(x)=(2k)=dx

dsin(x)(2k)

(−1)sin(x)k

sin(x)=(2k+1)(−1)cos(x)k

sin(x)x=0sin(x)x=0(−1)k

cis(x)=

e=cos(x)+isin(x)

ix

x=π

e=−1,e+1=0

iπiπ

e,i,π,0,1+,=

sin(x),cos(x),tan(x)

sin(x)=0

1.⾸先,有解

2.假设是多项式函数，由多项式有根的代数基本理论(thefundamentaltheoremofalgebra)即多项式逼近理论

(polynomialapproximationtheorem)

令，求极限得到

我们得到：

令,有

参见

因为,所以有偶函数的⽆穷级数乘积公式：

sin(x)=0{kπ，k=0,±1,±2,⋯.}

sin(x)

sin(x)=c∗(kπ−k=1

∏∞

x)(kπ+x)x

=

x

sin(x)

c∗(kπ−k=1

∏∞

x)(kπ+x)

x→0

c=k=1

∏∞

kπ22

1

sin(x)=x(1−)(1+)k=1

∏

kπ

x

kπ

x

=x1−k=1

∏∞

[(kπ)2

x2]

x=2

π

2

π

=⋅k=1

∏

2k−1

2k

2k+1

2k

=k=1

∏∞

1−4k2

1

1

=1+k=1

∏∞

[4k−12

1]

=⋅⋅⋅⋅⋅⋯

1

2

3

2

3

4

5

4

5

6

7

6

=⋅⋅⋅⋯

3

4

15

16

35

36

63

64

sin(π/2−x)=cos(x)cos(x)

fc(x)=cos(x)1≈(−x)1−1+

2

π

k=1

∏∞

(

kππ/2−x)(

kππ/2−x)

=(−x)1−+1+−

2

π

k=1

∏∞

(

2k1

kπ

x)(

2k1

kπ

x)

同理可得

所有正整数平⽅的倒数和

通过泰勒级数展开理论和多项式逼近理论这两种不同⽅法所得到的结果进⾏⽐较。

>>>ries(sin(x),n=10)

x-x**3/6+x**5/120-x**7/5040+x**9/362880+O(x**10)

1.泰勒展开理论得到

2.多项式逼近理论得到

3.⽐较的系数可得

所以

收敛速度还可以接受。

fromsympyimportsummation,Symbol,oo

k=Symbol('k',integer=True)

summation(1/k**2,(k,1,oo))

计算

格雷⼽⾥-莱布尼茨级数(Gregory-Leibnizries)

格雷⼽⾥-莱布尼茨级数于1671年2⽉15⽇提出。利⽤反正切函数的泰勒级数展开得到的级数算法。

>>>fromsympyimport*

>>>x,y,z=symbols('xyz')

>>>ries(atan(x),n=10)

x-x**3/3+x**5/5-x**7/7+x**9/9+O(x**10)

fc(x)=cos(x)2≈1−1+1−(

π

2x)k=1

∏

((2k−1)π

2x)((2k+1)π

2x)

=1−1+1−(

π

2x)k=1

∏∞

(k−1/2

x/π)(k+1/2

x/π)

sin(x)=x−++⋯

3!

x3

5!

x5

sin(x)=x(1−)k=1

∏∞

(kπ)2

x2

x3

3!

1

=k=1

∑∞

(kπ)2

1

6

π2

=k=1

∑∞

k2

1

=1+++⋯++⋯

4

1

9

1

n2

1

π=

6k=1

∑

k2

1

4

π

反正切函数的级数(riesfortheinvertangentfunction),也称为Gregory’sries:

具体推导过程如下：

莱布尼茨公式（Leibnizformula）

可以通过代⼊上述反正切函数的级数中得到.这⾥有⼀个疑问，既然要求,⼜怎么能取呢？这⾥存在⼀

个的问题。

正如你看到，这个级数的收敛极慢！不可取。但是，通过变形，得到很多计算的公式。

反正切函数的特殊⽤途

马青——JohnMachin(1680-1751)英国伦敦⼈,于1706年发现公式：

上式得益于Gregory有关反正切函数的幂级数展开式：

令，则有：

鉴于Gregory-Leibniz级数收敛速度极慢，⼈们想出很多办法来改进此算法。

arctan(x)=x−+−+⋯

3

x3

5

x5

7

x7

=(−1)k=0

∑∞

k2k+1

x2k+1

∵=

dx

d(arctanx)

1+x2

1

∵1−x=n(1−x)(1+x+x+2⋯+x)n−1

∴=1−x

1

x,∣x∣

∑∞

n1

=1+x2

1

=

1−(−x)2

1

(−x)=n=0

∑∞

2n(−1)xn=0

∑∞

n2n

∴arctan(x)=dx∫1+x2

1

=(−1)xdx∫n=0

∑∞

n2n

=(−1),∣x∣<1n=0

∑∞

n2n+1

x2n+1

4

π

x=1(TODO)∣x∣<1x=1

0×∞=0

=

4

π

=n=0

∑∞

2n+1

(−1)n

−

1

1

+

3

1

−

5

1

+

7

1

⋯.

π/4

π=4[4arctan(1/5)−arctan(1/239)]

arctan(x)=(−1)n=0

∑∞

n2n+1

x2n+1

x=1

=

4

π

4(1/5)−(1/239)n=0

∑∞

(2n+1)

(−1)n

(2n+12n+1)

(1)

上式是夏普在1699年发现的，计算π到72位⼩数位；

(2)

(3)

上式是英国伦敦⼈约翰马青JohnMachin于1706年发现的，计算π到101位⼩数位。经程序验证，收敛极快，精度提⾼也快。⼩数点位

数与迭代次数之⽐⼤约为1.4，即1000次迭代可以得到⼩数点后1400位。

(4)

(5)

欧拉于1755年发现，计算π到126位⼩数位。

(6)

计算π到136位⼩数位。

(7)

计算π到152位⼩数位。

(8)

计算π到200位⼩数位。

(9)

(10)

(11)

Gauss公式，收敛速度最快，⼩数点位数与迭代次数之⽐⼤约为2.5，即1000次迭代可以得到⼩数点后2500位。

NilakanthaSeries

Thisisthefasterconvergentmethodforpi.

这是计算圆周率的更快⼀点的收敛⽅法。

NilakantharieswithRationaliterationsforPAI

耗时(s),迭代次数,圆周率π

1.92592,4500,3.755

3.656277,6500,3.88833262423877236599992949

6.4071963,8500,3.93862988637989428921227038

9.2266,10000,3.95433134521397

59.0730983,20000,3.976

166.8244976,30000,3.978398941744860

Viete’sFormula

tan=−1

3

1

=6

π

(1−+−+...)3

1

3⋅3

1

5⋅32

1

7⋅33

1

tan1=−1tan+−1

2

1

tan=−1

3

1

(1−+−...)+2

1

3⋅22

1

5⋅24

1

(1−+−...)3

1

3⋅32

1

5⋅34

1

tan(1)=−14tan−−1

5

1

tan−1

239

1

tan(1)=−15tan+−1

7

1

2tan−1

79

3

tan(1)=−12tan+−1

3

1

tan−1

7

1

tan(1)=−14tan−−1

5

1

2tan+−1

408

1

tan−1

1393

1

tan(1)=−14tan−−1

5

1

tan+−1

70

1

tan−1

99

1

tan(1)=−1tan+−1

2

1

tan+−1

5

1

tan−1

8

1

tan(1)=−1tan+−1

2

1

tan−1

3

1

tan(1)=−12tan−−1

2

1

tan−1

7

1

tan(1)=−112tan+−1

18

1

8tan−−1

57

1

5tan−1

239

1

π=3+−+−+⋯

2⋅3⋅4

4

4⋅5⋅6

4

6⋅7⋅8

4

8⋅9⋅10

4

=3+−+−+⋯

1⋅3⋅2

1

2⋅5⋅3

1

3⋅7⋅4

1

4⋅9⋅5

1

=3+n=2

∑∞

(n−1)n(2n−1)

(−1)n

有关数学常数的⽆穷级数乘积形式的韦达公式:

发表于1593年，以FrançoisViète(1540–1603)的名字命名。

韦达公式可以⽤极限表⽰如下：

此处,初始值.

欧拉(LeonhardEuler)发现正弦函数的⼆倍⾓公式：

因为

且

所以

⽤替换：

然后,每⼀项⽤半⾓公式:

代⼊后也能得到韦达公式。

ItisalsopossibletoderivefromViète’sformulaarelatedformulaforthatstillinvolvesnestedsquarerootsoftwo,but

usonlyonemultiplication:

该算法效率最⾼，计算时间最短，精度也⾼。170次迭代报错，超出了迭代深度（RecursionError:maximumrecursiondepth

exceeded，缺省是5000）。

反三⾓函数

1676年，⽜顿给出的⽅法。令下列

π

=

π

2

⋅

2

2

⋅2

2+

2

⋯2

2+

2+

2

=n→∞

lim

i=1

∏n

2

ai

π

2

a=n

2+an−1a=1

2

sinx=2sincos=⋯=2sincos

2

x

2

x

n2n

x

k=1

∏n

2k

x

2sin=xn←∞

lim

n2n

x

=cos⋅cos⋅cos⋯

x

sinx

2

x

4

x

8

x

x=

2

π

=cos⋅cos⋅cos⋯

π

2

4

π

8

π

16

π

cos=

2

x

2

1+cosx

π

π=2k→∞

lim

k

ksquareroots

2−

2+

arcsin(x)

x=

2

1

(sinx)=

dx

d

−1⟹

1−x2

1

同理可以得到

Stirling’sApproximationforn!

所以对于很⼤的，Stirling斯特林逼近公式为

或

更精确的有:

Taylor级数和MaclaurinSeries的实现

在GeoGebra中可以轻松实现多项式之和逼近各种函数。

1.先定义次数Order的滑条

2.再定义要逼近的函数，如

3.调⽤函数

4.调⽤函数可以显⽰级数

sinx−1=(1−x)dx+c(constant)∫2−1/2

=(1++++⋯+c),∣x∣<1∫

2

x2

8

3x4

16

5x6

=x++++⋯+c

2

1

3

x3

2⋅4

1⋅3

5

x5

2⋅4⋅6

1⋅3⋅5

7

x7

∵sinx=−10atx=0⟹c=0

∴sin(x)=−1

x++

2

1

3

x3

+

2⋅4

1⋅3

5

x5

+

2⋅4⋅6

1⋅3⋅5

7

x7

⋯,∣x∣<1

=xn=0

∑∞

2(n!)(2n+1)2n2

(2n)!

2n+1

∴cos(x)=−1−

2

π

x++++⋯,∣x∣<(

2

1

3

x3

2⋅4

1⋅3

5

x5

2⋅4⋅6

1⋅3⋅5

7

x7)1

=−

2

π

xn=0

∑∞

2(n!)(2n+1)2n2

(2n)!

2n+1

lnn!=ln1+ln2+...+lnn

=lnkk=1

∑n

=lnxdx∫1

n

=(xlnx−x)∣1

n

=nlnn−n+1

≈nlnn−n

n

lnn!≈nlnn−nn!≈nen−n2πn

≤

2πn(

en)n

n!≤e

2πn(

en)n

12n

1

n=slider(1,20,1)

f(x)=arctan(x),g(x)=sin(x),h(x)=cos(x),⋯

g=TaylorPolinomial(f,x(A),n),A=(0,0)

FormulaText(g,true,true)

计算收敛速度极快算法

代数和⼏何结合的⽅法。

采⽤刘徽割圆术，⽤正边形逼近圆的⽅法，实现计算圆周率的⽬的。

记正边形的边长为,由正边形产⽣的正边形的边长为,则如图可有

单位圆的半径,则有正边形的边长,利⽤迭代关系式可以快速求出周长

即圆周率就是半周长，上述算法收敛性很快！

在GeoGebra中实现

1.先建⽴迭代次数滑动条

2.再建⽴迭代函数

3.然后⽤GGB的迭代命令，初始值取正6边形时的值1.

4.正边形的半周长为

这⼀算法的收敛也极快，正多边形的边数以指数级增长。边数

Python程序得到的结果

Python源码参见

采⽤符号运算得到前10个结论：

π

nπ

nLnn2nL2n

L=nAA,L=′

2nAB

L=2n

2

()+

2

Ln

2(1−OC)2

=L+

4

1

n

2

1−(

1−L

4

1

n

2)2

=2−2

1−L

4

1

n

2

L=2n

2−

4−Ln

2

r=16L=61

C=O2π=(n×n→∞

lim

L)n

π=×

2

1

(n×n→∞

lim

L)n

n=slider(1,10,1)

f(x)=sqrt(2−sqrt(4−xx))

value=Iteration(f,1,n)→L=2n,L=

2−

4−L×Lnn

61

2nπ=value×6×2n−1

=6×2n−1

12sqrt(2-sqrt(3))

24sqrt(2-sqrt(sqrt(3)+2))

48sqrt(2-sqrt(sqrt(sqrt(3)+2)+2))

96sqrt(2-sqrt(sqrt(sqrt(sqrt(3)+2)+2)+2))

192sqrt(2-sqrt(sqrt(sqrt(sqrt(sqrt(3)+2)+2)+2)+2))

384sqrt(2-sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(3)+2)+2)+2)+2)+2))

768sqrt(2-sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(3)+2)+2)+2)+2)+2)+2))

1536sqrt(2-sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(3)+2)+2)+2)+2)+2)+2)+2))

3072sqrt(2-sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(3)+2)+2)+2)+2)+2)+2)+2)+2))

6144sqrt(2-sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(sqrt(3)+2)+2)+2)+2)+2)+2)+2)+2)+2))

第151次到第159次的迭代结果：

3.97932384626433832795993758

856

3.97932384626433832795993758

765

3.97932384626433832795993758

992

3.97932384626433832795993758

799

3.97932384626433832795993758

251

3.97932384626433832795993758

864

3.97932384626433832795993758

767

3.97932384626433832795993758

993

3.97932384626433832795993758

799

第160次迭代结果，准确到⼩数点后第96位：

3.97932384626433832795993758

001

圆周率⼩数点后100位

97932384626433832795993751 592307816406

2862253421170679。

圆周率⽣成器

在⽹络上搜索到⼀个圆周率⽣成器，只要输⼊⼀个数：⽣成位数，就可以算出该圆周率。没有验证正确否。供⼤家参考。

在Python环境下，只要调⽤(n+1)就可以得到任意位⼩数的圆周率.

fromsympyimportpi

(1001)#得到1000位⼩数位的π

我输⼊了1000位数，计算结果如下：

π

nπ

3.97932384626433832795993758

0679

82211

294895493637867833460348610

454326648258725

9499

3956735973362441

394946395224737675238467481

4526356275778962146844125

892354247747798372978049951

45945534693526

528865875332776698738235378759375195778

959

有关迭代次数超限问题

在Python中，迭代次数超过时，会报错RecursionError:maximumrecursiondepthexceeded

改进办法：

importsys

#thetrecursionlimitfunctionis

#udtomodifythedefaultrecursion

#his,

#wecanincreatherecursionlimit

#tosatisfyourneeds

ursionlimit(10**6)

104