FFT学习笔记

本人蒟蒻，只能介绍FFT在OI中的应用，如有错误或不当之处还请指出。

首先先说一下那一堆什么什么 $TT$ 的都是什么

DFT：离散傅里叶变换
FFT ：快速傅里叶变换 用于求多项式乘法 $O(nlog(n))$
FNTT/NTT ：FTT的优化，常数及精度更优
FWT：快速沃尔什变换
MTT：任意模数NTT
FMT 快速莫比乌斯变换
SFT 稀疏傅里叶变换 百度百科说比FFT快10到100倍 但是在OI里真的有应用吗？

后面三个我都不会，这里也不介绍 NTT我也不会，仙姑了~

以下是一些前置知识（当然不只是前置知识🌚）

多项式

多项式有系数表示法和点值表示法两种，我们常见的 $1+x+x^2+x^3$ 这样的就叫做系数表示法，而用若干个点表示多项式的就叫做点值表示法。

把多项式看作函数，如果你学过一点插值的话，就可以知道给定坐标系中横坐标互不相同的 $n+1$ 个点，可以唯一确定一条 $n$ 次函数经过这 $n+1$ 个点。

用系数表示法来计算两个多项式之积，枚举各项系数相乘，时间复杂度是 $O(n^2)$ 的。但是点值表示法呢？我们只需要把在某一点处两个多项式的值相乘，就能得到结果多项式在该处的值，直接枚举 $n$ 个点即可，时间复杂度 $O(n)$ 。

那么我们就可以考虑将两个多项式全部转化为点值表示，然后计算乘积，再将结果多项式从点值表示转化为系数表示，这也就是FFT的大致流程。

那么点值表示和系数表示相互转化的时间复杂度呢？发现是 $O(n^2)$ 的，《这么一整直接把常数优化大了🌚》

那么我们就要引出快速傅里叶变换了，在这之前我们引出单位根这个概念。

单位根

如果有没学过复数运算的，请先去学复数运算。

首先对于方程 $x^n=1$ ，其在复数域中有 $n$ 个解，并且如果画出来的话，这 $n$ 个解恰好等分单位圆，我们按照逆时针顺序分别叫他们 ${\omega }_n^1、{\omega }_n^2、\cdots 、{\omega }_n^n$ ，不难发现其性质：

$${\omega }_{2n}^{2k}={\omega }_n^k$$

$${\omega }_{2n}^{k+n}=-{\omega }_{2n}^k$$

$$({\omega }_n^i{)}^j=({\omega }_n^j{)}^i$$

$$\sum _{i=0}^{n-1}({\omega }_n^k{)}^i=n\times [k==n]$$

前三条性质显然，我证明一下第四条：

首先看这个级数：$\sum _{i=0}^n{x}^i$ ，我们令 $f=\sum _{i=0}^n{x}^i$ ，那么 $xf=\sum _{i=1}^n{x}^i$ ，相减可以求出 $f=\frac{x^{n+1}-1}{x-1}$ ，当然在分母为 $0$ 时不能用该式求解。我们发现上式就是一个等比数列的形式，所以：

$$\begin{array}{rl}\sum _{i=0}^{n-1}({\omega }_n^k{)}^i& =\frac{({\omega }_n^k{)}^n}{{\omega }_n^k-1}=\frac{({\omega }_n^n{)}^k-1}{{\omega }_n^k-1}=\frac{1-1}{{\omega }_n^k-1}\end{array}$$

这里分子一定为零了，我们讨论分母。发现当且仅当 $k=n$ 时，分母为零，此时直接代入原式得到答案为 $n$ 。其余情况该式均为零。

快速傅里叶变换

那么介绍这个单位根干什么呢？我们可以取一个 $n-1$ 次多项式在这 $n$ 个单位根处的取值来作为这个多项式的点值表示，至于为什么这么取，看下面就知道了。

我们定义一个 $n-1$ 次多项式 $A(x)=\sum _{i=0}^{n-1}{a}_i{x}^i$ ，其各项系数为 $a_0,a_1,\cdots ,a_{n-1}$ 。

我们按照各项系数下标的奇偶性把系数分成两类，并各作为一个 $\frac{n}{2}-1$ 次多项式的系数，从而得到下面两个新的多项式：

$$A_1(x)=\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{x}^i$$

$$A_2(x)=\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{x}^i$$

并且有：$A(x)=A_1(x^2)+x{A}_2(x^2)$

再看一下我们要求什么，我们要求：$A({\omega }_n^1),A({\omega }_n^2),\cdots ,A({\omega }_n^n)$ 。
分别代入 $A({\omega }_n^k)$ 和 $A({\omega }_n^{k+\frac{n}{2}})$ ，我们可以得到以下两个式子：

$$\begin{array}{rl}A({\omega }_n^k)& =A_1(({\omega }_n^k{)}^2)+{\omega }_n^k{A}_2(({\omega }_n^k{)}^2)\\ & =A_1({\omega }_n^{2k})+{\omega }_n^k{A}_2({\omega }_n^{2k})\\ & =A_1({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k{A}_2({\omega }_{\frac{n}{2}}^k)\end{array}$$

$$\begin{array}{rl}A({\omega }_n^{k+\frac{n}{2}})& =A_1(({\omega }_n^{k+\frac{n}{2}}{)}^2)+{\omega }_n^{k+\frac{n}{2}}{A}_2(({\omega }_n^{k+\frac{n}{2}}{)}^2)\\ & =A_1({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}{A}_2({\omega }_n^{2k+n})\\ & =A_1({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k{A}_2({\omega }_{\frac{n}{2}}^k)\end{array}$$

发现这两个式子惊人的相似，那我们只需要求出 $A({\omega }_n^1),A({\omega }_n^2),\cdots ,A({\omega }_n^{\frac{n}{2}})$ ，然后就可以 $O(1)$ 地求出 $A({\omega }_n^{\frac{n}{2}+1}),A({\omega }_n^{\frac{n}{2}+2}),\cdots ,A({\omega }_n^n)$ 。

我们考虑等式右边的东西，它的形式和原来那个问题一样，于是我们就可以递归求解了。但是这里的除法均不是向下取整，所以在实现的时候我们需要把 $n$ 补成一个 $2$ 的整数次幂，也就是把原先的多项式都补成一个 $2^k$ 次的多项式，至于多出来的系数，直接置为零就好了，也不会影响答案。不难发现，这样的算法时间复杂度是 $O(nlog(n))$ 的。

放个码

struct complex{
    double x,y;
    complex operator+(const complex op)const{
        return complex{x+op.x,y+op.y};
    } 
    complex operator-(const complex op)const{
        return complex{x-op.x,y-op.y};
    }
    complex operator*(const complex op)const{
        return complex{x*op.x-y*op.y,x*op.y+y*op.x};
    }
};//当然你也可以直接用std里的complex
void FFT(complex*a,int n,int op){
    if(!n)return;
    complex a1[n],a2[n];
    for(int i=0;i<n;++i)a1[i]=a[i<<1],a2[i]=a[i<<1|1];
    FFT(a1,n>>1,op);FFT(a2,n>>1,op);
    complex W{cos(pi/n),sin(pi/n)*op},w{1,0};
    for(int i=0;i<n;++i,w=w*W)a[i]=a1[i]+w*a2[i],a[i+n]=a1[i]-w*a2[i];
}

于是我们有了 $O(nlog(n))$ 将一个系数表示的多项式转为点值表示的方法，然后呢？不要忘了我们还要将点值表示的结果多项式再转为系数表示，所以下面我们介绍快速傅里叶逆变换。

快速傅里叶逆变换

这里有线性代数的推法，如果你是会范德蒙德矩阵巨佬的话，就不要看这里的介绍了。🌚

我们再定义一个多项式 $B=\sum _{i=0}^{n-1}{y}_i{x}^i$ ，这里的 $y_i$ 是 $A$ 在 ${\omega }_n^i$ 处的值，即 $A({\omega }_n^i)$ （说了句废话🌚）。

上面的 ${\omega }_n^0$ 就等于 ${\omega }_n^n$ ，怕有人没想过来还是写一下🌚。

这是系数表示 $B$ ，我们考虑用点值表示法来表示 $B$ ，这次我们取 $b_i={\omega }_n^{-i}$ ，那为什么要这么取呢？我也不知道 可以去学一下用范德蒙德矩阵推导，然后就知道为什么是这个了🌚。

我们先试着把 $b_k$ 代入 $B$ 中：

$$\begin{array}{rl}B(b_k)& =\sum _{i=0}^{n-1}A({\omega }_n^i)({\omega }_n^{-k}{)}^i\\ & =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j({\omega }_n^i{)}^j({\omega }_n^{-k}{)}^i\\ & =\sum _{j=0}^{n-1}{a}_j\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i\end{array}$$

后面的部分是个等比数列，直接用上文中单位根的第四条性质可以得到只有当 $j=k$ 时，该部分值为 $n$ ，否则为零，不难得出：

$$B(b_k)=B({\omega }_n^{-k})=n{a}_k$$

发现我们只需要得到 $B$ 在 ${\omega }_n^{-1},{\omega }_n^{-2},\cdots ,{\omega }_n^{-n},$ 处的取值，就可以推出结果多项式 $A$ 的各项系数，也就有了我们想要的东西。

如果你有过人的观察能力，就发现这步和上一步相比多了个负号，那么我们是否还可以用之前的方法去做呢？当然是可以的，这里我就不再赘述了，自己推导即可。实现的时候多传一个参数区分开变换和逆变换就行了。

于是我们就得到了递归版的FFT：

点击查看代码

#include<bits/stdc++.h>
const int N=4e6+10;
const double pi=acos(-1);
struct complex{
    double x,y;
    complex operator+(const complex op)const{
        return complex{x+op.x,y+op.y};
    } 
    complex operator-(const complex op)const{
        return complex{x-op.x,y-op.y};
    }
    complex operator*(const complex op)const{
        return complex{x*op.x-y*op.y,x*op.y+y*op.x};
    }
}a[N],b[N];
int n,m;
inline int read(){
    int ans=0;char ch=getchar_unlocked();bool fl=0;
    while(ch<'0'||ch>'9'){if(ch=='-')fl=1;ch=getchar_unlocked();}
    while(ch>='0'&&ch<='9')ans=(ans<<3)+(ans<<1)+(ch^48),ch=getchar_unlocked();
    return fl?~ans+1:ans;
}
inline void print(int x){(x<0)&&(putchar_unlocked('-'),x=-x);(x>9)&&(print(x/10),0);putchar_unlocked(x%10|48);}
void FFT(complex*a,int n,int op){
    if(!n)return;
    complex a0[n],a1[n];
    for(int i=0;i<n;++i)a0[i]=a[i<<1],a1[i]=a[i<<1|1];
    FFT(a0,n>>1,op);FFT(a1,n>>1,op);
    complex W{cos(pi/n),sin(pi/n)*op},w{1,0};
    for(int i=0;i<n;++i,w=w*W)a[i]=a0[i]+w*a1[i],a[i+n]=a0[i]-w*a1[i];
}
int main(){
    // freopen("1.in","r",stdin);
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;++i)a[i]={1.0*read(),0};
    for(int i=0;i<=m;++i)b[i]={1.0*read(),0};
    m=m+n;n=1;
    while(n<=m)n<<=1;
    FFT(a,n>>1,1);FFT(b,n>>1,1);
    for(int i=0;i<n;++i)a[i]=a[i]*b[i];
    FFT(a,n>>1,-1);
    for(int i=0;i<=m;++i)printf("%.0f ",fabs(a[i].x)/n);
}

这里应该是过不了板子题的，但是它意外过了，没看懂，但是还是要引出优化。《就当上面那个码过不了吧》

于是我们引出一些优化，注意到上面的代码中有：a[i]=a1[i]+w*a2[i],a[i+n]=a1[i]-w*a2[i]; ，其中 w*a2[i] 我们计算了两次，而复数相乘常数会很大，所以我们拿个变量记一下就好了。于是我们完成了所有的优化

《显然不是的🌚》，发现递归常数大，我们考虑是否可以摆脱递归，就有了下面这个优化：

位逆序置换

我们观察每次分治后下标的位置变化，这里以七次多项式举例：

$[a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7]$
$[a_0,a_2,a_4,a_6][a_1,a_3,a_5,a_7]$
$[a_0,a_4][a_2,a_6][a_1,a_5][a_3,a_7]$
$[a_0][a_4][a_2][a_6][a_1][a_5][a_3][a_7]$

再用二进制表示出开始和结尾的位置：

$000,001,010,011,100,101,110,111$
$000,100,010,110,001,101,011,111$

如果我们知道每个系数递归到最后的位置，那么我们就可以自下而上的递推求解，也就摆脱了递归。

那么我们如何求呢？显然我们可以一位一位地翻转求解，这样的时间复杂度是 $O(nlog(n))$ 的，当然我们还有更好的方法。

记原先下标为 $n$ 的系数最后的位置为 $f(n)$ ，我们考虑递推求解：

我们可以先将 $n$ 右移一位，翻过来再右移一位，最后再加上 $n$ 最低位翻转的结果，即可得到 $f(n)$ ，写成数学式子是：$f(n)=\lfloor \frac{f(\lfloor \frac{n}{2}\rfloor )}{2}\rfloor +(nmod2)\times \frac{len}{2}$ 。证明我也不知道

非递归版FFT

#include<bits/stdc++.h>
const int N=4e6+10;
const double pi=acos(-1);
struct complex{
    double x,y;
    complex operator+(const complex op)const{
        return complex{x+op.x,y+op.y};
    } 
    complex operator-(const complex op)const{
        return complex{x-op.x,y-op.y};
    }
    complex operator*(const complex op)const{
        return complex{x*op.x-y*op.y,x*op.y+y*op.x};
    }
}a[N],b[N],buf[N];
int n,m,pos[N];
inline int read(){
    int ans=0;char ch=getchar_unlocked();bool fl=0;
    while(ch<'0'||ch>'9'){if(ch=='-')fl=1;ch=getchar_unlocked();}
    while(ch>='0'&&ch<='9')ans=(ans<<3)+(ans<<1)+(ch^48),ch=getchar_unlocked();
    return fl?~ans+1:ans;
}
inline void print(int x){(x<0)&&(putchar_unlocked('-'),x=-x);(x>9)&&(print(x/10),0);putchar_unlocked(x%10|48);}
void FFT(complex a[],int flag){
    for(int i=0;i<n;i++)if(i<pos[i])std::swap(a[i],a[pos[i]]);
    for(int mid=1;mid<n;mid<<=1){
        complex W={cos(pi/mid),flag*sin(pi/mid)};
        for(int r=mid<<1,j=0;j<n;j+=r){
            complex w={1,0};
            for(int k=0;k<mid;k++,w=w*W){
                complex op=w*a[j+k+mid];
                buf[j+k]=a[j+k]+op;
                buf[j+k+mid]=a[j+k]-op;
            }
        }
        for(int i=0;i<n;i++)a[i]=buf[i];
    }
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;++i)a[i]={1.0*read(),0};
    for(int i=0;i<=m;++i)b[i]={1.0*read(),0};
    m=m+n;n=1;
    while(n<=m)n<<=1;
    for(int i=1;i<=n;i++)pos[i]=(pos[i>>1]>>1)|((i&1)*(n>>1));
    FFT(a,1);FFT(b,1);
    for(int i=0;i<n;++i)a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=0;i<=m;++i)printf("%.0f ",fabs(a[i].x)/n);
}

发现递归改成非递归，直接慢了0.49s，这是怎么烩逝呢？我也不知道，但理论上是快了🌚，大概是我写的常数太大了吧🌚。

但是无所谓，我们认为它起到了一定的优化作用，但是我们发现多了一个数组，那我们可不可以不要这个数组呢？于是就有了下面这个优化：

蝶形运算优化

原本的方法可能要动脑子，我这里就说一种不用动脑子的方法。

观察这里：

for(int r=mid<<1,j=0;j<n;j+=r){
    complex w={1,0};
    for(int k=0;k<mid;k++,w=w*W){
        complex op=w*a[j+k+mid];
        buf[j+k]=a[j+k]+op;
        buf[j+k+mid]=a[j+k]-op;
    }
}
for(int i=0;i<n;i++)a[i]=buf[i];

发现两次赋值比较慢，直接拿两个变量临时存一下，直赋值给原数组即可，这样就少了个数组，并且不会有影响。

就是这样

#include<bits/stdc++.h>
const int N=4e6+10;
const double pi=acos(-1);
struct complex{
    double x,y;
    complex operator+(const complex op)const{
        return complex{x+op.x,y+op.y};
    } 
    complex operator-(const complex op)const{
        return complex{x-op.x,y-op.y};
    }
    complex operator*(const complex op)const{
        return complex{x*op.x-y*op.y,x*op.y+y*op.x};
    }
}a[N],b[N];
int n,m,pos[N];
inline int read(){
    int ans=0;char ch=getchar_unlocked();bool fl=0;
    while(ch<'0'||ch>'9'){if(ch=='-')fl=1;ch=getchar_unlocked();}
    while(ch>='0'&&ch<='9')ans=(ans<<3)+(ans<<1)+(ch^48),ch=getchar_unlocked();
    return fl?~ans+1:ans;
}
inline void print(int x){(x<0)&&(putchar_unlocked('-'),x=-x);(x>9)&&(print(x/10),0);putchar_unlocked(x%10|48);}
void FFT(complex a[],int flag){
    for(int i=0;i<n;i++)if(i<pos[i])std::swap(a[i],a[pos[i]]);
    for(int mid=1;mid<n;mid<<=1){
        complex W={cos(pi/mid),flag*sin(pi/mid)};
        for(int r=mid<<1,j=0;j<n;j+=r){
            complex w={1,0};
            for(int k=0;k<mid;k++,w=w*W){
                complex op=w*a[j+k+mid],opl=a[j+k];
                a[j+k]=opl+op;
                a[j+k+mid]=opl-op;
            }
        }
    }
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;++i)a[i]={1.0*read(),0};
    for(int i=0;i<=m;++i)b[i]={1.0*read(),0};
    m=m+n;n=1;
    while(n<=m)n<<=1;
    for(int i=1;i<=n;i++)pos[i]=(pos[i>>1]>>1)|((i&1)*(n>>1));
    FFT(a,1);FFT(b,1);
    for(int i=0;i<n;++i)a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=0;i<=m;++i)printf("%.0f ",fabs(a[i].x)/n);
}

一个挺简单也挺好的优化，比着上一版非递归直接少了1.07s，想再了解点的移步oi-wiki，当然还有一个优化

三步变两步优化

顾名思义，之前的代码调用了三次FFT函数，我们实际上是可以只调用两次的。

对于给定的两个多项式 $A(x)=\sum _{i=0}^{n-1}{a}_i{x}^i,B(x)=\sum _{i=0}^{m-1}{b}_i{x}^i$ ，我们令 $a、b$ 分别为新的多项式各项的虚部和实部（我这里用 $\mathrm{i}$ 来表示虚数单位，注意辨别），从而得到多项式 $H(x)=\sum _{i=0}^{len-1}(a_i+b_i\mathrm{i})x^i$

可以将其拆成两项相乘：

$$\begin{array}{rl}H(x)=(\sum _{i=0}^{len-1}{a}_i{x}^i+\sum _{i=0}^{len-1}{b}_i\mathrm{i}{x}^i{)}^2& =(\sum _{i=0}^{len-1}{a}_i{x}^i{)}^2-(\sum _{i=0}^{len-1}{b}_i{x}^i{)}^2+2\sum _{i=0}^{len}{a}_i{x}^i\sum _{j=0}^{len}{b}_j{x}^j\mathrm{i}\end{array}$$

这样我们只需要计算 $H^2$ ，再取其虚部的 $\frac{1}{2}$ 即可，但是这样做的误差会大一点。

到这里FFT的内容就结束了（其实还有一个DIF DIT来着，但我还没学会）。

$$ToBeContnued\dots \dots$$

后记

开始看FFT的时候感觉好难，甚至以为单位根是最难的东西，后来发现这比着后面的部分简单的要多。许多东西也是边学边写的，如果有不当之处还请各位大佬指出，本人将不胜感激。

最后再写个无关的东西放松大脑吧

三项式定理

$$(a+b+c{)}^n=\sum _{i=0}^n{\displaystyle (\genfrac{}{}{0ex}{}{n}{i})}(b+c{)}^i{a}^{n-i}=\sum _{i=1}^n\sum _{j=0}^i{\displaystyle (\genfrac{}{}{0ex}{}{n}{i})}{\displaystyle (\genfrac{}{}{0ex}{}{i}{j})}{a}^{n-i}{b}^{i-j}{c}^j$$

其实这个是可以推广为 $n$ 项式的。自己探索吧（其实是我觉得式子太麻烦不想写了）。