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日为“国际数学日”，旨在庆祝“数学在我们日常生活中的美丽与重要”。今天(2022.03.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$。该级数的收敛速度相当快，因此该公式是目前计算机求解圆周率上亿位有效数字的数学基础。

⚪ 计算机求解

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