_rqy's Blog

一只猫猫,想成为天才少女数学家!

0%

一种多项式多点求值算法

Problem

给定一个 m1m-1 次多项式 f(x)f(x) 以及 nn 个点 α0,,αn1\alpha_0,\dots,\alpha_{n-1},求 f(α0),f(α1),,f(αn1)f(\alpha_0),f(\alpha_1),\dots,f(\alpha_{n-1})

参考:《转置原理的简单介绍》rushcheyo, negiizhao, Created_Equal

参考:EI’s blog

参考:Tellegen’s Principle into Practice

Algorithm

大家可能都知道可以通过多项式取模来做多点求值(i.e.f(α)=f(x)mod(xα)f(\alpha)=f(x)\bmod(x-\alpha))。

我们考虑这个多点求值到底是个什么东西:

对于一个 m1m-1 次多项式 f(x)=i=0n1fixif(x)=\sum_{i=0}^{n-1}f_ix^i,我们将其看做一个 mm 维列向量 f=(f0,f1,,fm1)Tf=(f_0,f_1,\dots,f_{m-1})^T。(T-^T 表示转置)。那么我们求的是:对于每个 i[0,n)i\in[0,n)j=0n1fjαij\sum_{j=0}^{n-1}f_j\alpha_i^j。换句话说我们要求

(1α0α02α0m11α1α12α1m11α2α22α2m11αn1αn12αn1m1)f\begin{pmatrix} 1&\alpha_0&\alpha_0^2&\cdots&\alpha_0^{m-1}\\ 1&\alpha_1&\alpha_1^2&\cdots&\alpha_1^{m-1}\\ 1&\alpha_2&\alpha_2^2&\cdots&\alpha_2^{m-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&\alpha_{n-1}&\alpha_{n-1}^2&\cdots&\alpha_{n-1}^{m-1}\\ \end{pmatrix}f

固定 α0,αn1\alpha_0,\dots\alpha_{n-1},左边其实就是一个范德蒙矩阵(Vandermonde Matrix),记它为 V{\bf V}

所以我们要做的其实就是要求:将一个向量左乘一个范德蒙矩阵。

这个东西不太好搞;但是如果考虑左乘它的转置( uVTuu\mapsto {\bf V}^Tu),那么这个我们是知道如何去做的:把式子展开,我们要求的是:对每个 i[0,m)i\in[0,m),求出 j=0n1ujαji\sum_{j=0}^{n-1}u_j\alpha_j^i。换句话说,我们要求 j=0n1uj1αjx\sum_{j=0}^{n-1}\frac{u_j}{1-\alpha_jx} 的前 mm 项。这个东西我们很会做:这就是分治 FFT 上去。

那么我们已知了如何从 uu 去求 VTu{\bf V}^Tu。这对我们的目标(从 ff 求出 Vf{\bf V}f)有什么帮助?

Tellegen’s Principle

我们发现 uVTuu\mapsto {\bf V}^Tu 这个计算是线性的(假设 n,α0,,αn1n,\alpha_0,\dots,\alpha_{n-1} 给定,把 u0un1u_0\dots u_{n-1} 看作输入的话,可以发现整个算法中仅有加减以及乘除常数)。这也是很自然的,因为要计算的东西本来就正好是 uu 左乘一个矩阵。

既然这样,我们把 uVTuu\mapsto {\bf V}^Tu 的整个程序展开,就一定会得到这样的一个过程:每一步都是如下三种操作的一种:

  • x+=cy{\bf x}\mathrel{+}=c{\bf y}x,y\bf x, y 表示程序中的变量,cc 是某个和输入无关的常数,下同);
  • x=c{\bf x}\mathrel{*}=c
  • swap(x,y){\rm swap}({\bf x},{\bf y})(交换两个变量的值)。

注意这里我们一直认为 n,α0,,αn1n,\alpha_0,\dots,\alpha_{n-1}固定的,即不把它们看做输入。

这里其实我们就是将 VT{\bf V}^T 分解成了若干个初等行变换矩阵的乘积 VT=EkEk1E1{\bf V}^T=E_kE_{k-1}\dots E_1kk 表示把整个算法展开之后会有多少步)。因此我们就有 V=E1TE2TEkT{\bf V}=E_1^TE_2^T\dots E_k^T;换句话说,我们只需要把上面每一种操作也“转置”过来,并把整个过程“反过来”,那么我们就得到了一个 fVff\mapsto {\bf V}f 的算法;这里,第一种操作转置过来就会变成 y+=cx{\bf y}\mathrel{+}=c{\bf x};后两种操作转置过来是不变的。

这个操作(通过 uAuu\mapsto{\bf A}u 的算法来得到 fATff\mapsto{\bf A}^Tf 的算法)被称作“转置原理”,或“特勒根原理”(Tellegen’s Principle)。

举个例子:如果我们的输入是 (f0,f1,f2)(f_0,f_1,f_2),输出是 (f0,f0+f1,f0+f1+f2)(f_0,f_0+f_1,f_0+f_1+f_2),那么有这样一个算法:f1+=f0;  f2+=f1f_1\mathrel{+}=f_0;\;f_2\mathrel{+}=f_1。转置过来之后的输出是 (f0+f1+f2,f1+f2,f2)(f_0+f_1+f_2,f_1+f_2,f_2),对应的算法是 f1+=f2;  f0+=f1f_1\mathrel{+}=f_2;\;f_0\mathrel{+}=f_1

这样的话,我们至少可以“劝自己相信存在某个和原来一样复杂度一样常数的算法”,毕竟我们会做同样多次的加减法呢、同样多次的乘法。

实现

但是一个问题是,你显然不会蠢到在程序里面手动实现一个转置原理(比如,把你用到的加减乘除操作全都丢到一个队列里然后反过来做,虽然这是完全可行的)。当然你可以手动把程序反过来写,但是比较好的方法是手动分析这里用到的各种算法的”转置“是什么样的。

多项式乘法

首先肯定要考虑多项式乘法。固定一个 nn 次多项式 gg,然后把一个 mm 次多项式 ff 看做输入,那么我们应该得到一个 n+mn+m 次多项式 ansans。这个线性变换就是对每对 i,ji,j,把 gig_i 倍的 fjf_j 加到 ansi+jans_{i+j} 上。因此它的转置应该是:把一个 n+mn+m 次多项式 ff 作为输入,得到一个 nn 次多项式 ansans;然后对每对 i,ji,j,把 gig_i 倍的 fi+jf_{i+j} 加到 ansjans_j 上*(这个过程其实很简单,你只需要记住把“输入”到“输出”的系数反过来就好了)*。因此我们得到的式子就是 ansj=igifi+jans_j=\sum_ig_if_{i+j},其中 ansans 的次数不低于 f,gf,g 的次数之差(“不低于”是因为我们可能会把 ff 看做超过 n+mn+m 次多项式,因为我们正着 FFT 的时候会算 n+mn+m 次但只保留前若干项)。

这个过程可以通过把 gg 反过来 FFT 一次实现,我们把这样得到的结果 ansans 写作 g×nTfg\times^T_nf(上标 TT 表示转置,而下标标出了结果需要保留多少次)。

多点求值

考虑我们计算 VTf{\bf V}^Tf 的算法:要计算 i=0n1ui1αix\sum_{i=0}^{n-1}\frac{u_i}{1-\alpha_ix},记 g(x)=i=0n1(1αix)g(x)=\prod_{i=0}^{n-1}(1-\alpha_ix),然后我们计算 i=0n1fig(x)1αix\sum_{i=0}^{n-1}f_i\frac{g(x)}{1-\alpha_ix},最后再把它乘上 g1(x)g^{-1}(x)。首先我们先利用分治乘法来计算 g(x)g(x)。由于这部分并不依赖于输入,所以它不属于转置的范畴。(再强调一次我们只把 u0,,un1u_0,\dots,u_{n-1} 看作输入)

对于 i=0n1uip(x)1αix\sum_{i=0}^{n-1}u_i\frac{p(x)}{1-\alpha_ix},我们采用分治:

  • n=1n=1,则直接返回 f0f_0
  • n>1n>1,令 m=n/2m=\lfloor n/2\rfloor,递归计算左边(u0,,m1u_{0,\dots,m-1})的答案 ansLans_L 和右边(um,,n1u_{m,\dots,n-1})的答案 ansRans_R;然后返回 ansLgR+ansRgLans_Lg_R+ans_Rg_L;其中 gL,gRg_L,g_R 是两边对应的 gg,在最开始的分治里面已经计算过了。

若这一步得到的答案记为 hh,那么最终的答案(VTu{\bf V}^Tu)就是 g1×hg^{-1}\times h

将这个算法转置过来,不要忘记要把操作的顺序也反过来:

首先仍然是分治计算 gg,这一步是不变的;然后对于输入的 ff,我们计算 h=g1×nTfh=g^{-1}\times^T_nf(注意,考虑到前面我们的 hhnn 位,所以这里的 hh 也应该保留 nn 位。但是前面的 g1g^{-1} 要计算 mm 位,所以这里的 g1g^{-1} 也应该是计算完整的 mm 位)。

接下来把前面所说的分治转置过来:

  • n=1n=1,直接返回 h0h_0
  • n>1n>1,令 m=n/2m=\lfloor n/2\rfloor,然后记 hL=gR×mTh,hR=gL×nmThh_L=g_R\times^T_mh,h_R=g_L\times^T_{n-m}h,然后分别把 hL,hRh_L,h_R 向两边递归。

这样我们就可以得到 Vf{\bf V}f,即多点插值结果。

此外,在计算的时候不难发现 ×T\times^T 的操作所需要的 FFT 长度是和原来相同的(因为我们反过来卷积后只需要后面一半,所以 FFT 短一些,循环卷积溢出到前面那一半我们也用不到),所以这样的常数和前面的是相同的;我们只需要一次多项式求逆,别的都只是普通的fft。这相比于原来(每一步都要做一个多项式取模)快了不少。

Code

咕咕咕。