Ratio of circumference to diameter.

圆周率$\pi$，在几何学中定义为圆的周长与直径之比，在分析学中定义为满足$\sin x=0$的最小正实数，其数值约为$3.1415926$。

1988年3月14日，旧金山科学博物馆的物理学家Larry Shaw组织博物馆的员工围绕博物馆纪念碑$3\frac{1}{7}=\frac{22}{7}$圈，并一起吃水果派。之后每年旧金山科学博物馆会在当天举办庆祝活动。

2019年11月26日，联合国教科文组织宣布每年3月14日为“国际数学日”，旨在庆祝“数学在我们日常生活中的美丽与重要”。

2021年8月17日，瑞士研究人员使用超级计算机，历时108天将圆周率计算到小数点后$62.8$万亿位，创迄今最精确值记录。

借助计算机的力量，计算圆周率的精确近似已不再困难，对圆周率的数值计算逐渐成为检验计算机正确性和算力的测试范例。

本文回顾人类历史上对圆周率计算的主要尝试方法(包括几何法、分析法、概率法、计算机法)。至于计算圆周率的意义在哪里？正如《疑犯追踪(Person of Interest)》中Finch所说：

1. 几何法

圆周率$\pi$定义为圆的周长$C$与直径$d$之比$\pi=\frac{C}{d}$或圆的面积$S$与半径$r$平方之比$\pi=\frac{S}{r^2}$。若指定半径为$1$的单位圆，求得其周长或面积即可得到圆周率的值。

圆的周长或面积无法精确求得，可以通过内接或外切正多边形近似求得。由于内接(外切)多边形的周长和面积总是小于(大于)圆，此时可求得圆周率的下界(上界)。当正多边形的边数足够多时，近似结果也越精确。

⚪ 内接多边形（割圆术）

下面给出内接正$2^{n+1}$边形($n=1,2,3,…$)的边长$x_n$的递推公式，假设单位圆半径$r=1$。

正$2^{n+1}$边形的边长$x_n=AB$，则正$2^{n+2}$边形的边长$AD=x_{n+1}$计算为：

\[x_{n+1} = \sqrt{AC^2+CD^2} = \sqrt{\frac{AB^2}{4}+(OA-OC)^2} \\ = \sqrt{\frac{AB^2}{4}+(OA-\sqrt{OA^2-AC^2})^2} \\ = \sqrt{\frac{x_n^2}{4}+(1-\sqrt{1-\frac{x_n^2}{4}})^2}\]

则圆周率的估计为：

\[\pi=\frac{2^{n+1}x_n}{2} =2^{n}x_n\]

当$n=1$时，对应正方形的边长为$x_1=\sqrt{2}$。通过递推计算，可以得到圆周率$\pi$的近似值。

⚪ 外切多边形

下面给出外切正$2^{n+1}$边形($n=1,2,3,…$)的边长$x_n$的递推公式，假设单位圆半径$r=1$。

正$2^{n+1}$边形的边长$x_n=EF$，则正$2^{n+2}$边形的边长$BD=x_{n+1}$满足关系式：

\[BC^2+CE^2=BE^2\]

或写作：

\[BC^2+(OE-OC)^2=(AE-AB)^2\] \[(\frac{BD}{2})^2+(\sqrt{OA^2+AE^2}-OC)^2=(\frac{EF}{2}-\frac{BD}{2})^2\]

带入$x_n=EF$与$x_{n+1}=BD$得到：

\[(\frac{x_{n+1}}{2})^2+(\sqrt{1+\frac{x_n^2}{4}}-1)^2=(\frac{x_n}{2}-\frac{x_{n+1}}{2})^2\]

解上述方程可得正多边形的边长递推式：

\[x_{n+1} = \frac{2\sqrt{4+x_n^2}-4}{x_n}\]

则圆周率的估计为：

\[\pi=\frac{2^{n+1}x_n}{2} =2^{n}x_n\]

当$n=1$时，对应正方形的边长为$x_1=2$。通过递推计算，可以得到圆周率$\pi$的近似值。

⚪ 几何法的发展

2. 分析法

分析法的引入使得数学家构造了圆周率$\pi$的无穷级数表达式，包括求和式、乘积式、连分数等。由于圆周率$\pi$的级数表达式比割圆术的计算效率更高，因此这一时期$\pi$的计算精度迅速增加。

⚪ Machin’s Equation

英国数学家Machin提出如下公式，将圆周率计算到小数点后$100$位：

\[\frac{\pi}{4} = 4 \arctan\frac{1}{5}-\arctan\frac{1}{239}\]

该公式的构造性证明使用了$\tan x = \frac{1}{5}$的相关结论。已知三角函数公式：

\[\tan(\alpha \pm \beta) = \frac{\tan \alpha \pm \tan \beta}{1 \mp \tan \alpha \tan \beta}\]

则有：

\[\tan(2\arctan\frac{1}{5}) =\frac{\frac{1}{5}+\frac{1}{5}}{1-(\frac{1}{5})^2} = \frac{5}{12}\] \[\tan(4\arctan\frac{1}{5}) =\frac{\frac{5}{12}+\frac{5}{12}}{1-(\frac{5}{12})^2} = \frac{120}{119}\]

能够找到以下等式：

\[\tan(4\arctan\frac{1}{5}-\arctan\frac{1}{239}) = \frac{\frac{120}{119} - \frac{1}{239}}{1 + \frac{120}{119}\cdot \frac{1}{239}} = 1\]

由于$\tan \frac{\pi}{4}=1$，即得到Machin公式。至于公式中的反三角函数，可以由Taylor展开式近似计算。

不难发现，可以构造一系列满足$\tan x = 1$的表达式$x$，从而构造Machin类公式。

英国的William耗费15年时间，于1874年通过Machin公式计算到圆周率的小数点后$707$位，并将其刻在墓碑上作为一生的荣誉。可惜，后人发现他从第$528$位开始算错了。

1948年英国的Ferguson使用这类方法计算到小数点后$808$位，创下人工计算圆周率值的最高纪录。

⚪ Ramanujan’s Equation

印度数学家Ramanujan不加证明地给出了如下公式：

\[\frac{1}{\pi} = \frac{2\sqrt{2}}{99^2}\sum_{k=0}^{∞} \frac{(4k)!}{(k!)^4}\frac{1103+26390k}{396^{4k}}\]

值得一提的是仅取$k=0$时，即可得到$\pi=3.1415927$。该级数的收敛速度相当快，因此该公式是目前计算机求解圆周率上亿位有效数字的数学基础。

⚪ Chudnovsky’s Equation

Chudnovsky圆周率公式由Chudnovsky兄弟于1988年提出，可认为是Ramanujan圆周率公式的变体，计算时每多一项，计算精度提升约$14$个数量级。

\[\frac{1}{\pi} = \frac{1}{53360\sqrt{640320}}\sum_{k=0}^{∞}(-1)^k \frac{(6k)!}{(k!)^3(3k)!}\frac{13591409+545140134k}{640320^{3k}}\]

⚪ BBP’s Equation

BBP公式由Baily,Borwein,Plouffe于$1996$年提出，且利用该公式证明了在十六进制下可直接计算出$\pi$的小数点后第$n$位数，而不需要计算出前$n-1$位数。

\[\pi = \sum_{n=0}^{∞}\left( \frac{4}{8n+1}-\frac{2}{8n+4}-\frac{1}{8n+5}-\frac{1}{8n+6} \right)\frac{1}{16^n}\]

下面给出证明：

\[\begin{aligned} & \sum_{n=0}^{∞}\left( \frac{4}{8n+1}-\frac{2}{8n+4}-\frac{1}{8n+5}-\frac{1}{8n+6} \right)\frac{1}{16^n} \\ = & \sum_{n=0}^{∞}\left( \int_0^14x^{8n}dx-\int_0^12x^{8n+3}dx-\int_0^1x^{8n+4}dx-\int_0^1x^{8n+5}dx \right)\frac{1}{16^n}\\ = & \int_0^1 \sum_{n=0}^{∞} (\frac{x^8}{16})^n \left( 4-2x^3-x^4-x^5 \right)dx \\ = & \int_0^1 \frac{1}{1-\frac{x^8}{16}} \left( 4-2x^3-x^4-x^5 \right)dx \\ = & \int_0^1 \frac{16(1-x)(x^2+2)(x^2+2x+2)}{(2-x^2)(2+x^2)(4+x^4)} dx \\ = & \int_0^1 \frac{16(1-x)(x^2+2x+2)}{(2-x^2)(4+x^4)} dx= \int_0^1 \frac{16(1-x)}{(2-x^2)(x^2+2-2x)} dx \\ = & \int_0^1\left( \frac{4x}{x^2-2}+\frac{4-4x}{x^2+2-2x}+\frac{4}{x^2+2-2x} \right) dx \\ = & \left( 2\ln(2-x^2)-2\ln(x^2+2-2x)++4\arctan(x-1) \right) |_0^1 \\ = & \pi \end{aligned}\]

3. 概率法

⚪ 蒲丰投针 Buffon’s Needle

1777年，法国数学家蒲丰和勒克莱尔提出了投针问题：

由于向平面投针是随机的，所以用二维随机变量（$X,Y$）来确定它在桌上的具体位置。设$X$表示针的中点到平行线的距离，$Y$表示针与平行线的夹角，如果$X<\frac{l}{2}\sin Y$时，针与直线相交。

由于$X$服从均匀分布$[0,\frac{d}{2}]$，$Y$服从均匀分布$[0,\pi]$，$XY$相互独立，因此可以写出（$X,Y$）的概率密度函数：

\[f(X,Y)=\begin{cases} \frac{2}{\pi d}, & 0<X<\frac{d}{2},0<Y<\pi \\ 0, & \text{otherwise} \end{cases}\]

因此针与平行线橡胶的概率计算为：

\[\begin{aligned} p(X<\frac{l}{2}\sin Y)&=\iint_{x<\frac{l}{2}\sin y}f(x,y)dxdy \\ &= \int_0^\pi \int_0^{\frac{l}{2}\sin y}\frac{2}{\pi d}dxdy \\ &= \frac{2 l}{\pi d} \end{aligned}\]

从而得到圆周率的近似计算公式：

\[\pi \approx \frac{2l}{pd}\]

基于此公式，可用概率方法得到圆周率的近似值。将投针试验重复进行多次，并记下相交的次数，从而得到$p$的值，即可算出$π$的近似值。

⚪ 蒙特卡罗方法 Monte Carlo Method

Monte Carlo方法的基本思想是首先建立一个概率模型，使所求问题的解正好是该模型的参数或其他有关的特征量。然后通过模拟多次随机抽样试验，统计出某事件发生的百分比。只要试验次数很大，该百分比便近似于事件发生的概率。这实际上就是概率的统计定义。

比如在单位正方形内随机撒点，统计落入单位圆的概率：

\[\frac{\pi}{4} \approx \frac{\text{圆内点数}}{\text{总点数}}\]

4. 计算机法

随着计算机的发明，圆周率的求解有了突飞猛进的发展，下面罗列一些里程碑事件。