如何实现 exp 函数

Posted on In Disqus: Word count in article: 2.7k Reading time ≈ 10 mins.

\(\newcommand{\d}{\mathrm{d}}\)说到解析方法数质数，由于精度原因，如果要计算 \(\pi(10^{18})\) 的话，64 位浮点数已经不能满足要求了，所以我不得不实现自己的高精度。当然可以选用现成的库，例如著名的 MPFR，但是我觉得这样 overhead 太大了，毕竟我只需要 100 bit 就行了。

既然要实现浮点数，自然得实现一堆函数，其中有四个函数尤为重要：\(\exp, \log, \sin, \cos\)。当然，在实现之前，我们需要参考下别人的实现，例如 libm 。我们这里就拿 \(\exp\) 来举例，并假设我们现在已经支持基本的加减乘除了。

第一步应该是判断数的 regularity。由于只给自己用，所以我就 assume 不会出现 NaN、Inf 这种东西，这些判断也就都省了。

第二步是 normalize 输入：对于输入 \(x\)，我们总能找到一个 \(k\) 和 \(r\)，使得 \[ e^x = 2^k e^r, \quad k \in \mathbb{N}, |r| \leq r_\max= \frac{\ln 2}{2}. \] 这里 \(2^k\) 是很好表示出来的，剩下的就是如何计算 \(e^r\) 了。接下来我们就来讨论一下如何求 \(e^{r_\max}\)。

Truncated Taylor series

既然 \(|r|\) 不大，我们就直接在 \(r = 0\) 出做 Taylor expansion 怎么样？ \[ e^r = 1 + r + \frac{1}{2!} r^2 + \frac{1}{3!} r^3 + o(r^3) \] 我们就取前四项，看看效果如何。于是我先在 Sage 里面初始化一下一些变量：

1
2
3
4
5
6
7
8
9
10
11
prec = 100
deg = 3
r = var('r')
F = RealField(prec)["r"]
r_max = F(log(2)) / 2

def plot_diff(expr, label):
    p = plot(exp(r) - expr, -r_max, r_max, legend_label=f"$e^r - {label}$")
    p.set_legend_options(font_size=14, shadow=False, loc="lower right")
    return p

然后再来画一下取前四项时的误差：

1
2
3
4
5
poly_taylor = exp(r).series(r, deg + 1).truncate()
print(F(poly_taylor))

plot_diff(poly_taylor, r"\operatorname{Taylor}(r)")

结果如下：最大误差在 \(r = r_\max\) 时取得，约为 \(6.45 \times 10^{-4}\)。

hmmm，看上去还凑合？

Chebyshev interpoloation

如你所愿！

我们先定义（第一类） Chebyshev 多项式如下： \[ \begin{aligned} T_0(x) & = 1, \\ T_1(x) & = x, \\ T_{n+1}(x) & = 2x^2 T_n(x) - T_{n-1}(x). \end{aligned} \] 另外一种理解方式是，把 \(\cos n\theta\) 看成是 \(\cos \theta\) 的一个多项式，则 \[ T_n(\cos \theta) = \cos n\theta. \] 从这个定义可以直接得到：\(T_n(x)\) 的第 \(1 \leq k \leq n\) 个根为 \(x_i = \frac{k - \frac{1}{2}}{n} \pi\)。我们考虑令 \(n = 4\)，然后用 \(\{(x_i, \exp(x_i))\}_{i \in [n]}\) 插值插出一个 3 次多项式出来：

1
2
3
4
5
6
7
cheby_nodes_a = [cos((k + 1/2) / (deg + 1) * pi) for k in range(deg + 1)]
values = [exp(r * r_max) for r in cheby_nodes_a]

poly_cheby_interp_a = F.lagrange_polynomial(list(zip(cheby_nodes_a, values)))
print(poly_cheby_interp_a)

plot_diff(poly_cheby_interp_a(r=r/r_max), r"\operatorname{ChebyInterpA}(r)")

结果如下：在 \(r = r_\max\) 时误差最大，为 \(8.10 \times 10^{-5}\)。

这个多项式的系数和 Taylor expansion 的系数有较大不同，但是误差却小很多！这说明了 Taylor expansion 得到的多项式不一定是最优的！

现在问题来了：Chebyshev interpolation 到底魔改了些啥，使得它能比 Taylor expansion 给出的多项式更优？从图上我们可以看出来：Taylor expansion 的误差集中在两端，而 Chebyshev interpolation 的误差则分布较为均匀，把 \(r_\max\) 对应的误差减小了，再把某些 \(r\) 上的误差增大了，但是由于这些点上原来的误差很小，增大了也没啥大事，相当于把误差从原来的两极分化均摊到每个 \(r\) 上了。

有人可能就要问了：为啥我们要选 Chebyshev 多项式的根作为 \(x\) 来 interpolate 呢？这就不得不提到著名的 Runge's Phenomenon 了，简单说来就是如果选的 \(x\) 不好的话（例如 \([-1, 1]\) 内均匀选点），那么有可能选的点越多拟合越差。这个链接稍微揭示了为什么 Chebyshev 多项式的根就不会出现这个问题。当然这里除了 Chebyshev 多项式的根外，还可以考虑 Legendre polynomial 的根。这两类多项式在数值计算中用的挺多的，特别是数值积分。最流行的两种数值积分方式，Newton-Legendre 方法和 Clenshaw-Curtis 方法，就分别基于 Legendre polynomial 和 Chebyshev polynomial。

另外我这里用的是 ChebyInterpA，第一类 Chebyshev nodes 也就是第一类 Chebyshev 多项式的根。事实上还有第二类 Chebyshev 多项式，对应的有第二类 Chebyshev nodes，分别为 \(x_k = \cos \frac{k}{n} \pi\) ，其中 \(0 \leq k \leq n\)。我试了一下，效果没有第一类 Chebyshev 多项式的好……

1
2
3
4
5
6
7
cheby_nodes_b = [cos(k / deg * pi) for k in range(deg + 1)]
values = [exp(r * r_max) for r in cheby_nodes_b]

poly_cheby_interp_b = F.lagrange_polynomial(list(zip(cheby_nodes_b, values)))
print(poly_cheby_interp_b)

plot_diff(poly_cheby_interp_b(r=r/r_max), r"\operatorname{ChebyInterpB}(r)")

结果如下：很显然，由于在 -1 和 1 处选了点，所以这两个点的误差均为 0，但是 0 处的误差就大很多了……

Truncated Chebyshev series

如你所愿！

刚才的两种方法都是基于 Chebyshev 多项式的。事实上，Chebyshev 多项式在 \(w(x) = \frac{1}{\sqrt{1 - x^2}}\) 下是正交的： \[ \int_{-1}^{1} T_n(x) T_m(x) w(x) \d x = \begin{cases} 0, & n \neq m, \\ \pi, & n = m = 0, \\ \frac{\pi}{2}, & n = m \neq 0. \\ \end{cases} \] 于是我们可以把 \(T_n(x)\) 当成是一组基来表达函数 \(f(x)\)： \[ f(x) = \frac{a_0}{2} + \sum_{k=1}^\infty a_k T_k(x), \quad a_k := \frac{2}{\pi} \int_{-1}^1 f(x) T_k(x) w(x) \d x. \] 但是求 \(a_k\) 似乎需要积分……但是没问题，我们变量代换 \(x = \cos \theta\) 后用均匀撒点即可： \[ a_k = \frac{2}{\pi}\int_{-1}^1 f(x) T_k(x) w(x) \d x = \frac{2}{\pi}\int_0^\pi f(\cos \theta) \cos k \theta \d \theta \approx \frac{2}{N} \sum_{j=1}^N f\left( \frac{j - \frac{1}{2}}{N} \pi \right) \cos\left( \frac{k (j - \frac{1}{2})}{N} \pi \right). \] 只要 \(N\) 足够大，\(a_k\) 就能足够好。实际上，随着 \(N\) 的增大，\(a_k\) 收敛地相当快，几乎是指数级别的。\(N = 100\) 时就能得到 \(a_2\) 的 100 位有效数字。

好了，你现在已经学会了 Chebyshev series 了，让我们来跑一下吧！

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
d = 50

poly_cheby_series = 0

for j in range(deg+1):
    c = 0
    for k in range(d):
        t = exp(cos(pi * (k + 1/2) / d) * r_max) * cos(pi * j * (k + 1/2) / d)
        c += F(t)
    c = c * 2 / d
    if j == 0:
        poly_cheby_series += c / 2
    else:
        poly_cheby_series += c * chebyshev_T(j, r)

print(poly_cheby_series)

plot_diff(poly_cheby_series(r=r/r_max), r"\operatorname{ChebySeries}(r)")


结果如下：

Q：咦，这和之前的图不是一模一样吗……

A：不会吧不会吧，不会真的有人看不出这个多项式在 \(r = r_\max\) 的时候误差仅有 \(7.83 \times 10^{-5}\) 吧？（手动狗头）

Minimax approximation 和 Remez algorithm

如你所愿！

这次我们一步到位，直接找出最优的多项式！首先我们祭出以下定理：

[Theorem 10.1, ATAP, informal] 若一个函数 \(f: [-1, 1] \to \mathbb{R}\) 连续，则其最优 \(n\) 次多项式逼近 \(p_n^* := \arg\min_{p: \text{deg}(p)\leq n} \|f - p\|_\infty\) 唯一，且 \(f - p^*_n\) 有至少 \(n+2\) 个极值点，这些极值点的值的绝对值相等。

代码我就不放了，因为是用 Julia 写的 Remez 算法……总之我们跑出来结果如下：

从这张图中我们也在来理解一下刚才的定理：多项式次数为 3，则 \(g(r) = e^r - \text{Remez}(r)\) 有 5 个极值点，分别为 \(x_1 \approx -r_\max, x_2 \approx -0.24, x_3 \approx 0, x_4 \approx 0.24, x_5 \approx r_\max\)，且 \(g(x_1) = -g(x_2) = g(x_3) = -g(x_4) = g(x_5)\)。之前用 Chebyshev 多项式搞出的多项式已经接近于这个性质。

虽然 Remez 算法给出来的结果很接近 truncated Chebyshev series，但是从系数上来说它更接近于 truncated Taylor series 的结果。

自然，一个很有趣的问题就是：用 Chebyshev polynomial 得到的多项式有多好？这里给出一个定理：

[Theorem 16.1, ATAP, informal] 函数 \(f(x)\) 在 \([-1, 1]\) 上连续。令 \(f_n(x)\) 为 \(f\) 的 \(n\) 次 truncated Chebyshev series，\(p_n(x)\) 为 \(f\) 的 \(n\) 次 Chebyshev interpolant，则有： \[ \begin{aligned} \|f - p_n\|_\infty & \leq \left( 2 + \frac{2}{\pi} \log (n+1) \right) \|f - p^*_n\|_\infty, \\ \|f - f_n\|_\infty & \leq \left( 4 + \frac{4}{\pi^2} \log (n+1) \right) \|f - p^*_n\|_\infty, \end{aligned} \]

也就是说，Chebyshev interpolation 和 truncated Chebyshev series 给出的结果几乎是最优的了……

变量代换

如你所愿！

刚才瞎扯扯了好久，改进都不大，这次我们整点有用的！

不知道是否还有人记得我们最开始的目标：我们是想求 \(e^r\) 来着。可是，我们真的想直接求 \(e^r\) 吗？

我们先求 \(u = r \frac{e^r + 1}{e^r - 1}\)，有了 \(u\) 后可以由 \(e^r = \frac{u + r}{u - r}\) 来求出 \(e^r\)。求 \(u\) 有什么好处呢？注意到 \(u\) 是一个偶函数，所以对 \(u\) 做 Taylor expansion 后，就只有偶数项了： \[ u = r \frac{e^r + 1}{e^r - 1} = 2 + \frac{r^2}{6} - \frac{r^4}{360} + \frac{r^6}{15120} + o(r^6). \] 虽然这个多项式次数比之前高，但是运算量却只大了一点点，只多了算出 \(r^2\) 的一步。让我们看看这个算法表现如何！

1
2
3
4
5
6
7
8
def u(x):
    return x * (exp(x) + 1) / (exp(x) - 1) if x != 0 else 2

poly_taylor2 = (r * (exp(r) + 1) / (exp(r) - 1)).series(r, 2 * (deg + 1)).truncate()
print(poly_taylor2)

plot_diff((poly_taylor2(r) + r) / (poly_taylor2(r) - r),
          r"\frac{\operatorname{Taylor2}(r) + r}{\operatorname{Taylor2}(r) - r}")

结果如下：

AAAmazing！手动微笑。之前 Chebyshev 多项式什么的都是渣渣……

变量代换 + Chebyshev interpolation

如你所愿！

虽然刚才喷了一发 Chebyshev 多项式，但是毕竟是另一个优化方向。我图方便就只试了 Chebyshev interpolation 了，代码如下：

1
2
3
4
5
6
7
8
cheby_nodes_a = [cos((k + 1/2) / (2 * deg + 1) * pi) for k in range(2 * deg + 1)]
values = [u(r * r_max) for r in cheby_nodes_a]

poly_cheby_interp_a = F.lagrange_polynomial(list(zip(cheby_nodes_a, values)))
print(poly_cheby_interp_a)

plot_diff((poly_cheby_interp_a(r=r/r_max) + r) / (poly_cheby_interp_a(r=r/r_max) - r),
          r"\frac{\operatorname{ChebyInterpA2}(r) + r}{\operatorname{ChebyInterpA2}(r) - r}")

结果为：

误差最大为 \(1.32 \times 10^{-12}\) 左右。图中可以看出误差的分布并不均匀，这是因为误差均匀的 \(u = r \frac{e^r + r}{e^r - r}\) 不代表 \(\frac{u + r}{u - r}\) 也误差均匀，但总之，这个方法比一开始的 Taylor expansion 展开 \(e^r\) 不知好了多少。libm 用了一个 10 次的多项式来近似 \(u\)，误差保证不超过 \(2^{-59}\) 了。它号称用了 Remez 算法来计算系数。

需要注意的一点是，无论是 Chebyshev 多项式还是 Remez 算法，优化的都是绝对误差，但是在浮点数计算过程中，相对误差可能更重要。但是可喜可贺的是，无论是直接计算 \(e^r\) 还是先计算 \(u\)，在 \([-r_\max, r_\max]\) 内 \(e^r\) 和 \(u\) 不小于某个常数，所以绝对误差也能 bound 住相对误差，但是对于另外的函数可能就不是这样了……

Reference

  • https://www.johndcook.com/blog/2020/03/11/chebyshev-approximation/
  • [ATAP] : Approximation Theory and Approximation Practice