计算机是如何计算对数函数 $\log_a{x}$ 的呢？

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

计算原理

泰勒展开

提起对数函数 $\log_a{x}$，大家最熟悉的对数函数应该是以自然底数 $e$ 为底的对数函数：$\ln{x}$，即 $\log_e{x}$。

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

对于不是以 $e$ 为底的对数函数 $\log_a{x}$ ，可以使用换底公式，将其转换为以 $e$ 为底的函数形式。

$换底公式：\log_a{b} = \frac{\log_c{b}}{\log_c{a}} \qquad\qquad\qquad$

根据换底公式：

$\log_a{x} = \frac{\log_e{x}}{\log_e{a}} = \frac{\ln{x}}{\ln{a}} \tag{1} \qquad\qquad\qquad$

于是，计算对数函数 $\log_a{x}$ 的问题就转换为计算 $\ln{x}$ 的问题。

对于计算 $\ln{x}$，可能我们第一时间想到的就是泰勒展开公式：

$\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}$ 展开是有条件的，即 $x$ 所属的区间为 $x \in (0, 2]$ 。而且，展开的项数一定是固定的，毕竟不能无限次的进行展开计算。

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

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

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

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

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

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

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

单精度存储

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

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

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

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

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

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

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

以 13.75 为例：

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

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

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

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

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

归约化

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

借助对数函数的另外两个性质：

$\log_a{(b \cdot c)} = \log_a{(b)} + \log_a{(c)} \\ \log_a{(b^j)} = j \cdot \log_a{(b)} \tag{3}$

可得：

$\ln{x} = \ln{[(1.m) \cdot 2^{j-127}]} = \ln{(2^{j-127})} + \ln{(m)} = (j-127) \cdot \ln{2} + \ln{m} \tag{4}$

在 (4) 式中，$(j-127)$ 为 $x$ 的阶码，可以想办法取；$\ln2 \approx 0.693147$ 为常数，因此，只需要计算 $\ln{m}$ 的值即可。

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

为了展开的精度更加精确，我们应该尽量让 $m$ 趋向于 $1$，根据 (3) 式：

$\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}$

因为 $m \in(1,2)$，所以 $\frac{\sqrt{2}}{2} m \in(\frac{\sqrt{2}}{2},\sqrt{2})$，即 $m$ 的大致范围为 $(0.707106, 1.414213)$，更加接近于 $1$。

于是，我们最终的式子就是：

$\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}$

对于 $\ln{(\frac{\sqrt{2}}{2}m)}$，使用泰勒公式将其展开：

$令 \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}$

实现思路

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

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

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

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

1 2	float x = 13.75; int i_x = (int )&x; // i_x = 0x415c0000

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

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

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

然后与 i_x 作与运算：

1	int i_j = i_x & 0x7f800000; // i_j = 0x41000000

将 i_j 右移 23 位：

1	i_j = i_j >> 23; // i_j = 0x00000082 = 130

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

1	i_j = i_j - 127; // i_j = 3

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

将其与 i_x 作与运算：

1	int i_m = i_x & 0x007fffff; // i_m = 0x005c0000

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

将其与 i_m 作或运算：

1	i_m = i_m \| 0x3f800000; //i_m = 0x7fdc0000

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

1	float x_m = (float )&i_m;

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

代码实现

源码

#include<stdio.h>
#include <math.h>

float my_log(float x){
	float ln_2 = 0.693147;		// ln(2) = 0.693147
	int i_x = *(int*)&x;
	int x_j = ((i_x & 0x7f800000) >> 23) - 127;
	int x_m = (i_x & 0x007fffff) | 0x3f800000;
	float m = *(float*)&x_m;
	m *= 0.707106;				// sqrt(2) / 2 = 0.707106
	m -= 1.0;
	float t_m = m * m;
	float k = 1.0;
	float result = m;
	for (int i = 2; i < 8; i++){
		k = -1.0 * k;
		result += k * t_m / (1.0 * i);
		t_m *= m;
	}
	return x_j * ln_2 + 0.5 * ln_2 + result;
}
int main(){
	for (int i = 1; i <= 30; i++){
		float x = i * 0.2;
		printf("x = %2.6f, my log(x) = %2.6f, log(x) = %2.6f, err = %2.6f\n", x, my_log(x), log(x), fabs(my_log(x) - log(x)));
	}
	return 0;
}

误差

运行结果如下