zoukankan      html  css  js  c++  java
  • 分治FFT&&多项式求逆学习笔记

    分治FFT:

    解决的是形似以下的问题:

    给定n次多项式(g(x)),求多项式(f(x)),其中(f)的第(i)项系数的表达式为(sum f_{j} imes g_{i-j})

    解法:

    不难发现式子也是卷积的形式,但是与普通多项式乘法不一样的是,每一项的系数依赖前面的项的系数,使得普通的FFT无法起作用。

    考虑分治,将区间([l,r])分为([l,mid],[mid+1,r])两个区间计算,计算完([l,mid])中的多项式的系数之后,可以很方便的将([l,mid])([mid+1,r])的贡献给算出来,贡献的计算同上面的递推式。

    由于贡献的计算也需要依赖前面的系数,所以分治时的进程是:
    divide(l,mid);
    calc([l,mid]->[mid+1,r]);
    divide(mid+1,r);

    不难发现每个点之前的部分会被拆成若干个区间累加到这个位置上,所以正确性显而易见。

    /*=======================================
     * Author : ylsoi
     * Time : 2019.1.29
     * Problem : luogu4721
     * E-mail : ylsoi@foxmail.com
     * ====================================*/
    #include<bits/stdc++.h>
    
    #define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
    #define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
    #define debug(x) cout<<#x<<"="<<x<<" "
    #define fi first
    #define se second
    #define mk make_pair
    #define pb push_back
    typedef long long ll;
    
    using namespace std;
    
    void File(){
        freopen("luogu4721.in","r",stdin);
        freopen("luogu4721.out","w",stdout);
    }
    
    template<typename T>void read(T &_){
        _=0; T fl=1; char ch=getchar();
        for(;!isdigit(ch);ch=getchar())if(ch=='-')fl=-1;
        for(;isdigit(ch);ch=getchar())_=(_<<1)+(_<<3)+(ch^'0');
        _*=fl;
    }
    
    const int maxn=1e5+10;
    const int mod=998244353;
    
    ll qpow(ll x,ll y){
        ll ret=1; x%=mod;
        while(y){
            if(y&1)ret=ret*x%mod;
            x=x*x%mod;
            y>>=1;
        }
        return ret;
    }
    
    struct ksslbh{
        int lim,cnt,dn[maxn<<2];
        ll g[maxn<<2],ig[maxn<<2];
        void ntt(ll *A,int ty){
            REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
            for(int len=1;len<lim;len<<=1){
                ll w= ty==1 ? g[len<<1] : ig[len<<1];
                for(int L=0;L<lim;L+=len<<1){
                    ll wk=1;
                    REP(i,L,L+len-1){
                        ll u=A[i],v=A[i+len]*wk%mod;
                        A[i]=(u+v)%mod;
                        A[i+len]=(u-v)%mod;
                        wk=wk*w%mod;
                    }
                }
            }
        }
        void work(int n,int m,ll *a,ll *b){
            lim=1,cnt=0;
            while(lim<=n+m)lim<<=1,++cnt;
            if(!cnt)cnt=1;
            REP(i,0,lim){
                dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
                if(i>n)a[i]=0;
                if(i>m)b[i]=0;
            }
            g[lim]=qpow(3,(mod-1)/lim);
            ig[lim]=qpow(g[lim],mod-2);
            for(int i=lim>>1;i;i>>=1){
                g[i]=g[i<<1]*g[i<<1]%mod;
                ig[i]=ig[i<<1]*ig[i<<1]%mod;
            }
            ntt(a,1),ntt(b,1);
            REP(i,0,lim-1)a[i]=a[i]*b[i]%mod;
            ntt(a,-1);
            ll inv=qpow(lim,mod-2);
            REP(i,0,lim-1)a[i]=a[i]*inv%mod;
        }
    }T;
    
    int n;
    ll f[maxn<<2],g[maxn<<2],a[maxn<<2],b[maxn<<2];
    
    void divide(int l,int r){
        if(l==r)return;
        int mid=(l+r)>>1;
        divide(l,mid);
        int nn=mid-l,mm=r-l-1;
        //nn means the pows of f/a
        //mm means the pows of g/b
        REP(i,0,nn)a[i]=f[l+i];
        REP(i,0,mm)b[i]=g[i+1];
        T.work(nn,mm,a,b);
        REP(i,mid+1,r)f[i]=(f[i]+a[i-1-l])%mod;
        divide(mid+1,r);
    }
    
    int main(){
        //File();
        read(n);
        REP(i,1,n-1)read(g[i]);
        f[0]=1;
        divide(0,n);
        REP(i,0,n-1)printf("%lld ",(f[i]+mod)%mod);
        return 0;
    }
    
    

    多项式求逆:

    对于一个多项式(A(x)),找到一个多项式(A^{-1}(x)),使得(A(x) imes A^{-1}(x))只有常数项的系数为1,其余项的系数全部为0,那么我们称(A^{-1}(x))为多项式(A(x))的逆。

    当常数项存在逆元的时候,利用递推的方法去求解多项式的逆,不难发现所得的多项式是一个无限项的级数,所以我们一般只取前面(n-1)次,即对(x^n)取模。

    暴力求法:

    从低次依次向高次递推,每次求解(n)次项的系数(a_n)时,将前(n)项带入,所得的和需为(0),直接解出来即可。

    优化:

    观察上面的递推过程,不难发现后面的项求解需要用到前面的项的信息,利用分治FFT可以优化到(nlog^2n)

    更加优秀的解法:

    直接一项一项地递推显然太浪费时间,考虑利用倍增的思想来利用前面的项推出后面的项。

    假设我们现在得到了前(frac{n}{2})项的系数,那么我们就要想办法得出前(n)项的系数。

    假设前(frac{n}{2})项的多项式为(B(x)),最终的多项式为(A^{-1}(x)),那么我们可以得到

    [A^{-1}(x)-B(x)equiv 0 mod x^{frac{n}{2}} ]

    为了使模数提升至(x^{n}),可以将等式两边平方得:

    [(A^{-1}(x)-B(x))^2equiv 0 mod x^{n} ]

    考虑证明这一步,任意一个(i (i < n))次项的系数的拆分必定有一个数(< frac{n}{2}),而任意(< frac{n}{2})次项的系数都为0,所以(i (i<n))次项的系数都为(0)

    同时注意这种变化只有在等号右边等于(0)时才成立,这也是为什么要拿两式相减的理由。

    然后将平方项拆开:

    [A^{-2}(x)-2 imes A^{-1}(x) imes B(x)+B^{2}(x)equiv 0 mod x^n\ ]

    两边同时乘以(A(x))

    [A^{-1}(x)-2 imes B(x)+A(x) imes B^2(x)equiv 0 mod x^n\ A^{-1}(x)equiv 2 imes B(x)-A(x) imes B^{2}(x) mod x^n\ ]

    然后我们就成功地得到了在模(x^n)意义下的多项式(A^{-1}(x))

    具体实现:

    考虑递归实现上述过程。

    实现过程中,理论上(A(x))(B(x))项数都很多,但为了保证复杂度,(A(x))每次取对计算有贡献的前(n)项,而(B(x))原本是在模(x^frac{n}{2})意义下的多项式,只取前(frac{n}{2})项即可,所以最后乘出来的多项式为(2n)次的,时间复杂度满足(T(n)=nlog n+T(frac{n}{2}))

    所以最终时间复杂度为(nlog n)

    
     /*=======================================
     * Author : ylsoi
     * Time : 2019.2.2
     * Problem : luogu4238
     * E-mail : ylsoi@foxmail.com
     * ====================================*/
    #include<bits/stdc++.h>
    
    #define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
    #define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
    #define debug(x) cout<<#x<<"="<<x<<" "
    #define fi first
    #define se second
    #define mk make_pair
    #define pb push_back
    typedef long long ll;
    
    using namespace std;
    
    void File(){
        freopen("luogu4238.in","r",stdin);
        freopen("luogu4238.out","w",stdout);
    }
    
    template<typename T>void read(T &_){
    	_=0; T fl=1; char ch=getchar();
    	for(;!isdigit(ch);ch=getchar())if(ch=='-')fl=-1;
    	for(;isdigit(ch);ch=getchar())_=(_<<1)+(_<<3)+(ch^'0');
    	_*=fl;
    }
    
    const int maxn=1e5+10;
    const int mod=998244353;
    
    int lim,cnt,dn[maxn<<2];
    ll g[maxn<<2],ig[maxn<<2];
    
    ll qpow(ll x,ll y){
    	ll ret=1; x%=mod;
    	while(y){
    		if(y&1)ret=ret*x%mod;
    		x=x*x%mod;
    		y>>=1;
    	}
    	return ret;
    }
    
    void ntt(ll *A,int ty){
    	REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
    	for(int len=1;len<lim;len<<=1){
    		ll w= ty==1 ? g[len<<1] : ig[len<<1];
    		for(int L=0;L<lim;L+=len<<1){
    			ll wk=1;
    			REP(i,L,L+len-1){
    				ll u=A[i],v=A[i+len]*wk%mod;
    				A[i]=(u+v)%mod;
    				A[i+len]=(u-v)%mod;
    				wk=wk*w%mod;
    			}
    		}
    	}
    }
    
    int n;
    ll f[maxn<<2],a[maxn<<2],b[maxn<<2];
    
    void init(){
    	read(n);
    	REP(i,0,n-1)read(f[i]);
    	lim=1;
    	while(lim<=n+n)lim<<=1;
    	g[lim]=qpow(3,(mod-1)/lim);
    	ig[lim]=qpow(g[lim],mod-2);
    	for(int i=lim>>1;i;i>>=1){
    		g[i]=g[i<<1]*g[i<<1]%mod;
    		ig[i]=ig[i<<1]*ig[i<<1]%mod;
    	}
    }
    
    void poly_inv(int k){
    	if(k==1){
    		b[0]=qpow(f[0],mod-2);
    		return;
    	}
    	poly_inv((k+1)>>1);
    	lim=1,cnt=0;
    	while(lim<=k+k)lim<<=1,++cnt;
    	if(!cnt)cnt=1;
    	REP(i,0,lim-1)dn[i]=dn[i>>1]>>1|((i&1)<<(cnt-1));
    	REP(i,0,k-1)a[i]=f[i];
    	REP(i,k,lim-1)a[i]=0;
    	ntt(a,1),ntt(b,1);
    	REP(i,0,lim-1)b[i]=b[i]*(2ll-a[i]*b[i]%mod)%mod;
    	ntt(b,-1);
    	ll inv=qpow(lim,mod-2);
    	REP(i,0,k-1)b[i]=b[i]*inv%mod;
    	REP(i,k,lim)b[i]=0;
    }
    
    int main(){
        File();
    	init();
    	poly_inv(n);
    	REP(i,0,n-1)printf("%lld ",(b[i]%mod+mod)%mod);
        return 0;
    }
    
    
    
  • 相关阅读:
    [Debug] Make python2.7 use the right version of opencv
    路线图 | 学习OpenCV路线图
    学习笔记 | Princeton
    书单 | 2017年阅读书单
    路线图 | 摄影师成长路线
    学习笔记 | Morvan
    如何在pycharm中进入shell脚本调试代码
    python3运行报错:TypeError: Object of type 'type' is not JSON serializable解决方法(详细)
    动态规划的引入 P1616 疯狂的采药【完全背包】
    动态规划的引入 P1048 采药【01背包】
  • 原文地址:https://www.cnblogs.com/ylsoi/p/10351114.html
Copyright © 2011-2022 走看看