浅谈 FFT （快速傅里叶变换）

发布于 6 天前  15 次阅读

直接开始吧。

前置芝士

  1. 简单数学（平面直角坐标系，多项式乘法）

  2. 码力（好吧其实并没有那么多代码要写）

  3. 简单分治（或者，用过线段树也可以）

数论的话，学 NTT 之前再看也可以。

如果你都会就直接开始了。

引入：多项式乘法

七年级课本数学，应该不用怎么讲吧。

例如：

(x2+5x+6)(3x3+x2−2x−4)(x^2+5x+6)(3x^3+x^2-2x-4)

拆括号：

3x5+x4−2x3−4x2+15x4+5x3−10x2−20x+18x3+6x2−12x−243x^5+x^4-2x^3-4x^2+15x^4+5x^3-10x^2-20x+18x^3+6x^2-12x-24

合并同类项：

3x5+16x4+21x3−8x2−6x−243x^5+16x^4+21x^3-8x^2-6x-24

上面的方法就相当于：

若 h(x)=f(x)×g(x)h(x)=f(x)\times g(x) ，则

hk=∑i+j=kfigjh_k=\sum_{i+j=k}f_ig_j

这样做复杂度是 O(n2)O(n^2) 的。

暴力代码：（模板题， 55pts55pts ）

#include<bits/stdc++.h>
using namespace std;
int n,m;
int f[1000005],g[1000005],h[1000005];
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++)scanf("%d",&f[i]);
    for(int i=0;i<=m;i++)scanf("%d",&g[i]);
    for(int i=0;i<=n;i++)
        for(int j=0;j<=m;j++)
            h[i+j]+=f[i]*g[j];
    for(int i=0;i<=n+m;i++)printf("%d ",h[i]);
    return 0;
}

点值转化

考虑优化这个暴力。

我们想，如果我们的乘法只要 O(n)O(n) ，那就很妙了。

根据拉格朗日插值，我们知道，在平面直角坐标系中， nn 个点就可以确定一个 n−1n-1 次多项式。

用 nn 个点来表示一个 n−1n-1 次多项式的方法，就是点值表示。

我们就会想到，能否将一个 n−1n-1 次多项式转化为 nn 个点呢？那么两个多项式相乘的点值表示，不就是这 2n2n 个点值的相乘吗？

例如，我们有两个多项式 f(x)=x2+5x+6f(x)=x^2+5x+6 和 g(x)=x2−2x−4g(x)=x^2-2x-4 。

那么 f(0)=6,f(1)=12,f(2)=20f(0)=6,f(1)=12,f(2)=20 （确定了 ff 这个多项式）

g(0)=−4,g(1)=−5,g(2)=−4g(0)=-4,g(1)=-5,g(2)=-4 （确定了 gg）

我们设 h(x)=f(x)×g(x)h(x)=f(x)\times g(x) 。

那么很显然， h(0)=−24,h(1)=−60,h(2)=−80h(0)=-24,h(1)=-60,h(2)=-80 。

好像有点不对啊， hh 的次数应该是 44 ，需要 55 个点啊。

那我们再找两个点不就好了……

好了现在我们知道 h(0)=−24,h(1)=−60,h(2)=−80,h(3)=−30,h(4)=168h(0)=-24,h(1)=-60,h(2)=-80,h(3)=-30,h(4)=168 。

根据拉格朗日插值公式， h(x)=∑0≤i≤4h(i)×∏1≤j≤4,i≠jx−ji−j=x4+3x3−8x2−32x−24h(x)=\sum_{0\le i\le 4}h(i)\times \prod_{1\le j\le 4,i\ne j}\dfrac{x-j}{i-j}=x^4+3x^3-8x^2-32x-24

嗯，就是这样。

对于两个次数为 nn 的多项式相乘，我们有一种策略：分别代入 2n+12n+1 个数，求出在这两个多项式上的点值，然后把这些点值相乘，最后再把点值转化为多项式。

我们把第一步——求点值叫做 DFT ，第二步——求多项式叫做 IDFT 。综合起来，就叫做 FFT 。

一个显而易见的事实是， IDFT 是 DFT 的CP逆运算。

有关复数

不会吧不会吧不会真的有人还没学会复数吧

一个复数的形式是 a+bia+bi ，这里 ii 是 −1-1 的平方根之一（另一个是 −i-i ，虽然他们两个哪个都不大于 00 （也不小于 00 （因为复数之间好像不能作大小比较来着）））。

复数的运算很简单，就跟多项式差不多（你可以认为 ii 是一个未知数）。

(a+bi)±(c+di)=(a±c)+(b±d)i(a+bi)\pm(c+di)=(a\pm c)+(b\pm d)i

(a+bi)(c+di)=(ac−bd)+(ad+bc)i(a+bi)(c+di)=(ac-bd)+(ad+bc)i

a+bic+di=(a+bi)(c−di)(c+di)(c−di)=(ac+bd)i+(bc−ad)ic2+d2\dfrac{a+bi}{c+di}=\dfrac{(a+bi)(c-di)}{(c+di)(c-di)}=\dfrac{(ac+bd)i+(bc-ad)i}{c^2+d^2}

而这些运算在几何上也有一些神奇的性质：

假设我们有一个复数 z=a+biz=a+bi ：

我们让它对应平面直角坐标系上的点 (a,b)(a,b) 。

定义一个复数的模长为 ∣z∣=a2+b2|z|=\sqrt{a^2+b^2} ，也就是它对应的点到原点 (0,0)(0,0) 的距离。

定义它的幅角为它对应的点到原点的这条线段与 xx 轴的夹角（带正负的角，如果不会自行百度），用式子就是 arctan⁡ba\arctan\dfrac{b}{a} 。

然后，我们看两个复数 z1=a+biz_1=a+bi 和 z2=c+diz_2=c+di 。

它们的乘积是 (ac−bd)+(ad+bc)i(ac-bd)+(ad+bc)i 。

我们把模长代进去： ∣z1×z2∣=(ac−bd)2+(ad+bc)2=a2c2+b2d2−2abcd+a2d2+b2c2+2abcd=a2c2+b2d2+a2d2+b2c2=(a2+b2)(c2+d2)=a2+b2×c2+d2=∣z1∣×∣z2∣|z_1\times z_2|=\sqrt{(ac-bd)^2+(ad+bc)^2}=\sqrt{a^2c^2+b^2d^2-2abcd+a^2d^2+b^2c^2+2abcd}=\sqrt{a^2c^2+b^2d^2+a^2d^2+b^2c^2}=\sqrt{(a^2+b^2)(c^2+d^2)}=\sqrt{a^2+b^2}\times\sqrt{c^2+d^2}=|z_1|\times|z_2|

把幅角的 tan⁡\tan 值代进去： ad+bcac−bd=ba+dc1−ba×dc\dfrac{ad+bc}{ac-bd}=\dfrac{\dfrac{b}{a}+\dfrac{d}{c}}{1-\dfrac{b}{a}\times\dfrac{d}{c}}

所以我们得到：复数相乘时，模长相乘，幅角相加。

单位根

学过三角函数的，一定都知道单位圆。

我们把单位圆搬到复平面上，那么在单位圆上的复数，模长一定为 11 。

Tips:Tips: 接下来讨论的复数，均是单位圆上的复数，也即模长为 11 的复数。

我们像切蛋糕一样，从 (1,0)(1,0) 开始，把这个圆切成 nn 块，切出的点就是 ωnk(0≤k<n)\omega_n^k(0\le k<n) 。

拿张图过来：

图中标出了 ω30,ω31,ω32\omega_3^0,\omega_3^1,\omega_3^2 。

显然， nn 次单位根是 xn=1x^n=1 的复数根。

因此 ωnk±n=ωnk\omega_n^{k\pm n}=\omega_n^k 。

同时，我们由图容易知道， ω2n2k=ωnk\omega_{2n}^{2k}=\omega_n^k （或者更一般地，上下标可以同时乘除任意一个正整数虽然没用）。

DFT 过程

我们把单位根代进多项式 F(x)=a0+a1x+⋯+an−1xn−1F(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1} ！

我们把 F(x)F(x) 拆成奇次项和偶次项，提出一个 xx ： F(x)=(a0+a2x2+⋯+an−2xn−2)+x(a1+a3x2+⋯+an−1xn−2)F(x)=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2}) ，设 FL(x)=(a0+a2x+⋯+an−2xn−22),FR(x)=x(a1+a3x+⋯+an−1xn−22)FL(x)=(a_0+a_2x+\cdots+a_{n-2}x^\frac{n-2}2),FR(x)=x(a_1+a_3x+\cdots+a_{n-1}x^\frac{n-2}2)

这里保证 nn 是 22 的幂，不会分得不均匀。

我们发现， FL(x)FL(x) 和 FR(x)FR(x) 都是关于 x2x^2 的多项式。

现在我们把单位根代进去。

当 2k<n2k<n 时，代入 ωnk\omega_n^k ：

F(ωnk)=FL(ωn2k)+ωnkFR(ωn2k)F(\omega_n^k)=FL(\omega_\frac{n}2^k)+\omega_n^kFR(\omega_\frac{n}2^k)

代入 ωnk+n2\omega_n^{k+\frac{n}2} ：

F(ωnk+n2)=FL(ωn2k+n2)+ωnk+n2FR(ωn2k+n2)=FL(ωn2k)−ωnkFR(ωn2k)F(\omega_n^{k+\frac{n}2})=FL(\omega_\frac{n}2^{k+\frac{n}2})+\omega_n^{k+\frac{n}2}FR(\omega_\frac{n}2^{k+\frac{n}2})=FL(\omega_\frac{n}2^k)-\omega_n^kFR(\omega_\frac{n}2^k)

哈！这两个等式不就只有一个正负号的区别吗？

辣么我们就可以用分治的思路，每次把 nn 减半，当然，到 n=1n=1 时，我们代入的只有 ω10=1\omega_1^0=1 ，此时什么也不要做。

IDFT 过程

上面说过， IDFT 是 DFT 的CP逆运算，它可以把点值表示转成系数表示。

其实它们就只有一句话的区别：只要把 DFT 中的 ω\omega 全部换成 ω−1\omega^{-1} ，做完之后除以 nn 即可。

证明只需证（使用矩阵求逆的证法）：

显然

我们假设

那么

验证：

所以 A×A−1A\times A^{-1} 确实是单位矩阵。

递归实现 FFT

#include<bits/stdc++.h>
using namespace std;
const double pi=acos(-1.0);
struct CP{
    double x,y;
    CP(double xx=0,double yy=0){
        x=xx,y=yy;
    }
    CP operator +(CP const &b)const{
        return CP(x+b.x,y+b.y);
    }
    CP operator -(CP const &b)const{
        return CP(x-b.x,y-b.y);
    }
    CP operator *(CP const &b)const{
        return CP(x*b.x-y*b.y,x*b.y+y*b.x);
    }
}f[2097160],g[2097160],sav[2097160];
int n,m;
void fft(CP *f,int len,bool flag){
    if(len==1)return;
    for(int k=0;k<len;k++)sav[k]=f[k];
    for(int k=0;k<len/2;k++)f[k]=sav[k<<1],f[k+len/2]=sav[k<<1|1];
    fft(f,len/2,flag);
    fft(f+len/2,len/2,flag);
    CP Ur(cos(2*pi/len),sin(2*pi/len));
    if(!flag)Ur.y*=-1;
    CP now(1,0);
    for(int k=0;k<len/2;k++){
        sav[k]=f[k]+now*f[k+len/2];
        sav[k+len/2]=f[k]-now*f[k+len/2];
        now=now*Ur;
    }
    for(int k=0;k<len;k++)f[k]=sav[k];
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++)scanf("%lf",&f[i].x);
    for(int i=0;i<=m;i++)scanf("%lf",&g[i].x);
    m+=n;
    n=1;
    while(n<=m)n<<=1;
    fft(f,n,1);fft(g,n,1);
    for(int i=0;i<n;i++)f[i]=f[i]*g[i];
    fft(f,n,0);
    for(int i=0;i<=m;i++)printf("%d ",(int)(f[i].x/n+0.499)); 
    return 0;
}

这样已经足以通过板子了，但众所周知递归非常慢，所以我们要卡卡常（

迭代实现 FFT

显然，递归可以直接改成迭代（手动从最后一层开始做）。可是我们没法知道哪个位置最后对应到哪个位置。

然后，我们发现：每个位置，经过原来的递归后，必定会到另一个位置等于没说。

关键在于：

0→00\to 0
1→41\to 4
2→22\to 2
3→63\to 6
4→14\to 1
5→55\to 5
6→36\to 3
7→77\to 7

哎！用二进制来表示，最后对应的位置就是原来位置的二进制反过来！

然后我们就可以递推每个数反过来是什么样子，然后就能做了。

#include<bits/stdc++.h>
using namespace std;
const double pi=acos(-1.0);
struct CP{
    CP (double xx=0,double yy=0){x=xx,y=yy;}
    double x,y;
    CP operator +(CP const &B)const{
        return CP(x+B.x,y+B.y); 
    }
    CP operator -(CP const &B)const{
        return CP(x-B.x,y-B.y); 
    }
    CP operator *(CP const &B)const{
        return CP(x*B.x-y*B.y,x*B.y+y*B.x); 
    }
}f[2097160],g[2097160];
int n,m;
int tr[2097160];
void fft(CP *f,int n,int flag){
    for(int i=0;i<n;i++)if(i<tr[i])swap(f[i],f[tr[i]]);
    for(int p=2;p<=n;p<<=1){
        int len=p>>1;
        CP Ur(cos(2*pi/p),sin(2*pi/p));
        if(!flag)Ur.y*=-1;
        for(int k=0;k<n;k+=p){
            CP now(1,0);
            for(int i=k;i<=k+len-1;i++){
                CP num=now*f[i+len];
                f[i+len]=f[i]-num;
                f[i]=f[i]+num;
                now=now*Ur;
            }
        }
    }
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++)scanf("%lf",&f[i].x);
    for(int i=0;i<=m;i++)scanf("%lf",&g[i].x);
    m+=n;
    n=1;
    while(n<=m)n<<=1;
    for(int i=0;i<n;i++)tr[i]=(tr[i>>1]>>1)|((i&1)?(n>>1):0);
    fft(f,n,1);fft(g,n,1);
    for(int i=0;i<n;i++)f[i]=f[i]*g[i];
    fft(f,n,0);
    for(int i=0;i<=m;i++)printf("%d ",(int)(f[i].x/n+0.499)); 
    return 0;
}

OIer and PVZer.（最近PVZ都不玩了

	</main>