FFT 学习笔记
\(\text{FFT}\) 学习笔记
建议先学习普通生成函数。
多项式
确定一个多项式，往往只需要知道每一次项前的系数是多少即可。众所周知，一个朴素的多项式往往可以被写成
的形式，在这种形式下的两个多项式 \(f,g\) 的乘积 \(h\) 往往可以按照
的方式计算。也就是说，得到的新多项式的每一次项前的系数可以写作
这种形式被叫做卷积，其特征是两个相乘的数所对应的下标相加是一个定值。不难发现，如此计算，对于两个次数为 \(n\) 的多项式来说，复杂度是 \(O(n^2)\)，然而这非常的慢，这个时候就需要利用 \(\text{FFT}\)，即快速傅里叶变换进行求解。
考虑 \(\text{FFT}\) 是如何优化的。对于一个 \(n\) 次多项式来说，只需要给定 \(n+1\) 个点，便可以确定一个唯一的 \(n\) 次多项式。那么我们在两个 \(n\) 次函数上取横坐标相同的 \(n+1\) 个点 \((x_i,y_i)\)，那么对于新函数来说，新得到的 \(n+1\) 个点的坐标一定为 \((x_i,{y_1}_i{y_2}_i)\)。不难发现，此时，我们只需要 \(O(n)\) 的时间复杂度即可确定新得到的多项式。那么，现在需要做得就是对于多项式的系数表示和点值表示进行一个快速的转换。
复数
即形如 \(a+bi\) 形式的数，其中 \(a,b\) 为实数，\(i\) 为虚数单位，即 \(\sqrt{-1}\)。
接着引入复平面，也就是以实数为 \(x\) 轴，虚数为 \(y\) 轴建立的平面直角坐标系，一个复数都在复平面上对应一个点。
复数的加减只需要将实数部分和虚数部分相加减即可。对于乘法则有
单位根
对于一个 \(n\) 次方程，一般来说都有 \(n\) 个根，而对于方程
的这 \(n\) 个根，对应到复平面上来看，我们以原点为圆心，做一个半径为 \(1\) 的圆，称为单位圆，并将其 \(n\) 等分，与 \(x\) 正半轴相交的点为一个根（即 \(1\)）。接着从 \(x\) 轴出发逆时针旋转到第 \(1\) 个点，这个点即为单位根，记作 \(w_n\)。接下来将第 \(i\) 个点记作 \(w_n^i\)，不难发现有以下性质：
- \(w_{n}^k=w_{n}^{n+k}\)
- \(w_n^k=w_{2n}^{2k}\)
- \(w_{2n}^{n+k}=-w_{2n}^{k}\)
- \((w_n^k)^a=w_n^{ak}\)
- \(w_n^k\) 对于 \(k\in[0,n-1]\) 互不相等
对于单位根的计算，可以利用三角函数简单求解，即
离散傅里叶变换
即 \(\text{DFT}\)，是一种将系数表示快速转换成点值表示的变换。
对于一个 \(2n-1\) 次多项式 \(f(x)\)：
如果我们分别设
那么显然有 \(f(x)=g(x^2)+xh(x^2)\)，不妨代入 \(x=w_{2n}^k\)，那么有
同理有
因此当我们求出 \(g(w_n^k)\) 和 \(h(w_n^k)\) 后便可以得到 \(f(w_{2n}^k)\) 和 \(f(w_{2n}^{n+k})\)。
由此我们考虑对于多项式 \(f\) 的 \(2n\) 个点分别代入 \(x=w_{2n}^i\)，那么每一个值均可以通过交换系数后求解 \(g,h\)，接着通过求解出的 \(g,h\) 更新所求的 \(f\)。为了让合并容易进行，不妨将次数补齐到 \(2^k-1\) 次。复杂度是 \(O(n\log n)\) 的。
但是这样写常数很大，我们考虑模拟交换系数的过程：
我们考虑这些数的二进制表示：
不难发现最初系数的位置和最终系数的位置满足下标进行了一次二进制下的翻转。
于是我们不难联想到与处理出位置是如何变化的，并存储在 \(r\) 数组中，接着考虑合并时本质上是将 \(2n\) 个位置中的 \(k,n+k\) 合并至 \(k,n+k\)，因此可以直接进行覆盖。这样我们就完成了常数上的优化。
代码
void DFT(int n){
Set(n-1);
for(int i=0;i<n;i++){
if(i<r[i])swap(a[i],a[r[i]]);//按照预处理交换系数
}
for(int len=1;len<n;len<<=1){
Complex step(cos(Pi/len),sin(Pi/len));//单位根
for(int l=0;l<n;l+=(len<<1)){
Complex w(1,0);
for(int k=0;k<len;k++,w=w*step){
Complex A=a[l+k],B=w*a[l+len+k];
a[l+k]=A+B;
a[l+len+k]=A-B;
}
}
}
return ;
}
逆离散傅里叶变换
即 \(\text{IDFT}\)，也就是将函数从点值表示转换成系数表示的变换。
考虑 \(\text{DFT}\) 在线性代数上的表示，即为
记点值矩阵为 \(E\)，中间的矩阵为 \(D\)，系数矩阵为 \(V\)，则 \(E=DV\)，那么有 \(V=ED^{-1}\)，也就是我们只需要求出 \(D\) 矩阵的逆矩阵即可。
考虑构造 \(A_{ij}=D_{ij}^{-1}\)，计算 \(AD\)：
不难发现 \(\frac{AD}{n}=I\)，所以有 \(D^{-1}=\frac{A}{n}\)。更具体的，因为 \(A_{ij}=D_{ij}^{-1}\)，所以每个位置代入的数就变成了 \(w_{2n}^{-i}\)，与 \(w_{2n}^{i}\) 关于 \(x\) 轴对称，也就是虚数部分相反。只需要对于 \(\text{DFT}\) 函数传于一个参数 \(\text{type}\) 表示如何转换即可。
代码
void DFT(int n,int type){
Set(n-1);
for(int i=0;i<n;i++){
if(i<r[i])swap(a[i],a[r[i]]);
}
for(int len=1;len<n;len<<=1){
Complex step(cos(Pi/len),type*sin(Pi/len));//是否需要取反
for(int l=0;l<n;l+=(len<<1)){
Complex w(1,0);
for(int k=0;k<len;k++,w=w*step){
Complex A=a[l+k],B=w*a[l+len+k];
a[l+k]=A+B;
a[l+len+k]=A-B;
}
}
}
if(type==-1){
for(int i=0;i<n;i++)a[i].a/=n;//除以 n 才是正确结果
}
return ;
}
下面可以给出一个带封装的总代码：
const int N=2097152;//注意需要取到最大次数向上的最大二次幂数
const double Pi=acos(-1.0);//处理 2pi
int n,m;
int r[N];//预处理数组
struct Complex{
double a,b;
Complex(double a=0,double b=0):a(a),b(b){}
Complex operator +(Complex x){return {a+x.a,b+x.b};}
Complex operator -(Complex x){return {a-x.a,b-x.b};}
Complex operator *(Complex x){return {a*x.a-b*x.b,a*x.b+b*x.a};}
};//手写复数
struct Polynomial{
int len;
vector<Complex>a;
Complex& operator [](int x){return a[x];}
void Set(int len){this->len=len,a.resize(len+1);}//节省空间的做法
void DFT(int n,int type){
Set(n-1);
for(int i=0;i<n;i++){
if(i<r[i])swap(a[i],a[r[i]]);
}
for(int len=1;len<n;len<<=1){
Complex step(cos(Pi/len),type*sin(Pi/len));
for(int l=0;l<n;l+=(len<<1)){
Complex w(1,0);
for(int k=0;k<len;k++,w=w*step){
Complex A=a[l+k],B=w*a[l+len+k];
a[l+k]=A+B;
a[l+len+k]=A-B;
}
}
}
if(type==-1){
for(int i=0;i<n;i++)a[i].a/=n;
}
return ;
}
Polynomial operator *(Polynomial x){
Polynomial y=*this,z;
int len=x.len+y.len;
int n=1;
while(n<=len)n<<=1;//向上取到二次幂
z.Set(n-1);
for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)*(n>>1));
x.DFT(n,1);
y.DFT(n,1);
for(int i=0;i<n;i++)z[i]=x[i]*y[i];
z.DFT(n,-1);
z.Set(len);
return z;
}
};
Polynomial f,g;
int main(){
ios::sync_with_stdio(0);
cin.tie(0);
cin>>n>>m;
f.Set(n);
g.Set(m);
for(int i=0;i<=n;i++)cin>>f[i].a;
for(int i=0;i<=m;i++)cin>>g[i].a;
f=f*g;
for(int i=0;i<=f.len;i++)cout<<(int)(f[i].a+0.5)<<" ";//精度问题
return 0;
}
快速数论变换
即 \(\text{NTT}\)，考虑到 \(\text{FFT}\) 存在精度问题，对于模意义下的数就需要用到 \(\text{NTT}\) 了。
考虑到此时单位根不再使用，我们考虑模意义下能够代替单位根的东西，也就是上面讲到的原根。考虑原根如何能够满足代替单位根进行运算，考虑一些性质。
-
\(w_n^k\) 对于 \(k\in[0,n-1]\) 互不相等。
若模数 \(p\) 是一个质数，那么对于原根 \(g\)，\(g^k\) 对于 \(k\in[0,p-1]\) 在模 \(p\) 意义下互不相等，那么 \((g^{\frac{p-1}{n}})^k\) 对于 \(k\in[0,n-1]\) 在模 \(p\) 意义下互不相等。
-
\(w_n^k=w_n^{n+k}\)。
\((g^{\frac{p-1}{n}})^{n+k}=g^{p-1}(g^{\frac{p-1}{n}})^k\equiv(g^{\frac{p-1}{n}})^k\pmod p\)
-
\(w_n^k=w_{2n}^{2k}\)。
\((g^{\frac{p-1}{n}})^k=g^{\frac{k(p-1)}{n}}=g^{\frac{2k(p-1)}{2n}}=(g^{\frac{p-1}{2n}})^{2k}\pmod p\)
-
\(w_{2n}^{n+k}=-w_{2n}^{k}\)。
\(\frac{(g^{\frac{p-1}{2n}})^{n+k}}{(g^{\frac{p-1}{2n}})^k}=(g^{\frac{p-1}{2n}})^n=g^{\frac{p-1}{2}}\equiv -1\pmod p\)，所以 \((g^{\frac{p-1}{2n}})^{n+k}\equiv-(g^{\frac{p-1}{2n}})^k\)。
-
\((w_n^k)^a=w_n^{ak}\)。
\(((g^{\frac{p-1}{n}})^k)^a=(g^{\frac{p-1}{n}})^{ak}\)。
因此，不难看出，单位根即为 \(g^{\frac{p-1}{n}}\)，同时，这也对 \(p\) 作出了相应的限制，假设运算过程中多项式的次数最多到达 \(2^t\)，那么 \(p=a2^\alpha+1,\alpha\ge t\)。下面是一些常见的 \(\text{NTT}\) 模数。
那么当题目给出的模数不是符合要求的 \(\text{NTT}\) 模数时应该怎么办呢——大骂出题人。
代码
#define ll long long
const int N=1e7+5,M=1e5+5,R=262144,P=1004535809,G=3,GI=334845270;
int n,m,t,s;
int w[M],r[R];
ll ans;
ll f[N],g[N];
namespace Math{
ll QuickPow(ll a,int b,const int p){
ll res=1;
while(b>0){
if((b&1)==1)res=res*a%p;
a=a*a%p;
b>>=1;
}
return res;
}
ll Inv(int x,int p){
return QuickPow(x,p-2,p);
}
};
struct Polynomial{
int len;
vector<ll>a;
ll& operator [](int x){return a[x];}
void Set(int len){this->len=len,a.resize(len+1);}
void NTT(int n,int type){
Set(n-1);
for(int i=0;i<n;i++){
if(i<r[i])swap(a[i],a[r[i]]);
}
for(int len=1;len<n;len<<=1){
ll step=Math::QuickPow(type==1?G:GI,(P-1)/(len<<1),P);
for(int l=0;l<n;l+=(len<<1)){
ll w=1;
for(int k=0;k<len;k++,w=w*step%P){
ll A=a[l+k],B=w*a[l+len+k]%P;
a[l+k]=(A+B)%P;
a[l+len+k]=(A-B)%P;
}
}
}
if(type==-1){
ll inv=Math::Inv(n,P);
for(int i=0;i<n;i++)a[i]=a[i]*inv%P;
}
return ;
}
Polynomial operator *(Polynomial x){
Polynomial y=*this,z;
int len=x.len+y.len;
int n=1;
while(n<=len)n<<=1;
z.Set(n-1);
for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)*(n>>1));
x.NTT(n,1);
y.NTT(n,1);
for(int i=0;i<n;i++)z[i]=x[i]*y[i]%P;
z.NTT(n,-1);
z.Set(len);
return z;
}
}A,B;
int main(){
ios::sync_with_stdio(0);
cin.tie(0);
cin>>n>>m;
A.Set(n);
B.Set(m);
for(int i=0;i<=n;i++)cin>>A[i];
for(int i=0;i<=m;i++)cin>>B[i];
A=A*B;
for(int i=0;i<=A.len;i++)cout<<A[i]<<" ";
return 0;
}
常见模型
事实上，所有的模型都离不开以下两种情况：
上面的就是朴素的卷积，而下面的称作差卷积，即两个数下标的差为定值。对于这种情况，我们通常选择讲一个序列倒转，即构造 \(b'_i=b_{n-i}\)，那么有
于是发现和始终为 \(n-k\)，可以求解。
字符串匹配
虽然很不想承认，但是 \(\text{FFT}\) 确实适用于字符串匹配的情景，考虑如何维护两个串是否相同，相同则意味着每一位都相同，因此我们构造
不难发现，当 \(c_k=0\) 时，意味着从模式串的第 \(k\) 位开始能够成功匹配，那么进行一些化简：
很显然，除了中间那一项，其他的项均可以通过预处理 \(O(1)\) 得出，而中间显然是一个差卷积形式，因此构造 \(b'_{i}=b_{m-1-i}\)，所以有
\(\text{FFT}\) 求解即可，复杂度是 \(O(n\log n)\) 的。而对于存在通配符的匹配，我们认为通配符对应值为 \(0\)，这意味着只要有通配符那么这一项就认为相等，因此构造
化简可得
不难看出只需要进行 \(3\) 次 \(\text{FFT}\) 即可。
背包问题
类似于求解给出 \(n\) 个数，选出 \(k\) 个数（可重复），求解和在某一范围内的方案数。不难想出构造 \(a_i\) 的值为有多少个值为 \(i\) 的数，所求即为询问范围内 \((a^k)_i\) 的和。需要注意，为了保证时间和空间的正确性，每次处理完后，对于值超过范围的数应当忽视，因为他们一定不符合要求。但是复杂度为 \(O(kn\log n)\) 依旧不优，考虑二进制优化，那么复杂度优化为 \(O(n\log n\log k)\)。
高精度乘法
据某不知名消息途径，\(\text{Python}\) 的高精度乘法就是用 \(\text{FFT}\) 实现的。
考虑一个十进制数的一种不常规写法
于是不妨假设 \(x=10\)，那么
这和多项式没有区别，于是两个数相乘就可以抽象成两个多项式相乘，因此可以利用 \(\text{FFT}\) 优化，复杂度从 \(O(n^2)\) 降至 \(O(n\log n)\)。
分治 \(\text{FFT}\)
针对 \(\prod f_i(x)\) 的形式，不难看出，如果一个一个相乘，复杂度是 \(O(n^2\log n)\) 的。不难结合乘法结合律通过分治将原式拆成长度相等的两段，分别计算后相乘即可，复杂度降至 \(O(n\log^2n)\)。
