计算机如何计算指数函数

计算机是如何计算指数函数 $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+\cdot \cdot \cdot +\frac{1}{n!}{x}^n(x\to 0，n\to \infty )\text{\hspace{1em}(1)}$$
对 $e^x$ 进行展开也是有条件的，那便是要求 $x\to 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}\text{\hspace{1em}(2)}$$
于是：
$$e^x=(2^{{\log }_{2}e}{)}^x=2^{x{\log }_{2}e}\text{\hspace{1em}(3)}$$
其中 ${\log }_{2}e\approx 1.442695$ 为常数，那么我们就成功地将计算 $e^x$ 的问题转换为计算 $2^{x{\log }_{2}e}$ 的问题了。

令 $t=x{\log }_{2}e$，则：
$$e^x=2^{x{\log }_{2}e}=2^t,(t=x^{{\log }_{2}e})\text{\hspace{1em}(4)}$$
但是问题又来了，通过构造阶码的形式来计算 2 的次幂，只能计算 2 的整数次幂，而 $t$ 为实数，那该怎么办呢？

这时要利用指数函数的一个性质：
$$a^{b+c}=a^b\cdot a^c\text{\hspace{1em}(5)}$$
我们将 $t$ 分为整数部分 $m$ 与 小数部分 $n$，则：
$$\begin{gathered} 令：t=[t]+(t-[t])=m+n,(m=[t],n=t-[t]) \\ 则：2^t=2^{m+n}=2^m\cdot 2^n\text{\hspace{1em}(6)} \end{gathered}$$
对于整数次幂 $2^m$ 可通过构造阶码的形式来实现，对于小数次幂 $2^n$ ，则可对其进行再次转化为 $e^x$ 的形式，对其进行展开：

因为 $2=e^{\ln 2}$ ，所以:
$$2^n=(e^{\ln 2}{)}^n=e^{n\ln 2}\text{\hspace{1em}(7)}$$
因为 $n\in (-1,1),\ln 2\approx 0.693147$ ，所以 $n\ln 2\in (-0.693147,+0.693147)$ ，已经足够接近 0 点，所以可以对其使用泰勒展开进行计算。

令 $k=n\ln 2,k\in (-0.693147,+0.693147)$，则：
$$2^n=e^k,(k=n\ln 2)\text{\hspace{1em}(8)}$$
那么：
$$\begin{array}{rl}{e}^x& =2^{x{\log }_{2}e}\\ & =2^t=2^{m+n}(t=xlo{g}_{2}e,m=[t],n=t-[t])\\ & =2^m\cdot 2^n\\ & =2^m\cdot e^{n\ln 2}=2^m\cdot e^k(k=n\ln 2,n\in (-\ln 2,+\ln 2))\end{array}$$
对于 $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$ 乘以 $\ln 2$ 得到 $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;
}

运行结果