FFT快速傅里叶变换
FFT快速傅里叶变换
some 多项式变换
DFT：离散傅里叶变换 $\longrightarrow $ \(O(n^2)\)计算多项式乘法
FFT：快速傅里叶变换 $\longrightarrow $ \(O(n\log n)\)计算多项式乘法
FNTT/NTT：快速傅里叶变换的优化版 $\longrightarrow $ 优化常数及误差
FWT：快速沃尔什变换 $\longrightarrow $ 利用类似FFT的东西解决一类卷积问题
MTT：毛爷爷的FFT $\longrightarrow $ 非常nb/任意模数
FMT 快速莫比乌斯变化
快速傅里叶变换
就 \(DFT\) 来说，它分治地来求当 \(x=\omega_n^k\) 的时候 \(f(x)\) 的值。
\(FFT\) 的分治思想体现在将多项式分为奇次项和偶次项处理。
对于一个 \(n-1\) 次的多项式：
按照次数的奇偶来分成两组，然后右边提出来一个 \(x\)：
分别用奇偶次次项数建立新的函数：
那么原来的 \(f(x)\) 用新函数表示为：
利用偶数次单位根的性质：\(\omega^i_n = -\omega^{i + \frac{n}{2}}_n\)
\(G\left(x^2\right)\) 和 \(H\left(x^2\right)\) 是偶函数
我们知道在复平面上 \(\omega^i_n\) 和 \(\omega^{i+\frac{n}{2}}_n\) 的 \(G(x^2)\) 的 \(H(x^2)\) 对应的值相同。
得到：
和：
因此我们求出了 \(G(\omega_{\frac{n}{2}}^k)\) 和 \(H(\omega_{\frac{n}{2}}^k)\) 后，就可以同时求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+\frac{n}{2}})\)。
于是对 \(G\) 和 \(H\) 分别递归 \(DFT\) 即可。
考虑到分治 \(DFT\) 能处理的多项式长度只能是 \(2^m(m \in \mathbf{N}^ \ast )\)，否则在分治的时候左右不一样长，右边就取不到系数了。
所以要在第一次 \(DFT\) 之前就把序列向上补成长度为 \(2^m(m \in \mathbf{N}^\ast )\)（高次系数补 \(0\)）、最高项次数为 \(2^m-1\) 的多项式。
位逆序置换
原序列的每个数用二进制表示，然后把二进制翻转对称一下，就是最终位置的下标。
可以在 \(O(n\log n)\) 或 \(O(n)\) 的时间内求出每个数变换后的结果
蝶形运算优化
已知 \(G(\omega_{\frac{n}{2}}^k)\) 和 \(H(\omega_{\frac{n}{2}}^k)\) 后，需要使用下面两个式子求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+\frac{n}{2}})\)：
使用位逆序置换后，对于给定的 \(n\), \(k\)：
\(G(\omega_{\frac{n}{2}}^k)\) 的值存储在数组下标为 \(k\) 的位置，\(H(\omega_{\frac{n}{2}}^k)\) 的值存储在数组下标为 \(k + \dfrac{n}{2}\) 的位置。
\(f(\omega_n^k)\) 的值将存储在数组下标为 \(k\) 的位置，\(f(\omega_n^{k+\frac{n}{2}})\) 的值将存储在数组下标为 \(k + \dfrac{n}{2}\) 的位置。
因此可以直接在数组下标为 \(k\) 和
\(k + \frac{n}{2}\) 的位置进行覆写，而不用开额外的数组保存值。
再详细说明一下如何借助蝶形运算完成所有段长度为 \(\frac{n}{2}\) 的合并操作：
令段长度为 \(s = \frac{n}{2}\)；
同时枚举序列 \(\{G(\omega_{\frac{n}{2}}^k)\}\) 的左端点 \(l_g = 0, 2s, 4s, \cdots, N-2s\) 和序列 \(\{H(\omega_{\frac{n}{2}}^k)\}\) 的左端点 \(l_h = s, 3s, 5s, \cdots, N-s\)；
合并两个段时，枚举 \(k = 0, 1, 2, \cdots, s-1\)，此时\(G(\omega_\frac{n}{2}^k)\) 存储在数组下标为 \(l_g + k\) 的位置，\(H(\omega_\frac{n}{2}^k)\) 存储在数组下标为 \(l_h + k\) 的位置；
使用蝶形运算求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+\frac{n}{2}})\)，然后直接在原位置覆写。
快速傅里叶逆变换
设 \(b_i=\omega_n^{-i}\)，则多项式 \(A\) 在 \(x=b_0,b_1,\cdots,b_{n-1}\) 处的点值表示法为 \(\left\{ A(b_0),A(b_1),\cdots,A(b_{n-1}) \right\}\)。
对 \(A(x)\) 的定义式做一下变换，可以将 \(A(b_k)\) 表示为
记
当 \(a=0 \pmod{n}\) 时，\(S\left(\omega_n^a\right)=n\)。
当 \(a\neq 0 \pmod{n}\) 时，我们错位相减
也就是说
那么代回原式
也就是说给定点 \(b_i=\omega_n^{-i}\)，则 \(A\) 的点值表示法为
综上所述，我们取单位根为其倒数，对 \(\{y_0,y_1,y_2,\cdots,y_{n-1}\}\) 跑一遍 \(FFT\)，然后除以 \(n\) 即可得到 \(f(x)\) 的系数表示。
递归写法
il void FFT(qwq *q,int lim,int k){
if(lim==1) return;
for(ri int i=0;i<lim;++i) o[i]=q[i];
ri int m=lim>>1;
for(ri int i=0;i<lim;i++){
if(i&1) q[m+(i>>1)]=o[i];
else q[i>>1]=o[i];
}
FFT(q,m,k),FFT(q+m,m,k);
wn={cos(Pi*2/lim),k*sin(Pi*2/lim)},w={1,0};
for(ri int i=0;i<m;++i,w=w*wn){
x=q[i],y=w*q[i+m];
q[i]=x+y,q[i+m]=x-y;
}
return;
}
code
#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
using namespace std;
namespace Q{
il int rd(){
ri int x=0;ri bool f=0;ri char c=getchar();
while(!isdigit(c)) f|=(c==45),c=getchar();
while(isdigit(c)) x=x*10+(c^48),c=getchar();
return f?-x:x;
}
il void wt(int x){
if(x<0) x=-x,putchar(45);
if(x>=10) wt(x/10);
return putchar(x%10|48),void();
}
} using namespace Q;
cs double Pi=3.141592653589793;
cs int N=2097152;
namespace FastFastTle{
struct qwq{
double r,i;
qwq operator +(cs qwq o)cs{
return (qwq){r+o.r,i+o.i};
}
qwq operator -(cs qwq o)cs{
return (qwq){r-o.r,i-o.i};
}
qwq operator *(cs qwq o)cs{
return (qwq){r*o.r-i*o.i,o.r*i+r*o.i};
}
}wn,w,f[N],g[N],o[N],x,y;
il void FFT(qwq *q,int lim,int k){
if(lim==1) return;
for(ri int i=0;i<lim;++i) o[i]=q[i];
ri int m=lim>>1;
for(ri int i=0;i<lim;i++){
if(i&1) q[m+(i>>1)]=o[i];
else q[i>>1]=o[i];
}
FFT(q,m,k),FFT(q+m,m,k);
wn={cos(Pi*2/lim),k*sin(Pi*2/lim)},w={1,0};
for(ri int i=0;i<m;++i,w=w*wn){
x=q[i],y=w*q[i+m];
q[i]=x+y,q[i+m]=x-y;
}
return;
}
} using namespace FastFastTle;
signed main(){
int n=rd(),m=rd(),lim=1;
while(lim<=n+m) lim<<=1;
for(ri int i=0;i<=n;++i) f[i].r=rd();
for(ri int i=0;i<=m;++i) g[i].r=rd();
FFT(f,lim,1),FFT(g,lim,1);
for(ri int i=0;i<lim;++i) f[i]=f[i]*g[i];
FFT(f,lim,-1);
for(ri int i=0,s;i<=n+m;++i){
wt(f[i].r/lim+0.5),putchar(32);
}
return 0;
}
迭代实现
il void FFT(qwq q[],int lim,int k){
for(ri int i=0;i<lim;++i){
if(i<r[i]) swap(q[i],q[r[i]]);
}
for(ri int mid=1;mid<lim;mid<<=1){
wn={cos(Pi/mid),k*sin(Pi/mid)};
for(ri int R=mid<<1,j=0;j<lim;j+=R){
w={1,0};
for(ri int t=0;t<mid;++t,w=w*wn){
x=q[j+t],y=w*q[j+mid+t];
q[j+t]=x+y,q[j+mid+t]=x-y;
}
}
}
return;
}
code
#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
using namespace std;
namespace Q{
il int rd(){
ri int x=0;ri bool f=0;ri char c=getchar();
while(!isdigit(c)) f|=(c==45),c=getchar();
while(isdigit(c)) x=x*10+(c^48),c=getchar();
return f?-x:x;
}
il void wt(int x){
if(x<0) x=-x,putchar(45);
if(x>=10) wt(x/10);
return putchar(x%10|48),void();
}
} using namespace Q;
cs double Pi=3.141592653589793;
cs int N=2097152;
namespace FastFastTle{
struct qwq{
double r,i;
qwq operator +(cs qwq o)cs{
return (qwq){r+o.r,i+o.i};
}
qwq operator -(cs qwq o)cs{
return (qwq){r-o.r,i-o.i};
}
qwq operator *(cs qwq o)cs{
return (qwq){r*o.r-i*o.i,o.r*i+r*o.i};
}
}wn,w,f[N],g[N],x,y;
int l,r[N];
il void FFT(qwq q[],int lim,int k){
for(ri int i=0;i<lim;++i){
if(i<r[i]) swap(q[i],q[r[i]]);
}
for(ri int mid=1;mid<lim;mid<<=1){
wn={cos(Pi/mid),k*sin(Pi/mid)};
for(ri int R=mid<<1,j=0;j<lim;j+=R){
w={1,0};
for(ri int t=0;t<mid;++t,w=w*wn){
x=q[j+t],y=w*q[j+mid+t];
q[j+t]=x+y,q[j+mid+t]=x-y;
}
}
}
return;
}
} using namespace FastFastTle;
signed main(){
int n=rd(),m=rd(),lim=1;
while(lim<=n+m) lim<<=1,l++;
for(ri int i=0;i<=n;++i) f[i].r=rd();
for(ri int i=0;i<=m;++i) g[i].r=rd();
for(ri int i=0;i<lim;++i){
r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
FFT(f,lim,1),FFT(g,lim,1);
for(ri int i=0;i<lim;++i) f[i]=f[i]*g[i];
FFT(f,lim,-1);
for(ri int i=0,s;i<=n+m;++i){
wt(f[i].r/lim+0.5),putchar(32);
}
return 0;
}
