计算机如何计算对数函数

计算机是如何计算对数函数 log ⁡ a x \log_a{x} loga​x 的呢？

本文并不是为了探究各种编程语言中数学计算函数库对 log ⁡ \log log 函数的具体实现，而是提供一种计算思路，来帮助理解 log ⁡ \log log 函数的计算原理，并使用 C 语言来实现这个算法。

提起对数函数 log ⁡ a x \log_a{x} loga​x，大家最熟悉的对数函数应该是以自然底数 e e e 为底的对数函数： ln ⁡ x \ln{x} lnx，即 log ⁡ e x \log_e{x} loge​x。

ln ⁡ x \ln{x} lnx 函数的图像如下所示：

对于不是以 e e e 为底的对数函数 log ⁡ a x \log_a{x} loga​x ，可以使用换底公式，将其转换为以 e e e 为底的函数形式。 换底公式： log ⁡ a b = log ⁡ c b log ⁡ c a 换底公式：\log_a{b} = \frac{\log_c{b}}{\log_c{a}} \qquad\qquad\qquad 换底公式：loga​b=logc​alogc​b​ 根据换底公式： log ⁡ a x = log ⁡ e x log ⁡ e a = ln ⁡ x ln ⁡ a (1) \log_a{x} = \frac{\log_e{x}}{\log_e{a}} = \frac{\ln{x}}{\ln{a}} \tag{1} \qquad\qquad\qquad loga​x=loge​aloge​x​=lnalnx​(1) 于是，计算对数函数 log ⁡ a x \log_a{x} loga​x 的问题就转换为计算 ln ⁡ x \ln{x} lnx 的问题。

对于计算 ln ⁡ x \ln{x} lnx，可能我们第一时间想到的就是泰勒展开公式： ln ⁡ ( x ) = ∑ n = 0 ∞ ( − 1 ) n n + 1 ( x − 1 ) n + 1 = ( x − 1 ) − 1 2 ( x − 1 ) 2 + 1 3 ( x − 1 ) 3 − 1 4 ( x − 1 ) 4 + 1 5 ( x − 1 ) 5 + ⋅ ⋅ ⋅ , x ∈ ( 0 , 2 ] (2) \ln{(x)} = \sum_{n = 0}^{\infty} \frac{(-1)^n}{n+1}{(x-1)}^{n+1} = (x - 1) - \frac{1}{2}{(x-1)}^2 + \frac{1}{3}{(x-1)}^3 - \frac{1}{4}{(x-1)}^4 + \frac{1}{5}{(x-1)}^5 + ···, \quad x \in (0, 2] \tag{2} ln(x)=n=0∑∞​n+1(−1)n​(x−1)n+1=(x−1)−21​(x−1)2+31​(x−1)3−41​(x−1)4+51​(x−1)5+⋅⋅⋅,x∈(0,2](2) 但是，使用泰勒公式对 ln ⁡ x \ln{x} lnx 展开是有条件的，即 x x x 所属的区间为 x ∈ ( 0 , 2 ] x \in (0, 2] x∈(0,2] 。而且，展开的项数一定是固定的，毕竟不能无限次的进行展开计算。

要知道的是，使用泰勒公式对 ln ⁡ x \ln{x} lnx 进行展开，始终是一种近似，而非绝对的等于。决定展开后的误差的因素有两点：

一： x x x 与 1.0 1.0 1.0 的距离。在展开项数相同的情况下， x x x 越接近 1.0 1.0 1.0，误差越小；

二：展开的项数。对于同一个 x x x，展开项数越多，误差越小。

上图就是 ln ⁡ x \ln{x} lnx 与不同项数的泰勒展开式的图像。

作为对比，下图列出了上述几个展开式与 ln ⁡ x \ln{x} lnx 之间的误差：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-LqjbDJGO-1676008929754)(计算机如何计算对数函数/ln与展开式误差.svg)]

结合上面两张图，可以看出， x x x 越趋向于 1.0 1.0 1.0 ，展开项数越多，精度越高。因此，要想办法把 x x x 转化到 1 1 1 的附近，就能够对 x x x 使用泰勒展开的方式进行计算了。

根据公式 ( 1 ) (1) (1) ，我们成功将 log ⁡ a x \log_a{x} loga​x 转化为求 ln ⁡ x \ln{x} lnx 的问题，对于 ln ⁡ x \ln{x} lnx ，我们希望能够使用泰勒公式对其进行展开。但是由于不能确定 x x x 的值，不能直接使用泰勒公式进行展开，那么我们就要想办法处理 x x x，让其满足泰勒公式展开的条件，并且尽可能的提高计算精度。

在处理 x x x 之前，我们要先弄清楚，数据在计算机中是如何存储的。

在计算机中，使用浮点数来表示一个带有小数的数（整数也可以使用浮点数来表示）。浮点数使用 IEEE 格式存储，分为单精度浮点数（float）、双精度浮点数（double）。

对于单精度浮点数而言，占 4 个字节，即 32 位。其中使用 1 位表示符号位，8 位表示二进制形式的阶码，剩余 23 位表示尾巴数。单精度浮点数的表示范围大约在 3. 4 − 38 − 3. 4 38 3.4^{-38} - 3.4^{38} 3.4−38−3.438 之间。

其中，最高位为符号位 f f f，符号位 f = 0 f = 0 f=0 代表该数是正数， f = 1 f = 1 f=1 代表该数是负数。

中间 8 位表示阶码 j j j ，阶码部分为纯整数，使用移码表示，真实值为 j − 127 j-127 j−127。

低 23 位表示尾数 m m m，在构造浮点数据的时候需要进行规格化，即要保证尾数的整数部分为 1，而这个 1 是不储存的，所以经过规格化后 m ∈ ( 1 , 2 ) m \in(1,2) m∈(1,2)。

因此采用 IEEE 格式保存的数据的值为： ( − 1 ) f ⋅ m ⋅ 2 j − 127 (-1)^f \cdot m \cdot 2^{j-127} (−1)f⋅m⋅2j−127 。

以 13.75 为例：

1、将十进制数转换为二进制数： 13.75 = 1101.11 B 13.75 = 1101.11B 13.75=1101.11B，即 13 = 8 + 4 + 1 = 2 3 + 2 2 + 2 0 = 1101 B 13 = 8 + 4 + 1 = 2^3 + 2^2 + 2^0 = 1101B 13=8+4+1=23+22+20=1101B， 0.75 = 0.5 + 0.25 = 2 − 1 + 2 − 2 = 0.11 B 0.75 = 0.5 + 0.25 = 2^{-1} + 2^{-2} = 0.11B 0.75=0.5+0.25=2−1+2−2=0.11B，因此 13.75 = 13 + 0.75 = 1101 B + 0.11 B = 1101.11 B 13.75 = 13 + 0.75 = 1101B + 0.11B= 1101.11B 13.75=13+0.75=1101B+0.11B=1101.11B；

2、将该二进制数规格化： 1101.11 = 1.10111 × 2 3 1101.11 = 1.10111 \times 2^3 1101.11=1.10111×23，尾数的整数位是隐藏的，因此尾数部分存储的是： 10111 10111 10111 。

3、将阶码转换为移码表示： 3 + 127 = 130 = 11 B + 01111111 B = 10000010 B 3 + 127 = 130 = 11B + 01111111B = 10000010B 3+127=130=11B+01111111B=10000010B，所以阶码中存储的数据为： 100000010 100000010 100000010。

因此，在计算机中， 13.75 存储为：

双精度浮点数（double）存储形式与单精度浮点数（float）相似，只不过双精度浮点数占 8 个字节，64 位，最高位代表符号位，中间 11 位为阶码，剩余 52 位表示尾数。

借助 float 数据的存储方式，我们将 x x x 进行转换： x = ( − 1 ) f ⋅ m ⋅ 2 j − 127 x = (-1)^f \cdot m \cdot 2^{j-127} x=(−1)f⋅m⋅2j−127，由于 ln ⁡ x \ln{x} lnx 函数的定义域为 ( 0 , + ∞ ) (0,+\infty) (0,+∞) ，因此阶码始终为 0，于是，化简为 x = m ⋅ 2 j − 127 x = m \cdot 2^{j-127} x=m⋅2j−127。

借助对数函数的另外两个性质： log ⁡ a ( b ⋅ c ) = log ⁡ a ( b ) + log ⁡ a ( c ) log ⁡ a ( b j ) = j ⋅ log ⁡ a ( b ) (3) \log_a{(b \cdot c)} = \log_a{(b)} + \log_a{(c)} \\ \log_a{(b^j)} = j \cdot \log_a{(b)} \tag{3} loga​(b⋅c)=loga​(b)+loga​(c)loga​(bj)=j⋅loga​(b)(3) 可得： ln ⁡ x = ln ⁡ [ ( 1. m ) ⋅ 2 j − 127 ] = ln ⁡ ( 2 j − 127 ) + ln ⁡ ( m ) = ( j − 127 ) ⋅ ln ⁡ 2 + ln ⁡ m (4) \ln{x} = \ln{[(1.m) \cdot 2^{j-127}]} = \ln{(2^{j-127})} + \ln{(m)} = (j-127) \cdot \ln{2} + \ln{m} \tag{4} lnx=ln[(1.m)⋅2j−127]=ln(2j−127)+ln(m)=(j−127)⋅ln2+lnm(4) 在 (4) 式中， ( j − 127 ) (j-127) (j−127) 为 x x x 的阶码，可以想办法取； ln ⁡ 2 ≈ 0.693147 \ln2 \approx 0.693147 ln2≈0.693147 为常数，因此，只需要计算 ln ⁡ m \ln{m} lnm 的值即可。

又因为 m m m 为 x x x 的尾数，根据规格化的规则， m ∈ ( 1 , 2 ) m \in(1,2) m∈(1,2)，符合对其进行泰勒展开的要求。

为了展开的精度更加精确，我们应该尽量让 m m m 趋向于 1 1 1，根据 (3) 式： ln ⁡ ( m ) = ln ⁡ ( 2 ⋅ 2 2 ⋅ m ) = ln ⁡ 2 + ln ⁡ ( 2 2 m ) = 1 2 ln ⁡ 2 + ln ⁡ ( 2 2 m ) (5) \ln{(m)} = \ln{(\sqrt{2} \cdot\frac{\sqrt{2}}{2} \cdot m)} = \ln{\sqrt{2}} + \ln{(\frac{\sqrt{2}}{2} m)} = \frac{1}{2}\ln{2} + \ln{(\frac{\sqrt{2}}{2} m)} \tag{5} ln(m)=ln(2 ​⋅22 ​​⋅m)=ln2 ​+ln(22 ​​m)=21​ln2+ln(22 ​​m)(5) 因为 m ∈ ( 1 , 2 ) m \in(1,2) m∈(1,2)，所以 2 2 m ∈ ( 2 2 , 2 ) \frac{\sqrt{2}}{2} m \in(\frac{\sqrt{2}}{2},\sqrt{2}) 22 ​​m∈(22 ​​,2 ​)，即 m m m 的大致范围为 ( 0.707106 , 1.414213 ) (0.707106, 1.414213) (0.707106,1.414213)，更加接近于 1 1 1。

于是，我们最终的式子就是： ln ⁡ x = ( j − 127 ) ⋅ ln ⁡ 2 + 1 2 ln ⁡ 2 + ln ⁡ ( 2 2 m ) x = m ⋅ 2 j − 127 (6) \ln{x} = (j-127) \cdot \ln{2} + \frac{1}{2}\ln{2} + \ln{(\frac{\sqrt{2}}{2} m)} \\ x = m \cdot 2^{j-127} \tag{6} lnx=(j−127)⋅ln2+21​ln2+ln(22 ​​m)x=m⋅2j−127(6) 对于 ln ⁡ ( 2 2 m ) \ln{(\frac{\sqrt{2}}{2}m)} ln(22 ​​m)，使用泰勒公式将其展开： 令 t = ( 2 2 m ) 则 ln ⁡ ( 2 2 m ) = ln ⁡ ( t ) = ( t − 1 ) − 1 2 ( t − 1 ) 2 + 1 3 ( t − 1 ) 3 − 1 4 ( t − 1 ) 4 + 1 5 ( t − 1 ) 5 − 1 6 ( t − 1 ) 6 + 1 7 ( t − 1 ) 7 + o ( t − 1 ) 7 (7) 令 \quad t = (\frac{\sqrt{2}}{2}m) \\ 则 \quad \ln{(\frac{\sqrt{2}}{2}m)} = \ln{(t)} = (t-1) - \frac{1}{2}(t-1)^2 + \frac{1}{3}(t-1)^3 -\frac{1}{4}(t-1)^4 + \frac{1}{5}(t-1)^5 - \frac{1}{6}(t-1)^6 + \frac{1}{7}(t-1)^7 + o(t-1)^7 \tag{7} 令t=(22 ​​m)则ln(22 ​​m)=ln(t)=(t−1)−21​(t−1)2+31​(t−1)3−41​(t−1)4+51​(t−1)5−61​(t−1)6+71​(t−1)7+o(t−1)7(7)

根据 (6) 式，要想完成 ln ⁡ x \ln{x} lnx 的计算，我们要先知道 x x x 的阶码 j j j、尾数 m m m 和 ln ⁡ 2 \ln2 ln2 的值。其中 ln ⁡ 2 \ln2 ln2 为常数，约等于 0.693147 0.693147 0.693147，那么剩下就是要通过 x x x 取得 j j j 和 m m m。

我们已经知道了浮点数 x x x 在内存中的二进制存放方式，因此只需要通过简单的位运算，就可以分别获取 j j j 和 m m m 了。

首先，计算机是不支持 float 型数据的位运算的，因此我们要将其内存中的数据以 int 型数据的存放方式读取出来，然后再对其进行位运算。

以 int 型数据格式读取 float 型二进制数据的方式如下：

其中 x 是我们想要计算的数，设其值为 13.75；

i_x 是以 int 型数据读取的出的 x 的二进制数据，其十六进制的值为 0x415c0000，即：

为了取阶码，我们构造一个二进制 int 型数据 0x7f800000，其二进制表示为：

然后与 i_x 作 与 运算：

将 i_j 右移 23 位：

i_j 减去移码的 127，就得到 x 的阶码：

获取尾数的操作类似，我们先构造一个二进制 int 型数据 0x007fffff，其二进制表示为：

将其与 i_x 作 与 运算：

取到了尾数后，通过前面知道，浮点型的数据 x = m ⋅ 2 j − 127 x = m \cdot 2^{j-127} x=m⋅2j−127，而我们现在只取到尾数的二进制数据，需要将其构造出来，那么还需要将其阶码补齐。在此构造一个二进制 int 型数据 0x3f800000，其二进制表示如下：

将其与 i_m 作 或 运算：

再以 float 的格式将其读取出来：

后续就是将其进行缩放与泰勒展开，在此不再赘述。