多项式多点求值|快速插值 学习笔记

Date: Nov. 13th 2021

Tag: Code, Algorithm

简介：咕了，下次补。

多点求值 Algorithm 1：

传统做法，分治+多项式取模。

先咕了，下次补。

多点求值 Algorithm 2：

新科技，基于转置原理。

前置知识

转置原理

简述：

• 现有 $n\times n$ 的矩阵 $V$，和 $n$ 维向量 $f$。若可求 $f\to Vf$ ，则可以在相同复杂度内计算 $f\to V^Tf$。

不难发现，若一个问题可以转化为 $f\to Vf$ $\Leftrightarrow$ 求解过程只需对输入 $f$ 进行线性变换(将所有线性变换合并，表示成矩阵 $V$，最后右乘 $f$ 即可)。

转置原理的更多内容将在之后呈现，让我们先来学习一下相关概念和知识。

线性变换

对于一个向量 $f$ 可进行的简单线性变换有如下几种：()

1. $\operatorname{add}$ : $f_i = f_i + c * f_j$ ($c$ 必须是常数，下同)
2. $\operatorname{mul}$ : $f_i = f_i * c$
3. $\operatorname{swap}$ : $f_i, f_j = f_j, f_i$ (交换 $f_i, f_j$)

简单线性操作可以表示为 $f$ 左乘一个初等变换 $E$，其各其转置后的操作分别为：

1. $\operatorname{add}^T$ : $f_j = f_j + c * f_i$
2. $\operatorname{mul}^T$ : $f_i = f_i * c$
3. $\operatorname{swap}^T$ : $f_i, f_j = f_j, f_i$ (和 $\operatorname{swap}$ 没有区别)

不难发现 $E$ 对应的简单线性操作转置后对应初等变换 $E^T$。

线性操作转置后，相当于在原操作的基础上交换了输入和输出下标。

我们可以认为操作过程中存在一个列向量 $f^{memory}$ 作为内存空间，除原始输入外，允许使用一些额外位存储中间变量。这样理解可以更好的把抽象运算和我们的程序联系在一起。

另外，我们规定 $f_{1…n}$ 为程序开始时的输入位(有初始值)，也是程序结束时的输出位。再规定如下两个操作：

4. $\operatorname{input-mv}$ : 程序开始时执行，将 $f_{1…n}$ 通过 $\operatorname{swap}$ 移动至 $f_{i_1…i_n}$ 。
5. $\operatorname{output-mv}$ : 程序结束时执行，将 $f_{i_1…i_n}$ 通过 $\operatorname{swap}$ 移动至 $f_{1…n}$ 。

操作 $\operatorname{input-mv}$ 和 $\operatorname{output-mv}$ 本质都是多次 $\operatorname{swap}$, 且互为转置。

引入 $\operatorname{input/output-mv}$ 操作，是为了方便解释应用转置原理的算法中互逆的两个过程交换输入输出位置的原因。(详见后文的 “关于 process1 和 process2 中输入输出位置交换的解释。”)

转置原理的思想

考虑一个 ($n \times n$) 的操作矩阵 $A$，和一个 ($1 \times n$) 的输入向量 $X$，输出一个 ($1 \times n$) 的向量 $Y$，其中：

$$
AX =
\left ( \begin{matrix}
a_{1,1} & a_{1,2} & \cdots & a_{1,n}\\
a_{2,1} & a_{2,2} & \cdots & a_{2,n}\\
\vdots & \vdots & \ddots &\vdots \\
a_{n,1} & a_{n,2} & \cdots & a_{n, n}\\
\end{matrix} \right )
\times
\left (
\begin{matrix}
x_1\\
x_2\\
\vdots\\
x_n
\end{matrix}
\right ) =
\left (
\begin{matrix}
y_1\\
y_2\\
\vdots\\
y_n
\end{matrix}
\right ) = Y
$$

若能求出 $AX=Y$，现要求出 $A^TX=Y’$，那么根据**特勒根原理(Tellegen’s Principle)**：

1. 将 $A$ 拆分成 $m$ 个初等变换 $E_{1…m}$，即 $A=\prod\limits_{i=1}^m E_i$ 。
2. 则 $A^T$ 也可以用这 $m$ 个初等变换表示，为 $A^T=\prod\limits_{i=1}^m E^T_{m-i+1}$ 。
3. 由此，可以得到 $Y’=A^TX = (\prod\limits_{i=1}^m E^T_{m-i+1})X$ 。

可以认为 $AX$ 是将 $E_{m…1}$ 以$\text{顺序}{(m\to1)}$作用在 $X$ 上得到 $Y$，那么 $A^TX$ 就是将 $E{m…1}$ 转置后$text{逆序}_{(1\to m)}$作用在 $X$ 上，得到 $Y’$。

(注 ：$AX$ 中 $X$ 是左乘 $A$，所以根据结合律，$E_{m…1}$ 对 $X$ 的作用顺序是 $m\to1$。 )

值得注意的是，求解 $AX$ 和 $A^TX$ 时，过程均为求 $m$ 个大小相同的初等变换 $E$ 和 $X$ 之积，复杂度是相同的。

因为显然有 $(A^T)^T = A$， 所以若能求出 $A^TX$ 那么也可以用此构造法求出 $AX$ 的结果。

关于卷积和线性变换

既然要将转置原理应用到多项式问题上，自然避免不了卷积，让我们分析一下卷积在线性变换下的表示。

现在有 $h = f \times g$ ，其中 $f$ 是输入，$g$ 是定值。(若 $f, g$ 都是不确定的，则卷积无法表示成线性变换。)

记 $n = \deg f, m = \deg g$，则：

$$
\begin{align}
h = f \times g \Leftrightarrow &
\text{ for }i\text{ in }[0, n),\ j \text{ in }[0, m): \text{do}\ h_{i+j}=h_{i+j} + g_j*f_i
\end{align}
$$

也就是进行了 $nm$ 次的 $\operatorname{add}$ 操作。

接下来考虑卷积的转置，记为 $\times^T$，还是使用上方 $h = f \times g$ 的例子($\deg h = n + m -1$)，可得：
$$
\begin{align}
f = h\times^Tg \Leftrightarrow &
\text{ for }i\text{ in }[0\to n),\ j \text{ in }[0\to m): \text{do}\ f_i=f_i+g_j*h_{i+j}\\
\Leftrightarrow & \text{ for } i \text{ in } [0\to n+m-1),\ j \text{ in } [0 \to \min(i,m-1)]: \\
& \text{ do } f_{i-j} = f_{i-j} + g_j * h_i
\end{align}
$$
也就是逆卷积。**可用 $g$ 反转后和 $h$ 卷积，并截取的第 $m$ 至 $m + n - 1$ 项来计算 $f$**。

更形式化的，记 $g$ 反转的多项式为 $g_r$ ，$\deg h = n + m - 1$，$\deg g = m$，

则 $f = h \times^T g$ 可表示为：
$$
f = h\times^Tg = \frac{h\times g_r - h \times g_r \bmod x^m}{x^m} \bmod x^n
$$

至此，我们已经学会了所以必要的前置知识。

算法思想

现在回归到多点求值的问题上，考虑将 $n$ 阶多项式 $f$ 看作是 $n$ 维列向量：
$$
f= \left ( \begin{matrix} [x^0]f\\ [x^1]f\\ [x^2]f\\ \vdots \\ [x^{n-1}]f\end{matrix} \right )
$$
为了方便起见，记 $f_i := [x^i]f$，表示 $f$ 中 $x^i$ 项的系数，并规定矩阵下标从 $0$ 开始。

记有 $n$ 个点坐标 $a_i$ 需要求出点值，构造其范德蒙德矩阵 $V$：
$$
V=\left ( \begin{matrix}
1 & a_0 & a_0^2 & \cdots & a_0^{n-1} \\
1 & a_1 & a_1^2 & \cdots & a_1^{n-1} \\
1 & a_2 & a_2^2 & \cdots & a_2^{n-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & a_{n-1} & a_{n-1}^2 & \cdots & a_{n-1}^{n-1}
\end{matrix} \right )
$$
则 $f$ 对 $a_i$ 求点值的结果为 $(V\times f)[i]$ ，**故 $f$ 对 $a_{0…n-1}$ 多点求值等价于求解 $Vf$**。

显然不能用暴力矩阵乘法求 $Vf$，复杂度为 $n^3$ （本质上和多次单点带入求值没有区别）。

考虑利用转置原理，先考虑计算 $V^Tf$ 的过程(不必实际求解)，再构造求解 $Vf$。

计算 $V^Tf$ 的方法

观察 $V^T$ 的结构：
$$
V^T=\left ( \begin{matrix}
1 & 1 & 1 & \cdots & 1\\
a_1 & a_2 & a_3 & \cdots & a_n \\
a_1^2 & a_2^2 & a_3^2 & \cdots & a_n^{2} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
a_1^{n-1} & a_2^{n-1} & a_3^{n-1} & \cdots & a_n^{n-1}
\end{matrix} \right )
$$

发现 $V^Tf$ 不难求，即对于 $\forall d \in[0, n)$, 求出 $(V^Tf)[d]=\sum\limits_{i=0}^{n-1} f_ia_i^d$。

运用一点点生成函数知识可得：$(V^Tf)[d] = [x^d]\sum\limits_{i=0}^{n-1} \frac{f_i}{1-a_ix}$，所以我们只要求出 $G(x) = \sum\limits_{i=0}^{n-1} \frac{f_i}{1-a_ix}$的前 $n$ 项即可。

考虑分治求解 $V^Tf$ （即求解 G）的过程，记为 process1：

1. 对于 $n=1$ 的情况，返回 $u=f_i$，$v = 1-a_ix$ ($u$ 为分母， $v$ 为分子)。
2. 对于 $n > 1$ 的情况分治处理，将 $\le\lfloor \frac n2\rfloor$ 部分的结果记为 $u_0, v_0$，将 $\gt \lfloor \frac n2\rfloor$ 部分的结果记为 $u_1, v_1$。
   • 因为 $\frac{u_0}{v_0}+\frac{u_1}{v_1} = \frac{u_0v_1 + u_1v_0}{v_0v_1}$，所以返回 $u=u_0v_1+u_1v_0,v=v_0v_1$ 。
3. 最后的结果是 $G = u_{final}\times v_{final}^{-1}$。

观察可以发现， process1 中的分母 $v$ 与输入多项式 $f$ 无关(可视为定值)，分子 $u$ 则与 $f$ 有关。

故分治过程中的 $u=u_0v_1+u_1v_0$ 可以看作两次分步操作，即 $u = u + v_1 \times u_0, u = u + v_0 \times u_1$。

也就是说，$u$ 接受了两次子分治结果的上传操作，这样的理解对接下来运用转置原理会有帮助。

构造求解 $Vf$ 的方法

首先有个很简单并且正确(并且相当暴力)的做法：

• 执行 process1, 记下过程中所有的加法乘法等操作，最后转置逆序执行所有操作,得到结果。

当然这么做很不优美，实现起来也很[难以描述]。尝试简化这个过程。

相信你已经理解了卷积的转置，这在接下来的算法中非常重要，如果没有理解请回看一下。

我们即将进入算法实现的核心部分。

考虑构造求解 $Vf$ 的过程，该过程可视为 process1 的反向，记为 process2：

1. 先用分治 FFT 单独预处理出所有 $v$，并记录保留。
   • 在 $n=1$ 处初始化 $v=1-a_ix$，合并时 $v = v_0v_1$。
2. 初始化 $G = f$， $u_{final} = G \times^T v_{final}^{-1}$，以此为起点向下“反向分治”递归。
3. 对于 $n>1$ 的情况，当前点的 $u, v$ 已知，将 $\le\lfloor \frac n2\rfloor$ 部分初始化值记为 $u_0, v_0$，将 $\gt \lfloor \frac n2\rfloor$ 部记初始化值为 $u_1, v_1$。 其中 $v_0, v_1$ 已在预处理过程中求出。
   • 令 $u_0$ = $u \times^T v_1$, $u_1 = u \times^T v_0$ ，分别往左右两侧继续递归。
     • 这就 process1 中的上传操作在 process2 中的逆操作，下传操作。
     • 把一次上传操作 $u=u+u_0 \times v_1$ 分步表示为 $h = u_0 \times v_1, u = u + h$ ，不难发现其逆操作是 $h = h + u, u_0 = u \times^T v_1$。其中 $h$ 是用于存储卷积结果的临时空间，逆操作中的 $h = h + u$ 可以忽略。
     • 所以，$(u=u+u_0\times v_1)^T$ 可化简为 $u_0 = u\times^Tv_1$。
     • 同理，$(u=u+u_1\times v_0)^T$ 可化简为 $u_1 = u\times^Tv_0$。
4. 对于 $n=1$ 的情况，记其对应下标为 $i$。$[x^0]u$ 就是 $a_i$ 的点值，记录 $(Vf)[i] = [x^0]u$ ，返回。
5. 递归处理完所有 $n=1$ 的情况后，$Vf$ 已求得（即叶子节点处 $[x^0]u$ 拼接而成的向量）。

至此，我们已经找到了求解 $Vf$ 的方法，代码只需要实现 process2 即可，process1 是我们用来构造 process2 的思维过程。

时间复杂度为 $T(n) = 2T(\frac n2) + O(n\log n) = O(n\log^2 n)$。

理论复杂度和传统做法一样，但只需要用到常规 FFT 和一次多项式求逆，不需要多项数取模，常数要小很多。经测试，新做法的实际运行效率大概是老做法的两倍，非常恐怖。

一些 trick 和细节补充：

• 若要求点值的点数 $m$ 和多项式的长度 $n$ 不同，需各自补齐至 $\max(n, m)$，不然 $V^Tf$ 显然不满足矩乘的要求，没有意义。对于补齐点，随便补什么都可以，毕竟多求出几个其它点值也没什么坏处。对于补齐多项式，高位补 $0$ 就好了，常规操作。

  同时，建议像 DFT 那样补齐至 $2^k$。这样分治/递归时每段的长度都是 $2$ 的整数次，且左右儿子的长度相同，写起来方便，不用思考太多长度细节问题。

• 由于 $v$ 之后的大部分使用都需要反转，不如一开始就在叶子节点处初始化成 $a_i -x$ 而不是原本的 $1-a_ix$。也就是预处理出 $v_r$ 代替 $v$ 。(当然如果这样初始化，要用到 $v$ 不反转的状态时，记的反转。)

  可能有一定的常数级别优化，但是影响不大。

• 程序需要用到 $\times^T$ 的地方非常多，可以实现一个函数处理 $h\times^Tg$。而 $\times^T$ 中要用到 $g_r$，所以需要保证 $g$ 的长度正确(不同长度 $g$ 反转后的 $g_r$ 不同 ，会影响 $\times^T$ 的结果 )。长度差异问题主要来自卷积 DFT 时将 $\deg$ 补齐至 $2^k$，将卷积结果调整为实际长度即可(e.g. $h = f \times g$，那么 $\deg h =\deg f + \deg g - 1$)。

  该算法中所有的 $\times^T$ 运算中，$g$ 都来自预处理的分母 $v$。故处理 $v = v_0v_1$ 时，注意将 $\deg v$ 调整为 $v$ 实际长度(通常是 $\text{分治区间长度}+1$ ，也可以用 $\deg v_0 + \deg v_1 - 1$ 表示)。

代码实现：

先咕了，下次补。暂时贴个多项式板子在这里。

多点求值 Algorighm：

先咕了，下次补。

参考链接：

• x义x - 多点求值的新做法
• _rqy - 一种多项式多点求值算法
• rushcheyo - 转置原理及其应用