分治FFT与多项式点值插值

一.多项式乘法.

作为最基础的多项式运算，这里不再介绍，具体参见快速傅里叶变换FFT与快速数论变换NTT入门.

这里的多项式乘法卡常技巧不会过于深入，且只适用于FFT.


二.分治FFT.

分治FFT其实并不是一个算法，它只是分治里面用FFT处理信息，不过这个东西有一些经典的模型.

比如给定$g_i$，要求：
$$f_i=\sum _{j=1}^i{f}_{i-j}{g}_j$$

初始$f_0=1$.

这种类卷积的式子一般都是用cdq分治+FFT来做的（当然还可以多项式求逆）.

具体来说，每次我们分治区间$[l,r]$中点为$m$时，先递归计算$[l,m]$，然后计算$[l,m]$对$[m+1,r]$的贡献，再递归计算$[m+1,r]$.

时间复杂度$O(n{\log }^{2}n)$.

顺带一提，分治到区间长度足够小的时候，可以选择变成暴力求解，很多情况下能够显著减小常数.

代码如下：

int ta[N+9],tb[N+9],ans[N+9];

void Divide_ans(int l,int r){
  if (l==r) {ans[l]+=!l;return;}
  int mid=l+r>>1;
  Divide_ans(l,mid);
  for (int i=l;i<=mid;++i) ta[i-l]=ans[i];
  for (int i=l;i<=r;++i) tb[i-l]=a[i-l];
  Get_len(r-l);
  NTT(ta,len,0);
  NTT(tb,len,0);
  for (int i=0;i<len;++i) smul(ta[i],tb[i]),tb[i]=0;
  NTT(ta,len,1);
  for (int i=mid+1;i<=r;++i) sadd(ans[i],ta[i-l]);
  for (int i=0;i<len;++i) ta[i]=0;
  Divide_ans(mid+1,r);
}

三.多项式求逆.

多项式的逆：对于一个多项式$F(x)$，其逆$G(x)$定义为：
$$G(x)F(x)=1$$

但很多情况下$G(x)$是一个无穷级数的形式，所以一般我们只需要取其最后$n$项（次数低于$n$次），记为：
$$G(x)F(x)\equiv 1(mod{x}^n)$$

当然可以暴力求解，但是有一种倍增的思路：我们假设已经求出了模$x^{\lceil \frac{n}{2}\rceil }$意义下$F(x)$的逆$H(x)$，推出$G(x)$与$H(x)$之间的关系.

现在有两个条件：
$$\begin{gathered} G(x)F(x)\equiv 1(mod{x}^{\lceil \frac{n}{2}\rceil }) \\ H(x)F(x)\equiv 1(mod{x}^{\lceil \frac{n}{2}\rceil }) \end{gathered}$$

将两式相减可得：
$$\begin{gathered} (G(x)-H(x))F(x)\equiv 0(mod{x}^{\lceil \frac{n}{2}\rceil }) \\ G(x)-H(x)\equiv 0(mod{x}^{\lceil \frac{n}{2}\rceil }) \end{gathered}$$

再平方一下可以得到：
$$\begin{gathered} (G(x)-H(x){)}^2\equiv 0(mod{x}^n) \\ {G}^2(x)+H^2(x)-2G(x)H(x)\equiv 0(mod{x}^n) \end{gathered}$$

两边同时乘上$F(x)$得到：
$$\begin{gathered} F(x)G^2(x)+F(x)H^2(x)-2F(x)G(x)H(x)\equiv 0(mod{x}^n) \\ G(x)+F(x)H^2(x)-2H(x)\equiv 0(mod{x}^n) \\ G(x)\equiv 2H(x)-F(x)H^2(x)(mod{x}^n) \\ G(x)\equiv H(x)(2-F(x)H(x))(mod{x}^n) \end{gathered}$$

然后就可以快乐倍增了.

时间复杂度$O(n\log n)$.

代码如下：

int ta[N+9];

void Poly_inv(int *a,int n,int *res){
  if (!n) {res[0]=Get_inv(a[0]);return;}
  Poly_inv(a,n>>1,res);
  for (int i=0;i<=n;++i) ta[i]=a[i];
  Get_len(n<<1);
  NTT(res,len,0);
  NTT(ta,len,0);
  for (int i=0;i<len;++i) smul(res[i],sub(2,mul(res[i],ta[i]))),ta[i]=0;
  NTT(res,len,1);
  for (int i=n+1;i<=n<<1;++i) res[i]=0;
}

四.多项式除法.

多项式除法：一个多项式$n$次多项式$F(x)$和$m$次多项式$G(x)$的商$P(x)$和余数$Q(x)$定义为：
$$F(x)=G(x)P(x)+Q(x)$$

其中$P(x)$是一个$n-m$次多项式，$Q(x)$是一个$m-1$次多项式.

我们定义一个$n$次多项式$F(x)$的翻转$F^R(x)$为：
$$F^R(x)=x^{n}F\left(\frac{1}{x}\right)$$

那么有：
$$\begin{gathered} x^{n}F\left(\frac{1}{x}\right)=x^{n}G\left(\frac{1}{x}\right)P\left(\frac{1}{x}\right)+x^{n}Q\left(\frac{1}{x}\right) \\ {F}^R(x)=x^{m}G\left(\frac{1}{x}\right)x^{n-m}P\left(\frac{1}{x}\right)+x^{n-m+1}{x}^{m-1}Q\left(\frac{1}{x}\right) \\ {F}^R(x)=G^R(x)P^R(x)+x^{n-m+1}{Q}^R(x) \\ {F}^R(x)\equiv G^R(x)P^R(x)(mod{x}^{n-m+1}) \\ {P}^R(x)\equiv \frac{F^R(x)}{G^R(x)}(mod{x}^{n-m+1}) \end{gathered}$$

求出了$P(x)$后就很容易求出$Q(x)$了.

一般情况下$n<m$，故时间复杂度为$O(n\log n)$.

由于接下来我们只会用到多项式取模，所以只给出多项式取模的代码.

代码如下：

int tb[N+9];

void Poly_mod(int *a,int *b,int n,int m,int *res){
  if (n<m){
	for (int i=0;i<=n;++i) res[i]=a[i];
	return;
  }
  reverse(b,b+m+1);
  Poly_inv(b,n-m,tb);
  reverse(b,b+m+1);
  for (int i=0;i<=n-m;++i) ta[i]=a[n-i];
  Get_len(n-m<<1);
  NTT(ta,len,0);
  NTT(tb,len,0);
  for (int i=0;i<len;++i) smul(ta[i],tb[i]),tb[i]=0;
  NTT(ta,len,1);
  for (int i=n-m+1;i<=n-m<<1;++i) ta[i]=0;
  reverse(ta,ta+n-m+1);
  for (int i=0;i<=m;++i) tb[i]=b[i];
  Get_len(n);
  NTT(ta,len,0);
  NTT(tb,len,0);
  for (int i=0;i<len;++i) smul(ta[i],tb[i]),tb[i]=0;
  NTT(ta,len,1);
  for (int i=0;i<m;++i) res[i]=sub(a[i],ta[i]);
  for (int i=0;i<=n;++i) ta[i]=0;
}

五.多项式多点求值.

给定一个多项式$F(x)$和$n$个横坐标$x_i$，要求所有$x_i$对应的$y_i=F(x_i)$.

用秦九韶算法暴力代入最多只能做到$O(n^2)$，考虑将所有$y_i$一起求.

定理8.1：对于一个$n$次多项式$F(x)$和一个常数$a$，有：
$$F(a)=F(x)mod(x-a)$$

定理8.2：对于三个多项式$F(x),A(x),B(x)$，有：
$$F(x)mod(A(x)B(x))modA(x)=F(x)modA(x)$$

然后就可以快乐分治了.

具体来说，我们先用类似于构建线段树的方式预处理出每一个在分治求值时会用到的区间$[l,r]$对应的$(x-a_i)$之积，然后分治下去.

时间复杂度$O(n{\log }^{2}n)$，空间复杂度$O(n\log n)$.

代码如下：

int *prod[N+9];

void Divide_prod(int *x,int l,int r,int k){
  prod[k]=new int[r-l+2];
  if (l==r) {prod[k][0]=sub(0,x[l]);prod[k][1]=1;return;}
  int mid=l+r>>1,ls=k<<1,rs=k<<1|1;
  Divide_prod(x,l,mid,ls);
  Divide_prod(x,mid+1,r,rs);
  for (int i=0;i<=mid-l+1;++i) ta[i]=prod[ls][i];
  for (int i=0;i<=r-mid;++i) tb[i]=prod[rs][i];
  Get_len(r-l+1);
  NTT(ta,len,0);
  NTT(tb,len,0);
  for (int i=0;i<len;++i) smul(ta[i],tb[i]),tb[i]=0;
  NTT(ta,len,1);
  for (int i=0;i<=r-l+1;++i) prod[k][i]=ta[i],ta[i]=0;
}

int *polynow[N+9];

void Divide_polynow(int l,int r,int k,int *res){
  if (l==r) {res[l]=polynow[k][0];return;}
  int mid=l+r>>1,ls=k<<1,rs=k<<1|1;
  Poly_mod(polynow[k],prod[ls],r-l,mid-l+1,polynow[ls]=new int[mid-l+1]);
  Divide_polynow(l,mid,ls,res);
  Poly_mod(polynow[k],prod[rs],r-l,r-mid,polynow[rs]=new int[r-mid]);
  Divide_polynow(mid+1,r,rs,res);
}

void Poly_val(int *a,int *x,int n,int m,int *res){
  Divide_prod(x,0,m,1);
  Poly_mod(a,prod[1],n,m+1,polynow[1]=new int[m]);
  Divide_polynow(0,m,1,res);
}

六.多项式快速插值.

给定$n+1$个点$(x_i,y_i)$，求一个$n$次多项式$F(x)$使得$F(x)$经过所有$n+1$个点.

一个$O(n^2)$的做法是，利用拉格朗日插值：
$$F(x)=\sum _{i=0}^n{y}_i\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$$

我们现在要在与多点求值同样的复杂度内完成这个过程.

设$\pi (x)={\prod }_{i=0}^n(x-x_i)$.

那么有：
$$F(x)=\sum _{i=0}^n{y}_i\frac{\pi (x)}{x-x_i}\frac{x_i-x_i}{\pi (x_i)}$$

这个式子中有一个$\frac{0}{0}$的形式，并不严谨，于是将这部分换成极限，并使用洛必达法则：
$$\underset{t\to x_i}{\lim }\frac{t-x_i}{\pi (t)}=\frac{1}{{\pi }^{\prime }(x_i)}$$

那么原式相当于：
$$F(x)=\sum _{i=0}^n\frac{y_i}{{\pi }^{\prime }(x_i)}\frac{\pi (x)}{x-x_i}$$

可以用多点求值求得$a_i=\frac{y_i}{{\pi }^{\prime }(x_i)}$，那么剩下的式子变成：
$$F(x)=\sum _{i=0}^n{a}_i\frac{\pi (x)}{x-x_i}$$

考虑再一次使用分治FFT，设${\pi }_{l,r}(x)$和$F_{l,r}(x)$：
$$\begin{gathered} {\pi }_{l,r}(x)=\prod _{i=l}^r(x-x_i) \\ {F}_{l,r}(x)=\sum _{i=l}^r{a}_i\frac{{\pi }_{l,r}(x)}{x-x_i} \end{gathered}$$

合并${\pi }_{l,r}(x)$可以直接多项式乘法，考虑合并$F_{l,r}(x)$.

设区间$[l,r]$中点为$m$，则：
$$\begin{gathered} F_{l,r}(x)=\sum _{i=l}^r{a}_i\frac{{\pi }_{l,r}(x)}{x-x_i} \\ =\sum _{i=l}^m{a}_i\frac{{\pi }_{l,r}(x)}{x-x_i}+\sum _{i=m+1}^r{a}_i\frac{{\pi }_{l,r}(x)}{x-x_i} \\ =\sum _{i=l}^m{a}_i\frac{{\pi }_{l,m}(x){\pi }_{m+1,r}(x)}{x-x_i}+\sum _{i=m+1}^r{a}_i\frac{{\pi }_{l,m}(x){\pi }_{m+1,r}(x)}{x-x_i} \\ ={\pi }_{m+1,r}(x)\sum _{i=l}^m{a}_i\frac{{\pi }_{l,m}(x)}{x-x_i}+{\pi }_{l,m}(x)\sum _{i=m+1}^r{a}_i\frac{{\pi }_{m+1,r}(x)}{x-x_i} \\ ={\pi }_{m+1,r}(x)F_{l,m}(x)+{\pi }_{l,m}(x)F_{m+1,r}(x) \end{gathered}$$

时间复杂度$O(n{\log }^{2}n)$，空间复杂度$O(n\log n)$.

代码如下：

int tc[N+9],te[N+9];

void Divide_inter(int *res,int l,int r,int k){
  if (l==r) return;
  int mid=l+r>>1,ls=k<<1,rs=k<<1|1;
  Divide_inter(res,l,mid,ls);
  Divide_inter(res,mid+1,r,rs);
  for (int i=0;i<=mid-l;++i) ta[i]=res[i+l],tc[i]=prod[ls][i];
  tc[mid-l+1]=prod[ls][mid-l+1];
  for (int i=0;i<r-mid;++i) tb[i]=res[i+mid+1],te[i]=prod[rs][i];
  te[r-mid]=prod[rs][r-mid];
  Get_len(r-l);
  NTT(ta,len,0);
  NTT(tb,len,0);
  NTT(tc,len,0);
  NTT(te,len,0);
  for (int i=0;i<len;++i) ta[i]=add(mul(ta[i],te[i]),mul(tb[i],tc[i])),tb[i]=tc[i]=te[i]=0;
  NTT(ta,len,1);
  for (int i=l;i<=r;++i) res[i]=ta[i-l],ta[i-l]=0;
}

void Poly_inter(int *x,int *y,int n,int *res){
  Divide_prod(x,0,n,1);
  polynow[1]=new int[n+1];
  for (int i=1;i<=n+1;++i) polynow[1][i-1]=mul(prod[1][i],i);
  Divide_polynow(0,n,1,res);
  for (int i=0;i<=n;++i) res[i]=mul(y[i],Get_inv(res[i]));
  Divide_inter(res,0,n,1);
}