快速傅里叶变换 FFT

更新日志
2025/01/08：开工。
2025/01/09：添加整合模板，包括加减操作。优化代码细节。

2025/06/24：更新模板以适配缺省源。

2025/06/25：优化模板细节。提供单次 DFT 优化。

前言

文中所有代码块，只有最后的整合封装模板是我实时更新的，其余代码仅供方便理解，实现可能不优雅或更劣，当然正确性没有问题。

整合封装模板可能需要搭配缺省源使用，可以去我首页置顶的缺省源归档复制一份。

简介

\(\text{FFT}\)，用于优化多项式乘法。

所以，我们将会从多项式乘法开始，然后再进入优化过程。

鉴于优化中会使用到复数，所以作为前置知识进行讲解。

多项式

系数表示法

令 \(A(i)\) 表示一个 \(n\) 次多项式，则 \(A(i)=\sum\limits_{i=0}^{n}a_ix^i\)。

众所周知的，这样计算多项式乘法的复杂度为 \(O(n^2)\)。

所以，我们有必要对其进行优化！

点值表示法

根据某一定理，我们可以代入 \(n\) 个值，得到 \(n\) 组数对 \((x,y)\)，唯一确定一个 \(n-1\) 次多项式。

两个点值表示法表示的多项式相乘是 \(O(n)\) 的，具体的，\(C(x_i)=A(x_i)\times B(x_i)\)

但是，转换系数表示法和点值表示法的复杂度还是 \(O(n^2)\) 的。所以，我们就可以从这里下手优化了！

复数

概念与表示

一个复数可以用下列格式表示：

\[a+bi \]

其中 \(i\) 为虚数单位，即为 \(\sqrt{-1}\)。

表现在复平面上，横坐标为实部，纵轴为虚部。这样，我们就可以把复数表现在一个平面中，表现为一个向量。

运算法则

加减法

和向量一样，满足平行四边形法则。实部与虚部分别加减即可。

乘法

表现在几何上，模长相乘，幅角相加。描述的简单一些，到原点的距离相乘，与正横轴的夹角相加（任意角）。

在代数中表示，可得：

\[\begin{align*} &(a+bi)(c+di)\\ =&ac+adi+bci+bdi^2\\ =&ac+adi+bci-bd\\ =&(ac-bd)+(bc+ad)i \end{align*} \]

单位根

在复平面上，我们将单位圆 \(n\) 等分，起始点位于 \((1,0)\)。

我们从起始点开始，从 \(0\) 开始逆时针标号至 \(n-1\)，将编号为 \(k\) 分割点记作 \(\omega_n^k\)。

根据上面的乘法几何定义，我们可以得到：

\[\omega_n^k=(\omega_n^1)^k \]

简单解释的话，模长都是 \(1\)，\(\omega_n^k\) 的角度为 \(\omega_n^1\) 的 \(k\) 倍。

同时，我们可以通过三角函数得到它们的值：

\[\omega_n^k=\cos(\frac{k}{n}2\pi)+\sin(\frac{k}{n}2\pi)i \]

借助几何图像，我们可以简单得到一些性质（\(n\) 为 \(2\) 的正整数次幂）：

\[\omega_n^k=\omega_{2n}^{2k}\\ \omega_n^{k+\frac{n}{2}}=-\omega_n^k\\ \omega_n^0=\omega_n^n=1 \]

代码

如果不想用 \(\text{STL}\) 的话，可以用下方模板：

struct clx{
    flt x,y;
    clx(flt a=0,flt b=0){x=a,y=b;}
    clx operator+(clx b){return clx(x+b.x,y+b.y);}
    clx operator-(clx b){return clx(x-b.x,y-b.y);}
    clx operator*(clx b){return clx(x*b.x-y*b.y,x*b.y+y*b.x);}
};

FFT

快速傅里叶变换 DFT

首先，我们强制规定多项式系数为 \(2^k,k\in Z^*\)。如果不满足，在多项式后面补系数为 \(0\) 的项即可。

我们将代入的值固定为 \(w_n^0,w_n^1,\dots,w_n^{n-1}\)。同时，我们设 \(A(x)\) 的系数为 \(a_0,a_1,\dots,a_{n-1}\)。

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

带入之后，我们将得到 \(A(x)\) 的离散傅里叶变换 \(y_0,y_1,\dots,y_{n-1}\)，其中 \(y_i=A(\omega_n^i)\)。

我们新建两个多项式 \(A_1,A_2\)：

\[A_1(x)=a_0x^0+a_2x^1+\dots+a_{n-2}x^{\frac{n}{2}-1}\\ A_2(x)=a_1x^0+a_3x^1+\dots+a_{n-1}x^{\frac{n}{2}-1} \]

可得：

\[A(x)=A_1(x^2)+xA_2(x^2) \]

下面我们代入 \(\omega_n^k(k< \frac{n}{2})\)：

\[\begin{align*} A(\omega_n^k)=&A_1(\omega_n^{2k})+\omega_n^{k}A_2(\omega_n^{2k})\\ =&A_1(\omega_\frac{n}{2}^k)+\omega_n^{k}A_2(\omega_\frac{n}{2}^k) \end{align*} \]

下面我们代入 \(\omega_n^{k+\frac{n}{2}}(k< \frac{n}{2})\)：

\[\begin{align*} A(\omega_n^{k+\frac{n}{2}})=&A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})\\ =&A_1(\omega_\frac{n}{2}^k\omega_n^n)+\omega_n^{k+\frac{n}{2}}A_2(\omega_\frac{n}{2}^k\omega_n^n)\\ =&A_1(\omega_\frac{n}{2}^k)-\omega_n^kA_2(\omega_\frac{n}{2}^k) \end{align*} \]

不难发现，二者只有中间的常数不同，也就是说，我们可以在求前半部分的时候，\(O(1)\) 求出右半部分！

容易看出这是一种分治思想，且复杂度为 \(O(n\log n)\)。

当多项式只剩下一个常数项时，就可以返回了。

快速傅里叶逆变换 IDFT

别忘了我们还需要把点值表示法转换回系数表示法。

设 \((y_0,y_1,\dots,y_{n-1})\) 为多项式 \(A(x)\) 的离散傅里叶变换，我们以它们为一个新多项式的系数，即为：

\[B(x)=\sum_{i=0}^{n-1}y_ix^i \]

然后我们把上面那些单位根的倒数代入 \(B(x)\)，也就是代入 \(\omega_n^{0},\omega_n^{-1},\dots,\omega_n^{-(n-1)}\)，得到新的离散傅里叶变换为：

\[\begin{align*} z_k=&B(\omega_n^{-k})\\ =&\sum_{i=0}^{n-1}y_i\omega_n^{-ki}\\ =&\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j\omega_n^{ij})\omega_n^{-ki}\\ =&\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i) \end{align*} \]

令人振奋的是，\(\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i\) 可求。当 \(j=k\) 时，它为 \(n\)。否则，可以推出它为 \(0\)。

总之我们终于得到了下列公式：

\[z_i=na_i\\ a_i=\frac{z_i}{n} \]

我们成功获得了系数！

分治版模板

struct clx{
    flt x,y;
    clx(flt a=0,flt b=0){x=a,y=b;}
    clx operator+(clx b){return clx(x+b.x,y+b.y);}
    clx operator-(clx b){return clx(x-b.x,y-b.y);}
    clx operator*(clx b){return clx(x*b.x-y*b.y,x*b.y+y*b.x);}
};
typedef vec<clx> poly;
const flt Pi=acos(-1);

void fft(int lim,poly &a,int type){
    if(lim==1)return;
    a.resize(lim);
    poly a1(lim>>1),a2(lim>>1);
    repl(i,0,lim>>1)a1[i]=a[i<<1],a2[i]=((i<<1|1)<lim?a[i<<1|1]:0);
    fft(lim>>1,a1,type);fft(lim>>1,a2,type);
    clx w1=clx(cos(Pi*2/lim),type*sin(Pi*2/lim)),w=clx(1,0);
    repl(i,0,lim>>1){
        a[i]=a1[i]+w*a2[i];
        a[i+(lim>>1)]=a1[i]-w*a2[i];
        w=w*w1;
    }
}
poly operator*(poly a,poly b){
    int lim=1,len=a.size()+b.size()-1;
    while(lim<a.size()+b.size())lim<<=1;
    fft(lim,a,1);fft(lim,b,1);
    poly c(lim);
    repl(i,0,lim)c[i]=a[i]*b[i];
    fft(lim,c,-1);
    c.resize(len);
    repl(i,0,c.size())c[i].x/=lim;
    return c;
}

如果是逆变换，\(type\) 传入 \(-1\)。否则，\(type\) 传入 \(1\)。

例题代码

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef __int128 i128;
typedef double flt;
typedef long double ld;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
typedef pair<int,ll> pil;
typedef pair<ll,int> pli;
template <typename Type>
using vec=vector<Type>;
template <typename Type>
using grheap=priority_queue<Type>;
template <typename Type>
using lrheap=priority_queue<Type,vector<Type>,greater<Type> >;
#define fir first
#define sec second
#define pub push_back
#define pob pop_back
#define puf push_front
#define pof pop_front
#define chmax(a,b) a=max(a,b)
#define chmin(a,b) a=min(a,b)
#define rep(i,x,y) for(int i=(x);i<=(y);i++)
#define repl(i,x,y) for(int i=(x);i<(y);i++)
#define per(i,x,y) for(int i=(x);i>=(y);i--)

const int inf=0x3f3f3f3f;
const ll INF=0x3f3f3f3f3f3f3f3f;
const int mod=998244353;

const int N=4e6+5;

int n,m;

struct clx{
    flt x,y;
    clx(flt a=0,flt b=0){x=a,y=b;}
    clx operator+(clx b){return clx(x+b.x,y+b.y);}
    clx operator-(clx b){return clx(x-b.x,y-b.y);}
    clx operator*(clx b){return clx(x*b.x-y*b.y,x*b.y+y*b.x);}
};
typedef vec<clx> poly;
const flt Pi=acos(-1);

void fft(int lim,poly &a,int type){
    if(lim==1)return;
    a.resize(lim);
    poly a1(lim>>1),a2(lim>>1);
    repl(i,0,lim>>1)a1[i]=a[i<<1],a2[i]=((i<<1|1)<lim?a[i<<1|1]:0);
    fft(lim>>1,a1,type);fft(lim>>1,a2,type);
    clx w1=clx(cos(Pi*2/lim),type*sin(Pi*2/lim)),w=clx(1,0);
    repl(i,0,lim>>1){
        a[i]=a1[i]+w*a2[i];
        a[i+(lim>>1)]=a1[i]-w*a2[i];
        w=w*w1;
    }
}
poly operator*(poly a,poly b){
    int lim=1,len=a.size()+b.size()-1;
    while(lim<a.size()+b.size())lim<<=1;
    fft(lim,a,1);fft(lim,b,1);
    poly c(lim);
    repl(i,0,lim)c[i]=a[i]*b[i];
    fft(lim,c,-1);
    c.resize(len);
    repl(i,0,c.size())c[i].x/=lim;
    return c;
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    cin>>n>>m;
    poly a(n+1),b(m+1);
    rep(i,0,n)cin>>a[i].x;
    rep(j,0,m)cin>>b[j].x;
    poly c=a*b;
    repl(i,0,c.size())cout<<int(c[i].x+0.1)<<" ";//精度问题
    return 0;
}

迭代优化

借图一张：

观察最后序列的下标，发现是原序列下标二进制反转！

所以我们可以不用按奇偶排序，可以借此节省很大空间。同时我们也无须递归处理，可以节省很多时间。

具体的，我们先把序列推到最后的样子，然后从下往上进行迭代，同时求出点值表示。

至于如何推到最后的样子，就不得不给出米奇妙妙代码了！

r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));

迭代版模板

struct clx{
    flt x,y;
    clx(flt a=0,flt b=0){x=a,y=b;}
    clx operator+(clx b){return clx(x+b.x,y+b.y);}
    clx operator-(clx b){return clx(x-b.x,y-b.y);}
    clx operator*(clx b){return clx(x*b.x-y*b.y,x*b.y+y*b.x);}
};
typedef vec<clx> poly;
const flt Pi=acos(-1);

vec<int> r;
void fft(int lim,poly &a,int type){
    a.resize(lim);
    repl(i,0,lim)if(i<r[i])swap(a[i],a[r[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        clx w1(cos(Pi/mid),type*sin(Pi/mid));
        for(int j=0;j<lim;j+=(mid<<1)){
            clx w(1,0);
            repl(k,0,mid){
                clx x=a[j+k],y=w*a[j+mid+k];
                a[j+k]=x+y;
                a[j+mid+k]=x-y;
                w=w*w1;
            }
        }
    }
}
poly operator*(poly a,poly b){
    int lim=1,l=0,len=a.size()+b.size()-1;
    while(lim<a.size()+b.size())lim<<=1,l++;
    r.resize(lim);
    repl(i,0,lim)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    fft(lim,a,1);fft(lim,b,1);
    poly c(lim);
    repl(i,0,lim)c[i]=a[i]*b[i];
    fft(lim,c,-1);
    c.resize(len);
    repl(i,0,c.size())c[i].x/=lim;
    return c;
}

例题代码

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef __int128 i128;
typedef double flt;
typedef long double ld;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
typedef pair<int,ll> pil;
typedef pair<ll,int> pli;
template <typename Type>
using vec=vector<Type>;
template <typename Type>
using grheap=priority_queue<Type>;
template <typename Type>
using lrheap=priority_queue<Type,vector<Type>,greater<Type> >;
#define fir first
#define sec second
#define pub push_back
#define pob pop_back
#define puf push_front
#define pof pop_front
#define chmax(a,b) a=max(a,b)
#define chmin(a,b) a=min(a,b)
#define rep(i,x,y) for(int i=(x);i<=(y);i++)
#define repl(i,x,y) for(int i=(x);i<(y);i++)
#define per(i,x,y) for(int i=(x);i>=(y);i--)


const int inf=0x3f3f3f3f;
const ll INF=0x3f3f3f3f3f3f3f3f;
const int mod=998244353;

const int N=4e6+5;

int n,m;

struct clx{
    flt x,y;
    clx(flt a=0,flt b=0){x=a,y=b;}
    clx operator+(clx b){return clx(x+b.x,y+b.y);}
    clx operator-(clx b){return clx(x-b.x,y-b.y);}
    clx operator*(clx b){return clx(x*b.x-y*b.y,x*b.y+y*b.x);}
};
typedef vec<clx> poly;
const flt Pi=acos(-1);

vec<int> r;
void fft(int lim,poly &a,int type){
    a.resize(lim);
    repl(i,0,lim)if(i<r[i])swap(a[i],a[r[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        clx w1(cos(Pi/mid),type*sin(Pi/mid));
        for(int j=0;j<lim;j+=(mid<<1)){
            clx w(1,0);
            repl(k,0,mid){
                clx x=a[j+k],y=w*a[j+mid+k];
                a[j+k]=x+y;
                a[j+mid+k]=x-y;
                w=w*w1;
            }
        }
    }
}
poly operator*(poly a,poly b){
    int lim=1,l=0,len=a.size()+b.size()-1;
    while(lim<a.size()+b.size())lim<<=1,l++;
    r.resize(lim);
    repl(i,0,lim)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    fft(lim,a,1);fft(lim,b,1);
    poly c(lim);
    repl(i,0,lim)c[i]=a[i]*b[i];
    fft(lim,c,-1);
    c.resize(len);
    repl(i,0,c.size())c[i].x/=lim;
    return c;
}

int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    cin>>n>>m;
    poly a(n+1),b(m+1);
    rep(i,0,n)cin>>a[i].x;
    rep(j,0,m)cin>>b[j].x;
    poly c=a*b;
    repl(i,0,c.size())cout<<int(c[i].x+0.1)<<" ";
    return 0;
}

单次 DFT 优化

优雅地利用复数的虚数部分减少一次 DFT 操作。

记做乘法的两个多项式为 \(A,B\)，额外记录两个多项式：

\[P(x)=A(x)+iB(x)\\ Q(x)=A(x)-iB(x) \]

令对两多项式 DFT 的结果为 \(\mathrm{DFT}(P),\mathrm{DFT}(Q)\)，会有一性质（证明可以参考这篇博客）：

\[\mathrm{DFT}(Q)\left[k\right]=\overline{\mathrm{DFT}(P)\left[N-k\right]} \]

那么我们对 \(P\) 进行一次 DFT 即可直接得到 \(Q\) 的值。

\[\mathrm{DFT}(A)\left[k\right]=\frac{\mathrm{DFT}(P)\left[k\right]+\mathrm{DFT}(Q)\left[k\right]}{2}\\ \mathrm{DFT}(B)\left[k\right]=-i\frac{\mathrm{DFT}(P)\left[k\right]-\mathrm{DFT}(Q)\left[k\right]}{2} \]

相乘的结果为：

\[-i\frac{\mathrm{DFT}(P)\left[k\right]^2-\mathrm{DFT}(Q)\left[k\right]^2}{4} \]

然后对得到的序列执行 IDFT 即可。这样就只调用了两次 FFT 操作。

代码实现

inline void operator*=(poly &a,poly b){
	int lim=1,l=0,len=a.size()+b.size()-1;
	while(lim<len)lim<<=1,l++;
	Rt.resize(lim);
	repl(i,0,lim)Rt[i]=(Rt[i>>1]>>1)|((i&1)<<(l-1));
	a.resize(lim),b.resize(lim);poly c(lim);
	repl(i,0,lim)c[i]=clx(a[i].x,b[i].x);
	fft(lim,c,1);
	repl(i,0,lim)a[i]=(c[i]*c[i]-!c[i?lim-i:i]*!c[i?lim-i:i])*clx(0,-0.25);
	fft(lim,a,-1);
	a.resize(len);
	repl(i,0,a.size())a[i].x/=lim;
}

整合封装模板

namespace FFT{
    struct clx{
        flt x,y;
        clx(flt a=0,flt b=0){x=a,y=b;}
        inline clx operator+=(const clx &b){return x+=b.x,y+=b.y,*this;}
        friend inline clx operator+(clx a,clx b){return a+=b;}
        inline clx operator-=(const clx &b){return x-=b.x,y-=b.y,*this;}
        friend inline clx operator-(clx a,clx b){return a-=b;}
        inline clx operator*=(const clx &b){return *this=clx(x*b.x-y*b.y,x*b.y+y*b.x);}
        friend inline clx operator*(clx a,clx b){return a*=b;}
        inline clx operator!(){return clx(x,-y);}
    };
    typedef vec<clx> poly;
    const flt Pi=acos(-1);
    vec<int> Rt;
    inline void fft(int lim,poly &a,int type){
        repl(i,0,lim)if(i<Rt[i])swap(a[i],a[Rt[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            clx w1(cos(Pi/mid),type*sin(Pi/mid));
            for(int j=0;j<lim;j+=(mid<<1)){
                clx w(1,0);
                repl(k,0,mid){
                    clx x=a[j+k],y=w*a[j+mid+k];
                    a[j+k]=x+y;
                    a[j+mid+k]=x-y;
                    w=w*w1;
                }
            }
        }
    }
    inline void operator*=(poly &a,poly b){
        int lim=1,l=0,len=a.size()+b.size()-1;
        while(lim<len)lim<<=1,l++;
        Rt.resize(lim);
        repl(i,0,lim)Rt[i]=(Rt[i>>1]>>1)|((i&1)<<(l-1));
        a.resize(lim),b.resize(lim);poly c(lim);
        repl(i,0,lim)c[i]=clx(a[i].x,b[i].x);
        fft(lim,c,1);
        repl(i,0,lim)a[i]=(c[i]*c[i]-!c[i?lim-i:i]*!c[i?lim-i:i])*clx(0,-0.25);
        fft(lim,a,-1);
        a.resize(len);
        repl(i,0,a.size())a[i].x/=lim;
    }
    inline void operator+=(poly &a,poly b){
        if(a.size()<b.size())a.resize(b.size());
        repl(i,0,b.size())a[i]=a[i]+b[i];
    }
    inline void operator-=(poly &a,poly b){
        if(a.size()<b.size())a.resize(b.size());
        repl(i,0,b.size())a[i]=a[i]-b[i];
    }
    poly operator*(poly a,poly b){return a*=b,a;}
    poly operator+(poly a,poly b){return a+=b,a;}
    poly operator-(poly a,poly b){return a-=b,a;}
}using namespace FFT;

posted @ 2025-01-08 13:17 LastKismet 阅读(104) 评论(2) 收藏举报

刷新页面返回顶部

LastKismet

愿群星将你我连结、愿此刻的流光永镌

快速傅里叶变换 FFT

前言

简介

多项式

系数表示法

点值表示法

复数

概念与表示

运算法则

加减法

乘法

单位根

代码

FFT

快速傅里叶变换 DFT

快速傅里叶逆变换 IDFT

分治版模板

例题代码

迭代优化

迭代版模板

例题代码

单次 DFT 优化

代码实现

整合封装模板

公告