Skip to content
Lex Blog

计算机如何计算对数函数

从数值分析的角度，分析计算机中对数函数的底层实现。

glibc 实现方法简述：

对于任意底数的对数，可以利用换底公式：

log⁡kx=ln⁡xln⁡k(1)\log_k{x} = \frac{\ln{x}}{\ln k} \tag{1}

方程 (1)(1) 中的 ln⁡k\ln k 是常数，因此对于任意底数的对数，我们只需计算 ln⁡x\ln x。

计算 ln⁡x\ln x

在高等数学中，我们学过下面的泰勒级数：

ln⁡(1+x)=x−12x2+13x3−14x4+⋯=∑k=1∞(−1)k−1xk(2)\ln(1+x) = x - \frac{1}{2}x^2 + \frac{1}{3}x^3 - \frac{1}{4}x^4 + \dots = \sum_{k=1}^{\infty}\frac{(-1)^{k-1}x}{k} \tag{2}

假设要求 ln⁡x\ln x，自然可以令 t=x−1t=x-1，于是：

ln⁡(x)=ln⁡(1+t)=∑k=1∞(−1)kktk(2)\ln(x) = \ln(1+t) = \sum_{k=1}^{\infty}\frac{(-1)^k}{k}t^k \tag{2}

然而，在实际计算中通常不可能展开到无穷项。若展开到第 nn 项，带皮亚诺余项的泰勒级数可以表示为：

ln⁡(1+t)=t−12t2+13t3−14t4+⋯=∑k=1n(−1)kktk+o(tn)(3)\ln(1+t) = t - \frac{1}{2}t^2 + \frac{1}{3}t^3 - \frac{1}{4}t^4 + \dots = \sum_{k=1}^{n}\frac{(-1)^k}{k}t^k + o(t^n) \tag{3}

注意这里的皮亚诺余项是 o(tn)o(t^n)，代表误差项。若 t>1t>1（即 x>2x > 2），随着 tt 或 xx 增大，上面公式的误差可能会变得相当大。

不过，还有另一个公式可用：

ln⁡x=2(x−1x+1+13(x−1x+1)3+15(x−1x+1)5+… )=2∑k=0∞12k+1(x−1x+1)2k+1(4)\ln x = 2\left(\frac{x-1}{x+1} + \frac{1}{3}\left(\frac{x-1}{x+1}\right)^3 + \frac{1}{5}\left(\frac{x-1}{x+1}\right)^5 + \dots\right) = 2\sum_{k=0}^{\infty}\frac{1}{2k+1}\left(\frac{x-1}{x+1}\right)^{2k+1} \tag{4}

以下是对方程 (4)(4) 的非正式证明：

在方程 (2)(2) 中分别用 1x\frac{1}{x} 与 −1x-\frac{1}{x} 替换 xx：

ln⁡(1+1x)=ln⁡(x+1x)=∑k=1∞(−1)kk1xkln⁡(1−1x)=−∑k=1∞1kxk\begin{aligned} \ln\left(1+\frac{1}{x}\right) &= \ln\left(\frac{x+1}{x}\right) = \sum_{k=1}^{\infty}\frac{(-1)^k}{k}\frac{1}{x^k} \\ \ln\left(1-\frac{1}{x}\right) &= -\sum_{k=1}^{\infty}\frac{1}{kx^k} \end{aligned}

注意：

ln⁡(x−1x)=−∑k=1∞1kxk⟹−ln⁡(x−1x)=∑k=1∞1kxk⟹ln⁡(xx−1)=∑k=1∞1kxk\begin{aligned} \ln\left(\frac{x-1}{x}\right) &= -\sum_{k=1}^{\infty}\frac{1}{kx^k} \\ \Longrightarrow -\ln\left(\frac{x-1}{x}\right) &= \sum_{k=1}^{\infty}\frac{1}{kx^k} \\ \Longrightarrow \ln\left(\frac{x}{x-1}\right) &= \sum_{k=1}^{\infty}\frac{1}{kx^k} \end{aligned}

接着：

ln⁡(x+1x)+ln⁡(xx−1)=ln⁡(x+1x−1)∑k=1n(−1)kkxk+∑k=1n1kxk=2∑k=0n12k+11x2k+1\begin{aligned} &\ln\left(\frac{x+1}{x}\right) + \ln\left(\frac{x}{x-1}\right) = \ln\left(\frac{x+1}{x-1}\right) \\ &\sum_{k=1}^n\frac{(-1)^k}{kx^k} + \sum_{k=1}^n\frac{1}{kx^k} = 2\sum_{k=0}^n{\frac{1}{2k+1}\frac{1}{x^{2k+1}}} \end{aligned}

由此可得：

ln⁡(x+1x−1)=2∑k=0n12k+11x2k+1\ln\left(\frac{x+1}{x-1}\right) = 2\sum_{k=0}^n{\frac{1}{2k+1}\frac{1}{x^{2k+1}}}

令 z=(x+1)(x−1)z = \frac{(x+1)}{(x-1)}，则 x=(z−1)(z+1)x = \frac{(z-1)}{(z+1)}。将其代入上式即得方程 (4)。

这个公式的好处在于，若展开到第 nn 项：

ln⁡x=2∑k=0n12k+1(x−1x+1)2k+1+o((x−1x+1)2k+1)(5)\ln x = 2\sum_{k=0}^{n}\frac{1}{2k+1}\left(\frac{x-1}{x+1}\right)^{2k+1} + o\left(\left(\frac{x-1}{x+1}\right)^{2k+1}\right) \tag{5}

其余项为 o((x−1)(x+1)2k+1)o\left(\frac{(x-1)}{(x+1)^{2k+1}}\right)，由于 lim⁡x→+∞(x−1x+1)2k+1=1\lim_{x\rightarrow+\infty}\left(\frac{x-1}{x+1}\right)^{2k+1} = 1，该余项保证误差保持在 11 的范围之内。

然而这仍然不够，因为误差为 11 还是太大了。根据上面余项的分析，若要更高精度，xx 的值应当更接近 11。

考虑到 IEEE 754 浮点标准，我们知道任何浮点数在计算机中都被表示为：

x=1.xxx×2expx = 1.xxx \times 2^{exp}

这意味着我们可以通过位运算直接获得指数（exp\text{exp}）与小数部分（1.xxx1.xxx）。

由于浮点数已经是这种形式，我们可以直接利用小数部分 1.xxx1.xxx。若用前面的公式计算 ln⁡(1.xxx)\ln(1.xxx)，就等价于计算 ln⁡(x2exp)\ln\left(\frac{x}{2^{exp}}\right)。根据对数公式：

ln⁡(x2exp)=ln⁡x−exp×ln⁡2\ln\left(\frac{x}{2^{exp}}\right) = \ln x - \text{exp}\times \ln 2

因此：

ln⁡x=ln⁡(x2exp)+exp×ln⁡2\ln x = \ln\left(\frac{x}{2^{exp}}\right) + \text{exp} \times \ln 2

至此就完成了 ln⁡x\ln x 的计算，等价于计算任意底数的对数。

不过这只是一个总体概述，glibc 在实现中采用了各种提升精度的技巧。