快速傅里叶变换

定义

快速傅里叶变换(Fast Fourier Transform,FFT)用于在\(O(n\log n)\)时间内求解多项式乘法。

有关概念

  1. 虚数初步:人教数学A版选修2-2第三章;

  2. 单位根:

    \(n\)次单位根指满足\(z^n=1\)的复数,这样的数总共有\(n\)个,分别为\(\omega_n^k=e^{\frac{2\pi ik}{n}},k=0,1,\ldots n-1\)

    根据欧拉公式\(e^{ik}=\cos(k)+i\sin(k)\),故\(\omega_n^k=\cos(\frac{2\pi ik}{n})+i\sin(\frac{2\pi ik}{n})\)

    几何意义上,\(n\)次单位根均匀分布在复平面的单位圆上,在\(n=8\)时分布如下:

  3. 多项式的表示方法:

    系数表达:如多项式\(2x^4+x^2+3x\)可以被表示为向量\(\vec{a}=(3,1,0,2)\)

    点值表达:对于给定的\(n+1\)个不同的点,可以唯一确定一个\(n\)次多项式函数,一个多项式有无数种点值表达,将点值表达转换为系数表达的过程叫做插值,常用的有拉格朗日插值法

  4. 多项式乘法:

    给定两个多项式\(A(x),B(x)\)

    \[A(x)=\sum_{i=0}^n a_ix^i\]

    \[B(x)=\sum_{i=0}^n b_ix^i\]

    相乘得到\(C(x)\)

    \[C(x)=\sum_{j+k=i,0\leq j,k\leq n} a_jb_kx^i\]

    系数表达求解的复杂度为\(O(n^2)\),但是点值表达如果采样点相同,即如果\(A(x)\)\(B(x)\)点值表达是横坐标对应相同的\(n+1\)个点,就可以直接相乘,求解复杂度为\(O(n)\)

DFT和IDFT

离散傅里叶变换(Discrete Fourier Transform, DFT)可将系数表达转换为点值表达;

对于多项式\(A(x)=\sum_{i=0}^{n-1} a_ix^i\)的系数向量\(\vec{a}=(a_0,a_1,\ldots,a_{n-1}\),其DFT定义为点值向量\(\vec{y}=A(\omega_n^0),A(\omega_n^1),\ldots,A(\omega_n^{n-1})\)

IDFT(Inverse Discrete Fourier Transform)是DFT的逆,用于将点值表达转换为系数表达,即求解系数向量\(\vec{a}\),给出方程组:

\[\left \{ \begin{array}{ccccccccc}&a_0(\omega_n^0)^0&+&a_1(\omega_n^0)^1&+&\ldots&+&a_{n-1}(\omega_n^0)^{n-1}&=&A(\omega_n^0)&\\&a_0(\omega_n^1)^0&+&a_1(\omega_n^1)^1&+&\ldots&+&a_{n-1}(\omega_n^1)^{n-1}&=&A(\omega_n^1)&\\&\vdots&&\vdots&&\vdots&&\vdots&&\vdots&\\&a_0(\omega_n^{n-1})^0&+&a_1(\omega_n^{n-1})^1&+&\ldots&+&a_{n-1}(\omega_n^{n-1})^{n-1}&=&A(\omega_n^{n-1})&\end{array} \right.\]

矩阵形式:

\[\begin{bmatrix} &(\omega_n^0)^0&(\omega_n^0)^1&\ldots&(\omega_n^0)^{n-1}&\\ &(\omega_n^1)^0&(\omega_n^1)^1&\ldots&(\omega_n^1)^{n-1}&\\ &\vdots&\vdots&\vdots&\vdots&\\ &(\omega_n^{n-1})^0&(\omega_n^{n-1})^1&\ldots&(\omega_n^{n-1})^{n-1}&\\ \end{bmatrix} \begin{bmatrix} a_0\\a_1\\ \vdots\\a_{n-1}\end{bmatrix} =\begin{bmatrix} A(\omega_n^0)\\A(\omega_n^1)\\ \vdots\\A(\omega_n^{n-1})\end{bmatrix}\]

设:

\[V=\begin{bmatrix} &(\omega_n^0)^0&(\omega_n^0)^1&\ldots&(\omega_n^0)^{n-1}&\\ &(\omega_n^1)^0&(\omega_n^1)^1&\ldots&(\omega_n^1)^{n-1}&\\ &\vdots&\vdots&\vdots&\vdots&\\ &(\omega_n^{n-1})^0&(\omega_n^{n-1})^1&\ldots&(\omega_n^{n-1})^{n-1}&\\ \end{bmatrix}\quad D=\begin{bmatrix} &(\omega_n^0)^0&(\omega_n^0)^1&\ldots&(\omega_n^0)^{n-1}&\\ &(\omega_n^{-1})^0&(\omega_n^{-1})^1&\ldots&(\omega_n^{-1})^{n-1}&\\ &\vdots&\vdots&\vdots&\vdots&\\ &(\omega_n^{-(n-1)})^0&(\omega_n^{-(n-1)})^1&\ldots&(\omega_n^{-(n-1)})^{n-1}&\\ \end{bmatrix} \quad E=V\cdot D\]

\(e_{ij}=\sum_{k=0}^{n-1} v_{ik}d_{kj}=\sum_{k=0}^{n-1} \omega_n^{ik}\omega_n^{-kj}=\sum_{k=0}^{n-1} \omega_n^{(i-j)k}\)

\(i=j\)时,\(e_{ij}=n\)

\(i\not=j\)时,\(e_{ij}=\sum_{k=0}^{n-1} \omega_n^{(i-j)k}=\sum_{k=0}^{n-1} (\omega_n^{i-j})^k=\frac{1-(\omega_n^{i-j})^n}{1-\omega_n^{i-j}}=0\)

\(E\)矩阵在乘法中可看做数\(n\)

\(V^{-1}=\frac 1n D\),乘入原等式:

\[\begin{bmatrix} a_0\\a_1\\ \vdots\\a_{n-1}\end{bmatrix} =\frac 1n \begin{bmatrix} &(\omega_n^0)^0&(\omega_n^0)^1&\ldots&(\omega_n^0)^{n-1}&\\ &(\omega_n^{-1})^0&(\omega_n^{-1})^1&\ldots&(\omega_n^{-1})^{n-1}&\\ &\vdots&\vdots&\vdots&\vdots&\\ &(\omega_n^{-(n-1)})^0&(\omega_n^{-(n-1)})^1&\ldots&(\omega_n^{-(n-1)})^{n-1}&\\ \end{bmatrix} \begin{bmatrix} A(\omega_n^0)\\A(\omega_n^1)\\ \vdots\\A(\omega_n^{n-1})\end{bmatrix} \]

于是,IDFT和DFT的过程相同,只需要把\(\omega_n^i\)换成\(\omega_n^{-i}\),再乘\(\frac 1n\)即可。

FFT

递归实现

思路与简单证明

至此为止,求DFT和IDFT的复杂度还是\(O(n^2)\),于是我们可以利用单位根的几个性质来优化复杂度:

  1. 消去引理:对于非负整数\(n,d,k\),有\(\omega_{dn}^{dk}=\omega_n^k\)

    证明:\(\omega_{dn}^{dk}=e^{\frac{2\pi idk}{dn}}=e^{\frac{2\pi ik}{n}}=\omega_n^k\)

  2. 折半引理:对于正偶数\(n\),有\((\omega_n^{k+\frac{n}{2}})^2=(\omega_n^{k})^2\)

    证明:由于\((\omega_n^\frac{n}{2})^2=\omega_n^n=1\),故\(\omega_n^\frac{n}{2}=\pm 1\),又\(\omega_n^\frac{n}{2}\not=\omega_n^n\),故\(\omega_n^\frac{n}{2}=-1\)

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

FFT采用分治的思想,将原多项式\(A(x)\)分成两个多项式\(A^{[0]}(x),A^{[1]}(x)\)

\[ A^{[0]}(x)=a_0+a_2x+a_4x^2+\ldots+a_{n-2}x^{\frac n2-1}\\ A^{[1]}(x)=a_1+a_3x+a_5x^2+\ldots+a_{n-1}x^{\frac n2-1} \]

可得\(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)\)

这时,根据消去引理,\((\omega_n^k)^2=\omega_n^{2k}=\omega_{\frac n2}{k}\),故求\(A^{[0]}(x^2)\)\(A^{[1]}(x^2)\)只需要代入\(\frac n2\)次单位根即可,使得代入的值也可以减半;

具体来说,\(A(\omega_n^k)=A^{[0]}((\omega_n^k)^2)+\omega_n^kA^{[1]}((\omega_n^k)^2)=A^{[0]}(\omega_{\frac n2}^k)+\omega_n^kA^{[1]}(\omega_{\frac n2}^k)\)

又根据折半引理,\(A(\omega_n^{k+\frac{n}{2}})=A(-\omega_n^{k})=A^{[0]}(\omega_{\frac n2}^k)-\omega_n^kA^{[1]}(\omega_{\frac n2}^k)\)

于是可以递归执行,复杂度优化到\(O(n\log n)\),边界条件为\(n=1\)

注意执行FFT之前先将次数\(n\)补齐到\(2\)的整数次幂。

实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<complex>
#define INF 0x3f3f3f3f
using namespace std;
typedef long long LL;
typedef double db;
typedef complex<double> cp;
inline int read()
{
int x=0,f=1;
char ch=getchar();
while(ch<'0'||ch>'9') { if(ch=='-')f=-1; ch=getchar(); }
while(ch>='0'&&ch<='9') { x=(x<<1)+(x<<3)+ch-'0'; ch=getchar(); }
return x*f;
}
const int MAXN=3e5+10;
int n,m;
cp x[MAXN],y[MAXN];
void FFT(int n,cp *a,int flag)
{
if(n==1)return;
cp a0[(n>>1)+1],a1[(n>>1)+1];
cp omega_n(cos(2*M_PI/n),sin(flag*2*M_PI/n)),omega(1,0);
for(int i=0;i<(n>>1);++i)a0[i]=a[i<<1],a1[i]=a[i<<1|1];
FFT(n>>1,a0,flag);
FFT(n>>1,a1,flag);
for(int i=0;i<(n>>1);++i)
{
a[i]=a0[i]+omega*a1[i];
a[i+(n>>1)]=a0[i]-omega*a1[i];
omega*=omega_n;
}
if(flag==-1)
for(int i=0;i<n;++i)
x[i]=int(x[i].real()/n+0.5);
return;
}
int main()
{
n=read();m=read();
for(int i=0;i<=n;++i)x[i]=read();
for(int i=0;i<=m;++i)y[i]=read();
m+=n;n=1;
while(n<=m)n<<=1;
FFT(n,x,1);
FFT(n,y,1);
for(int i=0;i<n;++i)x[i]=x[i]*y[i];
FFT(n,x,-1);
for(int i=0;i<=m;++i)printf("%d ",(int)x[i].real());
puts("");
return 0;
}

迭代实现

思路与简单证明

考虑递归时的奇偶分组情况,如图(图中数字为系数下标):

观察叶子节点中下标的排列情况,用二进制表示它们:

0 1 2 3 4 5 6 7
000 001 010 011 100 101 110 111
0 4 2 6 1 5 3 7
000 100 010 110 001 101 011 111

可以发现改变顺序后的排列是原排列每个数的二进制位逆序置换,如果我们知道了当前位置应该填的数,可以用类似于加法进位的方法求出下一个位置应该填的数:

从最高位开始,循环执行:检查当前位,若为\(0\),将其置为\(1\)并退出循环;否则将其置为\(0\),并检查低位;

求得当前位置的应有的下标后,可以与对应位置的数交换,因为二进制位逆序置换的数是一一对应的,所以每个数至多只会被交换一次;

以上过程被称为雷德算法。在求得最终排列后,就可以简单地循环执行,从小到大把数组分段操作即可。

实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<complex>
#define INF 0x3f3f3f3f
using namespace std;
typedef long long LL;
typedef double db;
typedef complex<double> cp;
inline int read()
{
int x=0,f=1;
char ch=getchar();
while(ch<'0'||ch>'9') { if(ch=='-')f=-1; ch=getchar(); }
while(ch>='0'&&ch<='9') { x=(x<<1)+(x<<3)+ch-'0'; ch=getchar(); }
return x*f;
}
const int MAXN=3e5+10;
const db PI=acos(-1);
int n,m;
cp x[MAXN],y[MAXN];
void rader(int n,cp *x)
{
for(int i=0,j=0;i<n;++i)
{
if(i<j)swap(x[i],x[j]);
int k=n>>1;
while((j^=k)<k)k>>=1;
}
return;
}
void FFT(int n,cp *x,int flag)
{
rader(n,x);
for(int len=2;len<=n;len<<=1)
{
cp omega_n(cos(flag*2*PI/len),sin(flag*2*PI/len));
for(int i=0;i<n;i+=len)
{
cp omega(1,0);
for(int j=i;j<i+(len>>1);++j)
{
cp t1=x[j],t2=omega*x[j+(len>>1)];
x[j]=t1+t2;
x[j+(len>>1)]=t1-t2;
omega*=omega_n;
}
}
}
if(flag==-1)
for(int i=0;i<n;++i)
x[i]=int(x[i].real()/n+0.5);
return;
}
int main()
{
n=read();m=read();
for(int i=0;i<=n;++i)x[i]=cp(read(),0);
for(int i=0;i<=m;++i)y[i]=cp(read(),0);
m+=n;n=1;
while(n<=m)n<<=1;
FFT(n,x,1);
FFT(n,y,1);
for(int i=0;i<n;++i)x[i]=x[i]*y[i];
FFT(n,x,-1);
for(int i=0;i<=m;++i)printf("%d ",(int)x[i].real());
puts("");
return 0;
}
作者

xqmmcqs

发布于

2018-03-02

更新于

2021-08-18

许可协议

评论

Your browser is out-of-date!

Update your browser to view this website correctly.&npsb;Update my browser now

×