这东西对初中生挺友好的。
前置知识
-
复数
形如 \(a+bi(a,b\in \mathbb{R})\) 的数叫复数,其中 \(i^2=-1\)。
复数乘法:\((a+bi)(c+di)=ac-bd+(ad+bc)i\)。乘法分配律即可。 -
复平面
以 \(a\) 为 \(x\) 轴,\(b\) 为 \(y\)轴所组成的平面叫复平面。每个复数都对应复平面上一点。 -
单位根
以原点为圆心,\(1\) 为半径所形成的圆叫单位圆。以 \((1,0)\) 为起点,用 \(n\) 个点圆上的点将单位圆分成 \(n\) 等分。容易发现,这 \(n\) 个点正好是 \(w^n=1\) 的 \(n\) 个解。
如图为一个单位圆和一些单位根。
按角度大小记这 \(n\) 个点为 \(w_n^0,w_n^1\dots w_n^{n-1}\),当 \(n\) 为 \(2\) 的次幂时,有如下性质:
- \(w_n^k=-w_n^{k+\frac{n}{2}}\),也就是说他们两两互为相反数。
- \(w_n^k=w_n^{n+k}\),说明单位根有周期性。
- 当 \(n\) 个单位根平方后,他们依然是两两互为相反数。
这些性质决定了我们为什么选它来求值。
FFT
其实就是将两个多项式转化成点值表示相乘在转回来的过程。
设 \(F(x)=\sum\limits_{i=0}^m a_ix^i\),\(m\)是两个多项式相乘后的次数。默认 \(n\) 为最小的大于 \(m\) 的 \(2\) 的次幂(方便计算)。将 \(F\) 按次数奇偶分成两个函数。
则
显然 \(f_0,f_1\) 都是偶函数,也就是说 \(x\) 的值与 \(-x\) 是相等的。因此,我们还可以表示出
所以要取的点和多项式的次数都变为了原来的一半,可以用递归处理。但对于 \(f_0,f_1\) 而言所有的点值都是正的,不存在相反数,无法递归。所以我们尝试寻找一些更特殊的数,单位根正好满足这个性质,因为所有单位根平方后依然是两两互为相反数的。
举个例子:
这是当 \(n=4\) 时递归的 \(x\) 值。
这是当 \(n=8\) 时的 \(x\) 值。
当递归到只剩一个常数项时,递归结束。
我们可以用相同的方式将另一个多项式的点值求出来。
当我们求出两个多项式点值并相乘后,如何将点值重新转化为系数呢?
我个人认为这很人类智慧。
设相乘后的多项式为 \(G(x)=\sum\limits_{i=0}^{n-1} b_ix^i\),将 \(w_n^0,w_n^{-1},\dots w_n^{1-n}\) 带入得到的多项式中进行 \(\text{FFT}\),再将所得的多项式的每一项除以 \(n\) 就是每一项对应的系数。
证明 :
设 \(c_k=\sum\limits_{i=0}^{n-1}b_i(w_n^{-k})^i\),将 \(b_i\) 用多项式表示 \(b_i=\sum\limits_{j=0}^{n-1}a_j(w_n^i)^j\),带入上式中
当 \(j=k\) 时, \(c_k=na_j\)。
当 \(j \ne k\) 时
所以,\(c_k=na_j\)。\(\frac{c_k}{n}\) 就是我们要求的系数。
至此,\(\text{FFT}\) 的所有流程都已完成。
模板题
代码(递归)
bool _Start;
#include <bits/stdc++.h>
using namespace std;
#define il inline
#define Tp template<typename T>
#define Ts template<typename T,typename... _T>
Tp il void read(T& x) {x=0;bool f=0;char c=getchar();for(;!isdigit(c);c=getchar()) f|=c=='-';for(;isdigit(c);c=getchar()) x=(x<<1)+(x<<3)+(c^48);x=(f?-x:x);
}Ts il void read(T& x,_T&... y) {read(x),read(y...);}
using comp=complex<double>;
const double pi=acos(-1);
const int N=1<<21;
int n,m;
comp A[N],B[N];
void FFT(comp F[],int n,int tp) {if(n==1) return ;comp F0[(n>>1)+5],F1[(n>>1)+5];for(int i=0;i<n;i+=2)F0[i>>1]=F[i],F1[i>>1]=F[i+1];FFT(F0,n>>1,tp),FFT(F1,n>>1,tp);comp wn=exp(comp(0,2*pi*tp/n)),w=1;for(int i=0;i<(n>>1);i++,w*=wn) {F[i]=F0[i]+w*F1[i];F[i+(n>>1)]=F0[i]-w*F1[i];}
}
void mul(comp f[],comp g[]) {int lim=1;while(lim<=n+m) lim<<=1;FFT(f,lim,1),FFT(g,lim,1);for(int i=0;i<lim;i++) f[i]*=g[i];FFT(f,lim,-1);for(int i=0;i<lim;i++) f[i]/=lim;
}
bool _End;
int main() {fprintf(stderr,"Memory: %.4lf Mib\n",abs(&_End-&_Start)/1048576.0);read(n,m);for(int i=0,x;i<=n;i++) read(x),A[i]=x;for(int i=0,x;i<=m;i++) read(x),B[i]=x;mul(A,B);for(int i=0;i<=n+m;i++) printf("%d ",(int)(A[i].real()+0.5));return 0;
}
优化
balabalabala.