FFT快速傅里叶变换

FFT快速傅里叶变换

some 多项式变换

DFT：离散傅里叶变换 $⟶$ $O(n^2)$计算多项式乘法
FFT：快速傅里叶变换 $⟶$ $O(n\log n)$计算多项式乘法
FNTT/NTT：快速傅里叶变换的优化版 $⟶$ 优化常数及误差
FWT：快速沃尔什变换 $⟶$ 利用类似FFT的东西解决一类卷积问题
MTT：毛爷爷的FFT $⟶$ 非常nb/任意模数
FMT 快速莫比乌斯变化

快速傅里叶变换

快速傅里叶变换($FFT$) 算法的基本思想是分治。

就 $DFT$ 来说，它分治地来求当 $x={\omega }_n^k$ 的时候 $f(x)$ 的值。

$FFT$ 的分治思想体现在将多项式分为奇次项和偶次项处理。

对于一个 $n-1$ 次的多项式：

$$f(x)=a_0+a_1\cdot x+a_2\cdot x^2+...+a_{n-2}\cdot x^{n-2}+a_{n-1}\cdot x^{n-1}$$

按照次数的奇偶来分成两组，然后右边提出来一个 $x$：

$$\begin{array}{rl}f(x)& =(a_0+a_2\cdot x^2+a_4\cdot x^4+...+a_{n-2}\cdot x^{n-2})+(a_1\cdot x+a_3\cdot x^3+a_5\cdot x^5+...+a_{n-1}\cdot x^{n-1})\\ & =(a_0+a_2\cdot x^2+a_4\cdot x^4+...+a_{n-2}\cdot x^{n-2})+x\cdot (a_1+a_3\cdot x^2+a_5\cdot x^4+a_{n-1}{x}^{n-2})\end{array}$$

分别用奇偶次次项数建立新的函数：

$$\begin{array}{rl}G(x)& =a_0+a_{2}x+a_4{x}^2+...+a_{n-2}{x}^{\frac{n}{2}-1}\\ H(x)& =a_1+a_{3}x+a_5{x}^2+...+a_{n-1}{x}^{\frac{n}{2}-1}\end{array}$$

那么原来的 $f(x)$ 用新函数表示为：

$$f(x)=G\left(x^2\right)+x\times H\left(x^2\right)$$

利用偶数次单位根的性质：${\omega }_n^i=-{\omega }_n^{i+\frac{n}{2}}$

$G\left(x^2\right)$ 和 $H\left(x^2\right)$ 是偶函数

我们知道在复平面上 ${\omega }_n^i$ 和 ${\omega }_n^{i+\frac{n}{2}}$ 的 $G(x^2)$ 的 $H(x^2)$ 对应的值相同。

得到：

$$\begin{array}{rl}f({\omega }_n^k)& =G(({\omega }_n^k{)}^2)+{\omega }_n^k\times H(({\omega }_n^k{)}^2)\\ & =G({\omega }_n^{2k})+{\omega }_n^k\times H({\omega }_n^{2k})\\ & =G({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k\times H({\omega }_{\frac{n}{2}}^k)\end{array}$$

和：

$$\begin{array}{rl}f({\omega }_n^{k+\frac{n}{2}})& =G({\omega }_n^{2k+n})+{\omega }_n^{k+\frac{n}{2}}\times H({\omega }_n^{2k+n})\\ & =G({\omega }_n^{2k})-{\omega }_n^k\times H({\omega }_n^{2k})\\ & =G({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k\times H({\omega }_{\frac{n}{2}}^k)\end{array}$$

因此我们求出了 $G({\omega }_{\frac{n}{2}}^k)$ 和 $H({\omega }_{\frac{n}{2}}^k)$ 后，就可以同时求出 $f({\omega }_n^k)$ 和 $f({\omega }_n^{k+\frac{n}{2}})$。

于是对 $G$ 和 $H$ 分别递归 $DFT$ 即可。

考虑到分治 $DFT$ 能处理的多项式长度只能是 $2^m(m\in {\mathbf{N}}^{\ast })$，否则在分治的时候左右不一样长，右边就取不到系数了。

所以要在第一次 $DFT$ 之前就把序列向上补成长度为 $2^m(m\in {\mathbf{N}}^{\ast })$（高次系数补 $0$）、最高项次数为 $2^m-1$ 的多项式。

位逆序置换

原序列的每个数用二进制表示，然后把二进制翻转对称一下，就是最终位置的下标。
可以在 $O(n\log n)$ 或 $O(n)$ 的时间内求出每个数变换后的结果

蝶形运算优化

已知 $G({\omega }_{\frac{n}{2}}^k)$ 和 $H({\omega }_{\frac{n}{2}}^k)$ 后，需要使用下面两个式子求出 $f({\omega }_n^k)$ 和 $f({\omega }_n^{k+\frac{n}{2}})$：

$$\begin{array}{rl}f({\omega }_n^k)& =G({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k\times H({\omega }_{\frac{n}{2}}^k)\\ f({\omega }_n^{k+\frac{n}{2}})& =G({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k\times H({\omega }_{\frac{n}{2}}^k)\end{array}$$

使用位逆序置换后，对于给定的 $n$, $k$：

$G({\omega }_{\frac{n}{2}}^k)$ 的值存储在数组下标为 $k$ 的位置，$H({\omega }_{\frac{n}{2}}^k)$ 的值存储在数组下标为 $k+{\displaystyle \frac{n}{2}}$ 的位置。

$f({\omega }_n^k)$ 的值将存储在数组下标为 $k$ 的位置，$f({\omega }_n^{k+\frac{n}{2}})$ 的值将存储在数组下标为 $k+{\displaystyle \frac{n}{2}}$ 的位置。
因此可以直接在数组下标为 $k$ 和

$k+\frac{n}{2}$ 的位置进行覆写，而不用开额外的数组保存值。
再详细说明一下如何借助蝶形运算完成所有段长度为 $\frac{n}{2}$ 的合并操作：

令段长度为 $s=\frac{n}{2}$；

同时枚举序列 $\{G({\omega }_{\frac{n}{2}}^k)\}$ 的左端点 $l_g=0,2s,4s,\cdots ,N-2s$ 和序列 $\{H({\omega }_{\frac{n}{2}}^k)\}$ 的左端点 $l_h=s,3s,5s,\cdots ,N-s$；

合并两个段时，枚举 $k=0,1,2,\cdots ,s-1$，此时$G({\omega }_{\frac{n}{2}}^k)$ 存储在数组下标为 $l_g+k$ 的位置，$H({\omega }_{\frac{n}{2}}^k)$ 存储在数组下标为 $l_h+k$ 的位置；

使用蝶形运算求出 $f({\omega }_n^k)$ 和 $f({\omega }_n^{k+\frac{n}{2}})$，然后直接在原位置覆写。

快速傅里叶逆变换

设 $b_i={\omega }_n^{-i}$，则多项式 $A$ 在 $x=b_0,b_1,\cdots ,b_{n-1}$ 处的点值表示法为 $\{A(b_0),A(b_1),\cdots ,A(b_{n-1})\}$。

对 $A(x)$ 的定义式做一下变换，可以将 $A(b_k)$ 表示为

$$\begin{array}{rl}A(b_k)& =\sum _{i=0}^{n-1}f({\omega }_n^i){\omega }_n^{-ik}=\sum _{i=0}^{n-1}{\omega }_n^{-ik}\sum _{j=0}^{n-1}{a}_j({\omega }_n^i{)}^j\\ & =\sum _{i=0}^{n-1}\sum _{j=0}^{n-1}{a}_j{\omega }_n^{i(j-k)}=\sum _{j=0}^{n-1}{a}_j\sum _{i=0}^{n-1}{\left({\omega }_n^{j-k}\right)}^i\end{array}$$

记

$$S\left({\omega }_n^a\right)=\sum _{i=0}^{n-1}{\left({\omega }_n^a\right)}^i$$

当 $a=0(\mathrm{mod}n)$ 时，$S\left({\omega }_n^a\right)=n$。

当 $a\ne 0(\mathrm{mod}n)$ 时，我们错位相减

$$\begin{array}{rl}S\left({\omega }_n^a\right)& =\sum _{i=0}^{n-1}{\left({\omega }_n^a\right)}^i\\ {\omega }_n^{a}S\left({\omega }_n^a\right)& =\sum _{i=1}^n{\left({\omega }_n^a\right)}^i\\ S\left({\omega }_n^a\right)& =\frac{{\left({\omega }_n^a\right)}^n-{\left({\omega }_n^a\right)}^0}{{\omega }_n^a-1}=0\end{array}$$

也就是说

$$S({\omega }_n^a)=\{\begin{array}{l}n,a=0\\ 0,a\ne 0\end{array}$$

那么代回原式

$$A(b_k)=\sum _{j=0}^{n-1}{a}_{j}S\left({\omega }_n^{j-k}\right)=a_k\cdot n$$

也就是说给定点 $b_i={\omega }_n^{-i}$，则 $A$ 的点值表示法为

$$\begin{array}{rl}& \{(b_0,A(b_0)),(b_1,A(b_1)),\cdots ,(b_{n-1},A(b_{n-1}))\}\\ =& \{(b_0,a_0\cdot n),(b_1,a_1\cdot n),\cdots ,(b_{n-1},a_{n-1}\cdot n)\}\end{array}$$

综上所述，我们取单位根为其倒数，对 $\{y_0,y_1,y_2,\cdots ,y_{n-1}\}$ 跑一遍 $FFT$，然后除以 $n$ 即可得到 $f(x)$ 的系数表示。

递归写法

il void FFT(qwq *q,int lim,int k){
    if(lim==1) return;
    for(ri int i=0;i<lim;++i) o[i]=q[i];
    ri int m=lim>>1;
    for(ri int i=0;i<lim;i++){
        if(i&1) q[m+(i>>1)]=o[i];
        else q[i>>1]=o[i];
    }
    FFT(q,m,k),FFT(q+m,m,k);
    wn={cos(Pi*2/lim),k*sin(Pi*2/lim)},w={1,0};
    for(ri int i=0;i<m;++i,w=w*wn){
        x=q[i],y=w*q[i+m];
        q[i]=x+y,q[i+m]=x-y;
    } 
    return;
}

code

#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
using namespace std;

namespace Q{
    il int rd(){
        ri int x=0;ri bool f=0;ri char c=getchar();
        while(!isdigit(c)) f|=(c==45),c=getchar();
        while(isdigit(c)) x=x*10+(c^48),c=getchar();
        return f?-x:x;
    }
    il void wt(int x){
        if(x<0) x=-x,putchar(45);
        if(x>=10) wt(x/10);
        return putchar(x%10|48),void();
    }
} using namespace Q;

cs double Pi=3.141592653589793;
cs int N=2097152;

namespace FastFastTle{
    struct qwq{
        double r,i;
        qwq operator +(cs qwq o)cs{
            return (qwq){r+o.r,i+o.i};
        }
        qwq operator -(cs qwq o)cs{
            return (qwq){r-o.r,i-o.i};
        } 
        qwq operator *(cs qwq o)cs{
            return (qwq){r*o.r-i*o.i,o.r*i+r*o.i};
        }
    }wn,w,f[N],g[N],o[N],x,y; 
    il void FFT(qwq *q,int lim,int k){
        if(lim==1) return;
        for(ri int i=0;i<lim;++i) o[i]=q[i];
        ri int m=lim>>1;
        for(ri int i=0;i<lim;i++){
            if(i&1) q[m+(i>>1)]=o[i];
            else q[i>>1]=o[i];
        }
        FFT(q,m,k),FFT(q+m,m,k);
        wn={cos(Pi*2/lim),k*sin(Pi*2/lim)},w={1,0};
        for(ri int i=0;i<m;++i,w=w*wn){
            x=q[i],y=w*q[i+m];
            q[i]=x+y,q[i+m]=x-y;
        } 
        return;
    }
} using namespace FastFastTle;

signed main(){
    int n=rd(),m=rd(),lim=1;
    while(lim<=n+m) lim<<=1;
    for(ri int i=0;i<=n;++i) f[i].r=rd();
    for(ri int i=0;i<=m;++i) g[i].r=rd();
    FFT(f,lim,1),FFT(g,lim,1);
    for(ri int i=0;i<lim;++i) f[i]=f[i]*g[i];
    FFT(f,lim,-1);
    for(ri int i=0,s;i<=n+m;++i){
        wt(f[i].r/lim+0.5),putchar(32);
    }
    return 0;
}

迭代实现

il void FFT(qwq q[],int lim,int k){
    for(ri int i=0;i<lim;++i){
        if(i<r[i]) swap(q[i],q[r[i]]);
    }
    for(ri int mid=1;mid<lim;mid<<=1){
        wn={cos(Pi/mid),k*sin(Pi/mid)};
        for(ri int R=mid<<1,j=0;j<lim;j+=R){
            w={1,0};
            for(ri int t=0;t<mid;++t,w=w*wn){
                x=q[j+t],y=w*q[j+mid+t];
                q[j+t]=x+y,q[j+mid+t]=x-y;
            }
        }
    } 
    return;
}

code

#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
using namespace std;

namespace Q{
    il int rd(){
        ri int x=0;ri bool f=0;ri char c=getchar();
        while(!isdigit(c)) f|=(c==45),c=getchar();
        while(isdigit(c)) x=x*10+(c^48),c=getchar();
        return f?-x:x;
    }
    il void wt(int x){
        if(x<0) x=-x,putchar(45);
        if(x>=10) wt(x/10);
        return putchar(x%10|48),void();
    }
} using namespace Q;

cs double Pi=3.141592653589793;
cs int N=2097152;

namespace FastFastTle{
    struct qwq{
        double r,i;
        qwq operator +(cs qwq o)cs{
            return (qwq){r+o.r,i+o.i};
        }
        qwq operator -(cs qwq o)cs{
            return (qwq){r-o.r,i-o.i};
        } 
        qwq operator *(cs qwq o)cs{
            return (qwq){r*o.r-i*o.i,o.r*i+r*o.i};
        }
    }wn,w,f[N],g[N],x,y;
    int l,r[N];
    il void FFT(qwq q[],int lim,int k){
        for(ri int i=0;i<lim;++i){
            if(i<r[i]) swap(q[i],q[r[i]]);
        }
        for(ri int mid=1;mid<lim;mid<<=1){
            wn={cos(Pi/mid),k*sin(Pi/mid)};
            for(ri int R=mid<<1,j=0;j<lim;j+=R){
                w={1,0};
                for(ri int t=0;t<mid;++t,w=w*wn){
                    x=q[j+t],y=w*q[j+mid+t];
                    q[j+t]=x+y,q[j+mid+t]=x-y;
                }
            }
        } 
        return;
    }
} using namespace FastFastTle;

signed main(){
    int n=rd(),m=rd(),lim=1;
    while(lim<=n+m) lim<<=1,l++;
    for(ri int i=0;i<=n;++i) f[i].r=rd();
    for(ri int i=0;i<=m;++i) g[i].r=rd();
    for(ri int i=0;i<lim;++i){
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    }
    FFT(f,lim,1),FFT(g,lim,1);
    for(ri int i=0;i<lim;++i) f[i]=f[i]*g[i];
    FFT(f,lim,-1);
    for(ri int i=0,s;i<=n+m;++i){
        wt(f[i].r/lim+0.5),putchar(32);
    }
    return 0;
}

edit