一、引入
首先，定义多项式的形式为 (f(x)=sum_{i=0}^n a_ix^i)，其中 (a_i) 为系数，(n) 为次数，这种表示方法称为“系数表示法”，一个多项式是由其系数确定的。
可以证明，(n+1) 个点可以唯一确定一个 (n) 次多项式。对于 (f(x))，代入 (n+1) 个不同的 (x)，得到 (n+1) 个不同的 (y)。一个 (n) 次的多项式就可以等价地换成 (n+1) 个等式，相当于平面上的 (n+1) 组坐标 ((x_i,y_i))，这种表示方法称为“点值表示法”。
多项式乘法 (卷积)：设 (C(x)=A(x)cdot B(x))，(A,B,C) 的系数构成的数列分别为 (a,b,c)，则 (c_k=sum_{i=0}^ka_ib_{k-i})。理解：因为 (x^i imes x^{k-i}=x^k)，(a_i) 和 (b_{k-i}) 相乘后，它们后面的未知数就变成了 (x^k)，对 (c_k) 产生贡献。
暴力求解：两个 (n) 次多项式相乘，时间复杂度 (mathcal{O}(n^2))。
快速傅里叶变换（Fast Fourier Transform，简称 FFT）可以 (mathcal{O}(nlog n)) 求解。
二、基本步骤
利用点值表示法，分三步快速求出多项式乘积：
- 由系数表示法转换成点值表示法。（DFT）
- 利用点值表示法，求两个多项式的乘积。
- 再将点值表示法转化成系数表示法。（IDFT）
对于步骤二，给出两个 (n) 次多项式 (A) 和 (B) 的点值表达式，我们可以 (mathcal{O}(n)) 求出其乘积 (C) 的点值表达式。显然 (C) 是 (2n) 次的，取 (2n+1) 个不同的 (x_i)。
- 代入 (A)：({(x_0,y_0),(x_1,y_1),cdots,(x_{2n},y_{2n})})。
- 代入 (B)：({(x_0,y'_0),(x_1,y'_1),cdots,(x_{2n},y'_{2n})})。
- 则 (C) 的点值表达式为：({(x_0,y_0y'_0),(x_1,y_1y'_1),cdots,(x_{2n},y_{2n}y'_{2n})})。
下面重点分析步骤一。有一些前置概念。
三、复数
我们把形如 (z=a+bi)（(a,b) 为实数）的数称为 复数，其中，(a,b) 分别叫做复数 (z) 的实部与虚部，(i) 为虚数单位，(i^2=-1)。
我们把复数表示在 复平面 上（(x) 轴叫实轴，(y) 轴叫虚轴），就像把实数表示在数轴上。如图，复数 (z=a+bi) 与复平面内的点 (Z(a,b)) 一一对应。
连接 (OZ)。复数 (z=a+bi) 与平面向量 (vec{OZ}) 一一对应。
复数的四则运算：
- 加减法：((a+bi)pm (c+di)=(apm c)+(bpm d)i)。
- 乘法：((a+bi)(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(bc+ad)i)。
- 除法：(frac{a+bi}{c+di}=frac{(a+bi) imes (c-di)}{(c+di) imes(c-di)}=frac{(ac+bd)+(bc-ad)i}{c^2+d^2}=frac{ac+bd}{c^2+d^2}+frac{bc-ad}{c^2+d^2}i)。
complex<double>a; //STL。a.real() 返回复数 a 的实部，a.imag() 返回复数 a 的虚部 struct cp{ //手写，比 STL 快一点 double a,b; cp operator+(cp &x){return (cp){a+x.a,b+x.b};} cp operator-(cp &x){return (cp){a-x.a,b-x.b};} cp operator*(cp &x){return (cp){a*x.a-b*x.b,a*x.b+b*x.a};} cp operator/(cp &x){double v=x.a*x.a+x.b*x.b; return (cp){(a*x.a+b*x.b)/v,(b*x.a-a*x.b)/v};} //除法一般不用 qwq };
我们称 (overline{z}=a-bi) 为复数 (z=a+bi) 的 共轭复数。
对于复数 (z=a+bi)，模长 (|z|) 为点 (Z(a,b)) 到原点的距离，幅角 为对应的向量 (vec{OZ}) 与横轴正半轴的夹角。
复数的乘法可以表达为：模长相乘，辐角相加。
四、单位根
1. 定义
(n) 次单位根是满足 (x^n=1) 的复数 (x)。
首先，单位根的模长必然为 (1)。因为若 (|x|>1)，则 (|x^n|=|x|^n>1)；若 (|x|<1)，则 (|x^n|=|x|^n<1)。
所以单位根表示的点一定在 单位圆（圆心为原点，半径为 (1)）上。
其次，单位根的辐角 ( heta) 一定满足 (frac{n heta}{2pi}inmathbb{Z})。也就是一个向量从 ((1,0)) 开始，每次旋转 ( heta) 的角度，旋转 (n) 次后还落在 ((1,0)) 上，那么它一定旋转了整数圈。
然后发现 (n) 次单位根正好是模长为 (1)，辐角为 (frac{2kpi}{n}) 的向量对应的复数。
记模长为 (1)，辐角为 (frac{2kpi}{n}) 的向量对应的 (n) 次单位根为 (omega_n^k)，称为第 (k) 个 (n) 次单位根。
还能发现，(omega_n^k=omega_n^{kmod n})，所以一般情况下，我们认为 (n) 次单位根有 (n) 个，即 (omega_n^0,omega_n^1,cdots,omega_n^{n-1})。
2. 性质
单位根的性质：（这些性质在后文会被用到）
-
性质 1：(n) 次单位根对应的向量将单位圆 (n) 等分。
两个相邻的 (n) 次单位根对应的向量的夹角相等。单位根的辐角是周角的 (frac{1}{n})。 -
性质 2：(omega_n^k=omega_n^{kmod n})。
在弧度制下，任意弧度 ( heta) 与 ( heta+2kpi\,(kinmathbb{Z})) 表示相同的角。 -
性质 3：({(omega_n^k)}^p=omega_n^{kp})。
第 (k) 个 (n) 次单位根对应向量的辐角变为原来的 (p) 倍，相当于 (omega_n^{kp}) 对应的向量。 -
性质 4：(omega_{dn}^{dk}=omega_{n}^k)。
考虑两者对应向量的辐角，(frac{2dkpi}{dn}=frac{2kpi}{n})。也可以这样理解：(dn) 次单位根对应的向量将单位圆 (dn) 等分，取第 (dk) 个。(n) 次单位根对应的向量将单位圆 (n) 等分，取第 (k) 个。两者等价。 -
性质 5：(omega_n^{k+n/2}=-omega_n^k)。其中 (n) 为偶数。
相当于一个复数对应的向量进行一次中心对称，(a+bi) 变为 (-a-bi)。
根据性质 3 有，({(omega_n^k)}^2=omega_{n}^{2k})。根据性质 4 有，(omega_n^{2k}=omega_{n/2}^k)，其中 (n) 为偶数。
3. 求法
根据性质 3，有 (omega_n^k=(omega_n^1)^k)。也就是说，只要求出 (omega_n^1)，就能得到 (omega_n^0,omega_n^1,cdots,omega_n^{n-1})。
(omega_n^1) 所对应的向量模长为 (1)，辐角为 (frac{2pi}{n})，得到 (omega_n^1) 所对应的点为 ((cos(frac{2pi}{n}),sin(frac{2pi}{n})))。
求 (pi)：double pi=acos(-1)。
补充：(omega_n^k=e^{frac{2pi ik}{n}}=cos(frac{2pi k}{n})+icdot sin(frac{2pi k}{n}))。其中 (i) 为虚数单位。
五、DFT
将系数表示法转换成点值表示法。
1. 基本思路
对于 (n-1) 次多项式（也就是有 (n) 项） (f(x)=sum_{i=0}^{n-1} a_ix^i)，将奇偶次数分离。
(f(x)=(a_0+a_2x^2+a_4x^4+cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+cdots+a_{n-1}x^{n-1}))
定义两个新的多项式 (f_1(x)) 和 (f_2(x))：
- (f_1(x)=a_0+a_2x+a_4x^2+cdots+a_{n-2}x^{{n/2-1}})。
- (f_2(x)=a_1+a_3x+a_5x^2+cdots+a_{n-1}x^{n/2-1})。
于是有 (f(x)=f_1(x^2)+xf_2(x^2))。
将 (omega_n^k\,(k<frac{n}{2})) 代入得：(f(omega_n^k)=f_1(omega_n^{2k})+omega_n^kf_2(omega_n^{2k})=f_1(omega_{n/2}^k)+omega_n^kf_2(omega_{n/2}^k))。
同理，将 (omega_n^{k+n/2}\,(k<frac{n}{2})) 代入得：(f(omega_n^{k+n/2})=f_1(omega_n^{2k+n})+omega_n^{k+n/2}f_2(omega_n^{2k+n}))
(=f_1(omega_n^{2k})+omega_n^{k+n/2}f_2(omega_n^{2k})=f_1(omega_{n/2}^k)+omega_n^{k+n/2}f_2(omega_{n/2}^k)=f_1(omega_{n/2}^k)-omega_n^kf_2(omega_{n/2}^k))。
发现两者的右边只有正负号的区别。
第一个式子的 (k) 取遍 ([0,frac{n}{2}-1]) 时，(k+frac{n}{2}) 取遍 ([frac{n}{2},n-1])。
如果我们知道 (f_1(x),f_2(x)) 分别在 (x=omega_{n/2}^0,omega_{n/2}^1,cdots,omega_{n/2}^{n/2-1}) 的点值表示，就可以 (mathcal{O}(n)) 求出 (f(x)) 在 (x=omega_n^0,omega_n^1,cdots,omega_n^{n-1}) 的点值表示。
然后，发现 (f_1(x),f_2(x)) 和 (f(x)) 的性质完全相同，这样就把问题分成了两个子问题，对于这两个子问题再进行递归求解。这样就可以在 (mathcal{O}(nlog n)) 的时间复杂度内求出点值表达式。
2. 代码实现
DFT 是利用单位根的特殊性质进行分治。
考虑到它能处理的多项式长度只能为 (2^k)（(k) 为整数），否则在分治时左右的项数就会不同。我们可以在最高次补一些系数为 (0) 的项，把原来的 (n) 补到 (2^k\,(2^kgeq n))。这样不影响计算结果。
别忘了：(f(omega_n^k)=f_1(omega_{n/2}^k)+omega_n^kf_2(omega_{n/2}^k),f(omega_n^{k+n/2})=f_1(omega_{n/2}^k)-omega_n^kf_2(omega_{n/2}^k))。
递归实现 DFT：
#include<bits/stdc++.h> #define int long long using namespace std; const int N=1e6+5; int n,len; double x,pi=acos(-1); complex<double>a[N]; void DFT(complex<double>*a,int n){ if(n==1) return ; //边界 int m=n/2; complex<double>a1[m],a2[m]; for(int i=0;i<m;i++) a1[i]=a[i<<1],a2[i]=a[i<<1|1]; //按奇偶分类 DFT(a1,m),DFT(a2,m); //处理子问题 complex<double>x(cos(2*pi/n),sin(2*pi/n)),w(1,0); //x: w_n^1。w: 当前的 n 次单位根，初始值为 w_n^0，即 1。 for(int i=0;i<m;i++) a[i]=a1[i]+w*a2[i],a[i+m]=a1[i]-w*a2[i],w*=x; //w*=x: 得到下一个单位根。 } signed main(){ scanf("%lld",&n); for(int i=0;i<n;i++) scanf("%lf",&x),a[i]=x; //a[i] 的实部赋值为 x for(len=1;len<n;len<<=1); //把项数补到 2 的整幂，高次项的系数默认为 0 DFT(a,len); for(int i=0;i<len;i++) printf("(%.4lf,%.4lf) ",a[i].real(),a[i].imag()); return 0; }
六、IDFT
DFT 的逆运算。将点值表示法转化成系数表示法。
结论：把 DFT 中的 (omega_n^1) 换成它的共轭复数，即 ((cos(frac{2pi}{n}),-sin(frac{2pi}{n})))，得到的系数再除以 (n) 即可。证明略。
可以把 IDFT 和 DFT 放在一起写。
void FFT(complex<double>*a,int n,int opt){ //opt=1 为 DFT，opt=-1 为 IDFT if(n==1) return ; int m=n/2; complex<double>a1[m],a2[m]; for(int i=0;i<m;i++) a1[i]=a[i<<1],a2[i]=a[i<<1|1]; FFT(a1,m,opt),FFT(a2,m,opt); complex<double>x(cos(2*pi/n),sin(2*pi/n)*opt),w(1,0); for(int i=0;i<m;i++) a[i]=a1[i]+w*a2[i],a[i+m]=a1[i]-w*a2[i],w*=x; }
七、迭代实现
目前的代码如下。可以发现它的效率并不是很高。
//Luogu P3803 #include<bits/stdc++.h> using namespace std; const int N=3e6+5; int n,m,len; double x,pi=acos(-1); complex<double>a[N],b[N]; void FFT(complex<double>*a,int n,int opt){ //opt=1/-1: DFT/IDFT if(n==1) return ; int m=n/2; complex<double>a1[m],a2[m]; for(int i=0;i<m;i++) a1[i]=a[i<<1],a2[i]=a[i<<1|1]; FFT(a1,m,opt),FFT(a2,m,opt); complex<double>x(cos(2*pi/n),sin(2*pi/n)*opt),w(1,0); for(int i=0;i<m;i++) a[i]=a1[i]+w*a2[i],a[i+m]=a1[i]-w*a2[i],w*=x; //蝴蝶操作（只是一个名字 qwq）。这里 w*a2[i] 算了两次，先记录下来再算可以减小常数。 } signed main(){ scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) scanf("%lf",&x),a[i]=x; for(int i=0;i<=m;i++) scanf("%lf",&x),b[i]=x; n=n+m+1; //一个 n 次多项式和一个 m 次多项式的乘积是一个 n+m 次多项式 (有 n+m+1 项) for(len=1;len<n;len<<=1); FFT(a,len,1),FFT(b,len,1); for(int i=0;i<len;i++) a[i]*=b[i]; //点值直接乘 FFT(a,len,-1); for(int i=0;i<n;i++) printf("%d%c",(int)(a[i].real()/len+0.5),i==n-1?' ':' '); //注意这里是除以 len return 0; }
如图，考虑递归的结构：
求出第 (4) 层的状态，就能往上合并求出前 (3) 层的状态。
观察发现，最后的数组下标的序列是原序列的 二进制翻转。比如 (6=(110)_2)，反过来就是 ((011)_2=3)。而第 (4) 层 (a_3) 的位置就是原来 (a_6) 的位置。
//r[i] 表示 i 二进制翻转后的结果。求 0~len-1 在二进制位数为 log2(len)-1 意义下的二进制翻转。 for(int i=0;i<len;i++) //len 是 2 的幂次 r[i]=(r[i>>1]>>1)|((i&1)?len>>1:0);
理解：考虑 (i) 与 (frac{i}{2}) 在二进制下的关系。(i) 可以看作是 (frac{i}{2}) 在二进制下的每一位左移一位得到。翻转后，(i) 是 (frac{i}{2}) 在二进制下的每一位右移一位得到，然后判一下最后一位即可。
迭代实现：
//Luogu P3803 #include<bits/stdc++.h> using namespace std; const int N=3e6+5; int n,m,len,r[N]; double x,pi=acos(-1); complex<double>a[N],b[N]; void FFT(complex<double>*a,int n,int opt){ //opt=1/-1: DFT/IDFT for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]); //求出最后一层的序列 for(int k=2;k<=n;k<<=1){ //枚举区间长度 int m=k>>1; //待合并的长度 complex<double>x(cos(2*pi/k),sin(2*pi/k)*opt),w(1,0),v; for(int i=0;i<n;i+=k,w=1) //枚举起始点 for(int j=i;j<i+m;j++) //遍历区间 v=w*a[j+m],a[j+m]=a[j]-v,a[j]=a[j]+v,w*=x; //蝴蝶操作。注意先 a[j+m]=a[j]-v 再 a[j]=a[j]+v。 } } signed main(){ scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) scanf("%lf",&x),a[i]=x; for(int i=0;i<=m;i++) scanf("%lf",&x),b[i]=x; n=n+m+1; for(len=1;len<n;len<<=1); for(int i=0;i<len;i++) //二进制翻转 r[i]=(r[i>>1]>>1)|((i&1)?len>>1:0); FFT(a,len,1),FFT(b,len,1); for(int i=0;i<len;i++) a[i]*=b[i]; FFT(a,len,-1); for(int i=0;i<n;i++) printf("%d%c",(int)(a[i].real()/len+0.5),i==n-1?' ':' '); return 0; }
八、模板
以 P3803 【模板】多项式乘法（FFT） 为例。
递归实现：
#include<bits/stdc++.h> using namespace std; const int N=3e6+5; int n,m,len; double x,pi=acos(-1); complex<double>a[N],b[N]; void FFT(complex<double>*a,int n,int opt){ //opt=1/-1: DFT/IDFT if(n==1) return ; int m=n/2; complex<double>a1[m],a2[m]; for(int i=0;i<m;i++) a1[i]=a[i<<1],a2[i]=a[i<<1|1]; FFT(a1,m,opt),FFT(a2,m,opt); complex<double>x(cos(2*pi/n),sin(2*pi/n)*opt),w(1,0); for(int i=0;i<m;i++) a[i]=a1[i]+w*a2[i],a[i+m]=a1[i]-w*a2[i],w*=x; } signed main(){ scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) scanf("%lf",&x),a[i]=x; for(int i=0;i<=m;i++) scanf("%lf",&x),b[i]=x; n=n+m+1; for(len=1;len<n;len<<=1); FFT(a,len,1),FFT(b,len,1); for(int i=0;i<len;i++) a[i]*=b[i]; FFT(a,len,-1); for(int i=0;i<n;i++) printf("%d%c",(int)(a[i].real()/len+0.5),i==n-1?' ':' '); return 0; }
迭代实现：（比递归快）
#include<bits/stdc++.h> using namespace std; const int N=3e6+5; int n,m,len,r[N]; double x,pi=acos(-1); complex<double>a[N],b[N]; void FFT(complex<double>*a,int n,int opt){ //opt=1/-1: DFT/IDFT for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]); for(int k=2;k<=n;k<<=1){ int m=k>>1; complex<double>x(cos(2*pi/k),sin(2*pi/k)*opt),w(1,0),v; for(int i=0;i<n;i+=k,w=1) for(int j=i;j<i+m;j++) v=w*a[j+m],a[j+m]=a[j]-v,a[j]=a[j]+v,w*=x; } } signed main(){ scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) scanf("%lf",&x),a[i]=x; for(int i=0;i<=m;i++) scanf("%lf",&x),b[i]=x; n=n+m+1; for(len=1;len<n;len<<=1); for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)?len>>1:0); FFT(a,len,1),FFT(b,len,1); for(int i=0;i<len;i++) a[i]*=b[i]; FFT(a,len,-1); for(int i=0;i<n;i++) printf("%d%c",(int)(a[i].real()/len+0.5),i==n-1?' ':' '); return 0; }
记忆：
- (f(x)=f_1(x^2)+xf_2(x^2))。
- (f(omega_n^k)=f_1(omega_{n/2}^k)+omega_n^kf_2(omega_{n/2}^k))。
- (f(omega_n^{k+n/2})=f_1(omega_{n/2}^k)-omega_n^kf_2(omega_{n/2}^k))。