题解 P3803 【【模板】多项式乘法（FFT）】

· · 题解

前言

众所周知，这FFT是用来算多项式乘法的。

对于这种常人思维难以理解的东西，就少些理解，多背板子吧！

因此只总结一下思路和代码，什么概念和推式子就靠巨佬们吧

推荐自为风月马前卒巨佬的概念和定理都非常到位的总结

推荐ppl巨佬的简明易懂的总结

多项式乘法的蹊径——点值表示法

一般我们把两个长度为n的多项式乘起来，就类似于做竖式乘法，一位一位地乘再加到对应位上，是O(n^2)的

如何优化？直接看是没有思路的，只好另辟蹊径了。

多项式除了我们常用的系数表示法y=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}以外，还可以用点值表示法。

所谓点值表示法，就是相当于用函数图像上n个点的坐标(x_i,y_i)表示一个n次多项式，显然这个表示是唯一的。

我们可以把系数表示转化成点值表示。点值表示下的多项式怎么相乘呢？就是选相同的x_i，把对应的y_i相乘。

当然，两个长度为n的多项式相乘得到的是长度为2n-1的多项式，需要2n-1个点值才能唯一表示，因此一开始两个多项式也要选2n-1个点表示。

举一个例子

(x+1)(x+2)$→ → → → $x^2+3x+2

↓\qquad\qquad\qquad\qquad\qquad\qquad\ \ \ ↑

↓（点值）\qquad\qquad\qquad （系数）↑

↓\qquad\qquad\qquad\qquad\qquad\qquad\ \ \ ↑

$(0,2)(1,3)(2,4)$→ → → →$(0,2)(1,6)(2,12)

可是，随意选O(n)个点，计算y是O(n)的，总时间还是O(n^2)的，把它还原成系数表达式更不轻松。

所以，如果选的点很普通，这只是一条蹊径，并不是一条捷径。

神奇的单位根

介绍一个神奇的东西——n次单位根（记为\omega_n）。

它是n个复数的集合，每一个的n次方都等于1。其中的一个是e^{{2\pi i}\over n}，记为\omega_n。

普及一下欧拉公式，e^{\theta i}=\cos\theta+(\sin\theta)i，\theta就是这个复数向量的旋转角。显然满足(\omega_n)^n=e^{2\pi i}=1，那么它的k次方(\omega_n^k)^n=e^{2k\pi i}也等于1。

于是可以看出，满足n次方等于1的复数的取值只会有n个，为\omega_n^k(k\in[0,n-1])，因为会有\omega_n^{n+k}=\omega_n^k。

这n个复数向量在单位圆上呈放射状。下面是算导的图片：

介绍消去引理ω^{dk}_{dn}=ω^k_n，证明很容易的。

DFT&IDFT

单位根有什么用呢？

看看我们把\omega_n^k(k\in[0,n-1])分别带入多项式求点值会发生什么就知道了。这个过程叫DFT。

假设n为偶数，那么我们可以把它的奇偶项分开，用两个多项式表示它

A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{{\frac n 2}-1} A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{{\frac n 2}-1}

那么A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)。

注意，下面的变化用到了ω^{2k}_n=ω^k_{\frac n 2}，ω^n_n=1和ω^{\frac n 2}_n=-1。

首先带入单位根

\quad A(ω^k_n) =A^{[0]}(ω^{2k}_n)+ω^k_nA^{[1]}(ω^{2k}_n) =A^{[0]}(ω^k_{\frac n 2})+ω^k_nA^{[1]}(ω^k_{\frac n 2})(k<\frac n 2)

那k\geq\frac n 2时又会发生什么呢？把它变成\frac n 2+k

\quad A(ω^{\frac n 2+k}_n) =A^{[0]}(ω^{n+2k}_n)+ω^{\frac n 2+k}_nA^{[1]}(ω^{n+2k}_n) =A^{[0]}(ω^{2k}_n)+ω^{\frac n 2}_nω^k_nA^{[1]}(ω^{2k}_n) =A^{[0]}(ω^k_{\frac n 2})-ω^k_nA^{[1]}(ω^k_{\frac n 2})

也就是说，如果我们要求一个长度为n的多项式取\omega_n^k(k\in[0,n-1])的n个点值，我们只需要求出两个长度为n/2的多项式取\omega_{\frac n 2}^k(k\in[0,\frac n 2-1])的\frac n 2个点值，再通过上述两个式子合并结果。这完全就是个递归过程，很容易写一个函数来计算。

可以由算法写出DFT的复杂度T(n)=2T(\frac n 2)+O(n)=O(n\log n)

当然，别忘了还原成系数表示，这个过程叫做IDFT。

蒟蒻觉得这里理解太麻烦了，因此不再证明IDFT的过程，想了解的参见其它的总结。

只需要记住，IDFT的\omega_n是e^{-\frac{2\pi i}n}，最后的结果除以n，其它过程同DFT，可以写在一个函数里。具体看下面的代码：

void FFT(R complex<double>*a,R int n,R int op){//op=1为DFT，-1为IDFT
    if(!n)return;//为了方便，n的意义与上面不一样，这里的n是a0、a1的长度
    complex<double>a0[n],a1[n];
    for(R int k=0;k<n;++k)
        a0[k]=a[k<<1],a1[k]=a[k<<1|1];//奇偶项分离
    FFT(a0,n>>1,op);FFT(a1,n>>1,op);//递归处理
    R complex<double>wn(cos(PI/n),sin(PI/n)*op),w(1,0);//单位根
    for(R int k=0;k<n;++k,w*=wn)//k从到n/2
        a[k]=a0[k]+w*a1[k],a[k+n]=a0[k]-w*a1[k];
}

递归版过程

引入例题：洛谷P3803 【模板】多项式乘法（FFT）

由于系数很小，我们不必担心精度的问题了（当然float是使不得的）

递归版代码：

#include<cstdio>
#include<cmath>
#include<complex>
#include<iostream>
#define R register
#define G c=getchar()
using namespace std;
const int N=4.2e6;
const double PI=acos(-1);//自定义π
complex<double>f[N],g[N];
inline int in(){
    R char G;
    while(c<'-')G;
    return c&15;
}
void FFT(R complex<double>*a,R int n,R int op){//同上
    if(!n)return;
    complex<double>a0[n],a1[n];
    for(R int i=0;i<n;++i)
        a0[i]=a[i<<1],a1[i]=a[i<<1|1];
    FFT(a0,n>>1,op);FFT(a1,n>>1,op);
    R complex<double>W(cos(PI/n),sin(PI/n)*op),w(1,0);
    for(R int i=0;i<n;++i,w*=W)
        a[i]=a0[i]+w*a1[i],a[i+n]=a0[i]-w*a1[i];
}
int main(){
    R int n,m,i;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)f[i]=in();
    for(i=0;i<=m;++i)g[i]=in();
    for(m+=n,n=1;n<=m;n<<=1);//把长度补到2的幂，不必担心高次项的系数，因为默认为0
    FFT(f,n>>1,1);FFT(g,n>>1,1);//DFT
    for(i=0;i<n;++i)f[i]*=g[i];//点值直接乘
    FFT(f,n>>1,-1);//IDFT
    for(i=0;i<=m;++i)printf("%.0f ",fabs(f[i].real())/n);//注意结果除以n，小心“-0”
    puts("");
    return 0;
}

好记又好写的递归版结果怎样呢？

Fast Fast TLE!只有77分。

最主要的原因在于，空间调用太多了。

蝴蝶操作和Rader排序

针对效率太低的问题，我们继续观察FFT实现过程中的更多特点。

观察这一句代码a[k]=a0[k]+w*a1[k],a[k+n]=a0[k]-w*a1[k];，在操作过程中，取出了a0[k]和a1[k]的值，通过计算又求出了a[k]和a[k+n]的值。我们把这样的一次运算叫做“蝴蝶操作”。

这样的操作有什么特点呢？我们试着将a0和a1合并成一个数组a，每一次蝴蝶操作后，取出了a[k]和a[k+n]的值，又求出了a[k]和a[k+n]的值。最后，整个数组都完成了求值，而并没有用到两个数组！

以n=8为例，看看递归过程的结构

其实，我们完全可以递推求解！假设a数组已经变成了第四层的样子，先对a0和a4、a2和a6、a1和a5、a3和a7蝴蝶操作，变成第三层；再对a0和a2、a4和a6、a1和a3、a5和a7蝴蝶操作，变成第二层；最后对a0和a1、a2和a3、a4和a5、a6和a7蝴蝶操作，变成第一层，FFT就完成了，而空间复杂度仅为O(n)。这个过程可以用循环来控制。

剩下的问题就是把初始的数组变成最后一层的样子了。先别急着写一个递归函数暴力把位置换过去。来观察一下最后序列的编号的二进制表示000,100,010,110,001,101,011,111，是不是与原来000,001,010,011,100,101,110,111相比，每个位置上的二进制位都反过来了？这样的变化叫做Rader排序。

我们可以O(n)将Rader排序的映射关系求出。设iRader排序后的数为r_i，我们可以通过r_{\frac i 2}递推求出，具体方法可以看看代码&动动脑筋qwq

#include<cmath>
#include<cstdio>
#define R register
#define I inline
using namespace std;
const int N=4.2e6;
const double PI=acos(-1);
int n,r[N];
struct C{//手写complex，比STL快一点点
    double r,i;
    I C(){r=i=0;}
    I C(R double x,R double y){r=x;i=y;}
    I C operator+(R C&x){return C(r+x.r,i+x.i);}
    I C operator-(R C&x){return C(r-x.r,i-x.i);}
    I C operator*(R C&x){return C(r*x.r-i*x.i,r*x.i+i*x.r);}
    I void operator+=(R C&x){r+=x.r;i+=x.i;}
    I void operator*=(R C&x){R double t=r;r=r*x.r-i*x.i;i=t*x.i+i*x.r;}
}f[N],g[N];
I int in(){
    R char c=getchar();
    while(c<'-')c=getchar();
    return c&15;
}
I void FFT(R C*a,R int op){
    R C W,w,t,*a0,*a1;
    R int i,j,k;
    for(i=0;i<n;++i)//根据映射关系交换元素
        if(i<r[i])t=a[i],a[i]=a[r[i]],a[r[i]]=t;
    for(i=1;i<n;i<<=1)//控制层数
        for(W=C(cos(PI/i),sin(PI/i)*op),j=0;j<n;j+=i<<1)//控制一层中的子问题个数
            for(w=C(1,0),a1=i+(a0=a+j),k=0;k<i;++k,++a0,++a1,w*=W)
                t=*a1*w,*a1=*a0-t,*a0+=t;//蝴蝶操作
}
int main(){
    R int m,i,l=0;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)f[i].r=in();
    for(i=0;i<=m;++i)g[i].r=in();
    for(m+=n,n=1;n<=m;n<<=1,++l);
    for(i=0;i<n;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));//递推求r
    FFT(f,1);FFT(g,1);
    for(i=0;i<n;++i)f[i]*=g[i];
    FFT(f,-1);
    for(i=0;i<=m;++i)printf("%.0f ",fabs(f[i].r)/n);
    puts("");
    return 0;
}