隐藏目录
幂函数
迭代法
Newton迭代
对数函数
级数法
迭代法
指数函数
级数法
迭代法
矩形求和
三角函数与反三角函数
折半求和的应用
指数函数
余弦函数
其他函数

函数是数学中的核心概念,符号计算中许多数值问题都可以归结到函数求值上,因此一个友好的计算机代数系统应该提供大部分基本数学函数的计算方法.如果对精度要求不是太高,譬如说机器精度时,可以利用直接查表并结合函数的多项式逼近,有理分式(Pade)逼近或者将函数转化为数值积分等方法来计算函数值,这些都是数值计算中的重要技术.在这里我们主要讨论的是如何将函数值计算到任意精度,总的来看,与数学常数的情形类似,计算的方法大致分为级数法和迭代法两类,关于它们更多更详尽的讨论可以在Borwein的书[1]中找到.

回到顶部幂函数

正整数幂在基础数论算法中有详细的讨论,实用中一般采用简单的基于$$x^{2n}=(x^n)^2,\quad x^{2n+1}=x(x^n)^2$$的二分法,有理数幂$$x^{q\over p}=(\sqrt[ p]{x})^q$$归结到开方计算,更一般的幂函数$x^m$则利用恒等式$$x^m=e^{m\ln{x}}$$归结到对数函数与指数函数的计算,因此接下来主要讨论的是开方.

回到顶部迭代法

计算开方最有效的方法是Newton迭代法,设$f(x)=x^n-a$,$x_0$是$f(x)$的一个近似根且$f'(x_0)\neq 0$,那么迭代公式为$$x_{k+1}=(1-{1\over n})x_k+{a\over nx_k^{n-1}},\quad k=0,1,\ldots,$$其中$$\lim_{k\rightarrow\infty}x_k=\sqrt[ n]{a}.$$特别地,当$n=2$时迭代公式化为$$x_{k+1}={x_k\over 2}+{a\over 2x_k},$$每迭代一次只需要做一次除法.有时候为了避免除法也可以先求倒数的开方,设$f(x)={1\over x^n}-a$,对应的迭代公式为$$x_{k+1}=(1+{1\over n})x_k-{ax^{n+1}\over n},$$最后再利用$$\sqrt[ n]{a}=a({1\over\sqrt[ n]{a}})^{n-1}$$也可以求出$\sqrt[ n]{a}$.特别地,当$n=2$时迭代公式化为$$x_{k+1}=x_k+x_k{1-ax_k^2\over 2},$$最后一步只需要计算$\sqrt{a}=a{1\over\sqrt{a}}$.

Newton迭代法是二阶收敛并且自修正的,因此迭代过程中有效精度以2的幂次增长,这样一来,开方时大部分时间都是在计算最后一步迭代中的全精度乘(除)法,巧妙地将最后一步迭代和最后的乘法合成为一步可以加速开方的计算.以$n=2$为例,设$r$是倒数第二步迭代后得到的$1\over\sqrt{a}$的近似值,则$ar$是$\sqrt{a}$具有相同精度的近似值,将$\sqrt{a}$的迭代公式改写成$$x_{k+1}=x_k+{1\over x_k}{a-x_k^2\over 2}$$的形式,用$ar$代替$x_k$,$r$代替$1\over x_k$,迭代一步后得到的$x_{k+1}$就是$\sqrt{a}$的全精度近似值.

对于负幂的情形,只需要计算$x^{-1}$,即$x$的倒数.设$f(x)={1\over x}-a$,$x_0$是$f(x)$的一个近似根且$f'(x_0)\neq 0$,那么对应的Newton迭代公式为$$x_{k+1}=2x_k-ax_k^2,\quad k=0,1,\ldots,$$其中$$\lim_{k\rightarrow\infty}x_k={1\over a}.$$

数论中常常还会用到整数开方$[\sqrt[ n]{a}],a\in\mathbb{Z}$,这时候可以将Newton迭代公式中的除法全部换成整数带余除法.以$n=2$为例,新的迭代公式为$$x_{k+1}=\left[{x_k+[{a\over x_k}]\over 2}\right].$$可以证明,如果迭代初始值$x_0<[\sqrt[ n]{a}]$,那么$x_k$将单调上升,并且趋向于$[\sqrt[ n]{a}]$,因此一般取$$x_0=2^{[{\log_2{a}\over n}]}.$$

最后值得一提的是一种类似于除法而可以借助纸笔计算的开方方法,它可以从方根的最高位开始依次给出各位上的数字.以二进制下的开平方这种最简单的情形为例,设$N=[\log_2{a}]$,即$a$在二进制下共有$N$位,那么$\sqrt{a}$在二进制下将有$N/2$位,即第$N/2$上的数字是1.令$u=2^{N/2}$,$u$可以看成$\sqrt{a}$的1位近似值,再令$v=2^{N/2-1}$,比较$(u+v)^2$与$a$的大小,若$(u+v)^2\ge a$,那么$a$的第$N/2-1$位上的数字就是1,否则是0.接下来的每一步都让$u$等于$\sqrt{a}$已经求出的近似值,$ v$等于近似值最低位的下一位,这样就可以顺次求出$\sqrt{a}$各位上的数字.

回到顶部Newton迭代

从几何的角度来看,方程$f(x)=0$的根$x^*$是曲线$y=f(x)$与x轴的交点的横坐标,设$x_k$是根$x^*$的某个近似值且$f'(x_k)\neq 0$,过曲线$y=f(x)$上横坐标为$x_k$的点$P_k$引切线,并将该切线与x轴的交点的横坐标$x_{k+1}$作为$x^*$新的近似值,根据切线方程$$y=f(x_k)+f'(x_k)(x-x_k)$$可以解出$$x_{k+1}=x_k-{f(x_k)\over f'(x_k)},$$这就是Newton迭代,关于它的收敛性有$$\lim_{k\rightarrow\infty}{x_{k+1}-x^*\over(x_k-x^*)^2}={f''(x^*)\over 2f'(x^*)},$$因此有时也称之为二阶Newton迭代法或者切线法.设$$g(x)={f(x)\over \sqrt{|f'(x)|}},$$那么$g(x)$对应的Newton迭代公式为$$x_{k+1}=x_k-{2f(x_k)f'(x_k)\over 2(f'(x_k))^2-f(x_k)f''(x_k)},$$它是三阶收敛的,比$f(x)$对应的原始Newton迭代公式要快,使用这种迭代公式的迭代方法称为Halley迭代法.以开方为例,设$f(x)=x^n-a$,那么$f(x)$对应的Halley迭代公式为$$x_{k+1}=x_k\cdot{n(a+x_k^n)+(a-x_k^n)\over n(a+x_k^n)-(a-x_k^n)}.$$

Newton迭代法是数值稳定的,它常常被用来计算反函数的值.如果已知$f(x)=y_0$,改写成$g(x)=f(x)-y_0=0$的形式,那么$g(x)$对应的迭代公式为$$x_{k+1}=x_k-{f(x_k)-y_0\over f'(x_k)},$$这样就可以求出$x^*\approx f^{-1}(y_0)$. 实际计算中应该利用迭代公式的自修正性质适当安排中间精度,一般情况下每一步迭代计算时使用的精度可以只是下一步的一半,据此还可以证明如果使用FFT乘法等快速乘法算法,浮点数除法的计算复杂度与乘法相同.

如果在$P_k$处引$y=f(x)$的二次切线或三次切线,相应地可以得到三阶和四阶的Newton迭代法.以三次Newton迭代法为例,二次切线的抛物线方程对应于$y=f(x)$在$x_k$处的二阶Taylor近似$$y=f(x)\approx f(x_k)+f'(x_k)(x-x_k)+{f''(x_k)\over 2}(x-x_k)^2,$$可以解出$$x_{k+1}-x_k=-{f'(x_k)\over f''(x_k)}(1-\sqrt{1-{2f(x_k)f''(x_k)\over(f'(x_k))^2}}),$$利用近似关系$$1-\sqrt{1-\alpha}\approx{\alpha\over 2}+{\alpha^2\over 8}$$可以将迭代公式化为$$x_{k+1}-x_k=-{f(x_k)\over f'(x_k)}(1+{f(x_k)f''(x_k)\over 2(f'(x_k))^2}).$$

更高阶的迭代公式可以使用Shroeder方法或Householder方法来系统地构造,以$1\over\sqrt{a}$为例,设$y=1-ax^2$,那么Taylor展开式$${1\over\sqrt{a}}={x\over\sqrt{1-y}}=x(1+{y\over 2}+{3y^2\over 8}+\cdots)$$给出了迭代公式序列$$x_{k+1}=x_k(1+{1-ax_k^2\over 2}+{3(1-ax_k^2)^2\over 8}+\cdots).$$在介绍AGM迭代曾经提到过$m$阶迭代公式每迭代一步,有效位数变为原来的$m$倍,如果要求结果的有效位数是$ P$,总的迭代次数就近似为$$\ln{P}\over\ln{m}.$$迭代算法的复杂度在$m\approx 1$到$m=2$之间有一个飞跃,而当$m\ge 3$时不同阶的迭代公式所需迭代次数大致是同一个数量级,因此实用中很少使用超过4阶的迭代公式.

回到顶部对数函数

一般的对数函数$\log_a{x}$可以利用恒等式$$\log_a{x}={\ln{x}\over\ln{a}}$$归结到自然对数$\ln{x}$的计算,因此接下来所说的对数函数特指自然对数$\ln{x}$.在实际计算前可以利用一些函数关系,譬如说 \begin{align*} \ln{x}&=2^k\ln{x^{2^{-k}}},\\ \ln{x}&=k\ln{2}+\ln(2^{-k}x), \end{align*} 这样一来可以大大加快级数的收敛速度.

回到顶部级数法

$\ln(1-x)$在$x=0$处的Taylor展开式为$$-\ln(1-x)=x+{x^2\over 2}+{x^3\over 3}+\cdots,\quad 0<x<1,$$除非$x$很小,否则这个级数的收敛速度相当慢,只适用于较低精度的计算.由$$\ln{x+1\over x-1}=2({1\over x}+{1\over 3x^3}+{1\over 5x^5}+\cdots),\quad |x|>1,$$也可以得到$\ln{x}$的有理分式级数展开$$\ln{x}=2({x-1\over x+1}+{1\over 3}({x-1\over x+1})^3+{1\over 5}({x-1\over x+1})^5+\cdots),$$可以采用Taylor级数的矩形求和方法或超几何级数的折半求和方法来计算这些级数.

回到顶部迭代法

计算$\ln{x}$更正式与更快的方法总是基于迭代的,如果有快速计算$e^x$的方法,那么可以利用反函数关系$$e^y=x\Leftrightarrow y=\ln{x}$$使用Newton迭代来计算$\ln{x}$.二阶和三阶的Newton迭代公式分别是 \begin{align*} y_{k+1}&=y_k-1+{x\over e^{y_k}},\\ y_{k+1}&=y_k-2\frac{e^{y_k}-x}{e^{y_k}+x}, \end{align*} 迭代过程中主要计算量在$e^{y_k}$上,因此常常选择更快的三阶迭代公式,迭代的初始值取为$[\log_2{x}]\ln{2}$的近似值.

计算对数常数$\ln 2$时曾经用到过代数几何平均值(AGM)迭代,记$${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,10^{-N}x)$作为$\ln{x}$的近似值,依然借用计算$AGM(a,b)$时的迭代结构,那么在计算$AGM(a,b)$的同时就可以计算出$a_n,b_n$,进而可以计算出$R(a,b)$,这样就得到了计算$\ln{x}$的一个二阶收敛的迭代算法.这个算法用于计算$\ln{x}$时更常用的形式是以渐进关系$$\ln{x}=\pi x\frac{1+4x^{-2}(1-1/\ln{x})+O(x^{-4})}{2AGM(x,4)}$$给出的,当$|x|$很大时,譬如说$4x^{-2}<\epsilon$,其中$\epsilon$为需要计算的精度,渐进关系中的分子可以直接用1来代替,即$$\ln{x}=\frac{\pi x}{2AGM(x,4)},$$实际计算中根据精度要求可以利用函数关系 \begin{align*} \ln{x}&=2^{-k}\ln{x^{2^k}},\\ \ln{x}&=-k\ln{2}+\ln(2^kx), \end{align*} 扩大$x$以满足$x>{2\over\sqrt{\epsilon}}$的近似条件.

回到顶部指数函数

与对数函数类似,一般的指数函数$a^x$可以利用恒等式$$a^x=e^{x\ln{a}}$$归结到$e^x$的计算,因此接下来可以只讨论$e^x$.在实际计算前可以利用一些函数关系缩小$x$,譬如说 \begin{align*} e^x&=(e^{2^{-k}x})^{2^k},\\ e^x-1&=(e^{2^{-1}x}+1)(e^{2^{-2}x}+1)\cdots(e^{2^{-k}x}+1)(e^{2^{-k}x}-1), \end{align*} 使用第二个函数关系将产生更小的截断误差.

回到顶部级数法

和计算自然对数底$e$一样,可以直接利用$e^x$的Taylor展开式$$e^x=1+{1\over 1!}x+{1\over 2!}x^2+\cdots,$$取$$s_n(x)=\sum_{k=0}^n{x^n\over n!}$$作为$e^x$的近似值.除此之外,定义$$r_{m,n}(x)=\frac{\int_0^\infty t^n(t+z)^me^{-t}\,dt}{\int_0^\infty (t-z)^nt^me^{-t}\,dt},$$它的闭形式表达式为$$r_{m,n}(x)=\frac{\sum_{k=0}^m{C_m^k\over C_{n+m}^k}\cdot{x^k\over k!}}{\sum_{k=0}^n{C_n^k\over C_{m+n}^k}\cdot{(-x)^k\over k!}},$$则$r_{m,n}(x)$是$e^x$的有理分式级数展开.多项式逼近与有理分式逼近的误差估计分别为 \begin{align*} |e^x-s_n(x)|&\le{1\over(n+1)!}(1+{2\over n+1}),\\ |e^x-r_{n,n}(x)|&\le{8(n!)^2\over(2n+1)((2n)!)^2}. \end{align*} 类似的,也可以采用Taylor级数的矩形求和方法或超几何级数的折半求和方法来计算这些级数.

回到顶部迭代法

$e^x$没有很好的AGM迭代公式,如果有快速计算自然对数$\ln{x}$的方法,那么可以利用反函数关系$$y=e^x\Leftrightarrow x=\ln{y}$$使用Newton迭代法来计算$e^x$.其中二阶的Newton迭代公式为$$y_{n+1}=y_n(x+1-\ln{y_n}),$$三阶的Halley迭代公式为$$y_{n+1}=y_n\frac{1+(x+1-\ln{y_n})^2}{2}.$$进一步地,设$z=x-\ln{y}$,那么Taylor展开式$$e^x=ye^z=y(1+{z\over 1!}+{z^2\over 2!}+\cdots)$$给出了迭代公式序列$$y_{k+1}=y_k(1+(x-\ln{y_k})+{(x-\ln{y_k})^2\over 2!}+\cdots),$$实际计算中迭代公式的最佳阶数取决于自然对数$\ln{x}$的计算速度.

回到顶部矩形求和

截断后的Taylor级数一般形式为$$s_n(x)=\sum_{k=0}^na_kx^k,$$其中$a_k$常常可以表示成$k$的有理函数,因此可以事先写出$a_{k+1}/a_k$的通项公式.如果$n$不是很大,可以直接把它看成$n$次多项式来求值,例如使用类似于Horner法则的$$s_n(x)=a_0(1+{a_1\over a_0}x(1+{a_2\over a_1}x(1+{a_3\over a_2}x(\cdots+{a_n\over a_{n-1}}x)))),$$这个法则更常用的形式可以表示成如下算法.

算法1(Taylor级数求和)

输入:求值点$x_0$与级数$$s_n(x)=\sum_{k=0}^na_kx^k.$$

输出:$s_n(x_0)$.

  1. 令$t=a_0,s=a_0$.
  2. $k$从$0$到$n-1$,顺次计算$t={a_{k+1}\over a_k}x_0\cdot t$,$s=s+t$,若$t$达到指定精度,终止循环.
  3. 返回$s$.
这个算法用于Taylor级数求值一共需要$n$次短乘法.

更好的方法是矩形求和,设$n=r\cdot c-1$,那么多项式$s_n(x)$可以写成嵌套求和的形式 \begin{align*} s_n(x)&=(a_0+a_rx^r+\cdots+a_{(c-1)r}x^{(c-1)r})\\ &+x(a_1+a_{r+1}x^r+\cdots+a_{(c-1)r+1}x^{(c-1)r})\\ &+\cdots\\ &+x^{r-1}(a_{r-1}+a_{2r-1}x^r+\cdots+a_{cr-1}x^{(c-1)r}), \end{align*} 先使用快速求幂算法计算出$x^r$及$x^{2r},\ldots,x^{(c-1)r}$,然后再对外层的$r-1$次多项式使用Horner法则求值,计算各行时可以利用$a_k$与$a_{k+1}$之间的函数关系.

如果不计常系数乘法和计算$x^r$的过程,矩形求和的过程只需要大约$c+r<n$次乘法,这是因为它将大量乘以$x$的短乘法转化为大数之间的长乘法,这样就可以充分发挥快速乘法算法的威力.

现在将矩形的长和宽对换,那么多项式$s_n(x)$可以写成另一种嵌套求和的形式 \begin{align*} s_n(x)&=(a_0+a_1x+\cdots+a_{r-1}x^{r-1})\\ &+x^r(a_{r+1}+a_{r+2}x+\cdots+a_{2r-1}x^{r-1})\\ &+\cdots\\ &+x^{(c-1)r}(a_{(c-1)r}+a_{(c-1)r+1}x+\cdots+a_{cr-1}x^{r-1})\\ &=\sum_{k=0}^{c-1}x^{kr}p_{k,r}(x). \end{align*} 对于特别的Taylor级数,$p_{k,r}(x)$可以表示成仅关于$k,x$的函数,即$$p_{k,r}(x)=p_r(k,x).$$以$$-\ln(1-x)=x+{x^2\over 2}+{x^3\over 3}+\cdots$$为例,对应的$$s_n(x)=\sum_{k=1}^n{x^k\over k},$$那么$$s_{cr}(x)=\sum_{k=0}^{c-1}x^{kr}p_r(k),$$其中$$p_r(k)=\sum_{j=1}^r{x^j\over j+r\cdot k}.$$将$p_r(k)$看作关于$k$的有理分式函数$$p_r(k)={w_r(k)\over v_r(k)},$$利用多项式快速多点求值算法就可以同时求出$p_r(k)$在$0,1,\ldots,c-1$处的值,进而可以使用Horner法则计算出$$s_{cr}(x)=\sum_{k=0}^{c-1}p_r(k)(x^r)^k.$$事实上大部分初等函数的Taylor级数都满足上面的条件,这种基于多项式快速多点求值的算法的计算复杂度为$O(\sqrt{n}(\log n)^2M(n))$,其中$M(n)$表示两个$n$位整数相乘的计算复杂度.

回到顶部三角函数与反三角函数

从复分析中我们知道所有的初等函数都可以用$e^z$及其反函数$\ln{z}$表示出来,因此理论上我们只需要讨论对数函数和指数函数就可以解决所有的初等函数,关于这两个函数的计算已经涵盖了函数求值的主要方法,这里对三角函数与反三角函数的讨论将主要是一些技术性或描述性的内容.

除开具有直接倒数关系的函数,常用的三角函数和反三角函数共有3对,它们之间由简单的函数关系联系着,只需要计算出其中一个函数值就很容易求出其余两个函数值,通常选择三角函数中的$\cos{x}$和反三角函数中的$\tan^{-1}{x}$作为标准函数,这种计算体系对应的函数关系为 \begin{align*} \sin{x}&=\sqrt{1-\cos^2{x}},\\ \tan{x}&=\sqrt{{1\over\cos^2{x}}-1},\\ \sin^{-1}{x}&=\tan^{-1}{x\over\sqrt{1-x^2}},\\ \cos^{-1}{x}&=\tan^{-1}{\sqrt{1-x^2}\over x}. \end{align*}

$\cos{x}$一般直接通过Taylor展开式$$\cos{x}=1-{x^2\over 2!}+{x^4\over 4!}-\cdots$$来计算,先利用诱导公式将参数$x$限定到$(0,{\pi\over 2})$上,然后利用函数关系 \begin{align*} \cos{x}&=2\cos^2{x\over 2}-1,\\ 1-\cos{x}&=4(1-\cos{x\over 2})-2(1-\cos{x\over 2})^2, \end{align*} 缩小$x$,这样可以大大加快级数的收敛速度.

如果使用Taylor展开式$$\tan^{-1}{x}=x-{x^3\over 3}+{x^5\over 5}-\cdots$$来计算$\tan^{-1}{x}$,先根据$$\tan^{-1}{x}+\tan^{-1}{1\over x}={\pi\over 2},\quad x>0$$将参数$x$限定到$(0,1)$上,类似地,也可以利用函数关系$$\tan^{-1}{x}=2\tan^{-1}{x\over 1+\sqrt{1+x^2}}$$缩小$x$以加快级数的收敛速度.除此之外,由于$$\tan^{-1}{x}={1\over 2i}\ln{1+ix\over 1-ix},$$可以通过计算对数函数的快速AGM迭代来计算反正切函数,或者也可以直接据此构造出$\tan^{-1}{x}$的迭代算法,AGM迭代的计算复杂度比级数方法低.如果再利用反函数的Newton迭代计算出$\tan{x}$,就可以在$O(M(n)\log n)$的时间内计算出所有三角函数的$n$位有效数字.[2]

回到顶部折半求和的应用

折半求和[3]是用来计算超几何级数在有理点处的值的一项重要技术,大部分初等函数都具有超几何级数的展开式.

回到顶部指数函数

设求值点$x_0={q\over p}$,其中$0\le x_0<1,(p,q)=1$.如果$p,q$都不是很大,可以直接代入折半求和公式进行计算,否则先将$x_0$写成$$x_0=\sum_{k=0}^m2^{-2^k}u_k$$的形式,其中$0\le u_k<2^{2^{k-1}}$.记$r_k=2^{-2^k}u_k$,则$r_k$的分母是$2^{2^k}$,使用折半求和计算$e^{r_k}$时可以用移位来代替除法,最后再计算$$e^{x_0}=e^{r_0}e^{r_1}\cdots e^{r_m}$$就求出了$e^x$在有理点$x_0$处的值.

回到顶部余弦函数

类似的,先将$x_0$写成$$x_0=\sum_{k=0}^m2^{-2^k}u_k$$的形式,其中$0\le u_k<2^{2^{k-1}}$.记$r_k=2^{-2^k}u_k$,则$r_k$的分母是$2^{2^k}$,并令 \begin{align*} x_j&=\sum_{k=j}^mr_k,\\ S_j&=\sin{r_j},\\ C_j&=\cos{r_j}, \end{align*} 那么 \begin{align*} \sin x_j&=S_j\cos x_{j+1}+C_j\sin x_{j+1},\\ \cos x_j&=C_j\cos x_{j+1}-S_j\sin x_{j+1}. \end{align*} 使用折半求和计算$S_j$和$C_j$时可以用移位来代替除法,这样就可以递推地计算出$\cos{x_0}$和$\sin{x_0}$来.

回到顶部其他函数

下面给出了一些其他常用初等函数的定义或计算公式.

  • 双曲三角函数 \begin{align*} \cosh{x}&={1\over 2}(e^x+e^{-x})\\ \sinh{x}&={1\over 2}(e^x-e^{-x})\\ \tanh{x}&=1-{2\over e^{2x}+1}\\ \cosh^{-1}{x}&=\ln(x+\sqrt{x^2-1})\\ \sinh^{-1}{x}&=\ln(x+\sqrt{x^2+1})\\ \tanh^{-1}{x}&={1\over 2}\ln{1+x\over 1-x} \end{align*}

参考文献

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

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

[3]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年01月28日 22:41:14 返回主页 / 更新记录 / 许可协议 Copyright © 2007,2008 maTHμ Project Group