\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 里面初始化一下一些变量：

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

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

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 次多项式出来：

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 多项式的好……

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 了，让我们来跑一下吧！

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$ 的一步。让我们看看这个算法表现如何！

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 了，代码如下：

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