快速傅里叶变换（FFT）

Posted 2020-12-18 manmanjiangqwq

快速傅里叶变换，简称FFT，是一种可以O(nlogn)的时间内计算n次多项式乘法的算法。

写得很好的博客：自为风月马前卒大佬的博客。

大致步骤是：

1.先将两个多项式的系数表达法O(nlogn)都转化成点值表达法。

例如：y=A₀+A₁*x+A₂*x²+A₃*x³+...+A_n*xⁿ可以转化为n+1个点(x₁,y₁),(x₂,y₂),(x₃,y₃),...,(x_n,y_n),(x_n+1,y_n+1)。

这n+1个点也能确定这个多项式，就像两点确定一条直线，三点确定一条抛物线一样。

同理第二个多项式也转化为n+1个点，这n+1个点的x要与第一个多项式的点的x分别相等。

2.x相同的点所对应的y可以O(n)直接相乘，达到多项式相乘的效果。

3.再将点值表达法转化为系数表达法就行了。

然后发现暴力转化的话时间是O(n²)，那么如何优化呢，这里就要牵扯到复数的概念和性质了，前面的博客里写的非常好，我就不再赘述了其实就是懒。

下面是递归版的代码和非递归版的代码：（递归版的时间又长空间又大好像没什么用，但应该可以帮助理解吧w）

#include<bits/stdc++.h>
using namespace std;
const int N=4e6+1e5;
const double pai=acos(-1.0);
double Co[N],Si[N];
struct Complex{
    double x,y;
    Complex(double xx=0,double yy=0){
        x=xx; y=yy;
    }
}A[N],B[N];
Complex operator+(Complex X,Complex Y) { return Complex(X.x+Y.x,X.y+Y.y); }
Complex operator-(Complex X,Complex Y) { return Complex(X.x-Y.x,X.y-Y.y); }
Complex operator*(Complex X,Complex Y) { return Complex(X.x*Y.x-X.y*Y.y,X.x*Y.y+X.y*Y.x); }
int rd(){
    int s=0,ff=1;
    char w=getchar();
    while(!isdigit(w))
        ff^=w==‘-‘,w=getchar();
    while(isdigit(w))
        s=s*10+w-‘0‘,w=getchar();
    return ff?s:-s;
}
void FFt(int n,Complex *a,int Opt){
    if(n==1) return ;
    Complex c[n>>1],d[n>>1],k;
    for(int i=0;i<n;i+=2)//注意这里没有等于号 
        c[i>>1]=a[i],d[i>>1]=a[i+1];
    FFt(n>>1,c,Opt); FFt(n>>1,d,Opt);
    Complex w=Complex(Co[n],Opt*Si[n]),x=Complex(1,0); n>>=1;
    for(int i=0;i<n;i++,x=x*w)
        k=x*d[i],a[i]=c[i]+k,a[i+n]=c[i]-k;
}
int main(){
    int n=1,N,M;
    N=rd(); M=rd();
    while(n<=N+M) n<<=1;
    for(int i=1;i<=n;i++)
        Co[i]=cos(2.0*pai/i),Si[i]=sin(2.0*pai/i);
    for(int i=0;i<=N;i++) A[i].x=rd();
    for(int i=0;i<=M;i++) B[i].x=rd();
    FFt(n,A,1); FFt(n,B,1);
    for(int i=0;i<=n;i++) A[i]=A[i]*B[i];
    FFt(n,A,-1);
    for(int i=0;i<=N+M;i++)
        printf("%d ",(int)(A[i].x/n+0.5));
    return 0;
}

递归版

#include<bits/stdc++.h>
using namespace std;
const int N=4e6+1e5;
const double pai=acos(-1.0);
int r[N];
double Co[N],Si[N];
struct Complex{
    double x,y;
    Complex(double xx=0,double yy=0){
        x=xx; y=yy;
    }
}A[N],B[N];
Complex operator+(Complex X,Complex Y) { return Complex(X.x+Y.x,X.y+Y.y); }
Complex operator-(Complex X,Complex Y) { return Complex(X.x-Y.x,X.y-Y.y); }
Complex operator*(Complex X,Complex Y) { return Complex(X.x*Y.x-X.y*Y.y,X.x*Y.y+X.y*Y.x); }
int rd(){
    int s=0,ff=1;
    char w=getchar();
    while(!isdigit(w))
        ff^=w==‘-‘,w=getchar();
    while(isdigit(w))
        s=s*10+w-‘0‘,w=getchar();
    return ff?s:-s;
}
void FFt(int n,Complex *a,int Opt){
    Complex x,w,X,Y;
    for(int i=0;i<n;i++) //注意无等于号 
        if(i<r[i]) swap(a[i],a[r[i]]);
    for(int mid=1;mid<n;mid<<=1){
        for(int l=0,j=(mid<<1);l<n;l+=j){
            w=Complex(Co[mid<<1],Opt*Si[mid<<1]); x=Complex(1,0);
            for(int i=l;i<l+mid;i++,x=x*w){
                X=a[i],Y=x*a[i+mid];
                a[i]=X+Y,a[i+mid]=X-Y;
            }
        }
    }
}
int main(){
    int n=1,N,M,l=0; N=rd(); M=rd();
    while(n<=N+M) n<<=1,l++;
    for(int i=1;i<n;i++)
        r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
    for(int i=1;i<=n;i++)
        Co[i]=cos(2.0*pai/i),Si[i]=sin(2.0*pai/i);
    for(int i=0;i<=N;i++) A[i].x=rd();
    for(int i=0;i<=M;i++) B[i].x=rd();
    FFt(n,A,1); FFt(n,B,1);
    for(int i=0;i<=n;i++) A[i]=A[i]*B[i];
    FFt(n,A,-1);
    for(int i=0;i<=N+M;i++)
        printf("%d ",(int)(A[i].x/n+0.5));
    return 0;
}

非递归版