快速傅里叶变换（FFT）

写在前面的话：

感觉最近每次考试都是同样的一种情形：考试的时候绞尽脑汁地去思考正解无果，在Solution发下来以后，却发现这是一种我没有学习过的算法......并且感觉这些东西貌似就我一个人不会了，赶紧来学习了一下，首先，便是FFT

（一）预备知识

要学习FFT，首先得学习一些关于复数的东西
这一篇写得不错，理解起来很容易，具体看here

（二）开始学习

附上我在学习FFT的时候用的3篇博客，很多内容都是一样的，但是中间有些知识可以互补
首先学习博客1，这篇博客很多东西的证明写得很清楚，有利于初学者
然后再看一看博客2，这篇博客中有两副图有助于理解，并且也详细解释了使用矩阵来计算IDFT的过程
全部理解了以后，自然是要写代码的呀，看一看博客3，这篇代码写得还是蛮详细的。

（三）重要知识

为了方便我以后复习，先将3篇博客的重要部分浓缩至此，略去所有证明以及基础概念性的东西，为什么？因为我懒得打字啊~~~

（1）多项式乘法

定义$ A(x)=\sum_{i=0}^{n-1} a_{i}x^{i} $ 与 \(B(x)=\sum_{i=0}^{n-1} b_{i}x^{i}\) 相乘的结果为\(C(x)\)，那么\(C(x)\)的值为
\[C(x)=\sum_{k=0}^{2n-2}(\sum_{k=i+j} a_{i}b_{j})x^{k}\]
这样做的时间复杂度无疑是\(O(n^{2})\)，显然不行
但是如果我们知道了两个多项式在2n-1个点处所取得的点值表示，即将那2n-1个点作为自变量x带入A与B之后所得的结果，那么我们就可以在\(O(n)\)时间之内计算出\(C(x)\)
那么我们该如何做到这一点呢？这就是为什么要用FFT的原因了。

（2）单位根

在复平面上，以原点为圆心，1为半径作圆，所得的圆叫做单位圆。以原点为起点，单位圆的n等分点为终点，作n个向量。设所得的幅角为正且最小的向量对应的复数为\(\omega_{n}\)，称为n次单位根。

在代数中，若\(z^{n}=1\)，我们把z称为n次单位根

性质一 ：\(\omega_{2n}^{2k}=\omega_{n}^{k}\)

性质二 ：\(\omega_{n}^{k+ \frac{n}{2}}=-\omega_{n}^{k}\)

（3）快速傅里叶变换

考虑多项式\(A(x)\)的表示。将n次单位根的0到n-1次幂带入多项式的系数表示，所得点值向量\((A(\omega_{n}^{0}),...,A(\omega_{n}^{n-1}))\)称为其系数向量的离散傅里叶变换。

Cooley-Tukey 算法是FFT的最常用算法，具体见3个博客

（4）傅里叶逆变换

将点值表示的多项式转化为系数表示，同样可以使用快速傅里叶变换，这个过程叫做傅里叶逆变换。

设 \((y_{0},...,y{n-1})\) 为 \((a_{0},...,a{n-1})\) 的傅里叶变换。考虑另一个向量 \((c_{0},...,c_{n-1})\) 满足
\[c_{k}=\sum_{i=0}^{n-1} y_{i}(\omega_{n}^{-k})^{i}\]
即多项式 \(B(x)=y_{0}+y_{1}x+y_{2}x^{2}+...+y_{n-1}x^{n-1}\) 在 \(\omega_{n}^{0},\omega_{n}^{-1},...,\omega_{n}^{-(n-1)}\)处的点值表示。
。。。。。。（此处省略一大堆证明）
然后，我们可以得到：使用单位根的倒数代替单位根，做一次类似快速傅里叶变换的过程，再将结果每个数除以 n，即为傅里叶逆变换的结果。

（5）流程图

总得来说，FFT的程序流程图长这样

图还是蛮好的，很形象

（6）迭代实现

图像更清晰

但是细心的话可以发现，其实这个图右下角有些问题，画反了，反正不是我画的，所以呢，这锅我就不背了，哈哈。

（7）蝴蝶操作

懒得写了，看博客1，写得很详细

（8）代码

哈哈，终于可以上代码了！
原题，洛谷P3803

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
typedef long long ll;
typedef double dd;
#define For(i,j,k) for (int i=j;i<=k;++i)
#define Forr(i,j,k) for (int i=j;i>=k;--i)
#define Set(a,p) memset(a,p,sizeof(a))
using namespace std;

template<typename T>bool chkmax(T& a,T b) {return a<b?a=b,1:0;}
template<typename T>bool chkmin(T& a,T b) {return a>b?a=b,1:0;}

const int maxn=5e6+1e2;
const dd Pi=acos(-1.0);
struct Complex {
    dd x,y;
}a[maxn],b[maxn];
int n,m,N,cnt;
int p[maxn];

inline int read() {
    int x=0,p=1;
    char c=getchar();
    while (!isdigit(c)) {if (c=='-') p=-1; c=getchar();}
    while (isdigit(c)) {x=(x<<1)+(x<<3)+(c-'0'); c=getchar();}
    return x*p;
}

Complex operator + (const Complex &aa,const Complex &bb) {
    return (Complex){aa.x+bb.x,aa.y+bb.y};
}

Complex operator - (const Complex &aa,const Complex &bb) {
    return (Complex){aa.x-bb.x,aa.y-bb.y};
}

Complex operator * (const Complex &aa,const Complex &bb) {
    return (Complex){aa.x*bb.x-aa.y*bb.y,aa.x*bb.y+aa.y*bb.x};
}

inline void FFT(Complex *s,int type) {
    For (i,0,N-1)
        if (i<p[i]) swap(s[i],s[p[i]]);//得到新序列
    for (int mid=1;mid<N;mid<<=1) { //枚举合并区间的中点
        Complex Wn=(Complex){cos(Pi/mid),type*sin(Pi/mid)}; //单位根
        for (int len=mid<<1,j=0;j<N;j+=len) { 
                        //len是区间的长度，j是区间的左端点
            Complex w=(Complex){1,0}; //幂
            for (int k=0;k<mid;++k,w=w*Wn) {
                                  //枚举左区间
                Complex x=s[j+k],y=w*s[j+mid+k];
                s[j+k]=x+y; s[j+mid+k]=x-y;
                                 //蝴蝶操作
            }
        }
    }
}

int main() {
    n=read(); m=read();
    For (i,0,n) a[i].x=read();
    For (i,0,m) b[i].x=read();
    for (N=1;N<=n+m;N<<=1,++cnt) ;
    For (i,0,N-1)
        p[i]=(p[i>>1]>>1) | ((i&1)<<(cnt-1));
/*  这里就是翻转，对于原序列的第i个数，将i在二进制下翻转之后就是原来元素在新序列
    中的位置。这里预处理，因为要翻转i，其实就是先将i/2翻转之后，再看i的奇偶，判
    断最高位。*/
    FFT(a,1); FFT(b,1);
    For (i,0,N) a[i]=a[i]*b[i];
    FFT(a,-1);
    For (i,0,n+m) printf("%d ",(int)(a[i].x/N+0.5));
    return 0;
}

（8）注意

\(\Pi\) 是doule型的，不要习惯性的定义为int

转载于:https://www.cnblogs.com/Wuweizheng/p/8553081.html