FFT(快速傅里叶变换)

FFT F F T 在算法竞赛中的主要应用之一是加速多项式乘法的计算

信号处理和分析关系不大(表面上…)

多项式

系数表示

A(x)=i=0n1aixi=a0+a1x+a2x2+...+an1xn1 A ( x ) = ∑ i = 0 n − 1 a i x i = a 0 + a 1 x + a 2 x 2 + . . . + a n − 1 x n − 1

点值表示

对于上面的 n1 n − 1 次多项式,将一组互不相同的(n个) x x 带入进去得到对应的y,(即n个点)。

即这n个点可以唯一确定一个多项式

(由点值表示转系数表示,求法见多项式插值,朴素的插值算法时间复杂度为 O(n2) O ( n 2 )

多项式乘法

暂且称作“卷积”吧

有两个多项式, A(x)=n1i=0aixi A ( x ) = ∑ i = 0 n − 1 a i x i B(x)=n1i=0bixi B ( x ) = ∑ i = 0 n − 1 b i x i (假设两个多项式次数相同,若不同可在后面补零)


C(x)=A(x)B(x)=k=02n2(k=i+jaibj)xk C ( x ) = A ( x ) ∗ B ( x ) = ∑ k = 0 2 n − 2 ( ∑ k = i + j a i b j ) x k

两个 n1 n − 1 次多项式相乘,得到的是一个 2n2 2 n − 2 次多项式,时间复杂度为 O(n2) O ( n 2 )


另外,也可以用点值表达,即选择 2n1 2 n − 1 个互不相同的 xi x i 带入 A(x) A ( x ) B(x) B ( x ) 相乘,得到 2n1 2 n − 1

值,这 2n1 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 次单位根。

由复数乘法的定义(模长相乘,幅角相加)可知,其余的 n1 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 ,这为我们提供了一个计算单位根及其幂的公式

ωkn=cosk2πn+isink2πn ω n k = c o s k 2 π n + i s i n k 2 π 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到 n1 n − 1 次幂(共n个)带入多项式的系数表示,所得点值向量

(A(ω0n),A(ω1n),A(ω2n),...,A(ωn1n)) ( A ( ω n 0 ) , A ( ω n 1 ) , A ( ω n 2 ) , . . . , A ( ω n n − 1 ) )

(变换后的这n个点对,可以唯一确定这个多项式!)



如果从信号的角度看,这个过程就是时域到频域的转换, 即选择一段时间,取样n个点,即0时刻,1时刻,…,N-1时刻。转换的思路是这样的:我去“枚举”一些频率,对于一个固定的频率 f f ,去对时域图像(假设为x(t))做“缠绕操作”,而缠绕的时间取样就是上面的N个时刻,即X(f)=n=0N1x(n)wNfn (n是每个时间点)。简单来说,缠绕操作就是对一个无明显频率规律的时域图像,看它在指定(枚举)频率上是否有(这个频率),衡量的指标就是 X X 。(至于缠绕操作,就是一种方便的工具,使得x(n)分布在特定的角度上)。所以我们要枚举频率,一般也取 0N1 0 ∼ N − 1

我们对 X(f)=N1n=0x(n)wfnN X ( f ) = ∑ n = 0 N − 1 x ( n ) w N f n 简单变个形:

X(f)=n=0N1x(n)(wfN)n X ( f ) = ∑ n = 0 N − 1 x ( n ) ( w N f ) n

再看下朴素的多项式
A(x)=i=0n1aixi=a0+a1x+a2x2+...+an1xn1 A ( x ) = ∑ i = 0 n − 1 a i x i = a 0 + a 1 x + a 2 x 2 + . . . + a n − 1 x n − 1

f f 分别取0,1,2,...,n1,就是多项式 A(x) A ( x ) 带入 x=w0n,w1n,...,wn1n x = w n 0 , w n 1 , . . . , w n n − 1

而多项式的系数 a0,a1,...,an1 a 0 , a 1 , . . . , a n − 1 就对应 x(0),x(1),...,x(n1) x ( 0 ) , x ( 1 ) , . . . , x ( n − 1 ) (即信号在时域离散时间点上的n个强度值)

进一步可以参考

傅里叶级数与傅里叶变换

傅里叶变换和离散傅里叶变换



称这个结果 (A(ω0n),A(ω1n),A(ω2n),...,A(ωn1n)) ( A ( ω n 0 ) , A ( ω n 1 ) , A ( ω n 2 ) , . . . , A ( ω n n − 1 ) ) 为其系数向量 (a0,a1,a2,...,an1) ( a 0 , a 1 , a 2 , . . . , a n − 1 ) 离散傅里叶变换

按照朴素方法来求解原系数的离散傅里叶变换,时间复杂度仍为 O(n2) O ( n 2 )

而这个过程可以分治进行,因而可以优化,即FFT,这是算法竞赛的重点(但此时还不知道这n个点对求多项式乘法有什么用,继续看

至于时域频域的转换,只要理解思想就行了

快速傅里叶变换(FFT)

考虑将多项式按照系数下标的奇偶分为两部分

A(x)=(a0+a2x2+a4x4+...+an2xn2)+(a1x+a3x3+a5x5+...+an1xn1) A ( x ) = ( a 0 + a 2 x 2 + a 4 x 4 + . . . + a n − 2 x n − 2 ) + ( a 1 x + a 3 x 3 + a 5 x 5 + . . . + a n − 1 x n − 1 )


A1(x)=a0+a2x+a4x2+...+an2xn21A2(x)=a1+a3x+a5x2+...+an1xn21 A 1 ( x ) = a 0 + a 2 x + a 4 x 2 + . . . + a n − 2 x n 2 − 1 A 2 ( x ) = a 1 + a 3 x + a 5 x 2 + . . . + a n − 1 x n 2 − 1

则有
A(x)=A1(x2)+xA2(x2) A ( x ) = A 1 ( x 2 ) + x A 2 ( x 2 )

假设 k<n2 k < n 2 ,现在要求 A(ωkn) A ( ω n k )
A(ωkn)=A1(ωkn2)+ωknA2(ωkn2) A ( ω n k ) = A 1 ( ω n 2 k ) + ω n k A 2 ( ω n 2 k )

(使用了单位根的性质1)

由于前面 k<n2 k < n 2 ,对于 A(ωk+n2n) A ( ω n k + n 2 ) 的部分:

A(ωk+n2n)=A(wkn)=A1(ωkn2)ωknA2(ωkn2) A ( ω n k + n 2 ) = A ( − w n k ) = A 1 ( ω n 2 k ) − ω n k A 2 ( ω n 2 k )

(使用了单位根的性质1和性质2,这个指数范围的变化是精华)

k k 取遍 [0,n21] 时, k k k+n2取遍了 [0,n1] [ 0 , n − 1 ] ,即所求

那么,如果已知 A1(x) A 1 ( x ) A2(x) A 2 ( x ) ω0n2,ω1n2,...,ωn21n2 ω n 2 0 , ω n 2 1 , . . . , ω n 2 n 2 − 1 的值,就可以在 O(n) O ( n ) 时间内求得 A(x) A ( x ) ω0n,ω1n,ω2n,...,ωn1n ω 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 ) 的系数只有一项。

则该分治算法的时间复杂度为

T(n)=2T(n/2)+O(n)=O(nlogn) T ( n ) = 2 ∗ T ( n / 2 ) + O ( n ) = O ( n l o g n )

上述过程称为 Cooley-Tukey 算法(JW Cooley 和 John Tukey)

现在我们可以在 O(nlogn) O ( n l o g n ) 时间内求得 n个点对辣!

(ω0n,ω1n,ω2n,...,ωn1n) ( ω n 0 , ω n 1 , ω n 2 , . . . , ω n n − 1 ) (A(ω0n),A(ω1n),A(ω2n),...,A(ωn1n)) ( 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,...,yn1) ( y 0 , y 1 , y 2 , . . . , y n − 1 ) (a0,a1,a2,...,an1) ( a 0 , a 1 , a 2 , . . . , a n − 1 ) 的离散傅里叶变换。考虑另一个向量 (C0,C1,C2,...,Cn1) ( C 0 , C 1 , C 2 , . . . , C n − 1 ) ,满足

Ck=i=0n1yi(ωkn)i,k[0n1] C k = ∑ i = 0 n − 1 y i ( ω n − k ) i , k ∈ [ 0 , n − 1 ]

即多项式 B(x)=y0+y1x+y2x2+...+yn1xn1 B ( x ) = y 0 + y 1 x + y 2 x 2 + . . . + y n − 1 x n − 1 ω0n,ω1n,ω2n,...,ω(n1)n ω n 0 , ω n − 1 , ω n − 2 , . . . , ω n − ( n − 1 ) 处的点值表示

(B(x)以原系数的离散傅里叶变换作为新的系数)

Ck C k 展开

Ck=i=0n1yi(ωkn)i=i=0n1(j=0n1aj(ωin)j)(ωkn)i=i=0n1(j=0n1aj(ωjn)i)(ωkn)i=i=0n1(j=0n1aj(ωjn)i(ωkn)i)=i=0n1(j=0n1aj(ωjkn)i)=j=0n1aj(i=0n1(ωjkn)i) C k = ∑ i = 0 n − 1 y i ( ω n − k ) i = ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( ω n i ) j ) ( ω n − k ) i = ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( ω n j ) i ) ( ω n − k ) i = ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( ω n j ) i ( ω n − k ) i ) = ∑ i = 0 n − 1 ( ∑ j = 0 n − 1 a j ( ω n j − k ) i ) = ∑ j = 0 n − 1 a j ( ∑ i = 0 n − 1 ( ω n j − k ) i )

观察内层求和,为了解决这个求和式,先考虑一个式子先,令 (首项为1,公比为 ωkn ω n k 的等比数列前n项和为 S(ωkn) S ( ω n k )
S(ωkn)=1+ωkn+(ωkn)2+...+(ωkn)n1 S ( ω n k ) = 1 + ω n k + ( ω n k ) 2 + . . . + ( ω n k ) n − 1

k0 k ≠ 0 时,两边同时乘上 ωkn ω n k ,得
ωknS(ωkn)=ωkn+(ωkn)2+(ωkn)3+...+(ωkn)n ω n k S ( ω n k ) = ω n k + ( ω n k ) 2 + ( ω n k ) 3 + . . . + ( ω n k ) n

两式相减,整理后得
ωknS(ωkn)S(ωkn)=(ωkn)n1S(ωkn)=(ωkn)n1ωkn1 ω n k S ( ω n k ) − S ( ω n k ) = ( ω n k ) n − 1 即 S ( ω n k ) = ( ω n k ) n − 1 ω n k − 1

这个式子分子为0,分母一定不为0,因此
S(ωkn)=0 S ( ω n k ) = 0

当k=0时, S(ωkn)=n S ( ω n k ) = n

继续考虑 Ck C k

Ck=j=0n1aj(i=0n1(ωjkn)i)=j=0n1ajS(ωjkn) C k = ∑ j = 0 n − 1 a j ( ∑ i = 0 n − 1 ( ω n j − k ) i ) = ∑ j = 0 n − 1 a j S ( ω n j − k )

对每一个k,枚举j,只有当 j=k j = k 时, S(ωjkn)=n S ( ω n j − k ) = n ,否则 S(ωjkn)=0 S ( ω n j − k ) = 0 ,即
Ci=naiai=1nCi C i = n a i 即 a i = 1 n C i

使用单位根的倒数代替单位根原系数的离散傅里叶变换结果作为新的系数,做一次快速傅里叶变换的过程,再 将结果每个数除以 n,即为傅里叶逆变换的结果,即原系数 ai a i

这样就可以先做一次正变换再做一次逆变换求得系数啦

简单来说,求解多项式乘法的大致思路就是:

  • 确定结果多项式 C(x) C ( x ) 的次数,决定要取的点对的数量 N N N为2的次幂)
  • A(x) A ( x ) B(x) B ( x ) 在这 N N 个点的值 O(nlogn)
  • 做傅里叶逆变换求得系数值

小结&&实现细节

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 为结果多项式)

则有DFT(a()b)=DFT(a)DFT(b) 右边为点乘(很自然)

所以所求 a()b=DFT1(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

img

#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 ] 中,合并的操作可以表示为:

b[k]a[k]+wkna[n2+k]b[k+n2]a[k]wkna[n2+k] b [ k ] ← a [ k ] + w n k ∗ a [ n 2 + k ] b [ k + n 2 ] ← a [ k ] − w n k ∗ a [ n 2 + k ]

考虑加入一个临时变量t,使得这个过程可以在原地完成,而不需要数组b,
twkna[n2+k]a[k+n2]a[k]ta[k]a[k]+t t ← w n k ∗ a [ n 2 + k ] a [ k + n 2 ] ← a [ k ] − t a [ k ] ← a [ k ] + t

由于 k k k+n2 是对应的,所以不同的k之间不会相互影响

这一过程被称为蝴蝶操作

使用C++已经封装的Complex的迭代实现

img

#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;
}

推荐参考资料

https://oi.men.ci/fft-notes/

http://hzwer.com/3668.html

https://blog.csdn.net/iamzky/article/details/22712347

  • 3
    点赞
  • 14
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值