计算机如何计算指数函数

计算机是如何计算指数函数 $e^{x}$ 的呢？

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

计算原理

泰勒展开

函数 $\exp{(x)} = e^x$ 的图像如下所示：

与计算对数函数 $\log{x}$ 相似，指数函数的计算也是对 $x$ 进行处理后，使用泰勒公式进行展开计算。

指数函数 $e^x$ 的泰勒展开如下：

$$e^x = \sum_{n = 0}^{n = \infty} \frac{1}{n!} x^n = 1 + x + \frac{1}{2!} x^2 + \frac{1}{3!} x^3 + \frac{1}{4!} x^4 +···+ \frac{1}{n!} x^n \quad(x \rightarrow 0， n\rightarrow \infty) \tag{1}$$

对 $e^x$ 进行展开也是有条件的，那便是要求 $x \rightarrow 0$ ，即以上的展开式只对 $x$ 在 0 点附近有效。

同时，展开式的精度和 $x$ 与 0 点的距离和展开的项数有关，如下：

以上是对 $e^x$ 分别进行两项、四项、六项展开的图像。各自的误差如下：

可以看出，在相同的展开项数下，$x$ 越靠近 0 点，误差越小；$x$ 固定时，展开项数越多，误差越小。

参数归约化

为了能够更加精确计算出更大范围的 $x$ 值的 $e^x$ 结果，首先要对原式进行处理。

我们知道，在计算机中，实数使用 IEEE 浮点数来表示，那么计算机其实更加善于计算 $2^x$，只要对浮点数的阶码进行处理就能计算得出 $2^x$ 的结果。

那么我们就要想办法，将以 $e$ 为底的指数函数 $e^x$ 转化为以 2 为底的指数函数 $2^x$的计算，至少要部分转化。

我们知道：

$$e = 2^{\log_2{e}} \tag{2}$$

于是：

$$e^x = (2^{\log_2{e}})^x = 2^{x\log_2{e}} \tag{3}$$

其中 $\log_2e \approx 1.442695$ 为常数，那么我们就成功地将计算 $e^x$ 的问题转换为计算 $2^{x\log_2e}$ 的问题了。

令 $t = x\log_2{e}$，则：

$$e^x = 2^{x\log_2{e}} = 2^{t}, \qquad (t= x^{\log_2{e}}) \tag{4}$$

但是问题又来了，通过构造阶码的形式来计算 2 的次幂，只能计算 2 的整数次幂，而 $t$ 为实数，那该怎么办呢？

这时要利用指数函数的一个性质：

$$a^{b+c} = a^b \cdot a^c \tag{5}$$

我们将 $t$ 分为整数部分 $m$ 与 小数部分 $n$，则：

$$令：t = [t] + (t-[t]) = m + n, (m = [t], n = t-[t]) \\ 则：2^{t} = 2^{m+n} = 2^{m} \cdot 2^{n} \tag{6}$$

对于整数次幂 $2^m$ 可通过构造阶码的形式来实现，对于小数次幂 $2^n$ ，则可对其进行再次转化为 $e^x$ 的形式，对其进行展开：

因为 $2 = e^{\ln2}$ ，所以:

$$2^n = (e^{\ln2})^n = e^{n\ln2} \tag{7}$$

因为 $n \in (-1,1),\ln2 \approx 0.693147$ ，所以 $n\ln2 \in (-0.693147, +0.693147)$ ，已经足够接近 0 点，所以可以对其使用泰勒展开进行计算。

令 $k = n\ln2,\quad k \in (-0.693147, +0.693147)$，则：

$$2^n = e^k, \quad (k = n\ln2) \tag{8}$$

那么：

$$\begin{align} e^x &= 2^{x\log_2{e}} \\ & = 2^t = 2^{m+n} \qquad (t = xlog_2{e}, m = [t], n = t - [t]) \\ & = 2^m \cdot 2^n \\ & = 2^m \cdot e^{n\ln2} = 2^m \cdot e^k \qquad (k = n\ln2, n \in (-\ln2, +\ln2)) \end{align}$$

对于 $e^k$，根据 (1) 式：

$$e^k = \sum_{n = 0}^{\infty} \frac{1}{n!} k^n = 1 + k + \frac{1}{2!}k^2 + \frac{1}{3!}k^3 + \frac{1}{4!}k^4 + \frac{1}{5!}k^5 + o(k^5)$$

实现思路

先将 $x$ 乘以 $\log_2{e} \approx 1.442695$ 得到 $t$：

float x = 3.14;
float t = x * 1.442695;    // 1.442695 = log2(e)

分别取得 $t$ 的整数部分 $m$ 和小数部分 $n$：

m = (round)((int)(t * 2.0) / 2);
float n = t - int(t);

通过 $m$ 来构造 $2^m$：

m = (m + 127) << 23; 

然后将 $n$ 乘以 $\ln2$ 得到 $k$：

float k = n * 0.693147;    // 0.693147 = ln(2)

再使用泰勒公式将其展开即可。

代码实现

源码

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

float my_exp(float x){
	float t = x * 1.442695;
	int m = 0;
	float n = 0;
	m = (round)((int)(t * 2.0) / 2);
	n = t - m;

	m = (m + 127) << 23; 
	n = n * 0.693147;
	float tmp_n = n;
	float j = 1;
	float result = 1 + n;
	for (int i = 2; i < 7; i++){
		j = j * i;
		tmp_n = tmp_n * n;
		result += tmp_n / j;
	}
	return *(float *)&m * result;
}

int main(int argc, char *argv[]) {
	
	my_exp(3.14);
	for (int i = -10; i <= 30; i++){
		float x = i * 0.1;
		printf("x = %2.6f, my log(x) = %2.6f, log(x) = %2.6f, err = %2.6f\n", x, my_exp(x), exp(x), fabs(my_exp(x) - exp(x)));
	}
	return 0;
}

运行结果