poly 设计

本页解释为什么在对偶数上对 luna-poly 多项式求值就能得到其导数、稠密和稀疏两种求值顺序分别计算什么、结果的精度如何，以及为什么桥接仅限于单变量。

设计目标

复用 luna-poly 已有的求值算法，在不引入新的多项式表示、也不进行符号运算的情况下，为 luna-poly 用户提供多项式在某点处的值和导数。

数学背景

在对偶数上求值即得导数

对于交换半环 $R$ 上的 $p(x) = \sum_{k=0}^{n} c_k x^k$ ，dual 设计表明

p(a + b\varepsilon) = p(a) + p'(a)\,b\,\varepsilon, \qquad p'(x) = \sum_{k=1}^{n} k\,c_k x^{k-1},

其中 $k\,c_k$ 表示 $c_k$ 自加 $k$ 次。该恒等式是代数性的，因此对整数多项式精确成立，对浮点多项式在舍入意义下成立；它既不需要极限，也不需要除法。

对偶数上的 Horner 法则

DensePolynomial::eval 计算 $r_{n+1} = 0$ ，对 $k = n, \dots, 0$ 计算 $r_k = r_{k+1}\,x + c_k$ ，并返回 $r_0 = p(x)$ 。取 $x = a + b\varepsilon$ 、提升后的系数 $c_k + 0\varepsilon$ 以及 $r_k = p_k + q_k\varepsilon$ ，由对偶数的积与和得到

\begin{aligned} p_k &= p_{k+1}\,a + c_k, \\ q_k &= p_{k+1}\,b + a\,q_{k+1}, \qquad p_{n+1} = q_{n+1} = 0 . \end{aligned}

当 $b = 1$ 时，这就是同时计算多项式及其导数的经典方案。¹1 D. E. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., section 4.6.4. 由归纳可得 $p_k = \sum_{i \ge k} c_i a^{i-k}$ 以及 $q_k = b\sum_{i > k} (i - k)\,c_i a^{i-k-1}$ ，因此 $p_0 = p(a)$ ， $q_0 = p'(a)\,b$ 。

通过求幂计算稀疏项

SparsePolynomial::eval 对存储的各项求 $c\,x^e$ 之和，并用二进制快速幂计算 $x^e$ 。在对偶数上，每次乘法都应用乘积法则，因此 $x^e$ 变为 $a^e + e\,a^{e-1} b\,\varepsilon$ （即 dual 设计中的二项式恒等式），各项之和同样给出 $p(a) + p'(a)\,b\,\varepsilon$ 。

设计决策

求值而非符号求导

问题。 用户需要 $p'(x)$ 在若干点处的值。

选项。 用 DensePolynomial::derivative 构建导数多项式再求值；或在 Dual[T] 上对 $p$ 求值。

选择。 在 Dual[T] 上求值。它一遍就返回 $p(x)$ 和 $p'(x)$ ，只需要 Semiring（luna-poly 中的形式导数还需要 NatHomomorphism 来构成 $k\,c_k$ ），不会为每一阶导数额外分配一个多项式，并且可以通过 eval_dual 与其他对偶计算组合。当需要导数多项式本身时，形式导数仍是合适的工具；测试套件会检查两者是否一致。

复用 `luna-poly` 的求值

桥接用 Dual::constant 提升系数，用普通构造函数构建 DensePolynomial[Dual[T]] 或 SparsePolynomial[Dual[T]]，然后调用 eval。它并不重新实现 Horner 法则，因此 luna-poly 求值方式的任何改动（以及它对零系数的规范化，这也是需要 T : Eq 的原因）都会同样作用于这里。

稀疏多项式仅限单变量

问题。 SparsePolynomial 是多元的，但偏导数需要选定变量，而变量的身份存在于 luna-poly 的上下文中（VariableContext、ContextPolynomial）。

选择。 稀疏桥接用单一赋值 [x] 求值，并命名为 sparse_univariate_*。使用了其他变量的多项式会使 eval 中止。多元 API 需要接受变量上下文并返回梯度，这留待将来实现。持有上下文多项式的调用方可以先对其他变量做部分求值，正如本仓库的集成测试那样。

兼容名称

首个版本将这些函数命名为 derivative_at、value_and_derivative_at、eval_sparse_dual、sparse_derivative_at 和 sparse_value_and_derivative_at。dense_ 和 sparse_univariate_ 名称表明了所针对的表示形式和导数种类；旧名称作为普通别名保留，以便现有代码继续工作。

正确性与不变量

精确环中的精确性

对于具有精确算术的 T（Int、BigInt、精确有理数），由上述恒等式可知，结果恰好就是 $p(x)$ 和 $p'(x)$ 。

稠密求值的舍入误差

使用 $\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta)$ ， $|\delta| \le u$ ，以及乘积引理 $\prod (1 + \delta_i)^{\pm1} = 1 + \theta_k$ ， $|\theta_k| \le \gamma_k = ku/(1 - ku)$ 。²2 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, Lemma 3.1 and chapter 5.

值。根据投影同态， $\hat p_0$ 就是普通 Horner 法则的结果。每一步都把累加和乘以 $(1 + \delta_{\times})(1 + \delta_{+})$ ，并加上经过一次舍入的 $c_k$ ，因此当 $k < n$ 时 $c_k a^k$ 带有 $2k + 1$ 个因子；首项系数带有 $2n$ 个，因为第一步 $0 \cdot a + c_n$ 是精确的：

\hat p_0 = \sum_{k=0}^{n} c_k a^k (1 + \theta^{(k)}), \quad |\theta^{(k)}| \le \gamma_{2n}, \qquad |\hat p_0 - p(a)| \le \gamma_{2n} \sum_{k=0}^{n} |c_k|\,|a|^k .

导数（ $b = 1$ ）。 切向分量的每一步计算 $\mathrm{fl}\big(\hat p_{k+1} \cdot 1 + \mathrm{fl}(a\,\hat q_{k+1})\big)$ ，而加上常数 $c_k$ 只会给切向分量加上一个精确的零。导数包含 $k$ 份 $c_k a^{k-1}$ ，每一份对应一个步骤 $j < k$ ，在该步 $c_k$ 从值链转入切向链。沿这样一条路径， $c_k$ 在进入时舍入一次，在值的第 $k - 1, \dots, j + 1$ 步中每步舍入两次，转入切向链时舍入一次，在切向的第 $j - 1, \dots, 0$ 步中每步舍入两次：共 $1 + 2(k - 1 - j) + 1 + 2j = 2k$ 个因子。因此

|\hat q_0 - p'(a)| \le \gamma_{2n} \sum_{k=1}^{n} k\,|c_k|\,|a|^{k-1} ,

这与用 Horner 法则计算形式导数时的误差界相同。除非各项相互抵消，即除非 $p'$ 在 $a$ 处是病态的，否则相对误差都很小。

代价

稠密： $n + 1$ 步对偶乘加，即 T 上约 $3n$ 次乘法和 $3n$ 次加法，外加一个存放 $n + 1$ 个提升后系数的数组。稀疏：每一项一次二进制快速幂，即 $O(\log e)$ 次对偶乘法，外加一次对偶乘法和一次加法。

被否决的方案

先符号求导再求值。 需要两遍计算、一个额外的多项式，以及对 T 更强的约束；为需要导数多项式的用户，它以 luna-poly 自身的 DensePolynomial::derivative 形式保留。
私有的 Horner 循环。 速度只快一个常数倍，但会重复 luna-poly 的语义，并逐渐与之偏离。
为多元稀疏多项式猜测变量。 静默地对”第一个变量”求偏导很容易被误用；桥接会直接中止。

边界

仅支持一元：不支持多元多项式的偏导数或梯度，也不支持 ContextPolynomial。
仅支持在点处求一阶导数；不生成导数多项式（请使用 luna-poly），也不支持高阶导数。
没有带检查的变体：除了上面所述的稀疏中止情况外，稠密和稀疏求值都不会失败。
仅支持稠密和稀疏的 immut 表示。

D. E. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., section 4.6.4. ↩
N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, Lemma 3.1 and chapter 5. ↩