zoukankan      html  css  js  c++  java
  • 快速傅里叶变换FFT& 数论变换NTT

    相关知识

    时间域上的函数f(t)经过傅里叶变换(Fourier Transform)变成频率域上的F(w),也就是用一些不同频率正弦曲线的加 权叠加得到时间域上的信号。

    [F(omega)=mathcal{F}[f(t)]=intlimits_{-infty}^infty f(t)e^{-iwt}dt ]

    傅里叶逆变换是将频率域上的F(w)变成时间域上的函数f(t),一般称(f(t))为原函数,称(F(w))为象函数。原函数和象函数构成一个傅里叶变换对。

    [f(t)=mathcal{F^{-1}}[F(w)]=frac 1 {2pi}intlimits_{-infty}^infty F(w)e^{iwt}dw ]

    离散傅里叶变换(DFT)是傅里叶变换在时域和频域都呈现离散的形式。

    [x_n=sum_{k=0}^{N-1}X_ke^{frac{2pi i}{N} kn},n=0,...,N-1 ]

    当然也就有离散傅里叶逆变换(IDFT):

    [x_n=frac 1 Nsum_{k=0}^{N-1}X_ke^{-frac{2pi i}{N} kn},n=0,...,N-1 ]

    快速傅里叶变换FFT

    快速傅里叶变换(FFT)利用分治的思想,将 DFT 分解成更小的 DFT进行计算,可以在(O(nlog n))时间复杂度内完成 DFT 的算法。可用来加速卷积和多项式乘法。

    多项式乘法

    对于(C(x)=A(x)B(x)=sum_{i=0}^{n-1}sum_{j=0}^{i}a_jb_{i-j}x^i),(A的最高项次数加上B的最高项次数为n-1)直接计算是(O(n^2))的。我们将它们变成点值表达式

    [A(x):{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})}\\ B(x):{(x_0,y_0'),(x_1,y_1'),...,(x_{n-1},y_{n-1}')}\\ C(x):{(x_0,y_0y_0'),(x_1,y_1y_1'),...,(x_{n-1},y_{n-1}y_{n-1}')} ]

    那么就可以(O(n))计算出C(x)。

    将其转为点值表达式需要带入n个不同的x,利用霍纳法则,带入一个x需要(O(n))时间:

    [A(x_0)=a_0+x_0(a_1+x_0(a_2+...+x_0(a_{n-2}+x_0(a_{n-1}))...)) ]

    从点值表达式转换为系数表达式,用插值法,将n个不同的值带入A,唯一确定了一个系数表达式(证明可见《算法导论》)。

    利用 FFT 可以在(O(nlog n)) 内完成系数到点值,以及从点值到系数表达式。

    n次单位复数根

    n次单位复数根(w)是满足(w^n=1)的复数。 主n次单位根为(w_n=e^{frac {2pi } ni}),其它所有n次单位复数根都是(w_n)的幂次。

    我们知道复数的指数形式的定义是:

    [e^{ui}=cos(u)+isin(u) ]

    所以n个复数根均匀分布在以复平面原点为圆心的单位圆上。

    单位根的性质

    1. 消去引理 (w^{dk}_{dn}=w^k_n),带入定义即可证明。推论:(w^{n/2}_n=w_2=-1)

    2. 折半引理((w^{k+n/2}_n)^2=(w_n^k)^2),展开即可证明。推论:(w^{k+n/2}_n=-w^{k}_n)

    为了方便分治,我们将n补到2的次幂(多项式后面加0)。将n个n次单位根带入A(x)的过程就是对系数向量进行 DFT 。

    [A(w_n^k)=sum_{j=0}^{n-1}a_jw^{kj}_n,k=0,1,..,n-1 ]

    考虑按指数奇偶分类:

    [A_{even}(x)=a_0+a_2x+a_4x^2+...\ A_{odd}(x)=a_1+a_3x+a_5x^2+...\ ]

    那么

    [A(x)=A_{even}(x^2)+xA_{odd}(x^2)\ A(w_n^k)=A_{even}(w_{n}^{2k})+w_n^kA_{odd}(w_n^{2k})\ ]

    (k< n/2)时,由消去引理

    [A_{even}(w_n^{2k})+w_n^kA_{odd}(w_n^{2k})=A_{even}(w_{n/2}^k)+w_n^kA_{odd}(w_{n/2}^k) ]

    也就是相当于对两个长度为n/2的多项式(A_{even}(x)),(A_{odd}(x))进行 DFT。

    指数不小于(n/2)的部分与(k+n/2)一一对应,由折半引理

    [A_{even}((w_n^{k+n/2})^2)+w_n^{k+n/2}A_{odd}((w_n^{k+n/2})^2)=A_{even}(w_{n/2}^k)-w_n^kA_{odd}(w_{n/2}^k) ]

    所以这个过程是可以递归的,且(T(n)=2T(n/2)+O(n))

    递归的FFT(c++代码)

    typedef complex<double> CD;
    const double pi = acos(-1);
    CD tmp[N],epsilon[N];
    void init_epsilon(int n){
        for(int i = 0; i < n; ++i){
            epsilon[i] = CD(cos(2.0 * pi * i / n), sin(2.0 * pi * i / n)); 
            arti_epsilon[i] = conj(epsilon[i]);
        }
    }
    void recursive_fft(int n, CD* A,int offset, int step, CD* w){
        if(n==1)return;
        int m=n>>1;
        recursive_fft(m,A,offset,step<<1,w);
        recursive_fft(m,A,offset+step,step<<1,w);
        for(int k=0;k<m;++k){
            int pos=2*step*k;
            tmp[k]  =A[pos+offset]+w[k*step]*A[pos+offset+step];
            tmp[k+m]=A[pos+offset]-w[k*step]*A[pos+offset+step];
        }
        for(int i=0;i<n;++i)
            A[i*step+offset]=tmp[i];
    }
    

    迭代实现

    但是递归需要比较大的空间,如何实现迭代的写法?

    观察递归的过程,第一步:

    0(000)2(010)4(100)6(110),1(001)3(011)5(101)7(111)

    第二步:

    0 (000)4(100),2(010) 6(110),1(001)5(101)3(011)7(111)

    将下标的二进制翻转过来就是:

    000,001,010,011,100,101,110,111

    对应了0,1,2,3,...

    用reverse(i)表示i的二进制翻转后的数。就相当于二进制从高位到低位的+1,是从左往右遇到第一个0,就改为1,左边的1改为0。

    int reverse(int x){
        for(int l=1<<bit_length;(x^=l)<l;l>>=l);
        return x;
    }
    

    我们得到递归最后一步的数组:

    0,4,2,6,1,5,3,7

    就可以从下到上迭代了。

    void bit_reverse(CD* A,int n){
        for(int i=0,j=0;i<n;++i){
            if(i>j)swap(A[i],A[j]);
            for(int l=n>>1;(j^=l)<l;l>>=1);
        }
    }
    void fft(CD* A, int n, CD* w){
        bit_reverse(A,n);
        for(int i=2;i<=n;i<<=1)//自下到上,i为每一层的步长,或者说子问题长度
          for(int j=0,m=i>>1;j<n;j+=i)//j为偏移量,或者说这一层的每一个子问题的起点
            for(int k=0;k<m;++k){//k=0..i/2,计算出子问题的第k个和第k+i/2个的值
              CD b=w[n/i*k]*A[j+m+k];
              A[j+m+k]=A[j+k]-b;
              A[j+k]+=b;
            }
    }
    

    离散傅里叶逆变换(IDFT)

    将点值转回系数表达式,相当于解方程组(vec{y}=V_nvec{a}),其中(vec {a})是系数向量。(V_n)是如下范德蒙德矩阵:

    [egin{bmatrix} 1&1&1&1&cdots &1\ 1&w;_n&w;_n^2&w;_n^3&cdots &w;_n^{n-1}\ 1&w;_n^2&w;_n^4&w;_n^6&cdots &w;_n^{2(n-1)}\ vdots&vdots&vdots&vdots&ddots&vdots\ 1&w;_n^{n-1}&w;_n^{2(n-1)}&w;_n^{3(n-1)}&cdots &w;_n^{(n-1)(n-1)}\ end{bmatrix} ]

    那么(V_n^{-1}vec{y}=vec{a})

    可令([V_n^{-1}]_{kj}=w^{-kj}_n/n),只需证明(V_n^{-1}V_n=I_n)

    [[V_n^{-1}V_n]_{jj'}=sum_{k=0}^{n-1}(w_n^{-kj}/n)(w_n^{kj'})=sum_{k=0}^{n-1}w_n^{k(j'-j)}/n = left{ egin{aligned} 0,j' eq j \ 1,j'=j end{aligned} ight. ]

    所以只要用(w_n^{-1})替换(w_n),并将结果每个元素除以n,就可以计算出(DDF_n^{-1})

    FFT解决高精度乘法,c++代码

    51 Nod 1028 大数乘法 V2

    注意FFT因为用了cos和sin,以及是浮点数计算,所以会有精度误差,故加了0.5。

    #include <bits/stdc++.h>
    using namespace std;
    #define rep(i,l,r) for(int i=l;i<r;++i)
    #define per(i,l,r) for(int i=r-1;i>=l;--i)
    #define SZ(x) ((int)(x).size())
    
    typedef double dd;
    typedef complex<dd> CD;
    const dd PI=acos(-1.0);
    const int L=18,N=1<<L;
    
    CD eps[N],inv_eps[N],f[N],g[N];
    void init_eps(int p){
        rep(i,0,p)eps[i]=CD(cos(PI*i*2/p),sin(PI*i*2/p)),inv_eps[i]=conj(eps[i]);
    }
    void fft(CD p[], int n, CD w[]){
        for(int i=0,j=0;i<n;++i){
            if(i>j)swap(p[i],p[j]);
            for(int l=n>>1;(j^=l)<l;l>>=1);
        }
        for(int i=2;i<=n;i<<=1)
            for(int j=0,m=i>>1;j<n;j+=i)
                rep(k,0,m){
                    CD b=w[n/i*k]*p[j+m+k];
                    p[j+m+k]=p[j+k]-b;
                    p[j+k]+=b;
                }
    }
    int ans[N];
    int main(){
        string a,b;
        cin>>a>>b;
        int n=max(SZ(a),SZ(b)),p=1;
        while(p<n)p<<=1;p<<=1;
        rep(i,0,p)f[i]=g[i]=0;
        n=0;per(i,0,SZ(a))f[n++]=a[i]-'0';
        n=0;per(i,0,SZ(b))g[n++]=b[i]-'0';
    
        init_eps(p);
        fft(f,p,eps);fft(g,p,eps);
        rep(i,0,p)f[i]*=g[i];
        fft(f,p,inv_eps);
    
        int t=0;
        rep(i,0,p){
            ans[i]=t+(f[i].real()+0.5)/p;
            if(ans[i]>9){t=ans[i]/10;ans[i]%=10;}
            else t=0;
        }
        bool flag=0;
        per(i,0,p)if(ans[i]||flag){
            printf("%d",ans[i]);flag=1;
        }
        if(flag==0)puts("0");
        return 0;
    }
    

    数论变换(NTT)

    为了解决FFT产生的精度误差,就有了一种在模意义下的方法,运算过程只有整数。它就是NTT(Number-Theoretic Transform)。

    原根

    设m是正整数,a是整数,则使得同余式(a^r equiv 1 mod m) 成立的最小正整数 (r) 叫做 (a) 对于模 (m) 的指数。

    若a模m的指数等于φ(m)(欧拉函数),则称a为模m的一个原根。

    原根的性质

    1. 对于质数p来说原根 (g) 满足 (g^0,g^1,g^2,…,g^{p-1}) 构成模 (p) 的简化剩余系。

    我们令(p=ccdot 2^k+1) ,那么对于2的幂n有 (n|(p-1))。设

    [g_n=g^{frac {p-1}n} ]

    考虑FFT中,需要用到单位根的以下性质:

    1. (w_n^k(0le k < n))互不相同,保证点值表示的合法;

    2. (omega_{n} ^ {2k} = omega_{n/2} ^k),用于分治;

    3. (omega_n ^ { k + frac{n}{2} } = -omega_n ^ k),用于分治;

    4. 仅当(k eq 0)时,(sum_{j=0}^{n-1}(w_n^k)^j=0),用于逆变换。

    原根是否也有以上性质呢?

    (p=ccdot 2^k+1),取(n=2^m),则(g_n^k(0le k<n)=g^{frac {(p-1)k}="" n}="g0,g{ccdot" 2{k-m}},g{2ccdot="" 2{k-m}}…,g{(p-1){ccdot="" 2^{k-m}}})<="" span="">,互不相等。</n)=g^{frac>

    这里有个小知识点,若$a_1,a_2,..a_p$ 是模p的完全剩余系,那么$aa_1+b,aa_2+b,..aa_p+b$也是模p的完全剩余系,证明:若$aa_j+baequiv a_i+b mod p$则$a_i equiv a_j mod p$,矛盾。

    (g_n^{2k}=g^{frac {2k(p-1)} n}=g^{frac{k(p-1)}{n/2}}=g_{n/2}^k)

    (g_n^n=g^{p-1}equiv 1mod pRightarrow g^{frac {p-1}2}equivsqrt 1 mod p),因为(g^k)互不相同,所以(g^{frac {p-1}2}equiv -1mod p),于是(g^{k+frac n 2}_n=g_n^kcdot g_n^{frac n 2}=g_n^kcdot g^{frac {p-1} 2}=-g_n^k)

    1. (k eq 0)时,(sum_{j=0}^{n-1}(g_n^k)^j=frac{1-(g_n^k)^n}{1-g_n^k}),由3得,(g_n^n=1),所以(1-(g_n^k)^n=1-(g_n^n)^k=0)。否则等于n。

    因此在FFT中用原根代替单位根就是NTT了。如果要求模数任意,只要再用 中国剩余定理(CRT)合并即可。

    NTT的c++代码

    仍然是上面高精度乘法那题

    #include <bits/stdc++.h>
    using namespace std;
    #define rep(i,l,r) for(int i=l;i<r;++i)
    #define per(i,l,r) for(int i=r-1;i>=l;--i)
    #define SZ(x) ((int)(x).size())
    
    typedef long long LL;
    const int L=18,N=1<<L;
    
    const LL C = 479;
    const LL P = (C << 21) + 1;
    const LL G = 3;
    
    LL qpow(LL a, LL b, LL m){
        LL ans = 1;
        for(a%=m;b;b>>=1,a=a*a%m)if(b&amp;1)ans=ans*a%m;
        return ans;
    }
    
    LL eps[N],inv_eps[N],f[N],g[N];
    void init_eps(int n){
        LL t=(P-1)/n, invG=qpow(G,P-2,P);
        rep(i,0,n) eps[i]=qpow(G,t*i,P),inv_eps[i]=qpow(invG,t*i,P);
    }
    void fft(LL p[], int n, LL w[]){
        for(int i=0,j=0;i<n;++i){
            if(i>j)swap(p[i],p[j]);
            for(int l=n>>1;(j^=l)<l;l>>=1);
        }
        for(int i=2;i<=n;i<<=1)
            for(int j=0,m=i>>1;j<n;j+=i)
                rep(k,0,m){
                    LL b=w[n/i*k]*p[j+m+k]%P;
                    p[j+m+k]=(p[j+k]-b+P)%P;
                    p[j+k]=(p[j+k]+b)%P;
                }
    }
    LL ans[N];
    int main(){
        string a,b;
        cin>>a>>b;
        int n=max(SZ(a),SZ(b)),p=1;
        while(p<n)p<<=1;p<<=1;
        rep(i,0,p)f[i]=g[i]=0;
        n=0;per(i,0,SZ(a))f[n++]=a[i]-'0';
        n=0;per(i,0,SZ(b))g[n++]=b[i]-'0';
    
        init_eps(p);
        fft(f,p,eps);fft(g,p,eps);
        rep(i,0,p)f[i]=f[i]*g[i]%P;
        fft(f,p,inv_eps);
    
        int t=0;
        LL invp=qpow(p,P-2,P);
        rep(i,0,p){
            ans[i]=(t+f[i]*invp%P)%P;
            if(ans[i]>9){t=ans[i]/10;ans[i]%=10;}
            else t=0;
        }
        bool flag=0;
        per(i,0,p)if(ans[i]||flag){
            printf("%lld",ans[i]);flag=1;
        }
        if(flag==0)puts("0");
        return 0;
    }
    
  • 相关阅读:
    el-select下拉框选项太多导致卡顿,使用下拉框分页来解决
    vue+elementui前端添加数字千位分割
    Failed to check/redeclare auto-delete queue(s)
    周末啦,做几道面试题放松放松吧!
    idea快捷键
    解决flink运行过程中报错Could not allocate enough slots within timeout of 300000 ms to run the job. Please make sure that the cluster has enough resources.
    用.net平台实现websocket server
    MQTT实战3
    Oracle 查看当前用户下库里所有的表、存储过程、触发器、视图
    idea从svn拉取项目不识别svn
  • 原文地址:https://www.cnblogs.com/flipped/p/9669233.html
Copyright © 2011-2022 走看看