FFT快速傅立叶变换

学习数学真是一件赛艇的事
FFT是我到目前OI的数学相关学过最难的了
其实理解以后发现并不是很难,只是需要的基础知识比较多

---

前置技能

主要包括复数相关,线性代数相关,分治基础

---

复数相关,复平面向量

可汗学院讲的真不错,强烈推荐看一看,内容比较全面,不过时间有点长
这里对复数做一点简单的总结

定义:

$$i^2=-1$$

这里的\(i\)是虚数单位(imaginary number)
对于任何数\(z\)都可以写成形如:

$$z=a+b*i$$

对于我们之前说的实数,放在这个形式里就是b=0
其中a称为实部,\(bi\)称为虚部

---

复数的四则运算

在这里我们只需要掌握加减乘三则运算就好
设\(A=a+bi\), \(B=c+di\)

加法:

\(A+B=(a+c)+(b+d)i\)

减法

\(A+B=(a-c)+(b-d)i\)

乘法

\(A*B=(a+bi)(c+di)=ac+adi+bci+bdi^2 =ac-bd+adi+bci=(ac-bd)+(ad+bc)i\)

共轭复数

若\(z=a+bi,z'=a-bi\)则称\(z'是z的共轭复数\)

---

复平面

首先定义一个复平面,x轴表示实部,y轴表示虚部
那么任何一个复数都可以在这个平面上以一个向量的形式表示出来

然后考虑四则运算的几何表示
加法和减法依然满足平行四边形法则
乘法根据棣莫弗定理要记住模长相乘,幅角相加

---

单位根

在复平面上以原点为圆心,1为半径作圆得到\(单位圆\)
设\((\omega_n)^n=1\),称\(\omega_n\)为n次单位根
在复平面上我们可以想象成用n条射线,从实轴开始,把圆均分成n部分
第一条射线与单位圆的交点所形成的向量
假设n=16,如图

\((\omega_n)^n\)就是n个\(\omega_n\)相乘(n-1)得到的向量,
因为要满足模长相乘,幅角相加
模长都为1不变,幅角=\(\frac{2\pi*(n-1)}{n}\)
就相当于
整个平面被均分成了16份,所以旋转15次以后又回到(1,0),所以\(\omega_n^n=1\).
对于单位根我们要知道它的几个性质

1. \(\omega_{2n}^{2k}=\omega_n^k\)

把平面分成2n格,旋转2k格=把平面分成n格,旋转k格.

2. \(\omega_n^{n+k}=\omega_n^k\)

指数函数的性质.

3. \(\omega_n^{\frac{n}{2}+k}=-\omega_n^k\)

因为\(omega_n^{\frac{n}{2}}=-1\).

4. \(\omega_n=\cos\)\(\frac{2\pi}{n}+i\sin\)\({\frac{2\pi}{n}}\)

几何性质

---

线性代数相关

对于这一方面要求理解的并不是很多,只要了解矩阵乘法和矩阵的逆就好了
不了解也没有关系,接下来会详细说明

---

多项式乘法

例题:已知两个n次多项式\(A=a_0+a_1x+a_2x^2+...+a_nx^n\)和\(B=b_0+b_1x+b_2x^2+...+b_nx^n\),求A*B的各项系数.
首先看到这道题的做法就是\(O(n^2)\)的暴力,将A的每个系数与B的每个系数相乘
想一想有没有可以优化的地方?
然而并没有......
接下来就需要FFT的操作了

系数表示法与点值表示法

对于一个多项式,我们最常用的把他表示出来的方法是系数表示法
就是形如\(A=a_0+a_1x+a_2x^2+...+a_nx^n\)的式子
其实还有另外一种表示方法点值表示法
\((x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2))...(x_n,f(x_n)).\)
就像两点确定一条直线,三点确定一条抛物线一样,n+1个点能确定一个n次多项式
我们发现对于多项式A和B
\(A=(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2))...(x_n,f(x_n))\).
\(B=(x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2))...(x_n,g(x_n))\).
\(A*B=(x_0,f(x_0)\)\(*\)\(g(x_0)),(x_1,f(x_1)\)\(*\)\(g(x_1))...(x_n,f(x_n)\)\(*\)\(g(x_n)).\)

这个操作是\(O(n)\)的
这个\(O(n)\)给我们提供了一个很好的思路,我们可以通过某种方法将系数表示变成点值表示,
再O(n)计算A*B,最后再通过某种方法将点值表示变回系数表示
一张图

考虑第一个奇怪的方法,如果采用暴力赋值计算,复杂度还是\(O(n^2)\)的
快速幂?naive 更慢!\(O(n^2logn)\)
所以我们不得不采用一种特殊的方法
给你一个多项式

\(A(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+...+a_{n-1}x^{n-1}\)

我们设

\(A_0(x)=a_0+a_2x+a_4x^2+a_6x^3+...+a_{n-2}x^{\frac{n}{2}}\)

\(A_1(x)=a_1+a_3x+a_5x^2+a_7x^3+...+a_{n-1}x^{\frac{n}{2}}\)

可以发现

\(A(x)=A_0(x^2)+xA_1(x^2)\)

然后就是一步骚操作,令\(x=\omega_{n}^k\)

\[A(\omega_{n}^k)=A_0(\omega_{n}^{2k})+\omega_{n}^kA_1(\omega_{n}^{2k})$$$$=A_0(\omega_{\frac{n}{2}}^{k})+\omega_{n}^kA_1(\omega_{\frac{n}{2}}^{k})$$又令$x=\omega_n^{\frac{n}{2}+k}$$$A(\omega_n^{\frac{n}{2}+k})=A_0(\omega_{n}^{n+2k})+\omega_n^{\frac{n}{2}+k}A_1(\omega_{n}^{n+2k})$$$$=A_0(\omega_{n}^{2k})-\omega_{n}^kA_1(\omega_{n}^{2k})$$$$=A_0(\omega_{\frac{n}{2}}^{k})-\omega_{n}^kA_1(\omega_{\frac{n}{2}}^{k})$$比较上面两个式子,我们发现 假设我们已经知道了$A_0(\omega_{\frac{n}{2}}^{k})$和$A_1(\omega_{\frac{n}{2}}^{k})$,我们就能同时知道$A(\omega_{n}^k)$和$\omega_n^{\frac{n}{2}+k}$. 这样就把问题缩小了一半,对于每个问题都缩小一半,复杂度就从$O(n^2)$变成了$O(nlogn)$. 就可以利用一种类似于线段树的操作来计算,如果没有理解就看看这张盗来的图  每一层向上转移是$O(n)的$,因为树高只有$logn$层,总复杂度就变得很小 此时我们的第一步由系数到点值已经结束了,不过还要补充一点 观察树的最底层 序列为$0,4,2,6,1,5,3,7$ 原序列$0,1,2,3,4,5,6,7$ 转化成二进制来发现规律 新序列$000,100,010,110,101,011,111$ 原序列$000,001,010,011,101,110,111$ 我们发现新序列的每个数是原序列的二进制反转 现在我们已经知道了每个序列的下标,我们就可以利用下标实现 这种实现方法更像是倍增,而不是分治,虽然算法的本质还是分治 如果要看代码在最后面... --- ###IDFT 对于把点值表示转化成系数表示,我们已经知道了用分治的方法快速求. 抛开分治以及log算法不谈,我们可以从另一方面理解刚才的操作 #####线性代数 把所有系数放在一起组成一个系数列向量: \]

\left[
\begin{matrix}
a_0 \
a_1 \
a_2 \
\vdots \
a_{n-1} \
\end{matrix}
\right]

\[最终的n个点值也可以放在一起组成一个点值列向量 \]

\left[
\begin{matrix}
y_0 \
y_1 \
y_2 \
\vdots \
y_{n-1} \
\end{matrix}
\right]

\[系数列向量可以通过左乘一个矩阵得到点值列向量 \]

\left[
\begin{matrix}
(x_0)^0 & (x_0)^1 &(x_0)^2 &\cdots & (x_0)^{n-1} \
(x_1)^0 & (x_1)^1 &(x_1)^2 &\cdots & (x_1)^{n-1} \
(x_2)^0 & (x_2)^1 &(x_2)^2 &\cdots & (x_2)^{n-1} \
\vdots & \vdots &\vdots & \ddots & \vdots \
(x_{n-1})^0 & (x_{n-1})^1 &(x_{n-1})^2 &\cdots & (x_{n-1})^{n-1} \
\end{matrix}
\right]*\left[
\begin{matrix}
a_0 \
a_1 \
a_2 \
\vdots \
a_{n-1} \
\end{matrix}
\right]=\left[
\begin{matrix}
y_0 \
y_1 \
y_2 \
\vdots \
y_{n-1} \
\end{matrix}
\right]

\[因为我们代入的$x_0,x_1,x_2...x_{n-1}$分别是$\omega_n^0,\omega_n^1\omega_n^2...\omega_n^{n-1}$. 矩阵就变成了\]

\left[
\begin{matrix}
(\omega_n⁰⁾0 & (\omega_n⁰⁾1 &(\omega_n⁰⁾2 &\cdots & (\omega_n⁰⁾ \
(\omega_n¹⁾0 & (\omega_n¹⁾1 &(\omega_n¹⁾2 &\cdots & (\omega_n¹⁾ \
(\omega_n²⁾0 & (\omega_n²⁾1 &(\omega_n²⁾2 &\cdots & (\omega_n²⁾ \
\vdots & \vdots &\vdots & \ddots & \vdots \
(\omega_n^{n-1})0 & (\omega_n^{n-1})1 &(\omega_n^{n-1})2 &\cdots & (\omega_n^{n-1}) \
\end{matrix}
\right]*\left[
\begin{matrix}
a_0 \
a_1 \
a_2 \
\vdots \
a_{n-1} \
\end{matrix}
\right]=\left[
\begin{matrix}
y_0 \
y_1 \
y_2 \
\vdots \
y_{n-1} \
\end{matrix}
\right]

\[设这个矩阵为 ##$$A*B=C\]

在刚才的变换中,我们已知\(B\)向量,然后通过分治算法求得\(C\)向量.
现在我们通过\(O(n)\)时间的乘法得到新的\(C\)向量,也就是新的点值表示,我们要重新求回原来的\(B\)向量,怎么办
考虑逆矩阵.

$$A^{-1}*A*B=A*C$$

$$B=A^{-1}*C$$

我们只要求出\(A\)的逆矩阵就可以求出\(B\)向量了
\(A\)的逆矩阵并不好求,我们只需要知道它是什么就好了

\(A\)矩阵是一个特殊的范德蒙德矩阵,范德蒙德矩阵就是指每一行的元素为一个等比数列.
对于这个矩阵,我们有

$$A^{-1}=\frac{1}{n}*\overline{A}$$

其中\(\overline{A}\)是\(A的共轭矩阵\),共轭矩阵就是指矩阵内的所有元素都取共轭复数得到的矩阵.
具体证明最后再讲.

有了这个性质我们重新看看最开始的式子:

$$B=A^{-1}C\B=\frac{1}{n}\overline{A}*C$$

\[\left[ \begin{matrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \\ \end{matrix} \right]= \frac{1}{n} \left[ \begin{matrix} (\omega_n^0)^0 & (\omega_n^0)^1 &(\omega_n^0)^2 &\cdots & (\omega_n^0)^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 &(\omega_n^{-1})^2 &\cdots & (\omega_n^{-1})^{n-1} \\ (\omega_n^{-2})^0 & (\omega_n^{-2})^1 &(\omega_n^{-2})^2 &\cdots & (\omega_n^{-2})^{n-1} \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 &(\omega_n^{-(n-1)})^2 &\cdots & (\omega_n^{-(n-1)})^{n-1} \\ \end{matrix} \right]*\left[ \begin{matrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n-1} \\ \end{matrix} \right] \]

是不是很神奇?原来的只要将我们分治的时候代入的\(\omega_n^x\)换成\(\omega_n^{-x}\)就可以求出\(B\)向量,也就是系数表示了.

---

现在证明

$$A^{-1}=\frac{1}{n}*\overline{A}$$

设

$$\frac{1}{n}\overline{A}A=C$$

我们就是要证明\(C=U\),\(U\)是单位矩阵

\[\frac{1}{n}\left[ \begin{matrix} (\omega_n^0)^0 & (\omega_n^0)^1 &(\omega_n^0)^2 &\cdots & (\omega_n^0)^{n-1} \\ (\omega_n^1)^0 & (\omega_n^1)^1 &(\omega_n^1)^2 &\cdots & (\omega_n^1)^{n-1} \\ (\omega_n^2)^0 & (\omega_n^2)^1 &(\omega_n^2)^2 &\cdots & (\omega_n^2)^{n-1} \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 &(\omega_n^{n-1})^2 &\cdots & (\omega_n^{n-1})^{n-1} \\ \end{matrix} \right]\left[ \begin{matrix} (\omega_n^0)^0 & (\omega_n^0)^1 &(\omega_n^0)^2 &\cdots & (\omega_n^0)^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 &(\omega_n^{-1})^2 &\cdots & (\omega_n^{-1})^{n-1} \\ (\omega_n^{-2})^0 & (\omega_n^{-2})^1 &(\omega_n^{-2})^2 &\cdots & (\omega_n^{-2})^{n-1} \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 &(\omega_n^{-(n-1)})^2 &\cdots & (\omega_n^{-(n-1)})^{n-1} \\ \end{matrix} \right]=\left[ \begin{matrix} 1 & 0 &0 &\cdots & 0 \\ 0 & 1 &0 &\cdots & 0 \\ 0 & 0 &1 &\cdots & 0 \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ 0 & 0 &0 &\cdots & 1 \\ \end{matrix} \right]\]

\(C_{i,j}=\frac{1}{n}\sum_{k=0}^{n-1}A_{i,k}*\overline{A}_{k,j}\)

当\(i=j\)时

\(A_{i,k}*\overline{A}_{k,i}=1\),因为模长相乘,幅角相加,模长为一的两个共轭复数相乘以后等于1.

\(C_{i,j}=\frac{1}{n}\sum_{k=0}^{n-1}1=1\)

当\(i\neq j\)时

\(C_{i,j}=\frac{1}{n}\sum_{k=0}^{n-1}A_{i,k}*\overline{A}_{k,j}\)

\(=\frac{1}{n}\sum_{k=0}^{n-1}(\omega_n^i)^k*(\omega_n^{-k})^j\)

\(=\frac{1}{n}\sum_{k=0}^{n-1}(\omega_n^{i-j})^k\)

\(\omega_n^{i-j}\)可以视为常数项,就变成了等比数列求和公式.

\(=\frac{1}{n}\frac{\omega_n^{i-j}[1-(\omega_n^{i-j})^n]}{1-\omega^{i-j}}\)

因为\((\omega_n^{i-j})^n=1,i!=j\).

\(C_{i,j}=0\).

所以\(C=U\)就是我们要证的单位矩阵.

---

对我们刚才的那些操作取个名字,第一步系数->点值的操作叫DFT,
第三步点值->系数的操作叫IDFT
所有操作合起来叫做FFT.
最后看一看代码
洛谷上有模板题.

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
#define maxn 5000010
const double Pi=acos(-1.0);
inline int read(){
    int ret=0,ff=1;
    char ch=getchar();
    while(ch<'0'||ch>'9'){
        if(ch=='-') ff=-ff;
        ch=getchar();
    }
    while(ch>='0'&&ch<='9'){
        ret=ret*10+ch-'0';
        ch=getchar();
    }
    return ret*ff;
}
int lim=1,l=0;
int r[maxn];
struct Complex{
    double x,y;
    Complex(double _x=0,double _y=0){x=_x,y=_y;}
    Complex operator+(Complex t){ return Complex(x+t.x,y+t.y);}
    Complex operator-(Complex t){ return Complex(x-t.x,y-t.y);}
    Complex operator*(Complex t){ return Complex(x*t.x-y*t.y,x*t.y+y*t.x);}
}a[maxn],b[maxn];
void fft(Complex *A,int op){
    for(int i=0;i<lim;++i) if(i<r[i]) swap(A[i],A[r[i]]);
    //用之前的r[i]存的顺序对A重新排列
    for(int Mid=1;Mid<lim;Mid<<=1){
    //披着倍增外套的分治
        Complex Wn(cos(Pi/Mid),op*sin(Pi/Mid));//因为Mid是我们枚举的中点,本来就是要求区间的2倍,所以把2约掉
        for(int R=Mid<<1,j=0;j<lim;j+=R){
            Complex w(1,0);
            for(int k=j;k<Mid+j;++k,w=w*Wn){
                Complex t1=A[k],t2=w*A[k+Mid];
                A[k]=t1+t2,A[k+Mid]=t1-t2;
            }
        }
    }
}
int main(){
//	freopen("mod.in","r",stdin);
    int n=read(),m=read();
    for(int i=0;i<=n;++i) a[i].x=read();
    for(int j=0;j<=m;++j) b[j].x=read();
    while(lim<=n+m) lim<<=1,++l;
    for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    //这一步就是二进制的转置部分,r[i]表示i的二进制转置,比如说r[6(110)]=3(011)
    //r[i]由r[i/2]递推得来,对比i和i/2的二进制规律,我们发现i=(i>>2)<<1+(i&1)
    //因为r[i]是i的倒序,所以也应该是倒序递推
    fft(a,1),fft(b,1);
    for(int i=0;i<lim;++i) a[i]=a[i]*b[i];
    fft(a,-1);
    for(int i=0;i<=n+m;++i) printf("%d ",int(a[i].x/lim+0.5));
    return 0;
}