0%

如何实现 exp 函数

\(\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
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
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
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