前言 关于为什么又开始在这个博客上写算法

FFT是我高中学竞赛阶段接触到的最后一个算法，也是我高中一直没有学会的算法之一，在高中毕业的时候认为自己大学本科一定学的是CS，所以想要把这些自己之前没有学会的算法都再琢磨一下，可是基于当时的知识储备我还完全不知道复数域和复平面以及多项式的各种定理，于是半途而废了，FFT也成了在算法学习上我的一个心结之一。
但是巧合的是，我大学本科并没有学CS，而是学了数学，这些我一直未曾学会的算法也就被一直搁置了两年，直到昨天我在上数据结构课时，数学院数据结构课所学内容跟OI比起来简直不值一提，于是无聊的我突发奇想翻开了自己之前洛谷的提交记录，发现FFT的版子还一直放在我的题单，于是就想利用这些时间把这些之前没学懂的算法解决掉好了，也为明年考研跨考CS打下基础。
这次在数学院历练了两年多的我拿着高等代数数学分析复变函数的知识再来看FFT，顿时觉得这些知识难度与我两年前看的时候相比要简单了许多，在昨天晚上解决了自己所有的疑难点实现了代码拿到了AC。
看到我的这个博客之前发布的线段树网络流等等算法讲解总浏览量都接近一万了，那就继续吧，让这个断更了4年的博客重新活起来！

奈芙莲封面是这个博客的灵魂！！

问题背景

多项式乘法
给定一个次多项式和次多项式，求出与的卷积。

输入格式

第一行两个整数。
接下来一行个数字，从低到高表示的系数。
接下来一行个数字，从低到高表示的系数。

输出格式

一行个数字，从低到高表示的系数。

输入输出样例

输入样例

1 2
1 2
1 2 1

输出样例

1 4 5 2

问题分析

假设两个多项式和，两个多项式可以写作：

传统方法是利用两个多项式的系数进行卷积运算，得到一个次多项式：

这种卷积运算的时间复杂度为，显然在数据范围较大的情况下难以承受，而利用快速傅里叶变换可将时间复杂度降为。

FFT介绍

快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。
FFT（Fast Fourier Transformation） 是离散傅氏变换（DFT）的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。
——百度百科

要解决的问题是多项式的乘法，而一个多项式的表示方法并不唯一，传统意义上多项式的表示利用的是系数表示法，即对于一个多项式可由一个系数向量唯一表示。

而除了系数表示法之外，多项式也可以利用点值表示，对于多项式，选定个值带入多项式进行计算，得到个点值这个点值也可唯一表示多项式。

因此当我们利用同一个向量得到了两个两个多项式的点值表示法后，用对应点值相乘，得到即为两个多项式的乘积多项式的点值表示，这个过程的时间复杂度为。

需要注意的一点是，当两个多项式次数为时，他们的乘积多项式次数为，因此利用点值表示计算时，计算两个多项式的点值表示时应选用个变量，才能使得到的结果唯一表示乘积多项式。

然而对于一个多项式来说，代入任意选定的变量计算他的点值表示法时间复杂度依然是，并没有起到优化的效果，而快速傅里叶变换解决了这个问题，使得系数表示法转化为点值表示法的时间复杂度降低为。

快速傅里叶变换FFT

FFT在计算多项式的系数表示法变换为点值表示法时，选定复平面上单位圆上的单位复根作为变量计算多项式的点值，在这里单位根满足以下的一些性质（如果有不理解的可以自行查阅复数的一些相关知识）：

以上的这些性质都可以由Euler公式得到，推导过程并非FFT重点这里就省略了。

对于一个多项式，我们可以对其进行划分，将偶数次项与奇数次项分开，在这里假设为奇数，得到：

我们分别定义两个多项式：

那么原多项式就可以表示为：

将代入上式：

将代入上式：

由此可以发现，和在计算的过程中只有一个符号不同，因此在进行枚举计算时即可直接得到的值，利用这种方法进行分治，便可以将复杂度降至。

逆傅里叶变换IFFT

利用上述的方法得到了乘积多项式的点值表示法，那么现在需要解决的问题时如何将点值表示法再转换回系数表示法。

假设得到多项式的FFT点值表示为，其系数表示为，根据FFT原理，可如下表示：

取的共轭复数，如下定义向量：

那么由定义可以推导出如下的公式：

对于复平面上的单位根，有如下的性质：

因此可以得到：

因此，利用FFT得到了多项式的点值表示后只需要将变量换为原本选定的单位根的共轭复数再进行一次FFT就能得到多项式的系数表示。

这里需要说明一个问题，我们以上的讨论都是建立在为的幂次的条件下的，那么当不是的幂次时需要将扩大为大于的最小的的幂次，在进行逆傅里叶变换时，通过上述的推导可以发现，当我们计算的中的值大于原本的时，就不存在的情况了，因此求得的为，表示该多项式的次项系数为。

代码实现

下面是利用递归实现FFT的代码，具体的解释见代码注释好了。

#include<bits/stdc++.h>
#define rep(i,a,b) for(register int i=a;i<=b;++i)
using namespace std;

const int N=1e6+10;
const double Pi=acos(-1.0);
int n,m,l=1;

inline int read(){
    int x=0,f=1;
    char c=getchar();
    while(c<'0' || c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c<='9' && c>='0'){x=(x<<3)+(x<<1)+(int)(c-'0');c=getchar();}
    return x*f;
}

struct Complex{
    double re,im;
    Complex (double xx=0,double yy=0){re=xx,im=yy;}
};

Complex operator + (Complex a,Complex b){
    return Complex(a.re+b.re,a.im+b.im);
}

Complex operator - (Complex a,Complex b){
    return Complex(a.re-b.re,a.im-b.im);
}

Complex operator * (Complex a,Complex b){
    return Complex(a.re*b.re-a.im*b.im,a.im*b.re+a.re*b.im);
}

Complex a[N<<1],b[N<<1];

inline void FFT(Complex *a,int t,int inv){//inv表示当前进行FFT还是IFFT
    if(t==1)return;//计算长度为1时回溯
    int mid=t>>1;
    static Complex tmp[N];
    rep(i,0,mid-1)
        tmp[i]=a[2*i],tmp[i+mid]=a[2*i+1];//将多项式按照次数奇偶进行二分
    rep(i,0,t-1)a[i]=tmp[i];
    FFT(a,mid,inv);
    FFT(a+mid,mid,inv);
    Complex wn(cos(2.0*Pi/t),inv*sin(2.0*Pi/t)),w(1,0);//单位根
    rep(i,0,mid-1){
        tmp[i]=a[i]+w*a[i+mid];
        tmp[i+mid]=a[i]-w*a[i+mid];
        w=w*wn;
    }
    rep(i,0,t-1)a[i]=tmp[i];
}

int main(){
    n=read(),m=read();
    rep(i,0,n) a[i].re=read();
    rep(i,0,m) b[i].re=read();
    while(l<=(n+m))l<<=1;//找到大于等于n+m的最小2的幂次
    FFT(a,l,1);
    FFT(b,l,1);
    rep(i,0,l)
        a[i]=a[i]*b[i];
    FFT(a,l,-1);
    rep(i,0,n+m)
        printf("%d ",(int)(a[i].re/l+0.5));//输出整数，需要注意精度
    return 0;
}

但是问题到这里并没有结束

当有的毒瘤数据范围非常大的时候，用递归进行计算时，大量的递归会造成栈溢出，那么是否有不用递归的做法？

FFT的迭代实现

对于这样一个序列，我们观察对其进行二分的过程：

Rader算法

Rader算法即为实现上述操作的一种算法，对于个数，我们把递增自然数称为顺序数列；对顺序数列中的每一个数，将其二进制倒序后转化为十进制，称为倒序数列。
对于一个顺序数列，第个数的二进制可以视为将第（这里是整除）个数的二进制左移一位，再根据的奇偶性对其末尾加或者不加。
那么要得到它的倒序数列，只需要将这个操作反向进行即可，即第个数的二进制可以视为将第个数的二进制右移一位，再根据的奇偶性对其最高加或者不加，这里最高位即为第位。

迭代进行FFT（蝴蝶变换）

利用Rader算法求得了递推序列，那么如何通过迭代得到最终的答案？
这其实跟迭代实现01背包的做法思路差不多，对于求在次单位根的各幂次的点值时，次单位根的各幂次在或处的点值已经被计算并且储存在了数组中，那么在下一层的迭代过程中直接使用数组存储的答案继续进行迭代计算即可。

迭代优化FFT代码实现

详细解释依旧见代码注释。

#include<bits/stdc++.h>
#define rep(i,a,b) for(register int i=a;i<=b;++i)
using namespace std;

const int N=1e7+10;
const double Pi=acos(-1.0);
int n,m,l=1,t=0;
int bin[N<<1];

inline int read(){
    int x=0,f=1;
    char c=getchar();
    while(c<'0' || c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c<='9' && c>='0'){x=(x<<3)+(x<<1)+(int)(c-'0');c=getchar();}
    return x*f;
}

struct Complex{
    double re,im;
    Complex (double xx=0,double yy=0){re=xx,im=yy;}
};

Complex operator + (Complex a,Complex b){
    return Complex(a.re+b.re,a.im+b.im);
}

Complex operator - (Complex a,Complex b){
    return Complex(a.re-b.re,a.im-b.im);
}

Complex operator * (Complex a,Complex b){
    return Complex(a.re*b.re-a.im*b.im,a.im*b.re+a.re*b.im);
}

Complex a[N<<1],b[N<<1];

inline void FFT(Complex *A,int inv){
    rep(i,0,l-1)
        if(i<bin[i]) swap(A[i],A[bin[i]]);//得到倒序数列
    for(int mid=1;mid<l;mid<<=1){//枚举当前迭代区间的中点
        Complex wn(cos(Pi/mid),inv*sin(Pi/mid));//单位根
        for(int r=mid<<1,j=0;j<l;j+=r){//j枚举迭代区间位置，用r来得到区间右端点
            Complex w(1,0);
            rep(k,0,mid-1){//k枚举当前迭代区间
                Complex x=A[k+j],y=w*A[k+j+mid];
                A[k+j]=x+y;
                A[k+j+mid]=x-y;
                w=w*wn;
            }
        }
    }
}

int main(){
    n=read(),m=read();
    rep(i,0,n) a[i].re=read();
    rep(i,0,m) b[i].re=read();
    while(l<=(n+m)){
        l<<=1,t++;
    }
    rep(i,0,l)//Rader算法
        bin[i]=(bin[i>>1]>>1)|((i&1)<<(t-1));
    FFT(a,1);FFT(b,1);
    rep(i,0,l)
        a[i]=a[i]*b[i];
    FFT(a,-1);
    rep(i,0,n+m)
        printf("%d ",(int)(a[i].re/l+0.5));
    return 0;
}