FFT学习笔记

2021/5/1 10:25:10

编程Tag： 笔记 学习 limit sum complex pi frac fft 2j

本文主要是介绍FFT学习笔记，对大家解决编程问题具有一定的参考价值，需要的程序猿们随着小编来一起学习吧！

这是多项式全家桶最基础的一个东西，后面的多项式操作大多都建立在 \(FFT\) 的基础上

这个比较麻烦和基础，咕了

最后还是因为自己的懒惰吃了亏啊，上面这一行咕了的记录就不删了，以便时不时地回来“嘲讽”自己两句。

前置知识：

  1. 点值表达式
  2. 复数
  3. 单位根

对于上面这三个前置知识就不细讲了，不是主体内容。

关于 \(FFT\) ，分为 \(DFT\) 和 \(IDFT\) 这样两个操作，其中 \(IDFT\) 是 \(DFT\) 的逆运算。

对于多项式 $ f(x) $ ，$ g(x) $ 如果能分别求出它们的点值表达式 \(F(x),G(x)\)，那么就意味着可以在 \(O(n)\) 的时间内得到 \(f(x)*g(x)\) 的点值表达式。

但是朴素的得到点值表达式的过程是 \(O(n^2)\) 的，如果无法解决这个最基本瓶颈那么接下来的操作都是空中楼阁。

一位叫傅立叶的法国人发现，如果把单位根代进去，可以在比较优秀的时间内得到由 \(f(x)\) 得到 \(F(x)\) ，这个过程称为 \(DFT\) 。

具体地说，设当前要得到 \(L\) 个点值（\(L=2^k\)），这 \(L\) 个点值分别为 \(w^{\tfrac{2i\pi}{L}},i\in[0,L)\)，其中 \(w\) 是单位根。

为了方便，接下来设 \(w^{i}_{L}=w^{\tfrac{2i\pi}{L}}\)

那么可以得到式子：

\[\begin{aligned} f(w^{i}_{L})&=\sum^{L-1}_{j=0}a_{j} w^{ij}_{L}\\ &=a_{0}w^{0}+a_{1}w^{1}_{L}+\dots+a_{L-1}w^{L-1}_{L}\\ &=(a_{0}w^{0}+a_{2}w^{2}_{L}+\dots+a_{L-2}w^{L-2)}_{L})+w^{1}_{L}(a_{1}w^{0}_{L}+a_{3}w^{2}_{L}+\dots+a_{L-1}w^{L-2}_{L})\\ &=A(w^{2}_{L})+w^{1}_{L}B(w^{2}_{L})\\ \end{aligned} \]

于是，可以考虑先求出 \(A(x)\) 和 \(B(x)\) ，再反推回 \(F(x)\) 。

这是一个分治的过程，每次以 \(O(n)\) 的代价实现了将一个问题变成两个规模减半的子问题。

也就是 \(T(n)=2T(\frac{n}{2})+O(n)\) ，根据主定理或者直接将 \(T(n)\) 暴力展开，都可以得到 \(T(n)=O(n\log n)\) 。

这样，我们就成功的以可以接受的时间复杂度实现了 \(f(x)\) 从系数表达式转换成了点值表达式。

这种分治的思想是非常巧妙的，如果继续往后学还会碰到很多类似的做法。

简单优化：

根据单位根的性质，有 \(w^{t+\tfrac{L}{2}}_{L}=-w^{t}_{L},w^{t+L}_{L}=w^{t}_{L}\)

所以有

\[\begin{aligned} f(w^{t+\frac{L}{2}}_{L})&=A(w^{2t+L}_{L})+w^{t+\frac{L}{2}}_{L}B(w^{2t+L}_{L})\\ &=A(w^{2t}_{L})-w^{t}_{L}B(w^{2t}_{L})\\ f(w^{t}_{L})&=A(w^{2t}_{L})+w^{t}B(w^{2t}_{L})\\ \end{aligned} \]

因此可以把求 \(f(w^{t+\frac{L}{2}}_{L})\) 和 \(f(w^{t}_{L})\) 的这个过程给合并起来。

这个操作的专业装B称呼为『蝴蝶操作』。

DFT过程的实现代码 
const double Pi=acos(-1.0);
struct complex
{
    double x,y;
    complex (double xx=0,double yy=0){x=xx,y=yy;}
    complex operator+(const complex &b){
        complex a;
        a.x=x+b.x;
        a.y=y+b.y;
        return a;
    }
    complex operator-(const complex &b){
        complex a;
        a.x=x-b.x;
        a.y=y-b.y;
        return a;
    }
    complex operator*(const complex &b){
        complex a;
        a.x=x*b.x-y*b.y;
        a.y=x*b.y+y*b.x;
        return a;
    }
};
void DFT(int limit,complex *a){
    if(limit==1) return;
    complex a1[limit>>1],a2[limit>>1];
    for(ri i=0;i<limit;i+=2){
        a1[i>>1]=a[i],a2[i>>1]=a[i+1];
    }
    DFT(limit>>1,a1);
    DFT(limit>>1,a2);
    complex Wn=complex(cos(2.0*Pi/limit),sin(2.0*Pi/limit)),w=complex(1,0);
    for(ri i=0;i<limit>>1;++i,w=w*Wn){
        a[i]=a1[i]+w*a2[i];
        a[i+(limit>>1)]=a1[i]-w*a2[i];
    }
}

但是不能光有 \(DFT\) 呀，一般多项式还是习惯使用系数表达的，如果没法转回来的话还是相当于啥都没干！

那么该怎么转回来呢？多项式快速插值

考虑继续利用一些单位根的性质，在上面已经求出了\(f(x)\) 的点值表达式，接下来尝试求得到的这个点值表达式的点值表达式听起来很绕，设之前的 \(F(x)\) 的点值表达式是 \(H(x)\) 。

\[\begin{aligned} H_{k}&=\sum^{L-1}_{i=0}F_{k}(w_{L}^{-k})^i\\ &=\sum^{L-1}_{i=0}(\sum^{L-1}_{j=0}f_{j}(w^{i}_{L})^{j})(w_{L}^{-k})^i\\ &=\sum^{L-1}_{i=0}(\sum^{L-1}_{j=0}f_{j}(w_{L}^{j})^{i})(w_{L}^{-k})^i\\ &=\sum^{L-1}_{i=0}\sum^{L-1}_{j=0}f_{j}(w_{L}^{j-k})^i\\ &=\sum^{L-1}_{j=0}f_{j}\sum^{L-1}_{i=0}(w_{L}^{j-k})^i\\ \end{aligned} \]

不妨令后面那一陀为 \(s(w_{L}^{j-k})\) 。

那么 \(s(x)=\sum\limits^{L-1}_{i=0}x^i\) ，根据等比数列求和可以得到 \(s(x)=\tfrac{1-x^L}{1-x}\)。

我们把 \(w^{j-k}_L\) 次代入，可以发现当 \(j\not =k\)时，分子显然为 \(0\) ，分母不为 \(0\) 。

而当 \(j=k\) ，即 \(w^{j-k}_{L}=1\) 的时候，可以发现 \(s(x)=L\) 。

所以有 \(H_{k}=Lf_{k},f_{k}=\frac{H_{k}}{L}\) 。

重新梳理一下上面的过程，如果我们使用了 \(DFT\) 得到了 \(f(x)\) 的点值表达式 \(F(x)\) ，然后再用一次 \(DFT\)，不过把单位根从 \(w^{t}_{L}\) 换成了\(w^{-t}_{L}\) 就可以得到 \(H(x)\) 。最后只要将 \(H(x)\) 的每一项除以 \(L\) 就可以得到 \(F(x)\) 了。

显然 \(DFT\) 和 \(IDFT\) 的代码可以合并。

IDFT过程的实现代码 
const double Pi=acos(-1.0);
struct complex
{
    double x,y;
    complex (double xx=0,double yy=0){x=xx,y=yy;}
    complex operator+(const complex &b){
        complex a;
        a.x=x+b.x;
        a.y=y+b.y;
        return a;
    }
    complex operator-(const complex &b){
        complex a;
        a.x=x-b.x;
        a.y=y-b.y;
        return a;
    }
    complex operator*(const complex &b){
        complex a;
        a.x=x*b.x-y*b.y;
        a.y=x*b.y+y*b.x;
        return a;
    }
};
void DFT(int limit,complex *a,int type){
    if(limit==1) return;
    complex a1[limit>>1],a2[limit>>1];
    for(ri i=0;i<limit;i+=2){
        a1[i>>1]=a[i],a2[i>>1]=a[i+1];
    }
    DFT(limit>>1,a1,type);
    DFT(limit>>1,a2,type);
    complex Wn=complex(cos(2.0*Pi/limit),sin(2.0*type*Pi/limit)),w=complex(1,0);
    for(ri i=0;i<limit>>1;++i,w=w*Wn){
        a[i]=a1[i]+w*a2[i];
        a[i+(limit>>1)]=a1[i]-w*a2[i];
    }
}

这样，就实现了 \(O(n\log n)\) 求得两个多项式的卷积了。

优化

当然，这样的代码可能并没有办法通过模板题，有很多细节还可以继续优化一下。

  1. 减少乘法次数

    复数的乘法运算次数比较多，注意到 \(w\times a_{2}\) 这个出现了两次，可以预先记录一下。

  2. 递归->迭代

    可以发现，递归版效率之所以低，很大原因是一直在开数组，如果能够用迭代实现的话就可以在原数组上进行滚动而不用额外开数组了。

    可以发现上面递归的过程实际上是先不断地改变 \(a\) 数组的位置，直到最后一层为止。然后不停地执行 \(a=a_{1}+w\times a_{2}\)

    \[\begin{bmatrix} 0,1,2,3,4,5,6,7 \end{bmatrix} \\ \downarrow\\ \begin{bmatrix} 0,2,4,6],[1,3,5,7 \end{bmatrix}\\ \downarrow\\ \begin{bmatrix} 0,4],[2,6],[1,5],[3,7 \end{bmatrix}\\ \downarrow\\ \begin{bmatrix} 0],[4],[2],[6],[1],[5],[3],[7 \end{bmatrix}\\ \]

    这样似乎看不出什么特点，但是转成二进制之后就有：

    \[ \begin{bmatrix} 000,001,010,011,100,101,110,111 \end{bmatrix}\\ \downarrow\\ \begin{bmatrix} 000,100,010,110,001,101,011,111 \end{bmatrix}\\ \]

    刚好每一位的二进制在操作之后都进行了一次翻转。

    因此，可以先记录一个 \(r\) 数组，记录一下下标之间的对应关系。

    而 \(r\) 数组的求法也比较简单，对于 \(r_{i}\) ，都是可以通过 \(r_{i>>1}\) 转移过来的，只需要对 \(i\) 分奇偶讨论就可以了。

    求r数组的过程 
    void rev(int limit,int ws){
		//ws=log2(limit)
		for(ri i=0;i<limit;i+=2){
			r[i]=r[i>>1]>>1;
			r[i|1]=(r[i>>1]>>1)|(1<<ws-1);
		}
	}

    然后，就是用三个for 来实现递归的过程。

    Code 
    void FFT(int limit,complex *a,int type){
    for(ri i=0;i<=limit;++i){
        if(i<r[i]) swap(a[i],a[r[i]]);
    }
    for(ri l=1;l<limit;l<<=1){
        complex Wn=complex(cos(Pi/l),sin(type*Pi/l));
        for(ri i=0;i<limit;i+=(l<<1)){
            complex w=complex(1,0);
            for(ri j=0;j<len;++j){
                complex x=a[i+j],y=w*a[i+j+len];
                a[i+j]=x+y;
                a[i+j+len]=x-y;
                w=w*Wn;
            }
        }
    }
}

  3. 预处理单位根

    就是把计算单位根的过程提前预处理出来。

    但是 \(IDFT\) 过程的单位根也要预处理吗？

    注意到 \(w_{L}^{-t}=w_{L}^{L-t}\)

    因此，并不需要在计算的过程中判断当前是在做 \(DFT\) 还是 \(IDFT\) ，只需要在最后判断是否需要reverse 就可以了。

    Code 
    vector<Complex> W[23];
const double Pi=acos(-1.0);
void init(){
  int d=log(N)/log(2)+0.5;
  for(ri t=1;t<=d;++t){
     W[t].resize(1<<t);
     W[t][0]=(Complex){1,0},W[t][1]=(Complex){cos(Pi/(1<<t-1)),sin(Pi/(1<<t-1))};
     for(ri i=2;i<(1<<t);++i) W[t][i]=W[t][i-1]*W[t][1];
  }
}
 void DFT(int limit,Complex *f,int flag){
     for(ri i=1;i<limit;++i)
         if(r[i]<i) swap(f[i],f[r[i]]);
     for(ri l=1,t=1;l<limit;++t,l<<=1){
         for(ri i=0;i<limit;i+=(l<<1)){
             Complex *w=&W[t][0];
             for(ri j=0;j<l;++j){
                 Complex tmp=f[i+j+l]*(*w++);
                 f[i+j+l]=f[i+j]-tmp;
                 f[i+j]=f[i+j]+tmp;
             }
         }
     }
     if(flag==-1){
         reverse(f+1,f+limit);
         for(ri i=0;i<limit;++i) f[i].x/=limit,f[i].y/=limit;
     }
 }
最终代码 
#include<cmath>
#include<cstdio>
#include<algorithm>
#include<vector>
using namespace std;
#define il inline
#define ri register int
#define ll long long
#define ui unsigned int
il ll read(){
    bool f=true;ll x=0;
    register char ch=getchar();
    while(ch<'0'||ch>'9') {if(ch=='-') f=false;ch=getchar();}
    while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    if(f) return x;
    return ~(--x);
}
il void write(const ll &x){if(x>9) write(x/10);putchar(x%10+'0');}
il void print(const ll &x) {x<0?putchar('-'),write(~(x-1)):write(x);putchar(' ');}
il ll max(const ll &a,const ll &b){return a>b?a:b;}
il ll min(const ll &a,const ll &b){return a<b?a:b;}
const int N=1<<21;
int r[N];
void rev(int limit,int ws){
    for(ri i=0;i<limit;i+=2){
        r[i]=r[i>>1]>>1;
        r[i|1]=(r[i>>1]>>1)|(1<<ws-1);
    }
}
struct Complex
{
    double x,y;
    Complex operator+(const Complex &p)const{
        return (Complex){x+p.x,y+p.y};
    }
    Complex operator-(const Complex &p)const{
        return (Complex){x-p.x,y-p.y};
    }
    Complex operator*(const Complex &p)const{
        return (Complex){x*p.x-y*p.y,x*p.y+y*p.x};
    }
};
vector<Complex> W[23];
const double Pi=acos(-1.0);
void init(){
    int d=log(N)/log(2)+0.5;
    for(ri t=1;t<=d;++t){
        W[t].resize(1<<t);
        W[t][0]=(Complex){1,0},W[t][1]=(Complex){cos(Pi/(1<<t-1)),sin(Pi/(1<<t-1))};
        for(ri i=2;i<(1<<t);++i) W[t][i]=W[t][i-1]*W[t][1];
    }
}
void DFT(int limit,Complex *f,int flag){
    for(ri i=1;i<limit;++i)
        if(r[i]<i) swap(f[i],f[r[i]]);
    for(ri l=1,t=1;l<limit;++t,l<<=1){
        for(ri i=0;i<limit;i+=(l<<1)){
            Complex *w=&W[t][0];
            for(ri j=0;j<l;++j){
                Complex tmp=f[i+j+l]*(*w++);
                f[i+j+l]=f[i+j]-tmp;
                f[i+j]=f[i+j]+tmp;
            }
        }
    }
    if(flag==-1){
        reverse(f+1,f+limit);
        for(ri i=0;i<limit;++i) f[i].x/=(limit<<1),f[i].y/=(limit<<1);
    }
}
Complex FFT_f[N];
void FFT(int n,ll *f,int m,ll *g,int lim=0){
    if(!lim) lim=n+m;
    int limit=1,ws=0;
    for(;limit<=n+m;limit<<=1,++ws);
    rev(limit,ws);
    for(ri i=0;i<limit;++i) FFT_f[i].x=FFT_f[i].y=0;
    for(ri i=0;i<=n;++i) FFT_f[i].x=f[i];
    for(ri i=0;i<=m;++i) FFT_f[i].y=g[i];
    DFT(limit,FFT_f,1);
    for(ri i=0;i<limit;++i) FFT_f[i]=FFT_f[i]*FFT_f[i];
    DFT(limit,FFT_f,-1);
    for(ri i=0;i<=lim;++i) f[i]=FFT_f[i].y+0.5;
}
ll n,m,f[N],g[N];
int main(){
    init();
    // freopen("rand.in","r",stdin);
    // freopen("1.out","w",stdout);
    n=read(),m=read();
    for(ri i=0;i<=n;++i) f[i]=read();
    for(ri i=0;i<=m;++i) g[i]=read();
    FFT(n,f,m,g);
    for(ri i=0;i<=n+m;++i) print(f[i]);
    return 0;
}

进一步优化

上面的写法用到了一些性质。

构造多项式 \(F(x)=A(x)+iB(x)\) ，那么有 \(F(x)^2=A(x)^2-B(x)^2+2iA(x)B(x)\)。

因此只需要做一次 \(DFT\) 和一次 \(IDFT\) 就可以实现卷积了。

其实上面这种做法还隐藏着更多的性质。

\[\begin{aligned} F(x)&=A(x)+iB(x)\\ G(x)&=A(x)-iB(x)\\ F(w^{t}_{L})&=\sum^{L-1}_{j=0}(a_{j}\sin(\frac{2j\pi}{L})-b_{j}\cos(\frac{2j\pi}{L}))+i(a_{j}\cos(\frac{2j\pi}{L})+b_{j}\sin(\frac{2j\pi}{L}))\\ G(w^{t}_{L})&=\sum^{L-1}_{j=0}(a_{j}\sin(\frac{2j\pi}{L})+b_{j}\cos(\frac{2j\pi}{L}))+i(a_{j}\cos(\frac{2j\pi}{L})-b_{j}\sin(\frac{2j\pi}{L}))\\ \end{aligned} \]

所以有

\[\begin{aligned} F(w^{t}_{L})+G(w^{t}_{L})&=\sum^{L-1}_{j=0}2a_{j}\sin(\frac{2j\pi}{L})+2ia_{j}\cos(\frac{2j\pi}{L})=2A(w^{t}_{L})\\ A(w^{t}_{L})&=\frac{F(w_{L}^{t})+G(w^{t}_{L})}{2}\\ \texttt{同理有}B(w^{t}_{L})&=-i\frac{F(w^{t}_{L})-G(w^{t}_{L})}{2}\\ \end{aligned} \]

但是这样似乎还是没有什么用呀，还是要做两次 \(DFT\) ，而且常数还变大了。

仔细观察 \(F(x)\) 与 \(G(x)\) ，发现它们有很多关联。

\[\begin{aligned} G(w^{t}_{L})&=\sum^{L-1}_{j=0}(a_{j}\sin(\frac{2j\pi}{L})+b_{j}\cos(\frac{2j\pi}{L}))+i(a_{j}\cos(\frac{2j\pi}{L})-b_{j}\sin(\frac{2j\pi}{L}))\\ &=\sum^{L-1}_{j=0}conj(a_{j}\sin(\frac{2j\pi}{L})+b_{j}\cos(\frac{2j\pi}{L}))-i(a_{j}\cos(\frac{2j\pi}{L})-b_{j}\sin(\frac{2j\pi}{L})))\\ &=conj(\sum^{L-1}_{j=0}a_{j}\sin(\frac{2j\pi}{L})+b_{j}\cos(\frac{2j\pi}{L}))-i(a_{j}\cos(\frac{2j\pi}{L})-b_{j}\sin(\frac{2j\pi}{L})))\\ &=conj(\sum^{L-1}_{j=0}a_{j}\sin(\frac{(2L-2j)\pi}{L})-b_{j}\cos(\frac{(2L-2j)\pi}{L}))+i(a_{j}\cos(\frac{(2L-2j)\pi}{L})+b_{j}\sin(\frac{(2L-2j)\pi}{L})))\\ &=conj(F(w^{L-j}_{L}))\\ \end{aligned} \]

也就是说，只要得到了 \(DFTF(x)\) ，就等价于得到 \(DFTG(x)\) ，同时也就得到了 \(DFTA(x)\) 和 \(DFTB(x)\) ,而这一切只需要做一次 \(DFT\) !

参考文献

国家集训队2016论文 - 再谈快速傅里叶变换



这篇关于FFT学习笔记的文章就介绍到这儿，希望我们推荐的文章对大家有所帮助，也希望大家多多支持为之网！


相关编程文章

更多>
本站官方QQ
官方QQ:
23002807

沪ICP备2023033072号-1

     

扫一扫关注最新编程教程