译自 halo2 book：background/polynomials.md

多项式(polynomial)

设 $A(X)$ 是 ${\mathbb{F}}_p$ 上以形式不定元(formal indeterminate) $X$ 表示的多项式。举例来说，

$$A(X)=a_0+a_{1}X+a_2{X}^2+a_3{X}^3$$

定义了一个 $3$ 次多项式。$a_0$ 称为常数项(constant term)。$n-1$ 次多项式有 $n$ 个系数。我们常常想要计算把形式不定元 $X$ 替换成某个具体值 $x$ 的结果，记作 $A(x)$。

在数学中这通常被称为"在点 $x$ 处求值 $A(X)$"。这里"点(point)"一词源自多项式 $y=A(x)$ 形式的几何用法，其中 $(x,y)$ 是二维空间中某点的坐标。然而，我们处理的多项式几乎总是被约束为等于零，并且 $x$ 将是 某个域的元素。这不应与 椭圆曲线(elliptic curve) 上的点混淆，我们也会用到椭圆曲线，但绝不会在多项式求值的语境中使用。

重要说明：

• 多项式相乘产生的乘积多项式，其次数是各因子次数之和。多项式除法则从次数中作减法。 $$\mathrm{deg}(A(X)B(X))=\mathrm{deg}(A(X))+\mathrm{deg}(B(X)),$$ $$\mathrm{deg}(A(X)/B(X))=\mathrm{deg}(A(X))-\mathrm{deg}(B(X)).$$
• 给定一个 $n-1$ 次多项式 $A(X)$，如果我们在不同的点上得到该多项式的 $n$ 个求值，那么这些求值就完美地确定了该多项式。换言之，给定这些求值，我们可以通过多项式插值(interpolation)得到唯一的 $n-1$ 次多项式 $A(X)$。
• $[a_0,a_1,\cdots ,a_{n-1}]$ 是多项式 $A(X)$ 的 系数表示(coefficient representation)。等价地，我们可以使用它在 $n$ 个不同点上的 求值表示(evaluation representation) $$[(x_0,A(x_0)),(x_1,A(x_1)),\cdots ,(x_{n-1},A(x_{n-1}))]$$ 两种表示都唯一地确定同一个多项式。

（旁注）霍纳法则(Horner's rule)

霍纳法则允许高效地求一个 $n-1$ 次多项式的值，只需 $n-1$ 次乘法和 $n-1$ 次加法。它就是下面这个恒等式： $$\begin{array}{rl}{a}_0& +a_{1}X+a_2{X}^2+\cdots +a_{n-1}{X}^{n-1}\\ & =a_0+X(a_1+X(a_2+\cdots +X(a_{n-2}+X{a}_{n-1})\left)\right),\end{array}$$

快速傅里叶变换(Fast Fourier Transform, FFT)

FFT 是在多项式的系数表示与求值表示之间转换的一种高效方法。它在 $n$ 次单位根 $\{{\omega }^0,{\omega }^1,\cdots ,{\omega }^{n-1}\}$ 处对多项式求值，其中 $\omega$ 是本原 $n$ 次单位根。通过利用单位根中的对称性，FFT 的每一轮都把求值问题缩减为一个只有原来一半大小的问题。最常见的是，我们使用长度为 2 的某个幂次 $n=2^k$ 的多项式，并递归地应用这个减半缩减。

动机：快速多项式乘法

在系数表示下，把两个多项式相乘 $A(X)\cdot B(X)=C(X)$ 需要 $O(n^2)$ 次运算：

$$\begin{array}{rl}A(X)& =a_0+a_{1}X+a_2{X}^2+\cdots +a_{n-1}{X}^{n-1},\\ B(X)& =b_0+b_{1}X+b_2{X}^2+\cdots +b_{n-1}{X}^{n-1},\\ C(X)& =a_0\cdot (b_0+b_{1}X+b_2{X}^2+\cdots +b_{n-1}{X}^{n-1})\\ & +a_{1}X\cdot (b_0+b_{1}X+b_2{X}^2+\cdots +b_{n-1}{X}^{n-1})\\ & +\cdots \\ & +a_{n-1}{X}^{n-1}\cdot (b_0+b_{1}X+b_2{X}^2+\cdots +b_{n-1}{X}^{n-1}),\end{array}$$

其中第一个多项式中的 $n$ 项每一项都必须乘以第二个多项式的 $n$ 项。

然而在求值表示下，多项式乘法只需要 $O(n)$ 次运算：

$$\begin{array}{rl}A& :\{(x_0,A(x_0)),(x_1,A(x_1)),\cdots ,(x_{n-1},A(x_{n-1}))\},\\ B& :\{(x_0,B(x_0)),(x_1,B(x_1)),\cdots ,(x_{n-1},B(x_{n-1}))\},\\ C& :\{(x_0,A(x_0)B(x_0)),(x_1,A(x_1)B(x_1)),\cdots ,(x_{n-1},A(x_{n-1})B(x_{n-1}))\},\end{array}$$

其中每个求值是逐点(pointwise)相乘的。

这提示了如下快速多项式乘法的策略：

1. 在全部 $n$ 个点处对多项式求值；
2. 在求值表示中执行快速逐点乘法（$O(n)$）；
3. 转换回系数表示。

现在的挑战是如何高效地 求值(evaluate) 和 插值(interpolate) 多项式。朴素地说，在 $n$ 个点处对一个多项式求值需要 $O(n^2)$ 次运算（我们在每个点处使用 $O(n)$ 的霍纳法则）：

$$\left[\begin{array}{c}A(1)\\ A(\omega )\\ A({\omega }^2)\\ ⋮\\ A({\omega }^{n-1})\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& \dots & 1\\ 1& \omega & {\omega }^2& \dots & {\omega }^{n-1}\\ 1& {\omega }^2& {\omega }^{2\cdot 2}& \dots & {\omega }^{2\cdot (n-1)}\\ ⋮& ⋮& ⋮& & ⋮\\ 1& {\omega }^{n-1}& {\omega }^{2(n-1)}& \cdots & {\omega }^{(n-1{)}^2}\end{array}\right]\cdot \left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_{n-1}\end{array}\right].$$

为方便起见，我们将上述矩阵记为： $$\hat{\mathbf{A}}={\mathbf{V}}_{\omega }\cdot \mathbf{A}.$$

（$\hat{\mathbf{A}}$ 被称为 $\mathbf{A}$ 的 离散傅里叶变换(Discrete Fourier Transform)；${\mathbf{V}}_{\omega }$ 也称为 范德蒙德矩阵(Vandermonde matrix)。）

（基-2）Cooley-Tukey 算法

我们的策略是把一个大小为 $n$ 的 DFT 分成两个交织(interleaved)的大小为 $n/2$ 的 DFT。给定多项式 $A(X)=a_0+a_{1}X+a_2{X}^2+\cdots +a_{n-1}{X}^{n-1},$ 我们把它拆分成偶次项和奇次项：

$$\begin{array}{rl}{A}_{\text{even}}& =a_0+a_{2}X+\cdots +a_{n-2}{X}^{\frac{n}{2}-1},\\ A_{\text{odd}}& =a_1+a_{3}X+\cdots +a_{n-1}{X}^{\frac{n}{2}-1}.\end{array}$$

为恢复原始多项式，我们做 $A(X)=A_{\text{even}}(X^2)+X{A}_{\text{odd}}(X^2).$

在点 ${\omega }_n^i$ 和 ${\omega }_n^{\frac{n}{2}+i}$（$i\in [\mathrm{0..}\frac{n}{2}-1]$）上试一下，我们开始注意到一些对称性：

$$\begin{array}{rl}A({\omega }_n^i)& =A_{\text{even}}(({\omega }_n^i{)}^2)+{\omega }_n^i{A}_{\text{odd}}(({\omega }_n^i{)}^2),\\ A({\omega }_n^{\frac{n}{2}+i})& =A_{\text{even}}(({\omega }_n^{\frac{n}{2}+i}{)}^2)+{\omega }_n^{\frac{n}{2}+i}{A}_{\text{odd}}(({\omega }_n^{\frac{n}{2}+i}{)}^2)\\ & =A_{\text{even}}((-{\omega }_n^i{)}^2)-{\omega }_n^i{A}_{\text{odd}}((-{\omega }_n^i{)}^2)\leftarrow \text{(negation lemma)}\\ & =A_{\text{even}}(({\omega }_n^i{)}^2)-{\omega }_n^i{A}_{\text{odd}}(({\omega }_n^i{)}^2).\end{array}$$

注意我们只在半个定义域 $\{({\omega }_n^0{)}^2,({\omega }_n{)}^2,\cdots ,({\omega }_n^{\frac{n}{2}-1}{)}^2\}=\{{\omega }_{n/2}^i\},i=[\mathrm{0..}\frac{n}{2}-1]$（减半引理）上对 $A_{\text{even}}(X)$ 和 $A_{\text{odd}}(X)$ 求值。这给了我们重构 $A(X)$ 在整个定义域 $\{{\omega }^0,\omega ,\cdots ,{\omega }^{n-1}\}$ 上所需的所有项：这意味着我们已经把一个长度为 $n$ 的 DFT 转换成了两个长度为 $\frac{n}{2}$ 的 DFT。

我们选择 $n=2^k$ 为 2 的某个幂次（必要时通过补零），并递归地应用这个分治策略。由主定理(Master Theorem)¹，这给了我们一个具有 $O(n{\log }_{2}n)$ 次运算的求值算法，也就是快速傅里叶变换(FFT)。

逆 FFT(Inverse FFT)

至此我们已经对多项式求了值并逐点相乘。剩下的就是把乘积从求值表示转换回系数表示。为此，我们只需对求值表示再调用一次 FFT。不过这一次我们还要：

• 在范德蒙德矩阵中把 ${\omega }^i$ 替换为 ${\omega }^{-i}$，并且
• 把最终结果乘以一个因子 $1/n$。

换言之： $$\mathbf{A}=\frac{1}{n}{\mathbf{V}}_{{\omega }^{-1}}\cdot \hat{\mathbf{A}}.$$

（要理解逆 FFT 为何与 FFT 形式相似，请参阅 ² 的 Slide 13-1。下图同样取自 ²。）

Schwartz-Zippel 引理

Schwartz-Zippel 引理非形式地指出"不同的多项式在大多数点上都不同"。形式上，它可以写成如下：

设 $p(x_1,x_2,\cdots ,x_n)$ 是一个 $n$ 元、次数为 $d$ 的非零多项式。设 $S$ 是一个至少含有 $d$ 个元素的有限数集。如果我们从 $S$ 中随机选取 ${\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_n$，那么 $$\text{Pr}[p({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_n)=0]\le \frac{d}{|S|}.$$

在我们熟悉的单变量情形 $p(X)$ 下，这归结为：一个 $d$ 次非零多项式至多有 $d$ 个根。

Schwartz-Zippel 引理用于多项式相等性测试。给定两个多变量多项式 $p_1(x_1,\cdots ,x_n)$ 和 $p_2(x_1,\cdots ,x_n)$，次数分别为 $d_1,d_2$，我们可以对随机的 ${\alpha }_1,\cdots ,{\alpha }_n\leftarrow S$ 测试是否 $p_1({\alpha }_1,\cdots ,{\alpha }_n)-p_2({\alpha }_1,\cdots ,{\alpha }_n)=0$，其中 $S$ 的大小至少为 $|S|\ge (d_1+d_2).$ 如果两个多项式相同，这将永远成立；而如果两个多项式不同，那么等式成立的概率至多为 $\frac{\max (d_1,d_2)}{|S|}$。

消失多项式(Vanishing polynomial)

考虑以本原单位根 $\omega$ 生成的 $n$ 阶乘法子群 $\mathcal{H}$。对所有 ${\omega }^i\in \mathcal{H},i\in [n-1],$ 我们有 $({\omega }^i{)}^n=({\omega }^n{)}^i=({\omega }^0{)}^i=1.$ 换言之，$\mathcal{H}$ 的每个元素都满足方程

$$\begin{array}{rl}{Z}_H(X)& =X^n-1\\ & =(X-{\omega }^0)(X-{\omega }^1)(X-{\omega }^2)\cdots (X-{\omega }^{n-1}),\end{array}$$

也就是说每个元素都是 $Z_H(X)$ 的根。我们称 $Z_H(X)$ 为 $\mathcal{H}$ 上的 消失多项式(vanishing polynomial)，因为它在 $\mathcal{H}$ 的所有元素上求值都为零。

这在检查多项式约束时尤其方便。例如，要检查 $A(X)+B(X)=C(X)$ 在 $\mathcal{H}$ 上成立，我们只需检查 $A(X)+B(X)-C(X)$ 是 $Z_H(X)$ 的某个倍数。换言之，如果把我们的约束除以消失多项式仍然得到某个多项式 $\frac{A(X)+B(X)-C(X)}{Z_H(X)}=H(X),$ 我们就确信 $A(X)+B(X)-C(X)=0$ 在 $\mathcal{H}$ 上成立。

拉格朗日基函数(Lagrange basis functions)

TODO：（简要地）解释一般意义上的基(basis)是什么。

多项式通常用单项式基(monomial basis)书写（例如 $X,X^2,...X^n$）。然而，当在 $n$ 阶乘法子群上工作时，我们发现用拉格朗日基(Lagrange basis)表达更为自然。

考虑以本原单位根 $\omega$ 生成的 $n$ 阶乘法子群 $\mathcal{H}$。对应于这个子群的拉格朗日基是一组函数 $\{{\mathcal{L}}_i{\}}_{i=0}^{n-1}$，其中

$${\mathcal{L}}_i({\omega }^j)=\{\begin{array}{ll}1& \text{if }i=j,\\ 0& \text{otherwise.}\end{array}$$

我们可以更紧凑地写成 ${\mathcal{L}}_i({\omega }^j)={\delta }_{ij},$ 其中 $\delta$ 是克罗内克 delta 函数(Kronecker delta function)。

现在，我们可以把多项式写成拉格朗日基函数的线性组合，

$$A(X)=\sum _{i=0}^{n-1}{a}_i{\mathcal{L}}_i(X),X\in \mathcal{H},$$

这等价于说 $A(X)$ 在 ${\omega }^0$ 处求值为 $a_0$，在 ${\omega }^1$ 处为 $a_1$，在 ${\omega }^2$ 处为 $a_2,\cdots ,$ 依此类推。

当在乘法子群上工作时，拉格朗日基函数有一种便利的稀疏表示，形如

$${\mathcal{L}}_i(X)=\frac{c_i\cdot (X^n-1)}{X-{\omega }^i},$$

其中 $c_i$ 是重心权重(barycentric weight)。（要理解这个形式是如何推导出来的，请参阅 ³。）对于 $i=0,$ 我们有 $c=1/n⟹{\mathcal{L}}_0(X)=\frac{1}{n}\frac{(X^n-1)}{X-1}$。

假设我们给定了一组求值点 $\{x_0,x_1,\cdots ,x_{n-1}\}$。由于我们不能假定这些 $x_i$ 构成一个乘法子群，我们也考虑一般情形下的拉格朗日多项式 ${\mathcal{L}}_i$。那么我们可以构造：

$${\mathcal{L}}_i(X)=\prod _{j\ne i}\frac{X-x_j}{x_i-x_j},i\in [\mathrm{0..}n-1].$$

这里，每个 $X=x_j\ne x_i$ 都会产生一个为零的分子项 $(x_j-x_j),$ 使得整个乘积求值为零。另一方面，$X=x_i$ 会在每一项处求值为 $\frac{x_i-x_j}{x_i-x_j}$，导致整体乘积为一。这就在集合 $\{x_0,x_1,\cdots ,x_{n-1}\}$ 上给出了所需的克罗内克 delta 行为 ${\mathcal{L}}_i(x_j)={\delta }_{ij}$。

拉格朗日插值(Lagrange interpolation)

给定一个处于求值表示下的多项式

$$A:\{(x_0,A(x_0)),(x_1,A(x_1)),\cdots ,(x_{n-1},A(x_{n-1}))\},$$

我们可以在拉格朗日基下重构它的系数形式：

$$A(X)=\sum _{i=0}^{n-1}A(x_i){\mathcal{L}}_i(X),$$

其中 $X\in \{x_0,x_1,\cdots ,x_{n-1}\}.$

参考文献

1. Dasgupta, S., Papadimitriou, C. H., & Vazirani, U. V. (2008). "Algorithms" (ch. 2). New York: McGraw-Hill Higher Education. ↩

2. Golin, M. (2016). "The Fast Fourier Transform and Polynomial Multiplication" [lecture notes], COMP 3711H Design and Analysis of Algorithms, Hong Kong University of Science and Technology. ↩ ↩2

3. Berrut, J. and Trefethen, L. (2004). "Barycentric Lagrange Interpolation." ↩