immut/dense 设计

设计目标

DensePolynomial[A] 是作为值的一元多项式：一个可以自由共享、用 == 比较并用环运算符组合的规范系数向量。它面向系数大多非零的多项式，此时把直到次数为止的每个系数都存储下来，既是最简单的布局，也是最快的布局。

数学背景

$R$ 上的一元多项式是一个有限支撑序列 $(c_0, c_1, \dots)$ ，记作 $f = \sum_i c_i x^i$ 。其次数为 $\deg f = \max\{ i \mid c_i \neq 0 \}$ ，并约定 $\deg 0 = -\infty$ （以 None 返回）。加法逐系数进行，乘法是 Cauchy 乘积

(fg)_k = \sum_{i + j = k} f_i \, g_j .

若 $R$ 是交换环，则 $R[x]$ 也是；若 $R$ 只是半环（如 UInt），则 $R[x]$ 是半环。次数满足

\deg(f + g) \le \max(\deg f, \deg g), \qquad \deg(fg) \le \deg f + \deg g ,

其中第二式取等号当且仅当首项系数之积非零，当 $R$ 没有零因子时这总是成立。

设计决策

规范形式：修剪尾部零

问题。 序列 $(1, 2)$ 和 $(1, 2, 0, 0)$ 描述同一个多项式。如果两者都能被存储，==、degree、length 和哈希都必须先进行规范化。

选择。 每个构造器和每个操作在返回前都修剪尾部零，零多项式是空向量。修剪在多项式与满足 $c_{n-1} \neq 0$ 的字 $(c_0, \dots, c_{n-1})$ （外加空字）之间定义了一个双射，因此

\texttt{p == q} \iff p = q \text{ in } R[x], \qquad \texttt{p.length()} = \deg p + 1 .

在每个操作之后都修剪，而不只是在构造之后修剪，这一点很重要，因为上面的次数不等式可能是严格的。在 Int（即 $\mathbb{Z}/2^{32}$ ）中， $(2^{16} x + 1)^2 = 2^{32} x^2 + 2^{17} x + 1 = 2^{17} x + 1$ ：首项系数消失，修剪后的结果正确地具有次数 $1$ 。

向量是持久化的 @immut/vector.Vector，并且 from_coefficients 会复制其输入，因此调用者之后修改自己的数组不会改变多项式。

每个操作的最小系数约束

问题。 在整个类型上施加诸如 A : Ring 的单一约束会排除有用的系数类型：UInt 没有取负，而某些类型有乘法却没有单位元。

选择。 每个函数都声明它所用到的最小 luna-generic 能力集合。加法需要 Eq + AddMonoid（Eq 和 Zero 用于修剪），乘法再加上 Mul 但不需要 One，取负需要 Neg 但不需要 Mul，只有 variable、one、pow 和 karatsuba 需要 One。类型本身没有约束，因此 zero()、length() 和 degree() 对任意 A 都可用。

默认使用教科书乘法

* 用两层嵌套循环遍历存储的系数来计算 Cauchy 乘积，对长度 $m$ 和 $n$ 进行 $mn$ 次系数乘法与加法，然后修剪。对于典型用法中占主导的短多项式，这比任何递归方法都快，不需要 Neg，并且对精确系数类型是精确的。

Karatsuba 作为显式方法

问题。 对于长操作数， $O(mn)$ 会成为瓶颈。

推导。 在 $k$ 处拆分每个操作数： $a = a_0 + a_1 y$ ， $b = b_0 + b_1 y$ ，其中 $y = x^k$ 且 $\deg a_0, \deg b_0 < k$ 。则

\begin{aligned} ab &= a_0 b_0 + (a_0 b_1 + a_1 b_0)\, y + a_1 b_1 \, y^2, \\ a_0 b_1 + a_1 b_0 &= (a_0 + a_1)(b_0 + b_1) - a_0 b_0 - a_1 b_1 , \end{aligned}

其中第二行仅依赖分配律，因此它在任何环上都成立。三个一半规模的乘积代替了四个。记 $T(n)$ 为长度 $n$ 的代价， $c\,n$ 为加法和移位的代价，

T(n) = 3\,T(n/2) + c\,n \;\Rightarrow\; T(n) = c\,n \sum_{j=0}^{\log_2 n} \left(\tfrac{3}{2}\right)^j = O\!\left(n \cdot \left(\tfrac32\right)^{\log_2 n}\right) = O\!\left(n^{\log_2 3}\right) \approx O(n^{1.585}).

选择。 karatsuba 在较长长度的一半处拆分，对三个乘积递归，并在较短操作数至多有 $32$ 个系数时退回到 *，此时递归开销超过节省的乘法。其中的减法需要 Neg，按 $y$ 的移位是 scale(k, One::one())，需要 One。它是一个单独的方法，而不是 * 的实现，这样 * 可以保持较弱的约束和可预测的代价；有测试检查两者在阈值之上结果一致。

Horner 求值

$p(a)$ 按如下方式计算

p(a) = c_0 + a\bigl(c_1 + a\bigl(c_2 + \cdots + a\,(c_{n-1})\bigr)\bigr),

每个存储的系数一次乘法和一次加法，从零开始。在不对系数进行预处理的情况下，Horner 法则在求值一般多项式时使用了尽可能少的乘法。¹1 Ostrowski 证明了次数至多为 4 时的最优性（1954），Pan 证明了任意次数的情形（1966）：任何对次数为 $n$ 的一般多项式求值的算法都至少需要 $n$ 次乘法和 $n$ 次加法。它只需要 AddMonoid + Mul，因此可以在任何系数类型上求值，并且它从不构造幂 $a^i$ ，这对定宽整数能使中间值保持较小，对浮点数则能保持良态。

复合即在多项式处求值

substitute(q) 以多项式算术运行同样的 Horner 循环，计算 $p(q) = \sum_i c_i q^i$ 。对交换的 $R$ ，这就是求值同态 $\mathrm{ev}_q : R[x] \to R[x]$ ， $x \mapsto q$ ，即固定 $R$ 并把 $x$ 送到 $q$ 的唯一环同态。在单项式上，

\mathrm{ev}_q(x^i \cdot x^j) = q^{i+j} = q^i q^j = \mathrm{ev}_q(x^i)\,\mathrm{ev}_q(x^j),

并且两边都是双线性的，因此 $\mathrm{ev}_q(fg) = \mathrm{ev}_q(f)\,\mathrm{ev}_q(g)$ 且 $\mathrm{ev}_q(f + g) = \mathrm{ev}_q(f) + \mathrm{ev}_q(g)$ 。复合满足结合律 $p(q(r)) = (p(q))(r)$ ，因为两边都是在 $x$ 上取值一致的同态。

设 $\deg p = n$ ， $\deg q = d$ ，Horner 的第 $j$ 步把一个次数为 $jd$ 的多项式乘以 $q$ ，因此教科书乘法的代价为 $\sum_{j=1}^{n} (jd + 1)(d + 1) = O(n^2 d^2)$ ，结果的次数至多为 $nd$ 。

经由 ℕ 的典范映射求形式导数

derivative 计算 $D\bigl(\sum_i c_i x^i\bigr) = \sum_{i \ge 1} (i \cdot 1_R)\, c_i\, x^{i-1}$ ，用 NatHomomorphism::from_nat 把整数 $i$ 映射到 $R$ 中。 $D$ 是可加的，并满足 Leibniz 法则。在单项式上，

\begin{aligned} D(x^i x^j) &= (i + j)\, x^{i+j-1} \\ &= i\,x^{i-1} x^j + x^i\, j\,x^{j-1} \\ &= D(x^i)\,x^j + x^i\,D(x^j), \end{aligned}

并且 $D(fg) = D(f)\,g + f\,D(g)$ 两边关于 $(f, g)$ 都是双线性的，因此该法则推广到所有多项式。由于 $i$ 按特征取模，在特征 $p$ 下 $D(x^p) = p\,x^{p-1} = 0$ ，与代数中完全一致。luna-generic 为 Float、Double 和 BigInt 提供了 NatHomomorphism；定宽整数系数在获得这样的实例之前没有导数（上游正在用 FromNat 取代该 trait）。

用二分幂求幂

pow(e) 委托给包内部的 pow_nat，它在把 exp 减半的同时保持不变量 $\text{state} \cdot \text{factor}^{\text{exp}} = p^e$ 。它执行 $\lfloor \log_2 e \rfloor$ 次平方，并且对 $e$ 的每个置位至多向 state 乘入一次，因此至多进行 $2\lfloor \log_2 e \rfloor + 1$ 次多项式乘法。由于操作数不断增长，最后一次平方占主导：使用教科书乘法时，总代价为 $O\bigl((e\,\deg p)^2\bigr)$ 。

pow(0) 对每个 $p$ （包括零）都返回 one()，因此 $0^0 = 1$ 。这是 @arithmetic.PowNatChecked 的约定，DensePolynomial 通过返回 Ok(self.pow(e)) 来实现它。

用于容器的结构序

派生的 Compare 先按存储长度排序，然后从常数项开始逐个系数比较。由于修剪，长度等于次数加一，因此次数较低的多项式排在前面。这个序的存在是为了让多项式可以作为有序映射和有序集合的键；它不是环序，与 + 或 * 不相容。

正确性 / 不变量

规范形式。 每个公开操作之后，最后一个存储的系数非零，零为空。因此相等即多项式相等。
环定律。 当系数类型满足交换环定律时，DensePolynomial[A] 也满足：各运算就是规范形式上的教科书公式。性质测试检查加法单位元和规范化的幂等性。
一致性。 karatsuba(a, b) == a * b；substitute 就是 $\mathrm{ev}_q$ ；derivative 满足 Leibniz 法则；pow(e) 等于 $e$ 重乘积。
值语义。 没有方法会修改其接收者或参数；构造器会复制输入数组。
复杂度（长度为 $m, n$ ）：+、-、neg、scale 为 $O(m + n)$ ；* 为 $O(mn)$ ；karatsuba 对阈值之上的平衡操作数为 $O(n^{\log_2 3})$ ；eval 为 $O(n)$ ；substitute 为 $O(n^2 d^2)$ ；pow(e) 为 $O(\log e)$ 次乘法。

被否决的替代方案

在 * 内部使用 Karatsuba。 这会迫使每次乘法都要求 Neg + One，并改变短乘积的代价特征。保持其显式，把选择权留给调用者。
FFT 或数论变换乘法。 这些需要系数类型中有单位根或合适的模数，而泛型 A 并不提供。
单独存储次数并允许尾部零。 这在某些操作上省去一次修剪，但使相等和哈希依赖于规范化；每个操作修剪一次代价很小，并且消除了这个问题。
稀疏一元类型。 诸如 $x^{1000} + 1$ 的多项式在这里浪费空间；对它们请使用带单元素指数向量的 SparsePolynomial。

边界

仅限一元。多元多项式见 TermPolynomial 和 SparsePolynomial。
不提供除法、GCD、因式分解、求根或插值。
精确修剪：只有当系数 == 零时才会被丢弃。对于 Float 或 Double 系数，本应相消但带有舍入误差的首项系数会被保留，因此存储的次数可能超过数学上的次数。
定宽整数系数会回绕，因此多项式是 $\mathbb{Z}/2^k$ 上的，存在零因子以及如上所示的次数下降。
pow 接受 UInt 指数；没有负幂。

Ostrowski 证明了次数至多为 4 时的最优性（1954），Pan 证明了任意次数的情形（1966）：任何对次数为 $n$ 的一般多项式求值的算法都至少需要 $n$ 次乘法和 $n$ 次加法。 ↩