FFT F F T 在算法竞赛中的主要应用之一是加速多项式乘法的计算
和信号处理和分析关系不大(表面上…)
多项式
系数表示
点值表示
对于上面的 n−1 n − 1 次多项式,将一组互不相同的(n个) x x 带入进去得到对应的,(即n个点)。
即这n个点可以唯一确定一个多项式。
(由点值表示转系数表示,求法见多项式插值,朴素的插值算法时间复杂度为 O(n2) O ( n 2 ) )
多项式乘法
暂且称作“卷积”吧
有两个多项式, A(x)=∑n−1i=0aixi A ( x ) = ∑ i = 0 n − 1 a i x i 和 B(x)=∑n−1i=0bixi B ( x ) = ∑ i = 0 n − 1 b i x i (假设两个多项式次数相同,若不同可在后面补零)
则
两个 n−1 n − 1 次多项式相乘,得到的是一个 2n−2 2 n − 2 次多项式,时间复杂度为 O(n2) O ( n 2 ) 。
另外,也可以用点值表达,即选择 2n−1 2 n − 1 个互不相同的 xi x i 带入 A(x) A ( x ) 和 B(x) B ( x ) 相乘,得到 2n−1 2 n − 1 个
值,这 2n−1 2 n − 1 个点值就唯一确定了这个多项式,时间复杂度 O(n) O ( n ) (注意只是得到点值表达式)
复数
设 a、b 为实数, i2=−1 i 2 = − 1 ,形如 a+bi a + b i 的数叫做复数,其中 i i 被称为虚数单位。复数域是已知最大的域。
复平面
在复平面中,x 轴代表实数、y 轴(除原点外的所有点)代表虚数。每一个复数 a+bi a + b i 对应复平面上一个从 (0, 0) 指向 (a, b) 的向量。
该向量的长度 a2+b2−−−−−−√ a 2 + b 2 叫做模长。
从 x 轴正半轴到该向量的转角的有向(以逆时针为正方向)角叫做幅角。
复数相加遵循平行四边形定则。
复数相乘时,模长相乘,幅角相加。
单位根
下文中,如不特殊指明,均取 n为 2 的正整数次幂。
在复平面上,以原点为圆心,1为半径作圆,所得的圆叫做单位圆。以原点为起点,单位圆的 n 等分点为终点,作 n个向量。设所得的幅角为正且最小的向量对应的复数为 ωn ω n ,称为 n 次单位根。
由复数乘法的定义(模长相乘,幅角相加)可知,其余的 n−1 n − 1 个向量对应的复数分别为 ω2n,ω3n,...,ωnn ω n 2 , ω n 3 , . . . , ω n n ,其中 ωnn=ω0n=1 ,ω1n=ωn ω n n = ω n 0 = 1 , ω n 1 = ω n
单位根的幅角为圆周角的
1n
1
n
,这为我们提供了一个计算单位根及其幂的公式
单位根的性质
性质1 ω2k2n=ωkn ω 2 n 2 k = ω n k 显然
性质2 ωk+n2n=−ωkn ω n k + n 2 = − ω n k 因为 ωn2n=−1 ω n n 2 = − 1 (带入即可验证)
离散傅里叶变换(DFT)
对于一个n个系数,n-1次的多项式
考虑多项式
A(x)
A
(
x
)
的表示。将n次单位根的0到
n−1
n
−
1
次幂(共n个)带入多项式的系数表示,所得点值向量
(变换后的这n个点对,可以唯一确定这个多项式!)
如果从信号的角度看,这个过程就是时域到频域的转换, 即选择一段时间,取样n个点,即0时刻,1时刻,…,N-1时刻。转换的思路是这样的:我去“枚举”一些频率,对于一个固定的频率 f f ,去对时域图像(假设为x(t))做“缠绕操作”,而缠绕的时间取样就是上面的N个时刻,即 (n是每个时间点)。简单来说,缠绕操作就是对一个无明显频率规律的时域图像,看它在指定(枚举)频率上是否有(这个频率),衡量的指标就是 X X 。(至于缠绕操作,就是一种方便的工具,使得分布在特定的角度上)。所以我们要枚举频率,一般也取 0∼N−1 0 ∼ N − 1 。
我们对
X(f)=∑N−1n=0x(n)wfnN
X
(
f
)
=
∑
n
=
0
N
−
1
x
(
n
)
w
N
f
n
简单变个形:
再看下朴素的多项式
即 f f 分别取,就是多项式 A(x) A ( x ) 带入 x=w0n,w1n,...,wn−1n x = w n 0 , w n 1 , . . . , w n n − 1
而多项式的系数 a0,a1,...,an−1 a 0 , a 1 , . . . , a n − 1 就对应 x(0),x(1),...,x(n−1) x ( 0 ) , x ( 1 ) , . . . , x ( n − 1 ) (即信号在时域离散时间点上的n个强度值)
进一步可以参考
称这个结果 (A(ω0n),A(ω1n),A(ω2n),...,A(ωn−1n)) ( A ( ω n 0 ) , A ( ω n 1 ) , A ( ω n 2 ) , . . . , A ( ω n n − 1 ) ) 为其系数向量 (a0,a1,a2,...,an−1) ( a 0 , a 1 , a 2 , . . . , a n − 1 ) 的离散傅里叶变换
按照朴素方法来求解原系数的离散傅里叶变换,时间复杂度仍为 O(n2) O ( n 2 )
而这个过程可以分治进行,因而可以优化,即FFT,这是算法竞赛的重点(但此时还不知道这n个点对求多项式乘法有什么用,继续看)
至于时域频域的转换,只要理解思想就行了
快速傅里叶变换(FFT)
考虑将多项式按照系数下标的奇偶分为两部分
令
则有
假设 k<n2 k < n 2 ,现在要求 A(ωkn) A ( ω n k )
(使用了单位根的性质1)
由于前面
k<n2
k
<
n
2
,对于
A(ωk+n2n)
A
(
ω
n
k
+
n
2
)
的部分:
(使用了单位根的性质1和性质2,这个指数范围的变化是精华)
当 k k 取遍 时, k k 和 取遍了 [0,n−1] [ 0 , n − 1 ] ,即所求
那么,如果已知 A1(x) A 1 ( x ) 和 A2(x) A 2 ( x ) 在 ω0n2,ω1n2,...,ωn2−1n2 ω n 2 0 , ω n 2 1 , . . . , ω n 2 n 2 − 1 的值,就可以在 O(n) O ( n ) 时间内求得 A(x) A ( x ) 在 ω0n,ω1n,ω2n,...,ωn−1n ω n 0 , ω n 1 , ω n 2 , . . . , ω n n − 1 处的取值。而关于 A1(x),A2(x) A 1 ( x ) , A 2 ( x ) 的问题正好是原问题规模缩小一半的子问题,分治的边界为一个常数项 a0 a 0 ,即 Aϕ(x) A ϕ ( x ) 的系数只有一项。
则该分治算法的时间复杂度为
上述过程称为 Cooley-Tukey 算法(JW Cooley 和 John Tukey)。
现在我们可以在 O(nlogn) O ( n l o g n ) 时间内求得 n个点对辣!
即 (ω0n,ω1n,ω2n,...,ωn−1n) ( ω n 0 , ω n 1 , ω n 2 , . . . , ω n n − 1 ) 和 (A(ω0n),A(ω1n),A(ω2n),...,A(ωn−1n)) ( A ( ω n 0 ) , A ( ω n 1 ) , A ( ω n 2 ) , . . . , A ( ω n n − 1 ) )
然而还是不知道和多项式乘法有什么关系……
事实上,这n对取值,将作为“原料”,帮助我们反过来去求多项式的系数
和一般的插值不同( O(n2) O ( n 2 ) )
这些特殊的点对可以在 O(nlogn) O ( n l o g n ) 内求出多项式的系数
即下面的傅里叶逆变换
仿佛离最终目标进了一大步……
傅里叶逆变换
将点值表示的多项式转化为系数表示,同样可以使用快速傅里叶变换,这个过程叫做傅里叶逆变换。
设
(y0,y1,y2,...,yn−1)
(
y
0
,
y
1
,
y
2
,
.
.
.
,
y
n
−
1
)
为
(a0,a1,a2,...,an−1)
(
a
0
,
a
1
,
a
2
,
.
.
.
,
a
n
−
1
)
的离散傅里叶变换。考虑另一个向量
(C0,C1,C2,...,Cn−1)
(
C
0
,
C
1
,
C
2
,
.
.
.
,
C
n
−
1
)
,满足
即多项式 B(x)=y0+y1x+y2x2+...+yn−1xn−1 B ( x ) = y 0 + y 1 x + y 2 x 2 + . . . + y n − 1 x n − 1 在 ω0n,ω−1n,ω−2n,...,ω−(n−1)n ω n 0 , ω n − 1 , ω n − 2 , . . . , ω n − ( n − 1 ) 处的点值表示
(B(x)以原系数的离散傅里叶变换作为新的系数)
将
Ck
C
k
展开
观察内层求和,为了解决这个求和式,先考虑一个式子先,令 (首项为1,公比为 ωkn ω n k 的等比数列前n项和为 S(ωkn) S ( ω n k ) )
①当 k≠0 k ≠ 0 时,两边同时乘上 ωkn ω n k ,得
两式相减,整理后得
这个式子分子为0,分母一定不为0,因此
②当k=0时, S(ωkn)=n S ( ω n k ) = n
继续考虑
Ck
C
k
,
对每一个k,枚举j,只有当 j=k j = k 时, S(ωj−kn)=n S ( ω n j − k ) = n ,否则 S(ωj−kn)=0 S ( ω n j − k ) = 0 ,即
即 使用单位根的倒数代替单位根, 原系数的离散傅里叶变换结果作为新的系数,做一次快速傅里叶变换的过程,再 将结果每个数除以 n,即为傅里叶逆变换的结果,即原系数 ai a i 。
这样就可以先做一次正变换再做一次逆变换求得系数啦
简单来说,求解多项式乘法的大致思路就是:
- 确定结果多项式 C(x) C ( x ) 的次数,决定要取的点对的数量 N N (为2的次幂)
- 求 A(x) A ( x ) 和 B(x) B ( x ) 在这 N N 个点的值
- 做傅里叶逆变换求得系数值
小结&&实现细节
DFT(离散傅里叶变换)是一种对n个元素的数组的变换,根据式子直接的方法是 O(n2) O ( n 2 ) 的。但是用分治的方法是 O(nlogn) O ( n l o g n ) ,即FFT(快速傅里叶变换)。
由于DFT满足cyclic convolution(循环卷积)的性质,即
定义 h=a(∗)b h = a ( ∗ ) b 为 hr=∑x+y=raxby h r = ∑ x + y = r a x b y (多项式乘法)( h h 为结果多项式)
则有 右边为点乘(很自然)
所以所求 a(∗)b=DFT−1(DFT(a)⋅DFT(b)) a ( ∗ ) b = D F T − 1 ( D F T ( a ) ⋅ D F T ( b ) )
即只要对a,b分别进行DFT变化之后点乘之后再逆变换就可以了
- 由于FFT本身算法的要求,n是2的次幂,注意补0
- DFT是定义在复数域上的,有与整数之间变换的要求
- 不要忘记最后除n
- C++ 的 STL 在头文件
complex
中提供一个复数的模板实现std::complex<T>
,其中T
为实数类型,一般取double
,在对精度要求较高的时候可以使用long double
或__float128
(不常用)。 - 单位根的倒数等于其共轭复数,使用
std::conj()
取得 IDFT 所需的单位根的倒数。 - 对时间要求苛刻时,可以自己实现复数类以提高速度
FFT的递归实现
直接按照上述思路实现即可
使用C++已经封装的Complex
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int debug_num=0;
#define debug cout<<"debug "<<++debug_num<<" :"
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
#define bit(a,b) ((a>>b)&1)
const int inf = 0x3f3f3f3f;
const ll inff = 0x3f3f3f3f3f3f3f3f;
const double pi=acos(-1.0);
const int maxn=1e5;
typedef complex<double> xu;
xu a[maxn<<2],b[maxn<<2];
int inver;
void FFT(xu *s ,int n){
if(n==1) return ;//n=1时值就是系数 因为只有w_1^0这一点
xu a1[n>>1],a2[n>>1];
for(int i=0;i<n;i+=2){
a1[i>>1]=s[i];
a2[i>>1]=s[i+1];
}
FFT(a1,n>>1); FFT(a2,n>>1);
xu w=xu(cos(2*pi/n),inver*sin(2*pi/n));
xu wn=xu(1,0);
for(int i=0;i<(n>>1);++i,wn=wn*w){
s[i]=a1[i]+wn*a2[i];
s[i+(n>>1)]=a1[i]-wn*a2[i];
}
}
int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt","r",stdin);
#endif
ios::sync_with_stdio(false);
int n,m;
cin>>n>>m;
int tp;
n++;
m++;
for(int i=0;i<n;++i){
cin>>tp;
a[i]=xu(tp,0);
}
for(int i=0;i<m;++i){
cin>>tp;
b[i]=xu(tp,0);
}
int N=1; while(N<n+m-1) N<<=1;
inver=1;
FFT(a,N); FFT(b,N);
for(int i=0;i<N;++i) a[i]=a[i]*b[i];
//现在a中的值是上面的傅里叶变化的结果啦 wn0,wn1,wn2,...,wn n-1
//下面作为逆变化的系数
inver=-1;
FFT(a,N);
for(int i=0;i<n+m-1;++i) cout<<ll(a[i].real()/N+0.5)<<' ';
return 0;
}
FFT的迭代实现
递归实现的 FFT 效率不高,因为有栈的调用和参数的传递,实际中一般采用迭代实现
二进制位翻转
对分治规律的总结
//以8项为例
000 001 010 011 100 101 110 111
0 1 2 3 4 5 6 7 //初始要求
0 2 4 6 / 1 3 5 7 //分治后的两个多项式的系数 对应于原多项式
0 4 / 2 6 / 1 5 / 3 7
0 / 4 / 2 / 6 / 1 / 5 / 3 / 7
000 100 010 110 001 101 011 111 //发现正好是翻转
这为迭代实现提供了理论基础
蝴蝶操作
考虑合并两个子问题的过程,假设
A1(wkn2)
A
1
(
w
n
2
k
)
和
A2(wkn2)
A
2
(
w
n
2
k
)
分别存放在
a[k]
a
[
k
]
和
a[n2+k]
a
[
n
2
+
k
]
中,
A(wkn)
A
(
w
n
k
)
和
A(wk+n2n)
A
(
w
n
k
+
n
2
)
将要存放在
b[k]
b
[
k
]
和
b[n2+k]
b
[
n
2
+
k
]
中,合并的操作可以表示为:
考虑加入一个临时变量t,使得这个过程可以在原地完成,而不需要数组b,
由于 k k 和 是对应的,所以不同的k之间不会相互影响
这一过程被称为蝴蝶操作
使用C++已经封装的Complex的迭代实现
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int debug_num=0;
#define debug cout<<"debug "<<++debug_num<<" :"
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
#define bit(a,b) ((a>>b)&1)
const int inf = 0x3f3f3f3f;
const ll inff = 0x3f3f3f3f3f3f3f3f;
const double pi=acos(-1.0);
const int maxn=1e5;
typedef complex<double> xu;
xu omega[maxn<<2],omegaInverse[maxn<<2];//辅助
void init(const int n){
for(int i=0;i<n;++i){
omega[i]=xu(cos(2*pi/n*i),sin(2*pi/n*i));
omegaInverse[i]=conj(omega[i]);
}
}
void trans(xu *a,const int n,const xu *omega){//a是系数 n是2的次幂 omega是要带入的点
int k=0;
while((1<<k)<n) k++;
//看N是2的多少次幂
for(int i=0;i<n;++i){
int t=0;
for(int j=0;j<k;++j) if(i&(1<<j)) t|=(1<<(k-j-1));
if(i<t) swap(a[i],a[t]);//原地翻转
}
for(int l=2;l<=n;l*=2){//l=1不用求 就是系数
int m=l/2;
for(xu *p=a;p!=a+n;p+=l){//l为一段要求 由两段l/2合并得到
//当前处理[p,p+l)
for(int i=0;i<m;++i){
//蝴蝶操作 omega_{l}^{i}=omega_{N}^{N/l*i}
xu t=omega[n/l*i]*p[i+m];
p[i+m]=p[i]-t;
p[i]+=t;
}
}
}
}
void dft(xu *a,const int n){
trans(a,n,omega);
}
void idft(xu *a,const int n){
trans(a,n,omegaInverse);
for(int i=0;i<n;++i) a[i]/=n;
}
xu a[maxn<<2];
xu b[maxn<<2];
int main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt","r",stdin);
#endif
ios::sync_with_stdio(false);
int n,m;
cin>>n>>m;
n++;
m++;
int tp;
for(int i=0;i<n;++i){
cin>>tp;
a[i]=xu(tp,0);
}
for(int i=0;i<m;++i){
cin>>tp;
b[i]=xu(tp,0);
}
int N=1;
while(N<m+n-1) N<<=1;
init(N);
dft(a,N);
dft(b,N);
for(int i=0;i<N;++i) a[i]*=b[i];
idft(a,N);
for(int i=0;i<n+m-1;++i){
cout<<(ll)(a[i].real()+0.5)<<" ";
}
return 0;
}
推荐参考资料