Problem
给定一个 m−1 次多项式 f(x) 以及 n 个点 α0,…,αn−1,求 f(α0),f(α1),…,f(αn−1)。
参考:《转置原理的简单介绍》rushcheyo, negiizhao, Created_Equal
参考:EI’s blog
参考:Tellegen’s Principle into Practice
Algorithm
大家可能都知道可以通过多项式取模来做多点求值(i.e.f(α)=f(x)mod(x−α))。
我们考虑这个多点求值到底是个什么东西:
对于一个 m−1 次多项式 f(x)=∑i=0n−1fixi,我们将其看做一个 m 维列向量 f=(f0,f1,…,fm−1)T。(−T 表示转置)。那么我们求的是:对于每个 i∈[0,n) 求 ∑j=0n−1fjαij。换句话说我们要求
⎝⎜⎜⎜⎜⎜⎜⎛111⋮1α0α1α2⋮αn−1α02α12α22⋮αn−12⋯⋯⋯⋱⋯α0m−1α1m−1α2m−1⋮αn−1m−1⎠⎟⎟⎟⎟⎟⎟⎞f
固定 α0,…αn−1,左边其实就是一个范德蒙矩阵(Vandermonde Matrix),记它为 V。
所以我们要做的其实就是要求:将一个向量左乘一个范德蒙矩阵。
这个东西不太好搞;但是如果考虑左乘它的转置( u↦VTu),那么这个我们是知道如何去做的:把式子展开,我们要求的是:对每个 i∈[0,m),求出 ∑j=0n−1ujαji。换句话说,我们要求 ∑j=0n−11−αjxuj 的前 m 项。这个东西我们很会做:这就是分治 FFT 上去。
那么我们已知了如何从 u 去求 VTu。这对我们的目标(从 f 求出 Vf)有什么帮助?
Tellegen’s Principle
我们发现 u↦VTu 这个计算是线性的(假设 n,α0,…,αn−1 给定,把 u0…un−1 看作输入的话,可以发现整个算法中仅有加减以及乘除常数)。这也是很自然的,因为要计算的东西本来就正好是 u 左乘一个矩阵。
既然这样,我们把 u↦VTu 的整个程序展开,就一定会得到这样的一个过程:每一步都是如下三种操作的一种:
-
x+=cy(x,y 表示程序中的变量,c 是某个和输入无关的常数,下同);
-
x∗=c;
-
swap(x,y)(交换两个变量的值)。
注意这里我们一直认为 n,α0,…,αn−1 是固定的,即不把它们看做输入。
这里其实我们就是将 VT 分解成了若干个初等行变换矩阵的乘积 VT=EkEk−1…E1(k 表示把整个算法展开之后会有多少步)。因此我们就有 V=E1TE2T…EkT;换句话说,我们只需要把上面每一种操作也“转置”过来,并把整个过程“反过来”,那么我们就得到了一个 f↦Vf 的算法;这里,第一种操作转置过来就会变成 y+=cx;后两种操作转置过来是不变的。
这个操作(通过 u↦Au 的算法来得到 f↦ATf 的算法)被称作“转置原理”,或“特勒根原理”(Tellegen’s Principle)。
举个例子:如果我们的输入是 (f0,f1,f2),输出是 (f0,f0+f1,f0+f1+f2),那么有这样一个算法:f1+=f0;f2+=f1。转置过来之后的输出是 (f0+f1+f2,f1+f2,f2),对应的算法是 f1+=f2;f0+=f1。
这样的话,我们至少可以“劝自己相信存在某个和原来一样复杂度一样常数的算法”,毕竟我们会做同样多次的加减法呢、同样多次的乘法。
实现
但是一个问题是,你显然不会蠢到在程序里面手动实现一个转置原理(比如,把你用到的加减乘除操作全都丢到一个队列里然后反过来做,虽然这是完全可行的)。当然你可以手动把程序反过来写,但是比较好的方法是手动分析这里用到的各种算法的”转置“是什么样的。
多项式乘法
首先肯定要考虑多项式乘法。固定一个 n 次多项式 g,然后把一个 m 次多项式 f 看做输入,那么我们应该得到一个 n+m 次多项式 ans。这个线性变换就是对每对 i,j,把 gi 倍的 fj 加到 ansi+j 上。因此它的转置应该是:把一个 n+m 次多项式 f 作为输入,得到一个 n 次多项式 ans;然后对每对 i,j,把 gi 倍的 fi+j 加到 ansj 上*(这个过程其实很简单,你只需要记住把“输入”到“输出”的系数反过来就好了)*。因此我们得到的式子就是 ansj=∑igifi+j,其中 ans 的次数不低于 f,g 的次数之差(“不低于”是因为我们可能会把 f 看做超过 n+m 次多项式,因为我们正着 FFT 的时候会算 n+m 次但只保留前若干项)。
这个过程可以通过把 g 反过来 FFT 一次实现,我们把这样得到的结果 ans 写作 g×nTf(上标 T 表示转置,而下标标出了结果需要保留多少次)。
多点求值
考虑我们计算 VTf 的算法:要计算 ∑i=0n−11−αixui,记 g(x)=∏i=0n−1(1−αix),然后我们计算 ∑i=0n−1fi1−αixg(x),最后再把它乘上 g−1(x)。首先我们先利用分治乘法来计算 g(x)。由于这部分并不依赖于输入,所以它不属于转置的范畴。(再强调一次我们只把 u0,…,un−1 看作输入)
对于 ∑i=0n−1ui1−αixp(x),我们采用分治:
- 若 n=1,则直接返回 f0。
- 若 n>1,令 m=⌊n/2⌋,递归计算左边(u0,…,m−1)的答案 ansL 和右边(um,…,n−1)的答案 ansR;然后返回 ansLgR+ansRgL;其中 gL,gR 是两边对应的 g,在最开始的分治里面已经计算过了。
若这一步得到的答案记为 h,那么最终的答案(VTu)就是 g−1×h。
将这个算法转置过来,不要忘记要把操作的顺序也反过来:
首先仍然是分治计算 g,这一步是不变的;然后对于输入的 f,我们计算 h=g−1×nTf(注意,考虑到前面我们的 h 有 n 位,所以这里的 h 也应该保留 n 位。但是前面的 g−1 要计算 m 位,所以这里的 g−1 也应该是计算完整的 m 位)。
接下来把前面所说的分治转置过来:
- 若 n=1,直接返回 h0;
- 若 n>1,令 m=⌊n/2⌋,然后记 hL=gR×mTh,hR=gL×n−mTh,然后分别把 hL,hR 向两边递归。
这样我们就可以得到 Vf,即多点插值结果。
此外,在计算的时候不难发现 ×T 的操作所需要的 FFT 长度是和原来相同的(因为我们反过来卷积后只需要后面一半,所以 FFT 短一些,循环卷积溢出到前面那一半我们也用不到),所以这样的常数和前面的是相同的;我们只需要一次多项式求逆,别的都只是普通的fft。这相比于原来(每一步都要做一个多项式取模)快了不少。
Code
咕咕咕。