圆周率（$\pi$）的计算方法

慈悲心待人，时时有平安；智慧心安己，处处有自在！
From 翔文数学 © kumath@outlook.com
\\[1em]

一、圆周率介绍

---

代数方面

从代数上来讲，圆周率 $\pi$ 就是一个无理数（irrational numbers, 即无限不循环小数），不能表示为两个整数之比；也是一个超越数(transcendental numbers, 即不能作为某个有理系数多项式方程的根的数)，与三角函数有密切关系，并且与无穷级数的和有很大关系。

几何方面

从几何上来说，就是一个比值，特殊图形圆的周长与其直径的比值（ratio of the circumference of a circle to its diameter），也就是单位圆的半周长，还可以是单位圆的面积或圆面积与其半径的平方之比；也有人提出：如果用圆的周长与其半径的比值作为圆周率的话，也不是不可以的。

用希腊字母 $\pi$ 表示圆周率，它的几何定义就是 $\pi =\frac{圆周长}{圆直径}=\frac{\text{circumference}}{\text{diameter}}$

由上可得：圆面积 $A=\pi r^2;$ 圆周长 $C=2\pi r=\pi d;$

权威网站

$\pi$ 是一个数学常数，互联网上有很多关于圆周率的介绍，比较权威的参见 http://www.numberworld.org/y-cruncher/

现代计算机已经计算出 $\pi$ 的小数位数

有很多计算圆周率的方法，比较常用的有：

• 测量方法，量出圆周长和直径，用几何方法计算；
• 用多边形逼近方法计算；
• 用无限级数求和方法计算；
• 用蒙特卡罗法计算；
• …

常用的 $\pi$ 的近似值有 $3.14,3.1415926,\frac{22}{7},\frac{355}{113},3.14<\pi <\frac{22}{7}$。

南北朝时期的祖冲之不满足刘徽的计算结果，成功计算出了圆周率小数点后7位数字，确定了圆周率介于 $3.1415926<\pi <3.1415927$ 之间；得到了两个近似分数（有理数）约率 $\frac{22}{7}$ 和祖率（也称为密率） $\frac{355}{113}$。

2022年6月，Google 云团队 Emma Haruka Iwao 算出了 $100$ 兆位为 $0$ (100-trillion decimal place is 0)，运行程序是 y-cruncher 程序，采用的算法是 Chudnovsky 算法。运行了 $157$ 天。一共用了 $515$ TB 的存储空间。验证用了公式 Bailey-Borwein-Plouffe formula。

2021年8月19日，由瑞士格力孙应用技术大学 (University of Applied Sciences of the Grisons—UAS Grisons) 打破记录，历时仅 114 天 (4月28日~8月19日)，计算到了 $62.8$ 万亿位(Trillion digits), 精确位数为 $62^{\prime }831^{\prime }853^{\prime }071^{\prime }750$ 。 最后10位数字为 $7817924264$。

打破了 2020年1月29日由美国阿拉巴马州的 Timothy Mullican 耗时 8 个月（303天）创造的 $50$ 万亿位!

谷歌工程师 Emma Haruka Iwao 于2019年1月宣布已经计算到了 $31.4$ 万亿 $=31.4\times 10^{12}=3.14\times 10^{13}$ 位，为了纪念 圆周率日-PAI Day（每年的3月14日）。

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

泰勒级数理论（Taylor’s theorem and Taylor series）

多项式逼近任意函数

多项式 (polynomial) 逼近 (approximate) 未知函数 $f$ 在给定点附近的值，实际上就是找到多项式的系数 (coefficients): $a_0,a_1,\cdots ,a_n$

$$\begin{array}{c}{P}_n(x)=a_0+a_1(x-a)+a_2(x-a{)}^2+\cdots +a_n(x-a{)}^n+R((x-a{)}^{n+1})\end{array}$$

靠近给定点 $x=a$ 逼近函数 $f(x)$, 此处假定 (assuming) 函数 $f(x)$ 在 $x=a$ 有 $nth$ 微分(differentiable). 从而有 $f(x)\approx P_n(x)$, 下面试着找出系数 $a_n,n=0,1,2,\cdots ,n$.

>>> f=Function('f')
>>> series(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(Derivative(f(_x), (_x, 4)), _x, 0)/24 + x**5*Subs(Derivative(f(_x), (_x, 5)), _x, 0)/120 + O(x**6)

泰勒级数(Taylor series) 实函数或复函数(a real or complex-valued function) $f(x)$ 在 $x=a$ 有无限阶微分，则可以表示为幂级数 (power series):

$$\begin{array}{c}f(x)\approx f(a)+\frac{f^{\prime }(a)}{1!}(x-a)+\frac{f^{\prime \prime }(a)}{2!}(x-a{)}^2+\frac{f^{\prime \prime \prime }(a)}{3!}(x-a{)}^3+\cdots \end{array}$$

此处 $n!$ 表示$n$的阶乘 (denotes the factorial), $f^{(n)}(a)$ 表示 $f$ 在$x=a$ 处的 $nth$ 微分 (derivative). 可以写成更加紧凑的求和式子：

泰勒级数 $f(x)=\sum _{n=0}^{\infty }\frac{f^{(n)}(a)}{n!}(x-a{)}^n.$

$f$ 的0阶导数就是它自己 $f$，且$(x-a{)}^0=0!=1$. 当 $a=0$, 此时的泰勒级数又叫做 麦克劳林级数 Maclaurin series.

举例 (For instance)

指数函数(exponential function)

• 指数函数是重要的基本初等函数之一。一般地，$y=a^x$ 函数 $(a$ 为常数且 $a>0，a\ne 1)$ 叫做 指数函数。
• 指数函数的定义域是 $x\in \mathbb{R}$;
• 值域 $\text{(Domain)}(0,+\infty )$;
• 单调递增： $a>1$ 时, 单调递减： $0<a<1$ 时。

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

指数函数应用到值 $a=e$ 上的这个函数写为 $\exp (x)$。还可以等价的写为 $e^x$，这里的 $e$ 是数学常数，就是自然对数的底数，近似等于 $2.718281828$，还称为欧拉数.它的严格定义是 $e=\underset{n\to +\infty }{\lim }{\left(1+\frac{1}{n}\right)}^n$ 。

作为实数变量 $x$ 的函数， $e^x$ 的图像总是正的(在 $x$ 轴之上)并递增(从左向右看)。它永不触及 $x$ 轴，尽管它可以无限程度地靠近 $x$ 轴(所以，$x$ 轴是这个图像的水平渐近线。它的反函数是 自然对数 $\ln (x)$，其定义域是所有正数 $x>0$ 。

指数函数 $e^x$ 在 $x_0=0$ 处的泰勒级数展开式为：

$\begin{array}{rl}\sum _{n=0}^{\infty }\frac{x^n}{n!}& =\frac{x^0}{0!}+\frac{x^1}{1!}+\frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\frac{x^5}{5!}+\cdots \\ & =1+x+\frac{x^2}{2}+\frac{x^3}{6}+\frac{x^4}{24}+\frac{x^5}{120}+\cdots \end{array}$

• $e^x$ 对 $x$ 求导，始终不变，即 $\frac{{\mathbf{d}}^k}{\mathbf{d}{x}^k}{e}^x=e^x$，
• 且 $e^0=1$.
• 剩下的分子 (numerator) 为 $(x-0{)}^n$, 分母 (denominator) 为 $n!$.

欧拉公式 (Euler formula): $e^{ix}=\cos (x)+\text{i}\sin (x)$

>>> series(sin(x),n=10)
x - x**3/6 + x**5/120 - x**7/5040 + x**9/362880 + O(x**10)
>>> series(cos(x),n=10)
1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320 + O(x**10)
>>> series(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)

指数函数 $e^z$ 在 $z_0=0$ 的泰勒级数展开式为

$\begin{array}{rl}{e}^z& =\sum _{n=0}^{\infty }\frac{z^n}{n!}\\ & =\frac{z^0}{0!}+\frac{z^1}{1!}+\frac{z^2}{2!}+\frac{z^3}{3!}+\frac{z^4}{4!}+\frac{z^5}{5!}+\cdots \\ & =1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}+\frac{z^5}{120}+\cdots \end{array}$

令 $z=ix$(纯虚数), 此处 $i$ 虚数单位(imaginary unit), 满足 $i^{4k}=1,i^{4k+1}=i,i^{4k+2}=-1,i^{4k+3}=-i,k\in {\mathbb{N}}^{+}$.

右边为 $\begin{array}{r}1+ix-\frac{x^2}{2!}-i\frac{x^3}{3!}+\frac{x^4}{4!}+i\frac{x^5}{5!}-\frac{x^6}{6!}-i\frac{x^7}{7!}+\cdots \end{array}$

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

$$\begin{gathered} \cos (x)=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\cdots （偶函数） \\ \sin (x)=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\cdots （奇函数） \end{gathered}$$

上面也是函数 $\cos (x)$ 和 $\sin (x)$ 的泰勒级数展开式. (请参照 Colin Maclaurin series, $\forall x$)

右边可以简写为 $\cos (x)+\text{i}\sin (x)$, 简写为 $\text{cis}(x)$, 三个字母的缩写，产生一个新的函数。

故有 $\boxed{e^{\text{i }x}=\cos (x)+\text{i}\sin (x)}$

写成紧凑的求和式为：
$\begin{array}{rl}\cos (x)& =\sum _{n=0}^{\infty }\frac{(-1{)}^n{x}^{2n}}{(2n)!}\\ & =1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\cdots .\end{array}$

$\begin{array}{rl}\sin (x)& =\sum _{n=0}^{\infty }\frac{(-1{)}^n{x}^{2n+1}}{(2n+1)!}\\ & =x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\cdots .\end{array}$

可见，$\cos (x)$ 是偶函数，$\sin (x)$ 是奇函数。

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

$$\begin{gathered} {\sin }^{(2k)}(x)=\frac{d^{(2k)}\sin (x)}{dx}=(-1{)}^k\sin (x) \\ {\sin }^{(2k+1)}(x)=(-1{)}^k\cos (x) \end{gathered}$$

$\sin (x)$ 在 $x=0$ 处的偶次导数为 0, $\sin (x)$ 在 $x=0$ 奇次导数为 $(-1{)}^k$

欧拉恒等式 (Euler identity) 的由来

只要在欧拉公式 $\text{cis}(x)=\boxed{e^{\text{i }x}=\cos (x)+\text{i }\sin (x)}$ 中令 $x=\pi$, 即可得到 上帝恒等式——欧拉恒等式

$$\begin{array}{c}\boxed{e^{\text{i }\pi }=-1,e^{\text{i }\pi }+1=0}\end{array}$$

神奇的是：这个欧拉恒等式包含了 $e,i,\pi ,0,1$ 这五个神奇的数, 还有 $+,=$ 这两个基本运算符，故称为上帝公式。

多项式逼近理论应用 (Polynomial approximation theorem)

下面应用多项式逼近理论，得到几个常用的三角函数（$\sin (x),\cos (x),\tan (x)$）的乘积展开式。主要还是代数学中的求根方法。

问题: 求方程 $\sin (x)=0$ 的解。

1. 首先, $\sin (x)=0$ 有解 { k π ， k = 0 , ± 1 , ± 2 , ⋯   . } \{kπ，k=0,\pm 1,\pm 2, \cdots .\} {kπ，k=0,±1,±2,⋯.}

2. 假设 $\sin (x)$ 是多项式函数，由多项式有根的代数基本理论 (the fundamental theorem of algebra) 即 多项式逼近理论 (polynomial approximation theorem)

   $$\begin{gathered} \sin (x)=c\ast \prod _{k=1}^{\infty }(k\pi -x)(k\pi +x)x \\ \frac{\sin (x)}{x}=c\ast \prod _{k=1}^{\infty }(k\pi -x)(k\pi +x) \end{gathered}$$

   令 $x\to 0$， 求极限得到
   $c=\prod _{k=1}^{\infty }\frac{1}{k^2{\pi }^2}$

   我们得到：

   $\begin{array}{rl}\sin (x)& =x\prod _{k=1}^{\infty }(1-\frac{x}{k\pi })(1+\frac{x}{k\pi })\\ & =x\prod _{k=1}^{\infty }\left[1-\frac{x^2}{(k\pi {)}^2}\right]\end{array}$

令 $x=\frac{\pi }{2}$, 有

$\begin{array}{rl}\frac{\pi }{2}& =\prod _{k=1}^{\infty }\frac{2k}{2k-1}\cdot \frac{2k}{2k+1}\\ & =\prod _{k=1}^{\infty }\frac{1}{1-\frac{1}{4k^2}}\\ & =\prod _{k=1}^{\infty }\left[1+\frac{1}{4k^2-1}\right]\\ & =\frac{2}{1}\cdot \frac{2}{3}\cdot \frac{4}{3}\cdot \frac{4}{5}\cdot \frac{6}{5}\cdot \frac{6}{7}\cdots \\ & =\frac{4}{3}\cdot \frac{16}{15}\cdot \frac{36}{35}\cdot \frac{64}{63}\cdots \end{array}$

参见 John Wallis’ product for $\pi$

因为$\sin (\pi /2-x)=\cos (x)$, 所以有偶函数$\cos (x)$ 的无穷级数乘积公式：
$\begin{array}{rl}f{c}_1(x)=\cos (x)& \approx (\frac{\pi }{2}-x)\prod _{k=1}^{\infty }\left(1-\frac{\pi /2-x}{k\pi }\right)\left(1+\frac{\pi /2-x}{k\pi }\right)\\ & =(\frac{\pi }{2}-x)\prod _{k=1}^{\infty }\left(1-\frac{1}{2k}+\frac{x}{k\pi }\right)\left(1+\frac{1}{2k}-\frac{x}{k\pi }\right)\end{array}$

同理可得

$\begin{array}{rl}f{c}_2(x)=\cos (x)& \approx \left(1-\frac{2x}{\pi }\right)\prod _{k=1}^{\infty }\left(1+\frac{2x}{(2k-1)\pi }\right)\left(1-\frac{2x}{(2k+1)\pi }\right)\\ & =\left(1-\frac{2x}{\pi }\right)\prod _{k=1}^{\infty }\left(1+\frac{x/\pi }{k-1/2}\right)\left(1-\frac{x/\pi }{k+1/2}\right)\end{array}$

所有正整数平方的倒数和

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

>>> series(sin(x),n=10)
x - x**3/6 + x**5/120 - x**7/5040 + x**9/362880 + O(x**10)

1. 泰勒展开理论得到
   $$\begin{array}{rl}\sin (x)& =x-\frac{x^3}{3!}+\frac{x^5}{5!}+\cdots \end{array}$$

2. 多项式逼近理论得到
   $$\begin{array}{rl}\sin (x)& =x\prod _{k=1}^{\infty }(1-\frac{x^2}{(k\pi {)}^2})\end{array}$$

3. 比较 $x^3$ 的系数可得
   $$\begin{array}{rl}\frac{1}{3!}& =\sum _{k=1}^{\infty }\frac{1}{(k\pi {)}^2}\end{array}$$

   所以
   $\begin{array}{rl}\frac{{\pi }^2}{6}& =\sum _{k=1}^{\infty }\frac{1}{k^2}\\ & =1+\frac{1}{4}+\frac{1}{9}+\cdots +\frac{1}{n^2}+\cdots \end{array}$
   收敛速度还可以接受。

   from sympy import summation, Symbol, oo
   k = Symbol('k',integer=True)
   summation(1/k**2,(k,1,oo))
   

   计算
   $$\pi =\sqrt{6\sum _{k=1}^{\infty }\frac{1}{k^2}}$$

格雷戈里-莱布尼茨级数(Gregory-Leibniz series)

格雷戈里-莱布尼茨级数于 1671年2月15日提出。利用反正切函数的泰勒级数展开得到 $\frac{\pi }{4}$ 的级数算法。

>>> from sympy import *
>>> x,y,z=symbols('x y z')
>>> series(atan(x),n=10)
x - x**3/3 + x**5/5 - x**7/7 + x**9/9 + O(x**10)

反正切函数的级数 (series for the inverse tangent function), 也称为 Gregory’s series:
$\begin{array}{rl}\arctan (x)& =x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\cdots \\ & =\sum _{k=0}^{\infty }(-1{)}^k\frac{x^{2k+1}}{2k+1}\end{array}$

具体推导过程如下：
$\because \frac{\mathbf{d}(\arctan x)}{\mathbf{d}x}=\frac{1}{1+x^2}$
$$\begin{gathered} \because 1-x^n=(1-x)(1+x+x^2+\cdots +x^{n-1}) \\ \therefore \frac{1}{1-x}=\sum _{n=0}^{\infty }{x}^n,|x|<1 \end{gathered}$$
$\frac{1}{1+x^2}=\frac{1}{1-(-x^2)}=\sum _{n=0}^{\infty }(-x^2{)}^n=\sum _{n=0}^{\infty }(-1{)}^n{x}^{2n}$

$\begin{array}{rl}\therefore \arctan (x)& =\int \frac{1}{1+x^2}\mathrm{d}x\\ & =\int \sum _{n=0}^{\infty }(-1{)}^n{x}^{2n}\mathrm{d}x\\ & =\sum _{n=0}^{\infty }(-1{)}^n\frac{x^{2n+1}}{2n+1},|x|<1\end{array}$

莱布尼茨公式 （Leibniz formula）

$\frac{\pi }{4}$ 可以通过 $x=1$ $(TODO)$ 代入上述反正切函数的级数中得到.这里有一个疑问，既然 要求 $|x|<1$, 又怎么能取 $x=1$ 呢？这里存在一个 $0\times \infty \stackrel{?}{=}0$ 的问题。

$\frac{\pi }{4}=\sum _{n=0}^{\infty }\frac{(-1{)}^n}{2n+1}=\frac{1}{1}-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots .$

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

反正切函数的特殊用途

马青——John Machin (1680-1751) 英国伦敦人, 于 1706 年发现公式：
$\pi =4[4\arctan (1/5)-\arctan (1/239)]$

上式得益于 Gregory 有关反正切函数的幂级数展开式：
$\arctan (x)=\sum _{n=0}^{\infty }(-1{)}^n\frac{x^{2n+1}}{2n+1}$
令$x=1$，则有：
$\frac{\pi }{4}=\sum _{n=0}^{\infty }\frac{(-1{)}^n}{(2n+1)}\left(4(1/5{)}^{2n+1}-(1/239{)}^{2n+1}\right)$

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

(1) ${\tan }^{-1}\frac{1}{\sqrt{3}}=\frac{\pi }{6}=\frac{1}{\sqrt{3}}(1-\frac{1}{3\cdot 3}+\frac{1}{5\cdot 3^2}-\frac{1}{7\cdot 3^3}+...)$
上式是夏普在1699年发现的，计算π到72位小数位；

(2) ${\tan }^{-1}1={\tan }^{-1}\frac{1}{2}+{\tan }^{-1}\frac{1}{3}=\frac{1}{2}(1-\frac{1}{3\cdot 2^2}+\frac{1}{5\cdot 2^4}-...)+\frac{1}{3}(1-\frac{1}{3\cdot 3^2}+\frac{1}{5\cdot 3^4}-...)$

(3) ${\tan }^{-1}(1)=4{\tan }^{-1}\frac{1}{5}-{\tan }^{-1}\frac{1}{239}$
上式是英国伦敦人约翰马青 John Machin 于1706年发现的，计算π到101位小数位。经程序验证，收敛极快，精度提高也快。小数点位数与迭代次数之比大约为1.4，即1000次迭代可以得到小数点后1400位。

(4) ${\tan }^{-1}(1)=5{\tan }^{-1}\frac{1}{7}+2{\tan }^{-1}\frac{3}{79}$

(5) ${\tan }^{-1}(1)=2{\tan }^{-1}\frac{1}{3}+{\tan }^{-1}\frac{1}{7}$
欧拉于1755年发现，计算π到126位小数位。

(6) ${\tan }^{-1}(1)=4{\tan }^{-1}\frac{1}{5}-2{\tan }^{-1}\frac{1}{408}+{\tan }^{-1}\frac{1}{1393}$
计算π到136位小数位。

(7) ${\tan }^{-1}(1)=4{\tan }^{-1}\frac{1}{5}-{\tan }^{-1}\frac{1}{70}+{\tan }^{-1}\frac{1}{99}$
计算π到152位小数位。

(8) ${\tan }^{-1}(1)={\tan }^{-1}\frac{1}{2}+{\tan }^{-1}\frac{1}{5}+{\tan }^{-1}\frac{1}{8}$
计算π到200位小数位。

(9) ${\tan }^{-1}(1)={\tan }^{-1}\frac{1}{2}+{\tan }^{-1}\frac{1}{3}$

(10) ${\tan }^{-1}(1)=2{\tan }^{-1}\frac{1}{2}-{\tan }^{-1}\frac{1}{7}$

(11) ${\tan }^{-1}(1)=12{\tan }^{-1}\frac{1}{18}+8{\tan }^{-1}\frac{1}{57}-5{\tan }^{-1}\frac{1}{239}$
Gauss公式，收敛速度最快，小数点位数与迭代次数之比大约为2.5，即1000次迭代可以得到小数点后2500位。

Nilakantha Series

$\begin{array}{rl}\pi & =3+\frac{4}{2\cdot 3\cdot 4}-\frac{4}{4\cdot 5\cdot 6}+\frac{4}{6\cdot 7\cdot 8}-\frac{4}{8\cdot 9\cdot 10}+\cdots \\ & =3+\frac{1}{1\cdot 3\cdot 2}-\frac{1}{2\cdot 5\cdot 3}+\frac{1}{3\cdot 7\cdot 4}-\frac{1}{4\cdot 9\cdot 5}+\cdots \\ & =3+\sum _{n=2}^{\infty }\frac{(-1{)}^n}{(n-1)n(2n-1)}\end{array}$

This is the faster convergent method for pi.
这是计算圆周率的更快一点的收敛方法。

Nilakantha series with Rational iterations for PAI

耗时(s), 迭代次数, 圆周率 π
1.9259201999999997, 4500, 3.1415926535870515825837405507560137150803433030977
3.6562770999999996, 6500, 3.1415926535888833262423877236599992905011247088849
6.4071963, 8500, 3.1415926535893862988637989428921227071088158696038
9.226604299999998, 10000, 3.1415926535895433134507693207779502215449257421397
59.0730983, 20000, 3.1415926535897619931497723041779282156055633567492
166.8244976, 30000, 3.1415926535897839801292611829194198667214211744860

Viete’s Formula

有关数学常数 $\pi$ 的无穷级数乘积形式的韦达公式:
$\frac{2}{\pi }=\frac{\sqrt{2}}{2}\cdot \frac{\sqrt{2+\sqrt{2}}}{2}\cdot \frac{\sqrt{2+\sqrt{2+\sqrt{2}}}}{2}\cdots$
发表于1593年，以 François Viète (1540–1603)的名字命名。

韦达公式可以用极限表示如下：
$\underset{n\to \infty }{\lim }\prod _{i=1}^n\frac{a_i}{2}=\frac{2}{\pi }$
此处 $a_n=\sqrt{2+a_{n-1}}$, 初始值 $a_1=\sqrt{2}$.

欧拉 (Leonhard Euler) 发现正弦函数的二倍角公式：

因为 $\sin x=2\sin \frac{x}{2}\cos \frac{x}{2}=\cdots =2^n\sin \frac{x}{2^n}\prod _{k=1}^n\cos \frac{x}{2^k}$
且 $\underset{n\leftarrow \infty }{\lim }{2}^n\sin \frac{x}{2^n}=x$

所以 $\frac{\sin x}{x}=\cos \frac{x}{2}\cdot \cos \frac{x}{4}\cdot \cos \frac{x}{8}\cdots$

用 $x=\frac{\pi }{2}$ 替换：

$\frac{2}{\pi }=\cos \frac{\pi }{4}\cdot \cos \frac{\pi }{8}\cdot \cos \frac{\pi }{16}\cdots$

然后, 每一项用半角公式:
$\cos \frac{x}{2}=\sqrt{\frac{1+\cos x}{2}}$
代入后也能得到韦达公式。

It is also possible to derive from Viète’s formula a related formula for $\pi$ that still involves nested square roots of two, but uses only one multiplication: ref Viete formula

$\pi =\underset{k\to \infty }{\lim }{2}^k\underset{k\mathrm{square}\mathrm{roots}}{\underbrace{\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\cdots +\sqrt{2}}}}}}}}$

该算法效率最高，计算时间最短，精度也高。170次迭代报错，超出了迭代深度（RecursionError: maximum recursion depth exceeded，缺省是5000）。

反三角函数 $\arcsin (x)$

1676年，牛顿给出的方法。令下列 $x=\frac{1}{2}$
$\frac{d}{dx}({\sin }^{-1}x)=\frac{1}{\sqrt{1-x^2}}⟹$

$\begin{array}{rl}{\sin }^{-1}x& =\int (1-x^2{)}^{-1/2}\mathrm{d}x+c(constant)\\ & =\int (1+\frac{x^2}{2}+\frac{3x^4}{8}+\frac{5x^6}{16}+\cdots +c),|x|<1\\ & =x+\frac{1}{2}\frac{x^3}{3}+\frac{1\cdot 3}{2\cdot 4}\frac{x^5}{5}+\frac{1\cdot 3\cdot 5}{2\cdot 4\cdot 6}\frac{x^7}{7}+\cdots +c\end{array}$

$\because {\sin }^{-1}x=0\text{at}x=0⟹c=0$

$$\begin{gathered} \therefore {\sin }^{-1}(x)=x+\frac{1}{2}\frac{x^3}{3}+\frac{1\cdot 3}{2\cdot 4}\frac{x^5}{5}+\frac{1\cdot 3\cdot 5}{2\cdot 4\cdot 6}\frac{x^7}{7}+\cdots ,|x|<1 \\ =\sum _{n=0}^{\infty }\frac{(2n)!}{2^{2n}(n!{)}^2(2n+1)}{x}^{2n+1} \end{gathered}$$

同理可以得到
$$\begin{gathered} \therefore {\cos }^{-1}(x)=\frac{\pi }{2}-\left(x+\frac{1}{2}\frac{x^3}{3}+\frac{1\cdot 3}{2\cdot 4}\frac{x^5}{5}+\frac{1\cdot 3\cdot 5}{2\cdot 4\cdot 6}\frac{x^7}{7}+\cdots \right),|x|<1 \\ =\frac{\pi }{2}-\sum _{n=0}^{\infty }\frac{(2n)!}{2^{2n}(n!{)}^2(2n+1)}{x}^{2n+1} \end{gathered}$$

Stirling’s Approximation for n!

$\begin{array}{rl}\ln n!& =\ln 1+\ln 2+...+\ln n\\ & =\sum _{k=1}^n\ln k\\ & ={\int }_1^n\ln x\mathbf{d}x\\ & =(x\ln x-x){|}_1^n\\ & =n\ln n-n+1\\ & \approx n\ln n-n\end{array}$

所以 对于很大的 $n$，Stirling斯特林逼近公式为
$\boxed{\ln n!\approx n\ln n-n}$ 或 $\boxed{n!\approx n^n{e}^{-n}\sqrt{2\pi n}}$

更精确的有:
$\sqrt{2\pi n}{\left(\frac{n}{e}\right)}^n\le n!\le \sqrt{2\pi n}{\left(\frac{n}{e}\right)}^n{e}^{\frac{1}{12n}}$

Taylor级数和Maclaurin Series的实现

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

1. 先定义次数Order的滑条 $n=slider(1,20,1)$
2. 再定义要逼近的函数，如 $f(x)=\arctan (x),g(x)=\sin (x),h(x)=\cos (x),\cdots$
3. 调用函数 $g=TaylorPolinomial(f,x(A),n),A=(0,0)$
4. 调用函数 $FormulaText(g,true,true)$ 可以显示级数

TaylorSeries的多项式逼近演示GGB

计算 $\pi$ 收敛速度极快算法

代数和几何结合的方法。

采用刘徽割圆术，用正 $n$ 边形逼近圆的方法，实现计算圆周率 $\pi$ 的目的。

记 正 $n$ 边形的边长为 $L_n$, 由正 $n$ 边形产生的正 $2n$ 边形的边长为 $L_{2n}$, 则如图可有

$$\begin{gathered} L_n=A{A}^{\prime },L_{2n}=AB \\ {L}_{2n}^2=(\frac{L_n}{2}{)}^2+(1-OC{)}^2 \\ =\frac{1}{4}{L}_n^2+{\left(1-\sqrt{1-\frac{1}{4}{L}_n^2}\right)}^2 \\ =2-2\sqrt{1-\frac{1}{4}{L}_n^2} \\ {L}_{2n}=\sqrt{2-\sqrt{4-L_n^2}} \end{gathered}$$

单位圆的半径 $r=1$, 则有正 $6$ 边形的边长 $L_6=1$, 利用迭代关系式可以快速求出 周长
$$\begin{gathered} C_O=2\pi =\underset{n\to \infty }{\lim }(n\times L_n) \\ \pi =\frac{1}{2}\times \underset{n\to \infty }{\lim }(n\times L_n) \end{gathered}$$

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

在GeoGebra中实现

1. 先建立迭代次数滑动条 $n=slider(1,10,1)$
2. 再建立迭代函数 $f(x)=sqrt(2-sqrt(4-xx))$
3. 然后用GGB的迭代命令 $value=Iteration(f,1,n)\to$ $L_{2n}=\sqrt{2-\sqrt{4-L_n\times L_n}},L_6=1$，初始值取正6边形时的值1.
4. 正 $2n$ 边形的半周长为 $\pi =value\times 6\times 2^{n-1}$

这一算法的收敛也极快，正多边形的边数以指数级增长。边数 $=6\times 2^{n-1}$

Python程序得到的结果

Python 源码 参见迭代计算圆周率π的Python源码

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

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.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825324499856
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825337712765
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825341015992
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825341841799
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342048251
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342099864
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342112767
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342115993
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342116799

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

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117001

圆周率小数点后100位

14 15 92 65 35 89 79 32 38 46 26 43 38 32 79 50 28 84 19 71 69 39 93 75 10 58 20 97 49 44　59 23 07 81 64 06 28 62 08 99 86 28 03 48 25 34 21 17 06 79 。

圆周率 $\pi$ 生成器

在网络上搜索到一个圆周率生成器，只要输入一个数：生成位数，就可以算出该圆周率。没有验证正确否。供大家参考。

圆周率$\pi$生成器

在Python环境下，只要调用 sympy.pi.evalf(n+1) 就可以得到任意 $n$ 位小数的 圆周率 $\pi$.

from sympy import pi

pi.evalf(1001)  # 得到1000位小数位的π

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

3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679 821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201988

有关迭代次数超限问题

在Python中，迭代次数超过 $10^4$ 时，会报错 RecursionError: maximum recursion depth exceeded
改进办法：

import sys
 
# the setrecursionlimit function is
# used to modify the default recursion
# limit set by python. Using this,
# we can increase the recursion limit
# to satisfy our needs
 
sys.setrecursionlimit(10**6)