隐藏目录
圆周率
级数法
迭代法
折半求和
自然对数底
级数法
对数常数
级数法
迭代法
Euler常数
级数法
其他常数

数学常数凝集着数学研究的精华,计算任意精度的数学常数是符号计算中一项重要的技术,许多涉及到微积分,特殊函数和函数求值的运算都会涉及到数学常数.特别地,如果运算结果关于这些数学常数是非线性的,那么数学常数的精度就显得尤为重要,下面主要以圆周率为例介绍一下计算数学常数时用到的主要技术.

回到顶部圆周率

人们对于计算圆周率$\pi$的兴趣似乎从远古时代就开始了,如果只从较为系统的研究看起,刘徽的割圆术大概是微积分被发现之前用来计算$\pi$的最佳技术,祖冲之利用它求出了著名的密率$355\over 113$,一共有8位有效数字.1706年Taylor的私人教师Machin利用反正切函数的Taylor级数展开式借助纸笔演算求出了$\pi$的小数点后100位数字,目前计算$\pi$的小数位数的世界纪录[1]保持者是日本科学家Yasumasa Kanada(金田康正),他于2002年在一台超级并行计算机上将$\pi$的小数位数推进到了惊人的12411亿位.

回到顶部级数法

回到顶部Machin型公式

John Machin用来计算$\pi$的公式是$${\pi\over 4}=4\tan^{-1}{1\over 5}-\tan^{-1}{1\over 239},$$利用$\tan^{-1}x$的Taylor展开式$$\tan^{-1}x=x-{x^3\over 3}+{x^5\over 5}+\cdots+(-1)^n{x^{2n+1}\over 2n+1}+O(x^{2n+3}),$$可以得到$${\pi\over 4}=4\sum_{k=0}^\infty\frac{(-1)^k}{(2k+1)5^{2k+1}}-\sum_{k=0}^\infty\frac{(-1)^k}{(2k+1)239^{2k+1}}.$$这个公式的优点是第二项收敛得很快,而第一项的分母中含有5的方幂,便于手工在十进制下计算.考虑恒等式$\tan({\pi\over 4}-x)=\frac{1-\tan{x}}{1+\tan{x}}$,两边同时取反正切得到$${\pi\over 4}=x+\tan^{-1}{\frac{1-\tan{x}}{1+\tan{x}}},$$再令$x=4\tan^{-1}{1\over 5}$就可以证明Machin公式.如果考察正切函数更复杂的恒等变换,可以得到一系列形如$${\pi\over 4}=\sum_{k=0}^{n-1}a_k\tan^{-1}{1\over b_k},\quad a_k,b_k\in\mathbb{N}_+,$$的公式,统称为Machin型公式,以$n=2$为例,其他三个2阶Machin型公式是 \begin{align*} {\pi\over 4}&=\tan^{-1}{1\over 2}+\tan^{-1}{1\over 3}\\ &=2\tan^{-1}{1\over 2}-\tan^{-1}{1\over 7}\\ &=2\tan^{-1}{1\over 3}+\tan^{-1}{1\over 7}.\\ \end{align*} 注意到$a_k,b_k\in\mathbb{N}_+$,把三角恒等式看成有理分式形式的递推公式,就可以将寻找Machin型公式的问题转化成寻找高阶不定方程整数解的问题,据此还可以证明$n=2$的Machin型公式只有以上四种,而$n$从1到21所有的Machin型公式共有1500种[2].不同Machin型公式的收敛速度可以用Lehmer数[3]$$\sum_{k=0}^{n-1}{1\over \log_{10}{b_k}}$$来衡量,Lemher数越小公式收敛得越快.目前已知收敛最快的Machin型公式是黄见利[4]发现的 \begin{align*} {\pi\over 4}&=183\tan^{-1}{1\over 239}+32\tan^{-1}{1\over 1023}-68\tan^{-1}{1\over 5832}\\ &+12\tan^{-1}{1\over 110443}-12\tan^{-1}{1\over 4841182}-100\tan^{-1}{1\over 6826318}, \end{align*} 对应的Lehmer数是1.51244.

与其他的计算方法相比,Machin型公式更适合于并行计算.目前$\pi$的小数位数世界纪录的保持者是日本科学家Yasumasa Kanada(金田康正),他于2002年在一台超级并行计算机上将$\pi$的小数位数推进到了惊人的12411亿位,计算时采用了两个$n=4$的Machin型公式 \begin{align*} {\pi\over 4}&=12\tan^{-1}{1\over 49}+32\tan^{-1}{1\over 57}-5\tan^{-1}{1\over 239}+12\tan^{-1}{1\over 110443}\\ &=44\tan^{-1}{1\over 57}+7\tan^{-1}{1\over 239}-12\tan^{-1}{1\over 682}+24\tan^{-1}{1\over 12943}. \end{align*}

回到顶部Ramanujan型公式

Ramanujan构造过一个$1\over\pi$的超几何级数展开式[5]$$\frac{1}{\pi}={{2\sqrt{2}\over 9801}\sum_{n=0}^\infty\frac{(4n)!}{(n!)^4}\cdot \frac{1103+26390n}{396^{4n}}},$$后来Chudnovsky兄弟[6]和Borwein兄弟[7]又分别给出形如$$\frac{1}{\pi}=12\sum_{n=0}^{\infty}\frac{(-1)^n(6n)!}{(3n)!(n!)^3}\cdot\frac{An+B}{C^{3n+\frac{3}{2}}}$$的三组公式,分别对应于$A=545140134,B=13591409,C=640320$, \begin{align*} A&=1657145277365+212175710912\sqrt{61},\\ B&=107578229802750+13773980892672\sqrt{61},\\ C&=5280(236674+30303\sqrt{61}) \end{align*} 和 \begin{align*} A &= 63365028312971999585426220 \\ &\quad + 28337702140800842046825600\sqrt{5} \\ &\quad + 384\sqrt{5} (10891728551171178200467436212395209160385656017 \\ &\qquad + 4870929086578810225077338534541688721351255040\sqrt{5})^{1/2},\\ B &= 7849910453496627210289749000 \\ &\quad + 3510586678260932028965606400\sqrt{5} \\ &\quad + 2515968\sqrt{3110}(6260208323789001636993322654444020882161 \\ &\qquad + 2799650273060444296577206890718825190235\sqrt{5})^{1/2},\\ C &= -214772995063512240 \\ &\quad - 96049403338648032\sqrt{5} \\ &\quad - 1296\sqrt{5}(10985234579463550323713318473 \\ &\qquad + 4912746253692362754607395912\sqrt{5})^{1/2}. \end{align*} 这四个公式统称为Ramanujan型公式,利用代数数论和二次域的知识还可以构造出更多这样的Ramanujan型公式[8],不同Ramanujan型公式的收敛速度可以用每计算一个级数项后结果所增加的十进制有效位数来衡量,以上这四个公式每向后计算一项,结果的有效位数分别大约增加8位,14位,31位和50位.

在个人计算机上计算$\pi$时,Ramanujan型公式是最常用的级数展开式,它可以利用折半求和的技术来计算(见算法2).目前个人计算机上$\pi$的小数位数世界纪录的保持者是由Steve Pagliarulo开发的QPI[9],它仅需数秒就能求出$\pi$的前100万位,计算时采用的就是Chudnovsky公式.

回到顶部BBP公式

最后值得一提的是Borwein兄弟和其合作者共同发现的BBP公式[10]$$\pi=\sum_{k=0}^{\infty} \frac{1}{16^k}(\frac{4}{8k+1}-\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6}),$$注意到通项前的系数是$1\over 16^k$,因此可以利用它直接求出$\pi$在十六进制下某一个指定位上的数字,而不必先求出所有之前位上的数字.在此基础上Fabrice Bellard又提出了一个公式[11] \begin{align*} \pi&={1\over 2^6}\sum_{n=0}^\infty{(-1)^n\over 2^{10n}}(-{2^5\over 4n+1}-{1\over 4n+3}+{2^8\over 10n+1}\\ &-{2^6\over 10n+3}-{2^2\over 10n+5}-{2^2\over 10n+7}+{1\over 10n+9}), \end{align*} 它比BBP公式要快43%,他还指出如果利用$$-\ln(1-x)=\sum_{n=1}^\infty{x^n\over n}$$以及$\tan^{-1}x$是$-\ln(1-ix)$的虚部,可以从Machin型公式构造出更多这样类似于BBP公式的公式.

回到顶部迭代法

顾名思义,迭代法就是利用递推公式来计算$\pi$,通过迭代逐步提高精度.将递推公式设法改写成$\pi_{n+1}=f(\pi_n)$的形式,则存在正整数$m$使得$|\pi_{n+1}-\pi|=|f(\pi_n)-\pi|\sim|\pi_n-\pi|^m$,于是每迭代一步后的误差将变为原来误差的$m$次幂,如果换一个角度来看,这就是说有效位数将变为原来有效位数的$m$倍,因此可以用$m$的大小作为衡量迭代法性能的指标.

回到顶部代数几何平均值

首先来看Brent-Salamin算法[12],这是一个二阶收敛的迭代算法,它用到了Gauss-Legendre的AGM(代数几何平均值)[13]迭代.

定义1(代数几何平均值) 设$a,b\in \mathbb{R}^+$,$a>b$.令$a_0=a$,$b_0=b$,$a_{n+1}=\frac{1}{2}(a_n+b_n)$,$b_{n+1}=\sqrt{a_nb_n}$,则$a$,$b$的代数几何平均值定义为$$AGM(a,b)=\lim\limits_{n\rightarrow\infty}a_n=\lim\limits_{n\rightarrow\infty}b_n.$$

如果将代数平均值或几何平均值换成更高阶的平均值,例如$$({a_n^3b_n+b_n^3a_n\over 2})^{1/4},$$就可以得到更高阶的迭代结构.

回到顶部完全椭圆积分

为了利用AGM来计算圆周率,还需要用到完全椭圆积分的概念.

定义2(完全椭圆积分) 形如$\int R(x,y)dx$,其中$R$是$x,y$的有理函数,$y^2=P(x)$是$x$的三次或四次多项式的积分称为椭圆积分. \begin{align*} K(k)&=\int_0^1\frac{dx}{\sqrt{(1-x^2)(1-k^2x^2)}}=\int_0^{\pi\over 2}\frac{d\phi}{\sqrt{1-k^2sin^2\phi}},\\ E(k)&=\int_0^1\sqrt{\frac{1-k^2x^2}{1-x^2}}\,dx=\int_0^{\pi\over 2}\sqrt{1-k^2sin^2\phi}\,d\phi, \end{align*} 这两个积分分别称为第一类和第二类完全椭圆积分,其中$|k|<1$.

记$I(a,b)={1\over a}K(\sqrt{1-{b^2\over a^2}})$,那么 \begin{align*} I(a,b)&=\int_0^{\pi\over 2}\frac{d\phi}{\sqrt{a^2-(a^2-b^2)\sin^2\phi}}\\ &=\int_0^{\pi\over 2}\frac{d\phi}{\sqrt{a^2\cos^2\phi+b^2\sin^2\phi}}. \end{align*}

通过变量替换可以得出$I(a,b)=I({a+b\over 2},\sqrt{ab})$,也可以写成$I(a_0,b_0)=I(a_n,b_n)=I(a_{n+1},b_{n+1})$,两边取极限得到 \begin{align*} I(a,b)&=I(AGM(a,b),AGM(a,b))\\ &={1\over AGM(a,b)}K(0)={\pi\over 2AGM(a,b)}. \end{align*} 注意到$I(1,{1\over\sqrt{2}})=K({1\over\sqrt{2}})$,因此$$K({1\over\sqrt{2}})={\pi\over 2AGM(1,{1\over\sqrt{2}})}.$$

类似的,记$J(a,b)=aE(\sqrt{1-{b^2\over a^2}})$,那么$$J(a,b)=\int_0^{\frac{\pi}{2}}\sqrt{a^2\cos^2\phi+b^2\sin^2\phi}\,d\phi.$$令$a=1,b=\cos\phi$,于是$K(\sin\phi)={\pi\over 2AGM(a,b)}$,再令$c_0=\sin\phi,c_{n+1}=a_n-a_{n+1}$,可以证明$$\sum_{n=0}^\infty2^{n-1}c_n^2=1-{E(\sin\phi)\over K(\sin\phi)}.$$因此$$E({1\over\sqrt{2}})=K({1\over\sqrt{2}})(1-\sum_{n=0}^\infty2^{n-1}c_n^2).$$

若$\phi+\psi={\pi\over 2}$,那么第一类与第二类完全椭圆积分之间存在勒让德关系$$K(\sin\phi)E(\sin\psi)+E(\sin\phi)K(\sin\psi)-K(\sin\phi)K(\sin\psi)={\pi\over 2}.$$令$\phi=\psi={\pi\over 4}$就可以得到$\pi=4K({1\over\sqrt{2}})E({1\over\sqrt{2}})-2K({1\over\sqrt{2}})^2$,即$$\pi=K^2(2-4\sum_{n=0}^\infty2^{n-1}c_n^2),$$也可以写成$$\pi=\frac{(2AGM(1,{1\over\sqrt{2}})^2}{2-4\sum_{n=0}^\infty2^{n-1}c_n^2}.$$

回到顶部Brent-Salamin算法

现在可以写出Brent-Salamin算法的详细过程了.

算法1(Brent-Salamin算法)

输入:迭代次数$m$.

输出:迭代$m$次后得到的$\pi$的近似值.

  1. 令$a_0=1,b_0={1\over\sqrt{2}},t_0={1\over 4},p_0=1$.
  2. $n$从0到$m-1$,顺次计算$a_{n+1}={a_n+b_n\over 2}$,$b_{n+1}=\sqrt{a_nb_n}$,$t_{n+1}=t_n-p_n(a_n-a_{n+1})^2$,$p_{n+1}=2p_n$.
  3. 返回$\pi_m={(a_m+b_m)^2\over 4t_m}$.
使用Brent-Salamin算法计算时,前三次迭代可以得到近似值 \begin{align*} 3.14,\\ 3.1415926,\\ 3.14159265358979.\\ \end{align*}

回到顶部高阶迭代公式

除了利用AGM之外,人们还发现了许多其他的高阶迭代公式,例如下面的三个著名的迭代公式[14]就是由Borwein兄弟发现的.

定理1(Borwein算法)

  1. 令$x_0=\sqrt{2}$,$y_0=\sqrt[ 4]{2}$,$p_0=2+\sqrt{2}$.递推公式为 \begin{align*} x_{n+1}&={1\over 2}(x_n^{1\over 2}+x_n^{-{1\over 2}}),\\ y_{n+1}&=\frac{x_n^{1\over 2}y_n+x_n^{-{1\over 2}}}{y_n+1},\\ p_{n+1}&=p_n{x_{n+1}+1\over y_{n+1}+1}. \end{align*} 那么 \begin{align*} |p_{n+1}-\pi|&<{2^{-n-1}\over\pi^2}|p_n-\pi|^2,\\ |p_n-\pi|&<\frac{\pi^22^{n+4}e^{-2^{n+1}\pi}}{AGM(1,{1\over\sqrt{2}})^2}\\ &<10^{-2^n}. \end{align*}
  2. 令$p_0=6-4\sqrt{2}$,$x_0=\sqrt{2}-1$.递推公式为 \begin{align*} x_{n+1}&=\frac{1-(1-x_n^4)^{1\over 4}}{1+(1-x_n^4)^{1\over 4}},\\ p_{n+1}&=p_n(1+x_{n+1})^4-2^{2n+3}x_{n+1}(1+x_{n+1}+x_{n+1}^2). \end{align*} 那么$|p_n-{1\over\pi}|<4^{n+2}e^{-2^{2n+1}\pi}$.
  3. 令$a_0={1\over 3}$,$r_0={\sqrt{3}-1\over 2}$,$s_0=(1-r_0^3)^{1\over 3}$.递推公式为 \begin{align*} t_{n+1} &= 1 + 2r_n,\\ u_{n+1} &= (9r_n (1 + r_n + r_n^2))^{1\over 3},\\ v_{n+1} &= t_{n+1}^2 + t_{n+1}u_{n+1} + u_{n+1}^2,\\ w_{n+1} &= \frac{27 (1 + s_n + s_n^2)}{v_{n+1}},\\ a_{n+1} &= w_{n+1}a_n + 3^{2n-1}(1-w_{n+1}),\\ s_{n+1} &= \frac{(1 - r_n)^3}{(t_{n+1} + 2u_{n+1})v_{n+1}},\\ r_{n+1} &= (1 - s_{n+1}^3)^{1\over 3}. \end{align*} 那么$|a_n-{1\over \pi}|$9阶收敛.

回到顶部Beeler公式

最后值得一提的是Beeler利用不动点定理和反正切函数的Taylor展开式给出的一系列收敛到$\pi$的递推公式[15] \begin{align*} p_{n+1}&=p_n-\tan p_n,\\ p_{n+1}&=p_n-\tan p_n+{1\over 3}\tan^3p_n,\\ p_{n+1}&=p_n-\tan p_n+{1\over 3}\tan^3p_n-{1\over 5}\tan^5p_n,\\ \cdots, \end{align*} 其中初始值$p_0$应该取为$\pi$的一个普通的近似值,例如${22\over 7},{355\over 113}$等.

回到顶部折半求和

折半求和[16]是用来计算超几何级数在有理点处的值的一项重要技术,在前面介绍级数法计算$\pi$时已经提到过超几何级数,简单地来看,超几何级数就是一类特殊的幂级数.

定义3(超几何级数) 称形如$$_pF_q(a_1,\cdots,a_p;b_1,\cdots,b_q;z)=\sum\limits_{n=0}^{\infty}\frac{(a_1)_n\cdots(a_p)_nz^n}{(b_1)_n\cdots(b_q)_nn!}$$的级数为超几何级数,其中$(a)_n,(b)_n$代表降阶乘(Falling Factorial),其定义为$$(x)_n=x(x-1)\cdots(x-n+1)=\frac{x!}{(x-n)!}.$$

许多函数都具有超几何级数的展开式,最简单的例子是$e^z={_0F_0}(z)$,除此之外,初等函数一般都具有$_2F_1$的展开式,高斯超几何级数定义为 \begin{align*} F(a,b,c,z)&={_2F_1}(a+1,b+1;c+1;z)\\ &=1+\frac{ab}{c}\cdot\frac{z}{1!}+\frac{a(a+1)b(b+1)}{c(c+1)}\cdot\frac{z^2}{2!}+\cdots, \end{align*} 那么 \begin{align*} \ln(1+z)&=zF(1,1,2,-z),\\ \tan^{-1}z&=zF(\frac{1}{2},1,\frac{3}{2},-z^2),\\ \tanh^{-1}z&=zF({1\over 2},1,{3\over 2},z^2),\\ \frac{1}{(1-z)^a}&=F(a,1,1,z). \end{align*} 高斯超几何级数的下面两条性质也是经常被用到的.
  • $F(a,b,c,z)$是关于$\omega$的微分方程$$z(1-z)\frac{d^2\omega}{dz^2}+(c-(a+b+1)z)\frac{d\omega}{dz}-ab\omega=0$$的解.
  • $F(a,b,c,z)=(1-z)^{-b}F(c-a,b,c,\frac{z}{z-1})$.

许多数学常数可以看作超几何级数在某个有理点处的值,一般来说,超几何级数都可以改写成如下形式$$S=\sum_{n=0}^{\infty}\frac{a(n)}{b(n)}\frac{p(0)\cdots p(n)}{q(0)\cdots q(n)},$$其中$a(n),b(n),p(n),q(n)$都是整系数多项式.以Chudnovsky公式为例,先将通项改写成$$s(n)=\frac{(-1)^n(An+B)}{1}\cdot \frac{(6n)!/(3n)!}{(n!)^3C^{3n}},$$于是 \begin{align*} a(n)&=(-1)^n(An+B),\\ b(n)&=1,\\ p(0)\cdots p(n)&={(6n)!\over (3n)!},\\ q(0)\cdots q(n)&=(n!)^3C^{3n} \end{align*} 相邻两项相除就可以解出$p(n)=(6n-1)(6n-3)(6n-5)$,$q(n)=n^3C^3$.

考虑部分和$$S=S(n_1,n_2)=\sum_{n_1\le n<n_2}\frac{a(n)}{b(n)}\frac{p(n_1)\cdots p(n)}{q(n_1)\cdots q(n)}.$$记 \begin{align*} P&=P(n_1,n_2)=p(n_1)\cdots p(n_2-1),\\ Q&=Q(n_1,n_2)=q(n_1)\cdots q(n_2-1),\\ B&=B(n_1,n_2)=b(n_1)\cdots b(n_2-1), \end{align*} 再设$T=T(n_1,n_2)=BQS$,那么$S=\frac{T}{BQ}$.折半求和使用了一种类似于通分的技术构造出递推公式来计算$P,Q,B,T$,设$n_m$满足$n_1<n_m<n_2$,将区间$n_1\le x<n_2$分成两段$n_1\le x<n_m$及$n_m\le x<n_2$,分别用下标$l,r$来表示对应这两个子区间的$P,Q,B,T$的值,例如$P_l=P(n_1,n_m),P_r=P(n_m,n_2)$,那么存在递推公式 \begin{align*} P&=P_l\cdot P_r,\\ Q&=Q_l\cdot Q_r,\\ B&=B_l\cdot B_r,\\ T&=T_l\cdot B_rQ_r+T_r\cdot B_lP_l. \end{align*}

现在可以写出折半求和算法的详细过程了.

算法2(折半求和)

输入:通项公式$a(n),b(n),p(n),q(n)$,区间$[n_1,n_2)$.

输出:部分和$P(n_1,n_2),Q(n_1,n_2),B(n_1,n_2),T(n_1,n_2)$.

  1. 如果$n_2-n_1<5$,直接根据定义计算出$(P,Q,B,T)$并返回.
  2. 设$n_m=\left[\frac{n_1+n_2}{2}\right]$,$(P_l,Q_l,B_l,T_l)$,$(P_r,Q_r,B_r,T_r)$分别是对应于子区间$[n_1,n_m)$和$[n_m,n_2)$上的$(P,Q,B,T)$,在两个子区间上分别应用折半求和计算出它们.
  3. 计算$P=P_lP_r$,$Q=Q_lQ_r$,$B=B_lB_r$和$T=T_lB_rQ_r+T_rB_lP_l$,返回$(P,Q,B,T)$.
使用折半求和计算出$[n_1,n_2)$上的部分和$(P,Q,B,T)$后,再利用$S={T\over BQ}$做一次除法就可以求出原级数的部分和$S$.注意到最终结果$S$只依赖于$T$和$BQ$的比值,因此在折半求和计算$(P,Q,B,T)$的过程中,最后一步应用递推公式之前可以将$P_l,Q_r$先"约分",即将它们同时除以其最大公因子$\gcd(P_l,Q_r)$,这样一来由递推公式求出的$P,Q,T$都变为原来的$1\over\gcd(P_l,Q_r)$倍,而比值$T\over BQ$则保持不变.为了快速地计算出$\gcd(P_l,Q_r)$,常用的技术是在计算过程中保留$P,Q$的部分素因子分解式,并在递推相乘时同步更新其部分素因子分解式.

从整体上来看,折半求和并没有减少运算次数,它只是让参与乘法的两个数的数量级更接近,这样就可以充分发挥Karatsuba,Toom-Cook和FFT等大整数快速乘法算法的优势.如果大整数乘法采用古典乘法算法,那么折半求和是不起作用的,甚至可能比普通的级数求和方法更慢.除此之外,折半求和实际上是将截断的级数先乘上分母的最小公倍数后再计算,这样可以避免对中间结果进行除法等不精确计算,而只需要在最后做一次除法.可以证明,采用折半求和算法计算出部分和$S$的前$N$位有效数字大约需要$O((\log N)^2M(N))$次基本运算,其中$M(N)$代表两个$N$位整数相乘所需要的基本运算的次数(如果采用有限域上的FFT乘法,那么$M(N)=O(N\log N\log\log N)$).折半求和利用两个独立的子区间的数据来合成整个区间的数据,因此也很适合于并行计算.

回到顶部自然对数底

回到顶部级数法

自然对数底常记为$e$,它的定义为$$e=\lim_{n\rightarrow\infty}e_n=\lim_{n\rightarrow\infty}(1+{1\over n})^n.$$根据Newton二项式公式可以得到 \begin{align*} e_n&=(1+{1\over n})^n=\sum_{k=0}^nC_n^k({1\over n})^k\\ &=\sum_{k=0}^n\frac{n(n-1)\cdots(n-k+1)}{k!}\cdot{1\over n^k}\\ &=\sum_{k=0}^n{1\over k!}(1-{1\over n})\cdots(1-{k-1\over n}), \end{align*} 记$$s_n=\sum_{k=0}^n{1\over k!},$$那么$e_n<s_n\le e$,即$$e=\lim_{n\rightarrow\infty}s_n=\lim_{n\rightarrow\infty} 1+\frac{1}{1!}+\frac{1}{2!}+\cdots+\frac{1}{n!}.$$进一步还可以推出$|e-s_n|<{1\over nn!}$,和计算圆周率$\pi$用到的很多级数的误差项不同,这个误差项很小,它关于$n$是指数级收敛的,因此我们可以直接利用级数$s_n$来计算自然对数底$e$.除此之外,另一个级数$$e^{-1}=\lim_{n\rightarrow\infty}1-{1\over 1!}+{1\over 2!}+\cdots+(-1)^n{1\over n!}$$也可以求出$e$,它常常被用来验证第一个级数的计算结果,这两个级数都可以使用折半求和方法来计算.通过级数$s_n$的计算公式还可以得到关于自然对数底$e$和圆周率$\pi$一个有趣的递推公式,令$u_0=v_0=0$,$u_1=v_1=1$,递推公式为 \begin{align*} u_{n+2}&=u_{n+1}+{1\over n}u_n,\\ v_{n+2}&={1\over n}v_{n+1}+v_n, \end{align*} 那么 \begin{align*} \lim_{n\rightarrow\infty}{n\over u_n}&=e,\\ \lim_{n\rightarrow\infty}{2n\over v_n^2}&=\pi. \end{align*}

回到顶部对数常数

回到顶部级数法

对数常数指的是$\ln 2$,在计算机被发明之前,科学家和工程师们只能借助于对数表和对数尺之类的工具来完成繁杂的乘除运算,注意到$$\ln x=\ln{x\over 2^n}+n\ln 2,$$可以适当的选取$n$使得$0<{x\over 2^n}<1$,这样就需要计算$\ln 2$来完成大数与标准$(0,1)$区间之间的转换.利用$\ln(1+x)$的Taylor展开式$$\ln(1+x)=x-{x^2\over 2}+{x^3\over 3}-\cdots,$$令$x=1$就可以得到$\ln 2$的形式化定义$$\ln 2=\lim_{n\rightarrow\infty}1-{1\over 2}+{1\over 3}+\cdots+(-1)^{n+1}{1\over n}.$$如果令$x=-{1\over 2}$,那么有$$\ln 2=-\ln{1\over 2}=\sum_{k=1}^\infty{1\over k2^k},$$进一步地,利用恒等式$\ln 2=-{1\over 2}\ln(1-{1\over 4})+\tanh^{-1}{1\over 2}$,还可以得到类似于计算圆周率时提到过的BBP公式的两个公式 \begin{align*} \ln 2&=\sum_{k=0}^\infty{1\over 4^k}({1\over 4k+2}+{1\over 8k+8}),\\ \ln 2&={2\over 3}+\sum_{k=1}^\infty{1\over 16^k}({1\over 16k+12}+{1\over 8k+4}+{1\over 4k+1}+{1\over 2k}). \end{align*}

与计算圆周率$\pi$类似,计算对数常数$\ln 2$也有相应的Machin型公式.首先来看反正切函数的Taylor展开式$$\tanh^{-1}x={1\over 2}\ln{1+x\over 1-x}=\sum_{k=0}^\infty{x^{2k+1}\over 2k+1},$$令$x={1\over 3}$,那么$$\ln 2=2\tanh^{-1}{1\over 3}={2\over 3}\sum_{k=0}^\infty{1\over (2k+1)9^k}.$$如果考虑类似于$\tanh^{-1}{1\over x}=\tanh^{-1}{1\over 2x-1}+\tanh^{-1}{1\over 2x+1}$的初等变换,还可以得到一系列$n=2$的Machin型公式,其中最有名的是Euler用来计算$\ln2$时用到的$$\ln 2=2\tanh^{-1}{1\over 5}+2\tanh^{-1}{1\over 7}.$$除此之外,人们最近还发现了一些高阶的Machin型公式,例如 \begin{align*} \ln 2 &=18\coth^{-1}26-2\coth^{-1}4801+8\coth^{-1}8749,\\ \ln 2&=144\coth^{-1}251+54\coth^{-1}449-38\coth^{-1}4801+62\coth^{-1}8749,\\ \ln 2&=72\coth^{-1}127+54\coth^{-1}449+34\coth^{-1}4801-10\coth^{-1}8749. \end{align*}

前面介绍超几何级数时提到过 \begin{align*} \tanh^{-1}x&=xF({1\over 2},1,{3\over 2},x^2)\\ &={x\over 1-x^2}F(1,1,{3\over 2},{x^2\over x^2-1}), \end{align*} 于是有$$\ln 2=2\tanh^{-1}{1\over 3}={3\over 4}(1-{1\over 4}\cdot{1\over 3}+{1\over 4^2}\cdot{1\cdot 2\over 3\cdot 5}-{1\over 4^3}\cdot{1\cdot 2\cdot 3\over 3\cdot 5\cdot 7}+\cdots),$$这个公式结构简单,很容易利用折半求和方法来计算它.

回到顶部迭代法

代数几何平均值(AGM)迭代同样适用于计算对数常数$\ln 2$[12],记$${1\over R(a,b)}=1-\sum_{n=0}^\infty2^{n-1}(a_n^2-b_n^2),$$ 可以看出$R(a,b)$实际上就是$\pi\over 2AGM(a,b)$.那么$$|\ln x-R(1,10^{-N})+R(1,10^{-N}x)|\le {N\over 10^{2(N-2)}}.$$取$R(1,10^{-N})-R(1,2\cdot 10^{-N})$作为$\ln 2$的近似值,依然借用计算$AGM(a,b)$时的迭代结构,那么在计算$AGM(a,b)$的同时就可以计算出$a_n,b_n$,进而可以计算出$R(a,b)$,这样就得到了计算$\ln 2$的一个二阶收敛的迭代算法.

回到顶部Euler常数

回到顶部级数法

Euler常数常记为$\gamma$,它的定义为$$\gamma=\lim_{n\rightarrow \infty} (1+\frac{1}{2}+\cdots+\frac{1}{n}-\ln n)=\lim_{n\rightarrow\infty}(H_n-\ln n).$$如果直接按照这个定义来计算$\gamma$将会发现它收敛得太慢了,利用Euler-Maclaurin积分和Bernoulli数的定义可以得到$\gamma$一个更好的级数展开式$$\gamma=H_n-\ln n-{1\over 2n}+\sum_{k=1}^\infty{B_{2k}\over 2k}\cdot{1\over n^{2k}},\quad B_{2k}\text{为Bernoulli数},$$它的误差项是$$\epsilon_k,n={B_{2k+2}\over(2k+2)n^{2k+2}}\approx{2(2k+2)!\over(2k+2)(2\pi n)^{2k+2}}.$$

将$\gamma=-\Gamma'(1)$分部积分后可以得到 \begin{align*} \gamma+\ln n&=I_n-R_n,\\ I_n&=\int_0^n{1-e^{-t}\over t}\,dt,\\ R_n&=\int_n^\infty{e^{-t}\over t}\,dt. \end{align*} 注意到$R_n=O(e^{-n})$,利用$1-e^{-t}\over t$的级数展开式可以得到$$I_n=\sum_{k=1}^\infty(-1)^{k-1}{n^k\over kk!},$$于是$$\gamma=\sum_{k=1}^{\alpha n}(-1)^{k-1}{n^k\over kk!}-\ln n+O(e^{-n}),\quad \alpha(\ln \alpha-1)=1.$$引入常数$\alpha$是为了使$n^{\alpha n}\over(\alpha n)!$是$e^{-n}$阶的,其近似值为$\alpha=3.5911$.如果考虑$R_n$的渐进展开式,还可以得到收敛更快的公式,不过公式的形式同时也会变得很复杂.使用Bessel函数做类似的工作,我们将得到$$\gamma={A_n\over B_n}-\ln n+O(e^{-4n}),$$其中 \begin{align*} A_n&=\sum_{k=0}^{\alpha n}({n^k\over k!})^2H_k,\\ B_n&=\sum_{k=0}^{\alpha n}({n^k\over k!})^2,\\ H_k&=1+{1\over 2}+{1\over 3}+\cdots+{1\over k}, \end{align*}

这里$\alpha$的定义与上面相同,这个级数收敛得很快,并且很容易计算.使用类似于Richard外推加速法的技术对误差项进行渐进展开,可以得到$$\gamma={A_n\over B_n}-{C_n\over B_n^2}-\ln n+O(e^{-8n}),$$其中$$C_n={1\over 4n}\sum_{k=0}^{2n}\frac{((2k)!)^3}{(k!)^4(16n)^{2k}},$$并且$\alpha$满足$\alpha(\ln\alpha-1)=3$,其近似值为$\alpha=4.970625759$,利用这个公式可以计算$\gamma$到小数点后1亿位.

回到顶部其他常数

下面是一些其他数学常数的定义或常用计算公式.

  • Catalan常数. \begin{align*} C&=1-\frac{1}{3^2}+\frac{1}{5^2}-\frac{1}{7^2}+\frac{1}{9^2}-\cdots.\\ C&=\int_0^1\frac{\tan^{-1}x}{x}\,dx.\\ C&=\frac{3}{8}\sum_{k=0}^{\infty}\frac{(k!)^2}{(2k)!(2k+1)^2}+\frac{\pi\ln(2+\sqrt 3)}{8}. \end{align*}
  • Brun常数.$$B=(\frac{1}{3}+\frac{1}{5})+(\frac{1}{5}+\frac{1}{7})+(\frac{1}{11}+\frac{1}{13})+\cdots.$$
  • Mertsen常数.$$\mu=\lim_{p\rightarrow\infty}(\frac{1}{2}+\frac{1}{3}+\frac{1}{5}+\cdots+\frac{1}{p}-\ln\ln p).$$

参考文献

[1]Current published world record of pi calculation, http://www.super-computing.org/pi_current.html.

[2]Weisstein, Eric W., Machin-Like Formulas. From MathWorld--A Wolfram Web Resource

[3]Lehmer, D. H., On Arccotangent Relations for ., Amer. Math. Monthly 45 (1938), 657-664.

[4]Hwang, C.-L., More Machin-Type Identities., Math. Gaz. 81 (1997), 120-121.

[5]Berndt, B. C., Ramanujan's Notebooks, Part IV, Springer-Verlag, New York, 1994.

[6]Chudnovsky, D. V. and Chudnovsky, G. V., Approximations and Complex Multiplication According to Ramanujan, Ramanujan Revisited: Proceedings of the Centenary Conference, University of Illinois at Urbana-Champaign, June 1-5, 1987, Academic Press, 1987, 375-472.

[7]Borwein, J. M. and Borwein, P. B., More Ramanujan-Type Series for 1/PI, Ramanujan Revisited: Proceedings of the Centenary Conference, University of Illinois at Urbana-Champaign, June 1-5, 1987, Academic Press, 1988, 359-374.

[8]Weisstein, Eric W., Pi Formulas. From MathWorld--A Wolfram Web Resource

[9]Steve Pagliarulo, Stu's pi page, http://home.istar.ca/~lyster/otherconstants.html.

[10]David H. Bailey, Peter B. Borwein and Simon Plouffe, On the computation of the n'th decimal digit of various transcendental numbers, Mathematics of Computation (1996).

[11]Fabrice Bellard, A new formula to compute the n'th binary digit of pi, Bellard's Website (1997).

[12]Brent, Richard, Multiple-precision zero-finding methods and the complexity of elementary function evaluation, Analytic Computational Complexity (1975), 151-176.

[13]Borwein, Pi and the AGM, Addison-Wesley Longman Publishing Co., Inc., 1987.

[14]Borwein, J. M.; Borwein, P. B.; and Bailey, D. H, Ramanujan, Modular Equations, and Approximations to Pi, or How to Compute One Billion Digits of Pi, Amer. Math. Monthly 96 (1989), 201-219.

[15]Beeler, M.; Gosper, R. W.; and Schroeppel, R., HAKMEM, MIT Artificial Intelligence Laboratory, Memo AIM-239, Cambridge, MA, 1972.

[16]David V. Chudnovsky and Gregory V. Chudnovsky, Computer algebra in the service of mathematical physics and number theory, Computers and Mathematics 09 (1990), 232.

最后更新： 2009年04月28日 14:30:24 返回主页 / 更新记录 / 许可协议 Copyright © 2007,2008 maTHμ Project Group