目录
「学习笔记」FFT 快速傅里叶变换
啥是 FFT 呀？它可以干什么？
必备芝士
点值表示
什么是点值表示
点值表示的乘法
复数
复数的定义
复数的乘法运算
单位复数
傅立叶正变换
傅里叶逆变换
FFT 的代码实现
还会有的 NTT 和三模数 NTT...
「学习笔记」FFT 快速傅里叶变换
几个星期之后，继 扩展欧拉定理 之后， \(lj\) 大佬又给我们来了一发数论...
虽然听得心态爆炸， 但是还好的是没有 \(ymx\) 大佬的飞机开得好...
至少我还没有坐飞机...
啥是 FFT 呀？它可以干什么？
首先，你需要知道 矩阵乘法 的相关知识。
通过 矩阵乘法 的知识，我们知道，对于一个 \(f(x)\) 与 \(g(x)\) ，如果我们用朴素的矩阵乘法计算 \(f(x)\cdot g(x)\) ，时间复杂度是 \(O(N^2M)\) ，但是对于这是函数，复杂度中的 \(M=1\) ，那么时间复杂度为 \(O(N^2)\) 。
似乎这是比较优秀的，但是看看这道 板题 ：
luoguOJ:P3803【模板】多项式乘法（FFT）
嗯？ \(n,m\le 10^6\) ，似乎过于太小了？（不要管这病句...）
而 FFT 就是用来解决多项式乘法的问题的，可以把它的时间复杂度优化到 \(O(nlog^n_2)\) 。
这里引用一句话：
实际上，FFT 并不是直接计算多项式乘法，而是把原来的多项式 \(f(x),g(x)\) 在 \(O(nlog^n_2)\) 的复杂度内转换为它的 点值表示（后面会讲），而点值表示的多项式相乘的时间复杂度是 \(O(n)\) 的。最后再用 \(O(nlog_2^n)\) 的时间复杂度把所得多项式的点值表示转化为一般形式。
必备芝士
在学习 FFT 之前，我们需要知道 点值表示 和 复数 。
点值表示
什么是点值表示
对于一个函数 \(f(x)\) ，我们可以用一个多项式表示它：
\[f(x)=a_nx^n+a_{n-1}x^{n-1}+...+a_{1}x+a_0
\]
这种表达方法我们叫做 一般形式 。
如果熟悉，那么可以确定一个东西，如果我们有在这个函数上的 \(n+1\) 个点的具体的坐标 \((x_i,y_i)\) ，那么我们就可以唯一确定这个函数的解析式，那么，我们也可以用这 \(n+1\) 个点的具体坐标来表示这个多项式。
当然，如果知道 拉格朗日插值法 的大佬，或许可以更好地理解这句话。
也就是说一个多项式与一个点值表示是一一对应的。
那么 FFT 完成的操作就是：
把已知的一个多项式转化成对应的点值表示；
把已知的点值表示转换成对应的多项式。复杂度都是 \(O(nlog_2^n)\) 。
点值表示的乘法
那么，假设我们用 \(\{x_0,x_1,x_2,...,x_n\}\) 来表示 \(f(x)\) 为 \(\{(x_0,f(x_0)),(x_1,f(x_1)),...,(x_n,f(x_n))\}\) ，表示 \(g(x)\) 为 \(\{(x_0,g(x_0)),(x_1,g(x_1)),...,(x_n,g(x_n))\}\) 。
那么我们就知道 \(f(x)\cdot g(x)\) 的点值表达式为 \(\{(x_0,f(x_0)\cdot g(x_0)),...,(x_n,f(x_n)\cdot g(x_n))\}\) ，其实就是对应项相乘，那么这个的时间复杂度为 \(O(n)\) ，似乎可以接受。
复数
复数的定义
\(Hint.\) 是 复数(complex) 而不是 负数(negative) 。
首先我们要知道 虚数 ，即有一个很特别的数 \(i\) ，满足 \(i^2=-1\) ，即 \(i=\sqrt{-1}\) 。
不用理解为什么，因为他本来就是特殊定义的，只需要记住就好。
然后一个复数 \(a\) 可以表示为实部 \(x\) 和虚部 \(y\times i\) 之和—— \(a=x+y\times i\) 。对于一个实数，我们可以把它放在“一维数轴”上，即数轴；那么对于一个复数，我们需要把它放在“二维数轴”，也就是直角坐标系上。
其实，我们可以把 \(a\) 表示成 向量(vector) 的形式，即将 \(a=x+y\times i\) 表示为向量 \((x,y)\) ，但是对于后面虚数的乘法，其实是不满足向量乘法的，但是这可以作为一个理解方式。
复数的乘法运算
对于两个虚数 \(a,b\) ，其中 \(a=x+y\times i,b=x'+y'\times i\) ，计算 \(ab\) 。
我们有两种方法：
- 可以直接展开相乘： \(ab=(x+yi)(x'+y'i)=(xx'-yy')+(xy'+x'y)i\) 。
可以从几何角度理解：两个复数相乘等于“幅角相加，模长相乘”（但并不满足向量的运算）
更推荐就用方法一来理解。
单位复数
\(Hint.\) 复数的模长 \(|a|\) （这个不是绝对值），表示把它表示成向量后的长度，即 \(|a|=\sqrt{x^2+y^2}\) 。
我们定义 \(n\) 次单位复根 \(\omega_n^i\) 为满足 \(\text{z}^n=1\) 的复数 \(\text{z}\) 。
本人的理解方法：单位复根 \(\omega_n^i\) 即为将一个单位圆分成 \(n\) 分，而 \(\omega_n^i\) 就是 \(i\) 个这样的角拼起来。特别地， \(\omega_n^0=1\) 。
其实，之后的理解都可以通过 三角函数 有更好的理解
可以发现 \(n\) 次单位复根满足以下性质（ \(\omega_n^i\) 表示第 \(i\) 个单位复根）：
\(|\omega_n^i|=1\) （都在单位圆上）；
\(\omega_n^i\) 的幅角为 \(\frac{2\pi i}{n}\) ；\(Tab.\) 这里的 \(i\) 是下标不是虚数，角度为弧度制（ \(360^{\circ}\) 弧度角为 \(2\pi\) ）
总共有 \(n\) 个单位复根，记为 \(\omega_n^0,\omega_n^1,\omega_n^2,...,\omega_n^{n-1}\)
另外，还有一些比较明显的性质：
\(\omega_n^i\cdot\omega_n^j=\omega_n^{i+j}\) ；
\(\omega_n^{i+n}=\omega_n^i\) ；
\(\omega_n^i=\omega_{kn}^{ki}\) ；
\(\omega_n^i=-\omega_n^{i+n/2}\) （这里的“负”表示大小相等、方向相反）
以上这些结论都可以通过在单位圆上画出单位根来证明。
单位复根有什么用呢？因为 \(n\) 次单位复根恰好有 \(n\) 个，就可以把这 \(n\) 个单位复根分别代入 \(n-1\) 次多项式 \(f(x)\) —— 得到点值表达 \(\{(\omega_n^0,f(\omega_n^0)),(\omega_n^1,f(\omega_n^1)),...,(\omega_n^{n-1},f(\omega_n^{n-1}))\}\) 。
傅立叶正变换
有了这些辅助知识 没错，虽然你可能已经晕了，但是他们真的只是辅助知识 ，我们终于可以进行正题了。
所谓变换，那么一定有正也有逆，现在我我们先来掌握它的正变换。
FFT 的正变换实现，是基于对多项式进行奇偶项分开递归再合并的分治进行的。
对于 n-1 次多项式，我们选择插入 n 次单位根求出其点值表达式。
记多项式 \(A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\) 。
再记 \(A_o(x)=a_1+a_3x+a_5x^2+...\) 。
再记 \(A_e(x)=a_0+a_2x+a_4x^2+...\) 。
有 \(A(x)=x*A_o(x^2)+A_e(x^2)\) 。
令 \(n = 2*p\) 。则有：
\(A(w_n^k)=w_n^k*A_o[(w_{n/2}^{k/2})^2]+A_e[(w_{n/2}^{k/2})^2]=w_n^k*A_o(w_p^k)+A_e(w_p^k)\) ；
\(A(w_n^{k+p})=w_n^{k+p}*A_o(w_p^{k+p})+A_e(w_p^{k+p})=-w_n^k*A_o(w_p^k)+A_e(w_p^k)\)
在已知 \(A_o(w_p^k)\) 与 \(A_e(w_p^k)\) 的前提下，可以 \(O(1)\) 算出 \(A(w_n^k)\) 与 \(A(w_n^{k+p})\) 。
因此，假如我们递归求解 \(A_o(x),A_e(x)\) 两个多项式 \(p\) 次单位根的插值，就可以在 \(O(n)\) 的时间内算出 \(A(x)\) n 次单位根的插值。
时间复杂度是经典的 \(T(n)=2*T(n/2)+O(n)=O(n\log n)\)。
傅里叶逆变换
刚刚研究完正的，现在我们来研究逆变换，其实也比较好理解。
观察我们刚刚的插值过程，实际上就是进行了如下的矩阵乘法。
\[\begin{bmatrix} (w_n^0)^0 & (w_n^0)^1 & \cdots & (w_n^0)^{n-1} \\ (w_n^1)^0 & (w_n^1)^1 & \cdots & (w_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (w_n^{n-1})^0 & (w_n^{n-1})^1 & \cdots & (w_n^{n-1})^{n-1} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \begin{bmatrix} A(w_n^0) \\ A(w_n^1) \\ \vdots \\ A(w_n^{n-1}) \end{bmatrix}
\]
我们记上面的系数矩阵为 \(V\) 。
对于下面定义的 \(D\) ：
\[D = \begin{bmatrix} (w_n^{-0})^0 & (w_n^{-0})^1 & \cdots & (w_n^{-0})^{n-1} \\ (w_n^{-1})^0 & (w_n^{-1})^1 & \cdots & (w_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (w_n^{-(n-1)})^0 & (w_n^{-(n-1)})^1 & \cdots & (w_n^{-(n-1)})^{n-1} \end{bmatrix}
\]
考虑 \(D*V\) 的结果：
\[(D*V)_{ij}
=\sum_{k=0}^{k<n}d_{ik}*v_{kj}
=\sum_{k=0}^{k<n}w_n^{-ik}*w_{n}^{kj}
=\sum_{k=0}^{k<n}w_n^{k(j-i)}
\]
当 \(i = j\) 时， \((D*V)_{ij}=n\) ；
当 \(i ≠ j\) 时， \((D*V)_{ij}=1+w_n^{j-i}+(w_n^{j-i})^2+...=\frac{1-(w_n^{j-i})^n}{1-w_n^{j-i}}=0\) ；
根据定义， \(n\) 次单位根的 \(n\) 次方都等于 \(1\)
所以：\(\frac1n*D=V^{-1}\)
因此将这个结果代入最上面那个公式里面，有：
\[\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \frac1n \begin{bmatrix} (w_n^{-0})^0 & (w_n^{-0})^1 & \cdots & (w_n^{-0})^{n-1} \\ (w_n^{-1})^0 & (w_n^{-1})^1 & \cdots & (w_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (w_n^{-(n-1)})^0 & (w_n^{-(n-1)})^1 & \cdots & (w_n^{-(n-1)})^{n-1} \end{bmatrix}\begin{bmatrix} A(w_n^0) \\ A(w_n^1) \\ \vdots \\ A(w_n^{n-1}) \end{bmatrix}
\]
“这样，逆变换 就相当于把 正变换 过程中的 \(w_n^k\) 换成 \(w_n^{-k}\)，之后结果除以 \(n\) 就可以了。”——摘自某博客。
……
还是有点难理解。比如为什么我们不直接把 \(w_n^k\) 换成 \(\frac1n*w_n^{-k}\) 算了。
实际上，因为 \(w_n^{-k}=w_n^{n-k}\) ，也就是说它 TM 还是一个 \(n\) 次单位根。所以我们插值还是正常的该怎么插怎么插。如果换成 \(\frac1n*w_n^{-k}\) 它就不是一个单位根，以上性质就不满足了。
FFT 的代码实现
我们有两个版本——递归、迭代，相信大家都也想到了吧？
毋庸置疑的，递归版本确实很好写，将 \(FFT\) 的思路全部放到上面即可，但是方便也有他的坏处——递归实现常数较大，这样对于一些题就会 \(T\) 飞，但是我还是写了，这里上代码了
#include<cstdio>
#include<cmath>
#define rep(i,__l,__r) for(register int i=__l,i##_end_=__r;i<=i##_end_;++i)
#define fep(i,__l,__r) for(register int i=__l,i##_end_=__r;i>=i##_end_;--i)
#define writc(a,b) fwrit(a),putchar(b)
#define mp(a,b) make_pair(a,b)
#define ft first
#define sd second
#define LL long long
#define ull unsigned long long
#define pii pair<int,int>
#define Endl putchar('\n')
// #define FILEOI
// #define int long long
#ifdef FILEOI
#define MAXBUFFERSIZE 500000
inline char fgetc(){
static char buf[MAXBUFFERSIZE+5],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,MAXBUFFERSIZE,stdin),p1==p2)?EOF:*p1++;
}
#undef MAXBUFFERSIZE
#define cg (c=fgetc())
#else
#define cg (c=getchar())
#endif
template<class T>inline void qread(T& x){
char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
if(f)x=-x;
}
inline int qread(){
int x=0;char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
return f?-x:x;
}
template<class T,class... Args>inline void qread(T& x,Args&... args){qread(x),qread(args...);}
template<class T>inline T Max(const T x,const T y){return x>y?x:y;}
template<class T>inline T Min(const T x,const T y){return x<y?x:y;}
template<class T>inline T fab(const T x){return x>0?x:-x;}
inline int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
inline void getInv(int inv[],const int lim,const int MOD){
inv[0]=inv[1]=1;for(int i=2;i<=lim;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
}
template<class T>void fwrit(const T x){
if(x<0)return (void)(putchar('-'),fwrit(-x));
if(x>9)fwrit(x/10);putchar(x%10^48);
}
inline LL mulMod(const LL a,const LL b,const LL mod){//long long multiplie_mod
return ((a*b-(LL)((long double)a/mod*b+1e-8)*mod)%mod+mod)%mod;
}
const int MAXN=3e6;
const double Pi=acos(-1.0);
class task{
private:
struct cplx{
double vr,vi;//实部和虚部
cplx(const double R=0,const double I=0):vr(R),vi(I){}//构造函数
//------------------overload----------------//
cplx operator + (const cplx a)const{return cplx(vr+a.vr,vi+a.vi);}//重载加法
cplx operator - (const cplx a)const{return cplx(vr-a.vr,vi-a.vi);}
cplx operator * (const cplx a)const{return cplx(vr*a.vr-vi*a.vi,vr*a.vi+a.vr*vi);}
cplx operator / (const double var)const{return cplx(vr/var,vi/var);}
};
int n,m;
cplx a[MAXN+5],b[MAXN+5];
void fft(cplx* f,const int len,const short opt=1){
//opt==-1 : FFT 的逆变换
if(!len)return;
cplx f0[len+5],f1[len+5];
for(int i=0;i<len;++i)
f0[i]=f[i<<1],f1[i]=f[i<<1|1];
fft(f0,len>>1,opt);
fft(f1,len>>1,opt);
cplx w=cplx(cos(Pi/len),opt*sin(Pi/len)),buf=cplx(1,0);
for(int i=0;i<len;++i,buf=buf*w){
f[i]=f0[i]+buf*f1[i];
f[i+len]=f0[i]-buf*f1[i];
}
}
public:
inline void launch(){
qread(n,m);
rep(i,0,n)scanf("%lf",&a[i].vr);
rep(i,0,m)scanf("%lf",&b[i].vr);
for(m+=n,n=1;n<=m;n<<=1);
fft(a,n>>1);
fft(b,n>>1);
rep(i,0,n-1)a[i]=a[i]*b[i];
fft(a,n>>1,-1);
rep(i,0,m)writc((int)((a[i].vr)/n+0.5),' ');
Endl;
}
}This;
signed main(){
#ifdef FILEOI
freopen("file.in","r",stdin);
freopen("file.out","w",stdout);
#endif
This.launch();
return 0;
}
似乎递归版本比较好写，现在我们来看一下迭代（递推）版本应该怎么做：
原序列： \(0,1,2,3,4,5,6,7\)
终序列： \(0,4,2,6,1,5,3,7\)
转换为二进制再来看看。
原序列： \(000,001,010,011,100,101,110,111\)
终序列： \(000,100,010,110,001,101,011,111\)
可以发现终序列是原序列每个元素的二进制翻转。
于是我们可以先把要变换的系数排在相邻位置，从下往上迭代。
这个二进制翻转过程可以自己脑补方法，只要保证时间复杂度 \(O(n\log n)\) ，代码简洁就可以了。
在这里给出一个参考的方法：
我们对于每个 \(i\) ，假设已知 \(i-1\) 的翻转为 \(j\) 。考虑不进行翻转的二进制加法怎么进行：从最低位开始，找到第一个为 \(0\) 的二进制位，将它之前的 \(1\) 变为 \(0\) ，将它自己变为 \(1\) 。因此我们可以从 \(j\) 的最高位开始，倒过来进行这个过程。
这是迭代版本：
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
#define rep(i,__l,__r) for(register int i=__l,i##_end_=__r;i<=i##_end_;++i)
#define fep(i,__l,__r) for(register int i=__l,i##_end_=__r;i>=i##_end_;--i)
#define writc(a,b) fwrit(a),putchar(b)
#define mp(a,b) make_pair(a,b)
#define ft first
#define sd second
#define LL long long
#define ull unsigned long long
#define pii pair<int,int>
#define Endl putchar('\n')
// #define FILEOI
// #define int long long
#ifdef FILEOI
#define MAXBUFFERSIZE 500000
inline char fgetc(){
static char buf[MAXBUFFERSIZE+5],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,MAXBUFFERSIZE,stdin),p1==p2)?EOF:*p1++;
}
#undef MAXBUFFERSIZE
#define cg (c=fgetc())
#else
#define cg (c=getchar())
#endif
template<class T>inline void qread(T& x){
char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
if(f)x=-x;
}
inline int qread(){
int x=0;char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
return f?-x:x;
}
template<class T,class... Args>inline void qread(T& x,Args&... args){qread(x),qread(args...);}
template<class T>inline T Max(const T x,const T y){return x>y?x:y;}
template<class T>inline T Min(const T x,const T y){return x<y?x:y;}
template<class T>inline T fab(const T x){return x>0?x:-x;}
inline int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
inline void getInv(int inv[],const int lim,const int MOD){
inv[0]=inv[1]=1;for(int i=2;i<=lim;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
}
template<class T>void fwrit(const T x){
if(x<0)return (void)(putchar('-'),fwrit(-x));
if(x>9)fwrit(x/10);putchar(x%10^48);
}
inline LL mulMod(const LL a,const LL b,const LL mod){//long long multiplie_mod
return ((a*b-(LL)((long double)a/mod*b+1e-8)*mod)%mod+mod)%mod;
}
const int MAXN=3e6;
const double Pi=acos(-1.0);
class task{
private:
struct cplx{
double vr,vi;//实部和虚部
cplx(const double R=0,const double I=0):vr(R),vi(I){}//构造函数
//------------------overload----------------//
cplx operator + (const cplx a)const{return cplx(vr+a.vr,vi+a.vi);}//重载加法
cplx operator - (const cplx a)const{return cplx(vr-a.vr,vi-a.vi);}
cplx operator * (const cplx a)const{return cplx(vr*a.vr-vi*a.vi,vr*a.vi+a.vr*vi);}
cplx operator / (const double var)const{return cplx(vr/var,vi/var);}
};
int n,m;
cplx a[MAXN+5],b[MAXN+5];
int revi[MAXN+5];
/*
f(w^x)
=f0(w^{2x})+w^x*f1(w^{2x})
=a0+a2+a4+...,a1+a3+a5+...
=a0+a4+a8+...+a2+a6+a10...+a1+a5+a9+...+a3+a7+a11...
=f00(w^{4x})+w^{2x}*f01(w^{4x})+w^x*f10(w^{4x})+w^3x*f11(w^{4x})
=a0+a4+a8+...,w^{2x}*(a2+a6+a10...),w^x*(a1+a5+a9+...),w^{3x}*(a3+a7+a11...)
f_s -> 下标将 s 反过来之后, a_i 的 i 的二进制反过来与 s 相同
.
.
.
s
000 001 010 011
|反过来
v
000 100 010 110
a0 a{k/2} a{3*k/4} a{k}
a0 a4 a2 a6
*/
/*
void fft(cplx* f,const int len,const short opt=1){
if(!len)return;
cplx f0[len+5],f1[len+5];
for(int i=0;i<len;++i)
f0[i]=f[i<<1],f1[i]=f[i<<1|1];
fft(f0,len>>1,opt);
fft(f1,len>>1,opt);
cplx w=cplx(cos(Pi/len),opt*sin(Pi/len)),buf=cplx(1,0);
for(int i=0;i<len;++i,buf=buf*w){
f[i]=f0[i]+buf*f1[i];
f[i+len]=f0[i]-buf*f1[i];
}
}
*/
inline void fft(cplx* f,const short opt=1){
for(int i=0;i<n;++i)if(i<revi[i])
swap(f[i],f[revi[i]]);
for(int p=2;p<=n;p<<=1){
//枚举层数
int len=p/2;//上一层的一半的长度
cplx tmp(cos(Pi/len),opt*sin(Pi/len));
//单位复根
for(int k=0;k<n;k+=p){
cplx buf(1,0);//记录 omega 的次方
for(int l=k;l<k+len;++l,buf=buf*tmp){
//每次 buf 累成单位 omega
cplx tt=buf*f[len+l];
f[len+l]=f[l]-tt;
f[l]=f[l]+tt;
/*
此处与递归版本的 FFT 中这一段是一样的:
f[i]=f0[i]+buf*f1[i];
f[i+len]=f0[i]-buf*f1[i];
*/
}
}
}
if(opt==-1)for(int i=0;i<n;++i)f[i]=f[i]/n;
//如果是 逆变换 , 那么需要全部 /n
}
public:
inline void launch(){
qread(n,m);
rep(i,0,n)scanf("%lf",&a[i].vr);
rep(i,0,m)scanf("%lf",&b[i].vr);
for(m+=n,n=1;n<=m;n<<=1);
rep(i,0,n-1)revi[i]=(revi[i>>1]>>1)|((i&1)?n>>1:0);//处理反转
fft(a);fft(b);
rep(i,0,n-1)a[i]=a[i]*b[i];
fft(a,-1);
rep(i,0,m)writc((int)(a[i].vr+0.5),' ');
Endl;
}
}This;
signed main(){
#ifdef FILEOI
freopen("file.in","r",stdin);
freopen("file.out","w",stdout);
#endif
This.launch();
return 0;
}
如果我可以回到过去，我一定会去当杀手。
为什么？因为我要去干翻欧某和傅某某...
