zoukankan      html  css  js  c++  java
  • 多项式全家桶学习笔记

    多项式全家桶学习笔记

    Upd:

    (2020/2/9),做一道题被卡常后发现有些limit可以改的更小......(详见代码)

    (2020/2/5),该了一下(exp)的实现,直接用上取整递归了,好像比原来快乐一倍

    终于肝完了...... 多项式三角函数就不学了吧qwq


    多项式求逆

    还是看板子:【模板】多项式乘法逆

    给一个(n-1)(n)项柿(F(x)),要你求一个(n-1)次多项式(G(x)),满足(F(x)G(x)equiv 1 (mod x^n))

    就是把(F(x)G(x))卷积起来忽略掉次数(ge n)的项后它(equiv 1)

    一个比较难的情况:(n = 1),即(F(x)G(x))的常数项为(1),答案就是(F[0]^{-1}),((F)的常数项的逆元),怎么样,难吧!

    好我们下面来看更一般的情况

    [F(x)G(x) equiv 1 (mod x^n) ]

    假设我们现在已经知道了

    [F(x)G'(x) equiv 1 (mod x^{leftlceil frac{n}{2} ight ceil}) ]

    那么由于(F(x)G(x) equiv 1 (mod x^n)),所以(F(x)G(x))必定(equiv 1 (mod x^{leftlceil frac{n}{2} ight ceil})),所以两式相减得

    [F(x)(G(x) - G'(x)) equiv 0 (mod x^{leftlceil frac{n}{2} ight ceil}) ]

    由于(F(x) ot= 0),所以

    [G(x) - G'(x) equiv 0 (mod x^{leftlceil frac{n}{2} ight ceil}) ]

    然后发现我们回不上去了23333...

    然而在这里我们可以直接平方一下

    [(G(x) - G'(x))^2 equiv 0 (mod x^n) ]

    为什么呢?

    分类讨论一下

    • 对于次数小于(leftlceil frac{n}{2} ight ceil)的项,它不管乘什么都是(0)
    • 对于次数大于(leftlceil frac{n}{2} ight ceil)的项,它只有乘一个次数小于(leftlceil frac{n}{2} ight ceil)的项才会对上面那个恒等式产生影响,显然这也是(0)

    我们继续化简,暴力展开

    [G(x)^2 + G'(x)^2 - 2G(x)G'(x) equiv 0 (mod x^n) ]

    因为我们知道(F(x)G(x) equiv 1 (mod x^n)),两边乘(F(x))

    [G(x) + G'(x)^2F(x) - 2G'(x) equiv 0 (mod x^n) ]

    移项得

    [G(x) equiv 2G'(x) - G'(x) ^ 2 F(x) (mod x^n) ]

    为了好看,我们可以更简单地提一个(G'(x))出来

    [G(x) equiv G'(x)(2 - G'(x)F(x)) (mod x^n) ]

    顺着上面那个柿子递归用(NTT)算就好了。

    复杂度:听别人说是 (O(n log n))

    #include <bits/stdc++.h>
    using namespace std;
    const int N=1e6+10,P=998244353,G=3,IG=(P+1)/G;
    inline int fpow(int x,int y){
        int ret=1; for (;y;y>>=1,x=1ll*x*x%P) if (y&1) ret=1ll*ret*x%P;
        return ret;	
    }
    inline int add(int x,int y){return x+y>=P?x+y-P:x+y;}
    inline int sub(int x,int y){return x-y<0?x-y+P:x-y;}
    int rev[N];
    void init(int len){
        for (int i=0;i<len;i++) rev[i]=rev[i>>1]>>1|((i&1)?len>>1:0);
    }
    void ntt(int *f,int n,int flg){
        for (int i=0;i<n;i++) if(rev[i]<i) swap(f[i],f[rev[i]]);
        for (int len=2,k=1;len<=n;len<<=1,k<<=1){
            int wn=fpow(flg==1?G:IG,(P-1)/len);
            for (int i=0;i<n;i+=len){
                for (int w=1,j=i;j<i+k;j++,w=1ll*w*wn%P){
                    int tmp=1ll*w*f[j+k]%P;
                    f[j+k]=sub(f[j],tmp),f[j]=add(f[j],tmp);
                }
            }
        }
    }
    int FF[N];
    void getinv(int *F,int *G,int n){
        if (n==1){G[0]=fpow(F[0],P-2);return;}
        getinv(F,G,(int)ceil(n/2.0));
        int limit=1; while (limit<=2*(n-1))limit<<=1;
        init(limit);
        for (int i=0;i<n;i++) FF[i]=F[i];
        for (int i=n;i<limit;i++) FF[i]=0;
        ntt(FF,limit,1),ntt(G,limit,1);
        for (int i=0;i<limit;i++) G[i]=1ll*sub(2,1ll*FF[i]*G[i]%P)*G[i]%P;
        ntt(G,limit,-1); int inv=fpow(limit,P-2);
        for (int i=0;i<limit;i++) G[i]=1ll*G[i]*inv%P;
        for (int i=n;i<limit;i++) G[i]=0;
    }
    int f[N],inv[N];
    int main(){
        int n;scanf("%d",&n);
        for (int i=0;i<n;i++)scanf("%d",&f[i]);
        getinv(f,inv,n);
        for (int i=0;i<n;i++)printf("%d ",inv[i]);
        return 0;
    }
    

    ps: 其实还有个迭代版的......尝试写了一下......绝对邪教......,总之这样也不慢。


    多项式开根

    【模板】多项式开根

    这个做法有点像多项式乘法逆,考虑倍增。

    我们要求的是一个多项式(G(x)),满足

    [G^2(x)equiv F(x) (mod x^n) ]

    假设我们现在已经知道了

    [H^2(x) equiv F(x) (mod x^{leftlceil frac{n}{2} ight ceil}) ]

    显然有

    [G^2(x) - H^2(x) equiv 0 (mod x^{leftlceil frac{n}{2} ight ceil}) ]

    平方差一下

    [(G(x) - H(x))(G(x) + H(x)) equiv 0 (mod x^leftlceilfrac{n}{2} ight ceil) ]

    这时候可以发现(G(x))应该有两个解,但在某些题目中我们并不希望(G(x))出现负数,所以我们不妨令

    [G(x)-H(x) equiv 0 (mod x^leftlceilfrac{n}{2} ight ceil) ]

    套路的平方一下

    [G^2(x) + H^2(x) - 2G(x)H(x) equiv 0 (mod x^n) ]

    发现(G^2(x))就是(F(x)),然后再移项

    [F(x) + H^2(x) equiv 2G(x)H(x) (mod x^n) ]

    那么

    [G(x) equiv frac{F(x) + H^2(x)}{2H(x)} (mod x^n) ]

    (2H(x))求逆后(NTT)就行了。

    (Code:)

    #include <bits/stdc++.h>
    using namespace std;
    const int N=3e5+10,P=998244353,g=3,ig=(P+1)/g;
    inline int add(int x,int y){return x+y>=P?x+y-P:x+y;}
    inline int sub(int x,int y){return x-y<0?x-y+P:x-y;}
    inline int sqr(int x){return 1ll*x*x%P;}
    inline int fpow(int x,int y){
        int ret=1; for (x%=P;y;y>>=1,x=1ll*x*x%P)
            if (y&1) ret=1ll*ret*x%P;
        return ret;
    }
    namespace Poly{
        int rev[N];
        void init(int limit){
            for (int i=0;i<limit;i++)
                rev[i]=rev[i>>1]>>1|((i&1)?limit>>1:0);
        }
        void ntt(int *f,int n,int flg){
            for (int i=0;i<n;i++)
                if (rev[i]<i) swap(f[i],f[rev[i]]);
            for (int k=1,len=2;len<=n;len<<=1,k<<=1){
                int wn=fpow(flg==1?g:ig,(P-1)/len);
                for (int i=0;i<n;i+=len){
                    for (int w=1,j=i;j<i+k;j++,w=1ll*w*wn%P){
                        int tmp=1ll*w*f[j+k]%P;
                        f[j+k]=sub(f[j],tmp),f[j]=add(f[j],tmp);
                    }
                }
            }
            if (flg!=1){
                int inv=fpow(n,P-2);
                for (int i=0;i<n;i++) f[i]=1ll*f[i]*inv%P;
            }
        }
        void getinv(int *f,int n,int *inv){
            if (n==1){inv[0]=fpow(f[0],P-2);return;}
            getinv(f,(n+1)>>1,inv);
            static F[N];
            int limit=1; while (limit<=2*(n-1))limit<<=1; init(limit);
            for (int i=0;i<n;i++) F[i]=f[i];
            for (int i=n;i<limit;i++) F[i]=inv[i]=0;
            ntt(F,limit,1),ntt(inv,limit,1);
            for (int i=0;i<limit;i++) inv[i]=1ll*inv[i]*sub(2,1ll*F[i]*inv[i]%P)%P;
            ntt(inv,limit,-1);
            for (int i=n;i<limit;i++) inv[i]=0;
        }
        void getsqrt(int *f,int n,int *sqt){
            if (n==1){sqt[0]=1;return;}
            getsqrt(f,(n+1)>>1,sqt);
            static int H[N],iH[N];
            int limit=1; while (limit<=2*(n-1)) limit<<=1;
            for (int i=0;i<limit;i++) H[i]=i>=n?0:2ll*sqt[i]%P;
            getinv(H,n,iH),init(limit);
            for (int i=0;i<limit;i++) F[i]=i>=n?0:f[i],sqt[i]=i>=n?0:sqt[i];
            ntt(F,limit,1),ntt(sqt,limit,1),ntt(iH,limit,1);
            for (int i=0;i<limit;i++) sqt[i]=1ll*add(F[i],sqr(sqt[i]))*iH[i]%P;
            ntt(sqt,limit,-1);
            for (int i=n;i<limit;i++) sqt[i]=0;
    //		cout<<n<<" wtf: "; for (int i=0;i<n;i++) cout<<sqt[i]<<" "; cout<<endl;
        }
    }
    int f[N],sqt[N];
    int main(){
        int n; scanf("%d",&n);
        for (int i=0;i<n;i++) scanf("%d",&f[i]);
        Poly::getsqrt(f,n,sqt);
        for (int i=0;i<n;i++) printf("%d ",sqt[i]);
        return 0;
    }
    

    多项式带余除法

    【模板】多项式除法

    给一个(n)次多项式(F(x)),和一个(m)次多项式(G(x)),求多项式(Q(x),R(x))满足

    • (Q(x))次数为(n-m),(R(x))的次数小于(m)
    • (F(x) = Q(x)G(x) + R(x))

    (998244353)取模。

    如果没有余数的话我们直接多项式求逆就完了,但这里有(R(x)),我们考虑把它消掉。

    之前求逆的时候会对(x^n)取模,这样我们可以消去一些高次项,但我们想的是要把后面(R(x))消掉,然后保留原来的柿子。

    下面有一些非常妙的做法。

    我们把(x^{-1})代入(F(x)),这样次数就都是负的,然后把他乘上(x^n),即

    [x^nF(x^{-1}) ]

    容易发现这样即翻转(F(x))的系数,把原来高次项放到后面去了,不放叫他(F^T(x) = x^nF(x^{-1}))

    下面来推一下

    (x^{-1})代替(x)

    [F(x^{-1}) = Q(x^{-1})G(x^{-1}) + R(x^{-1}) ]

    两边乘上(x^n)

    [x^nF(x^{-1}) = x^nQ(x^{-1})G(x^{-1}) + x^nR(x^{-1}) ]

    注意到(Q(x)G(x))的次数为(x^n)(R(x))的次数最大为(m-1),所以

    [F^T(x) = Q^T(x)G^T(x) + x^{n-m+1}R^T(x) ]

    发现了什么?(Q)的次数为(n-m),后面又有(x^{n-m+1}),所以我们机智的对(x^{n-m+1})取模,就消掉了(R(x))

    [F^T(x) equiv Q^T(x)G^T(x) (mod x^{n-m+1}) ]

    通过对(G^T)求逆元我们就能算出(Q^T),然后翻转就得到了(Q),带回去算(R)就好了。

    [Q^T(x) equiv frac{F^T(x)}{G^T(x)} (mod x^{n-m+1}) ]

    然后再回去用(F - Q)(R)就行了。

    /*
     *没有的部分参照上面
    */
    namespace Poly{
        void van(int *f,int n){
            for (int i=0,j=n;i<j;i++,j--)swap(f[i],f[j]);
        }
        void div(int *f,int n,int *g,int m,int *Q,int *R){
            static int F[N],G[N],iG[N],qwq[N];
            for (int i=0;i<=n;i++) F[i]=f[i];
            for (int i=0;i<=m;i++) G[i]=g[i];
            van(F,n),van(G,m); getinv(G,n-m+1,iG);
    /*		cout<<"F: "; for (int i=0;i<=n;i++) cout<<F[i]<<" "; cout<<endl;
            cout<<"G: "; for (int i=0;i<=m;i++) cout<<G[i]<<" "; cout<<endl;
            cout<<"iG: "; for (int i=0;i<=n-m+1;i++) cout<<iG[i]<<" "; cout<<endl; */
            int limit=1; while (limit<=2*n)limit<<=1; init(limit);
            for (int i=n+1;i<limit;i++) F[i]=0;
            for (int i=n-m+1;i<limit;i++) iG[i]=0;
            ntt(F,limit,1),ntt(iG,limit,1);
            for (int i=0;i<limit;i++) Q[i]=1ll*F[i]*iG[i]%P;
            ntt(Q,limit,-1);
            for (int i=n-m+1;i<limit;i++) Q[i]=0;
            van(Q,n-m);
            for (int i=0;i<=n-m;i++) qwq[i]=Q[i];
            for (int i=n-m+1;i<limit;i++) qwq[i]=0;
            for (int i=0;i<=m;i++) G[i]=g[i];
            for (int i=m+1;i<limit;i++) G[i]=0;
            limit=1; while(limit<=n)limit<<=1; init(limit);
            ntt(qwq,limit,1),ntt(G,limit,1);
            for (int i=0;i<limit;i++) G[i]=1ll*G[i]*qwq[i]%P;
            ntt(G,limit,-1);
            for (int i=n+1;i<limit;i++) G[i]=0;
            for (int i=0;i<m;i++) R[i]=sub(f[i],G[i]);
        }
    }
    int n,m,f[N],g[N],q[N],r[N];
    int main(){
        scanf("%d%d",&n,&m);
        for (int i=0;i<=n;i++) scanf("%d",&f[i]);
        for (int i=0;i<=m;i++) scanf("%d",&g[i]);
        Poly::div(f,n,g,m,q,r);
        for (int i=0;i<=n-m;i++) printf("%d ",q[i]); puts("");
        for (int i=0;i<m;i++) printf("%d ",r[i]);
        return 0;
    }
    

    多项式(ln)

    【模板】多项式对数函数

    ln是啥??

    哦~不就是自然对数吗。等等......(e)在同余系下是啥?咋就能取对数了?

    推荐个视频教程:【MIT公开课】单变量微积分

    Two Hours Later......


    多项式求导 & 积分

    首先我们知道可以对每一项分开来求导、积分。

    那我们就来看(a_nx^n),先(a^n)是常数,把他扔出去。

    导数就是((a_nx^n)' = a_n(x^n)' = a_nnx^{n-1})

    积分出来(int a_nx^ndx = a_nint x^n dx = a_nfrac{1}{n+1}x^{n+1})

    (Code)

    namespace Poly{
        void dao(int *f,int n,int *d){
            static int F[N]; for (int i=0;i<=n;i++) F[i]=f[i];
            for (int i=1;i<=n;i++) d[i-1]=1ll*F[i]*i%P;	d[n]=0;
        }
        void jifen(int *f,int n,int *jf){
            static int F[N]; for (int i=0;i<=n;i++) F[i]=f[i];
            for (int i=0;i<=n;i++) 
                jf[i+1]=1ll*F[i]*fpow(i+1,P-2)%P;
            jf[0]=0;
        }
    }
    

    这个(Ln)不是直接算出来的,是求导求出来的......

    首先我们知道((ln x)' = frac{1}{x}),现在我们要求的是

    [ln F(x) = G(x) ]

    咋办?两边同时对(x)求导!(右边要用一下链式法则)

    [F^{-1}(x)F'(x) = G'(x) ]

    然后两边同时积分

    [G(x) = int F^{-1}(x)F'(x)dx ]

    照着求导、求逆、乘法、再积回去就没了。

    #include <bits/stdc++.h>
    using namespace std;
    const int N=3e5+10,P=998244353,gen=3,igen=(P+1)/gen;
    inline int add(int x,int y){
        return x+y>=P?x+y-P:x+y;
    }
    inline int sub(int x,int y){
        return x-y<0?x-y+P:x-y;
    }
    inline int fpow(int x,int y){
        int ret=1; for (x%=P;y;y>>=1,x=1ll*x*x%P)
            if (y&1) ret=1ll*ret*x%P;
        return ret;
    }
    namespace Poly{
        int rev[N];
        void init(int limit){
            for (int i=0;i<limit;i++)
                rev[i]=rev[i>>1]>>1|((i&1)?limit>>1:0);
        }
        void dao(int *f,int n,int *d){
            static int F[N]; for (int i=0;i<=n;i++) F[i]=f[i];
            for (int i=1;i<=n;i++) d[i-1]=1ll*F[i]*i%P;	d[n]=0;
        }
        void jifen(int *f,int n,int *jf){
            static int F[N]; for (int i=0;i<=n;i++) F[i]=f[i];
            for (int i=0;i<=n;i++) 
                jf[i+1]=1ll*F[i]*fpow(i+1,P-2)%P;
            jf[0]=0;
        }
        void ntt(int *f,int n,int flg){
            for (int i=0;i<n;i++)
                if (rev[i]<i) swap(f[i],f[rev[i]]);
            for (int len=2,k=1;len<=n;len<<=1,k<<=1){
                int wn=fpow(flg==1?gen:igen,(P-1)/len);
                for (int i=0;i<n;i+=len){
                    for (int j=i,w=1;j<i+k;j++,w=1ll*w*wn%P){
                        int tmp=1ll*f[j+k]*w%P;
                        f[j+k]=sub(f[j],tmp),f[j]=add(f[j],tmp);
                    }
                }
            }
            if (flg==-1){
                int inv=fpow(n,P-2);
                for (int i=0;i<n;i++) f[i]=1ll*f[i]*inv%P;
            }
        }
        void getinv(int *f,int n,int *G){
            if (n==1){G[0]=fpow(f[0],P-2);return;}
            getinv(f,(n+1)>>1,G); static int F[N];
            int limit=1; while(limit<=2*(n-1))limit<<=1; init(limit);
            for (int i=0;i<limit;i++)
                F[i]=i>=n?0:f[i],G[i]=i>=n?0:G[i];
            ntt(F,limit,1),ntt(G,limit,1);
            for (int i=0;i<limit;i++) G[i]=1ll*G[i]*sub(2,1ll*F[i]*G[i]%P)%P;
            ntt(G,limit,-1);
            for (int i=n;i<limit;i++) G[i]=0;
        }
        void getln(int *f,int n,int *G){
            static int F[N],iF[N]; for (int i=0;i<n;i++) F[i]=f[i];
            getinv(F,n,iF),dao(F,n-1,F);
            int limit=1; while(limit<=2*(n-1))limit<<=1; init(limit);
            for (int i=n-1;i<limit;i++) F[i]=0;
            ntt(F,limit,1),ntt(iF,limit,1);
            for (int i=0;i<limit;i++) G[i]=1ll*F[i]*iF[i]%P;
            ntt(G,limit,-1);
            jifen(G,n-1,G); for (int i=n;i<limit;i++)G[i]=0;
        }
    }
    int f[N],n,ln[N];
    int main(){
        scanf("%d",&n);
        for (int i=0;i<n;i++)scanf("%d",&f[i]);
        Poly::getln(f,n,ln);
        for (int i=0;i<n;i++) printf("%d ",ln[i]);
        return 0;
    }
    

    多项式牛顿迭代

    这个好像是最有用的(...)

    首先您需要知道啥是牛顿迭代

    下面的一些内容感性理解一下就好了......并不会太严谨的证明qwq。

    好像这个结论也挺好记的。


    牛顿迭代

    它可以用来求(f(x) = 0)近似解

    一个简单的应用就是开方。

    我们还是直接来看具体怎么做吧

    首先我们一眼猜到(f(x) = 0)的一个粗略解(x_0),然后我们故意把(f(x))当成一个一次函数,即作一条切(f(x))((x_0,f(x_0)))的直线(l),根据点斜式,这条直线可以被描述为(y - f(x_0) = f'(x_0)(x-x_0)),然后我们在这条直线上找与(x)轴的交点,即(y = 0)的点,记作((x, 0)),那么有

    [x = x_0 - frac{f(x_0)}{f'(x_0)} ]

    重复以上步骤,我们会不断逼近(f(x)=0)的解。

    即先从一个(x_0)出发,然后按一下递推式得到后面更精确的(x)

    [x_{n+1} = x_n - frac{f(x_n)}{f'(x_n)} ]

    放张(wiki)上的图

    举个例子,求解平方根。我们现在要求(sqrt k),那么函数就是(f(x) = x^2-k),我们要求的就是(f(x)=0)的解。

    (f(x))求导得到(f'(x) = 2x)。按柿子牛顿迭代就行了

    #include <bits/stdc++.h>
    using namespace std;
    typedef double db;
    db calc(db k){
        #define F(x) (1.0*(x)*(x)-k)
        #define F1(x) (2.0*(x))
        db x=1;
        for (int it=0;it<100;it++)
            x=x-F(x)/F1(x);
        return x;
    }
    int main(){
        db k; scanf("%lf",&k);
        printf("%.9f
    ",calc(k));
        return 0;
    }
    

    值得注意的是每牛顿迭代一次,精度都会翻一倍,这和二分有点像

    另外如果知道泰勒展开的话,上面的迭代也应该很自然。

    首先我们知道一个粗略解(x_0),那么(f(x))(x_0)处的泰勒展开结果

    [f(x) = sumlimits_{n=0}^{infty} frac{f^{(n)}(x_0)}{n!}(x - x_0)^n ]

    其中(f^{(n)})表示(f)(n)阶导数。因为我们可以不断逼近解,所以之取前两项也就是我们就直接认为(f(x) = f(x_0)+f'(x_0)(x-x_0)),也能得到牛顿迭代公式。


    对多项式而言,是这样的问题:

    给函数(F(x)),求一个多项式(G(x)),满足(F(G(x)) = 0),一般来讲,都会让你(mod)一个(x^n)(这个类似实数的精度)。

    这个可以类似上面的牛顿迭代。

    具体过程(假设(n)(2)的整数次幂):

    我们要求的是(F(G(x))equiv 0 (mod x^n)),因为牛顿迭代一次后精度会翻倍

    假设我们已经知道了(F(G^*(x)) equiv 0 (mod x^{frac{n}{2}}))

    那么根据泰勒展开,(F(G(x)))(G^*(x))处的展开式

    [F(G(x)) = frac{F(G^*(x))}{0!} + frac{F'(G^*(x))}{1!}(G(x)-G^*(x)) + frac{F^{''}(G^*(x))}{2!}(G(x)-G^*(x))^2 + cdots ]

    2333,咋整?不要慌张,注意到我们是对(x^n)取模,所以后面(G(x)-G^*(x))的次数大于(1)的项全部报废了。(Nice!)

    这样就只剩下

    [F(G(x)) = F(G^*(x)) + F'(G^*(x))(G(x)-G^*(x)) ]

    而我们要的是(F(G(x)) equiv 0 (mod x^n)),即

    [F(G^*(x)) + F'(G^*(x))G(x)-F'(G^*(x))G^*(x) equiv 0 (mod x^n) ]

    移个项

    [F'(G^*(x))G(x) = F'(G^*(x))G^*(x) - F(G^*(x)) ]

    两边同时除以(F'(G^*(x)))

    [G(x) equiv G^*(x) - frac{F(G^*(x))}{F'(G^*(x))} (mod x^n) ]

    所以上面的柿子就是结论,非常有用!


    多项式(exp)

    【模板】多项式指数函数(多项式 exp)

    给一个多项式(n-1)次多项式(A(x)),让你求(mod x^n)意义下的(e^{A(x)}),系数对(998244353)取模。

    这个板子也是非常nb啊......

    (e)是个无理数,我们还不知道它长啥样,考虑牛顿迭代。

    我们要求的是(B(x) = e^{A(x)}),两边(ln)一下,然后移项得到

    [ln B(x) - A(x) = 0 ]

    (F(G(x)) = ln G(x) - A(x)),按上面的柿子牛顿迭代就好了。

    (A(x))是常数,所以(F)的导数(F'(G(x)) = frac{1}{G(x)})

    假设我们知道了(F(H(x)) equiv 0 (mod x^{leftlceilfrac{n}{2} ight ceil})),有

    [G(x) = H(x) - frac{F(H(x))}{F'(H(x))} = H(x) - H(x)(ln H(x)-A(x)) = H(x)(1-ln H(x)+A(x)) ]

    递归计算就好了。s

    /*
    其他的看上面......
    */
    namespace Poly{
        void getexp(int *f,int n,int *G){
            if (n==1){G[0]=1;return;}
            getexp(f,(n+1)>>1,G);
            static int F[N],lnG[N];
            int limit=1; while(limit<=2*(n-1))limit<<=1; init(limit);
            for (int i=0;i<limit;i++) F[i]=i>=n?0:f[i],G[i]=i>=n?0:G[i];
            getln(G,n,lnG);
            ntt(G,limit,1),ntt(lnG,limit,1),ntt(F,limit,1);
            for (int i=0;i<limit;i++) G[i]=1ll*G[i]*sub(add(1,F[i]),lnG[i])%P;
            ntt(G,limit,-1);
            for (int i=n;i<limit;i++) G[i]=0;
        }
    }
    int f[N],ans[N],n;
    int main(){
        scanf("%d",&n);
        for (int i=0;i<n;i++)scanf("%d",&f[i]);
        Poly::getexp(f,n,ans);
        for (int i=0;i<n;i++)printf("%d ",ans[i]);
        return 0;
    }
    

    多项式快速幂

    【模板】多项式快速幂

    给一个(n-1)次多项式(A(x)),和一个(k(k le 10^{10^5})),求(A^k(x) mod x^n),系数对(998244353)取模。

    这个有了(ln)(exp)应该挺好做了吧。

    直接写柿子

    [A^k(x) = e^{k ln A(x)} ]

    注意到(k)会乘到(ln A(x))的系数里,边读边取模就好了。

    inline int getnum(){
        int x=0,f=1; char ch=getchar(); 
        while(!isdigit(ch)){f=ch=='-'?-f:f;ch=getchar();}
        while(isdigit(ch)){x=(10ll*x+ch-48)%P;ch=getchar();}
        return normal(x*f);
    }
    namespace Poly{
        void getpow(int *f,int n,int k,int *G){
            static int F[N]; for (int i=0;i<n;i++) F[i]=f[i];
            getln(F,n,F);
            for (int i=0;i<n;i++) F[i]=1ll*F[i]*k%P;
            getexp(F,n,G);
        }
    }
    int f[N],ans[N],n,k;
    int main(){
        scanf("%d",&n),k=getnum();
        for (int i=0;i<n;i++)scanf("%d",&f[i]);
        Poly::getpow(f,n,k,ans);
        for (int i=0;i<n;i++)printf("%d ",ans[i]);
        return 0;
    }
    

    总结

    这里放一下自己丑的一比的板子吧

    #include <bits/stdc++.h>
    using namespace std;
    const int N=1e6,P=998244353,gen=3,igen=(P+1)/gen;
    inline int add(int x,int y){
        return x+y>=P?x+y-P:x+y;
    }
    inline int sub(int x,int y){
        return x-y<0?x-y+P:x-y;
    }
    inline int fpow(int x,int y){
        int ret=1; for (x%=P;y;y>>=1,x=1ll*x*x%P)
            if (y&1) ret=1ll*ret*x%P;
        return ret;
    }
    inline int sqr(int x){return 1ll*x*x%P;}
    inline int normal(int x){return (x%P+P)%P;}
    inline int getnum(){
        int x=0,f=1; char ch=getchar(); 
        while(!isdigit(ch)){f=ch=='-'?-f:f;ch=getchar();}
        while(isdigit(ch)){x=(10ll*x+ch-48)%P;ch=getchar();}
        return normal(x*f);
    }
    namespace Poly{
        int rev[N];
        void init(int n){
            for (int i=0;i<n;i++)
                rev[i]=rev[i>>1]>>1|((i&1)?n>>1:0);
        }
        void ntt(int *f,int n,int flg){
            for (int i=0;i<n;i++)
                if (rev[i]<i) swap(f[i],f[rev[i]]);
            for (int len=2,k=1;len<=n;len<<=1,k<<=1){
                int wn=fpow(flg==1?gen:igen,(P-1)/len);
                for (int i=0;i<n;i+=len){
                    for (int w=1,j=i;j<i+k;j++,w=1ll*w*wn%P){
                        int tmp=1ll*w*f[j+k]%P;
                        f[j+k]=sub(f[j],tmp),f[j]=add(f[j],tmp);
                    }
                }
            }
            if (flg==-1){
                int inv=fpow(n,P-2);
                for (int i=0;i<n;i++) f[i]=1ll*f[i]*inv%P;
            }
        }
        void dao(int *f,int n,int *g){
            static int F[N]; for (int i=0;i<=n;i++)F[i]=f[i];
            for (int i=1;i<=n;i++) g[i-1]=1ll*i*F[i]%P; g[n]=0;
        }
        void jifen(int *f,int n,int *g){
            static int F[N]; for (int i=0;i<=n;i++)F[i]=f[i];
            for (int i=0;i<=n;i++) g[i+1]=1ll*fpow(i+1,P-2)*F[i]%P; g[0]=0;
        }
        void getinv(int *f,int n,int *G){
            if (n==1){G[0]=fpow(f[0],P-2);return;}
            getinv(f,(n+1)>>1,G); static int F[N];
            int limit=1; while(limit<=2*(n-1))limit<<=1; init(limit);
            for (int i=0;i<limit;i++) F[i]=i>=n?0:f[i],G[i]=i>=n?0:G[i];
            ntt(F,limit,1),ntt(G,limit,1);
            for (int i=0;i<limit;i++) G[i]=1ll*G[i]*sub(2,1ll*F[i]*G[i]%P)%P;
            ntt(G,limit,-1);
            for (int i=n;i<limit;i++) G[i]=0;
        }
        void van(int *f,int n){
            for (int i=0,j=n;i<j;i++,j--)swap(f[i],f[j]);
        }
        void div(int *f,int n,int *g,int m,int *Q,int *R){
            static int F[N],G[N],iG[N],qwq[N];
            for (int i=0;i<=n;i++) F[i]=f[i];
            for (int i=0;i<=m;i++) G[i]=g[i];
            van(F,n),van(G,m); getinv(G,n-m+1,iG);
    /*		cout<<"F: "; for (int i=0;i<=n;i++) cout<<F[i]<<" "; cout<<endl;
            cout<<"G: "; for (int i=0;i<=m;i++) cout<<G[i]<<" "; cout<<endl;
            cout<<"iG: "; for (int i=0;i<=n-m+1;i++) cout<<iG[i]<<" "; cout<<endl; */
            int limit=1; while (limit<=2*n)limit<<=1; init(limit);
            for (int i=n+1;i<limit;i++) F[i]=0;
            for (int i=n-m+1;i<limit;i++) iG[i]=0;
            ntt(F,limit,1),ntt(iG,limit,1);
            for (int i=0;i<limit;i++) Q[i]=1ll*F[i]*iG[i]%P;
            ntt(Q,limit,-1);
            for (int i=n-m+1;i<limit;i++) Q[i]=0;
            van(Q,n-m);
            for (int i=0;i<=n-m;i++) qwq[i]=Q[i];
            for (int i=n-m+1;i<limit;i++) qwq[i]=0;
            for (int i=0;i<=m;i++) G[i]=g[i];
            for (int i=m+1;i<limit;i++) G[i]=0;
            limit=1; while(limit<=n)limit<<=1; init(limit);
            ntt(qwq,limit,1),ntt(G,limit,1);
            for (int i=0;i<limit;i++) G[i]=1ll*G[i]*qwq[i]%P;
            ntt(G,limit,-1);
            for (int i=n+1;i<limit;i++) G[i]=0;
            for (int i=0;i<m;i++) R[i]=sub(f[i],G[i]);
        }
        void getln(int *f,int n,int *G){
            static int F[N],iF[N]; for (int i=0;i<n;i++) F[i]=f[i];
            getinv(F,n,iF),dao(F,n-1,F);
            int limit=1; while(limit<=2*(n-1))limit<<=1; init(limit);
            for (int i=n-1;i<limit;i++) F[i]=0;
            ntt(F,limit,1),ntt(iF,limit,1);
            for (int i=0;i<limit;i++) G[i]=1ll*F[i]*iF[i]%P;
            ntt(G,limit,-1);
            jifen(G,n-1,G); for (int i=n;i<limit;i++)G[i]=0;
        }
        void getexp(int *f,int n,int *G){
            if (n==1){G[0]=1;return;}
            getexp(f,(n+1)>>1,G);
            static int F[N],lnG[N];
            int limit=1; while(limit<=2*(n-1))limit<<=1; init(limit);
            for (int i=0;i<limit;i++) F[i]=i>=n?0:f[i],G[i]=i>=n?0:G[i];
            getln(G,n,lnG);
            ntt(G,limit,1),ntt(lnG,limit,1),ntt(F,limit,1);
            for (int i=0;i<limit;i++) G[i]=1ll*G[i]*sub(add(1,F[i]),lnG[i])%P;
            ntt(G,limit,-1);
            for (int i=n;i<limit;i++) G[i]=0;
        }
        void getpow(int *f,int n,int k,int *G){
            static int F[N]; for (int i=0;i<n;i++) F[i]=f[i];
            getln(F,n,F);
            for (int i=0;i<n;i++) F[i]=1ll*F[i]*k%P;
            getexp(F,n,G);
        }
        void getsqrt(int *f,int n,int *sqt){
            if (n==1){sqt[0]=1;return;}
            getsqrt(f,(n+1)>>1,sqt);
            static int F[N],H[N],iH[N];
            int limit=1; while (limit<=2*(n-1)) limit<<=1;
            for (int i=0;i<limit;i++) H[i]=i>=n?0:2ll*sqt[i]%P;
            getinv(H,n,iH),init(limit);
            for (int i=0;i<limit;i++) F[i]=i>=n?0:f[i],sqt[i]=i>=n?0:sqt[i];
            ntt(F,limit,1),ntt(sqt,limit,1),ntt(iH,limit,1);
            for (int i=0;i<limit;i++) sqt[i]=1ll*add(F[i],sqr(sqt[i]))*iH[i]%P;
            ntt(sqt,limit,-1);
            for (int i=n;i<limit;i++) sqt[i]=0;
        }
    }
    int f[N],ans[N],n,k;
    int main(){
        return 0;
    }
    
  • 相关阅读:
    在同一asp.net website工程中引用app_code中定义的control
    perl中的ppm 安装包时,如果通过代理如何下载最新的包及离线安装的问题
    利用htc s900拨号上联通3G网
    转载一篇有关ramdisk 的使用的文章
    .net dataTable对象的处理性能
    TPlink 340G+设置外网访问内网机器
    C++中的const关键字(zz)
    存在 汪峰
    ETW (Event Tracing for Windows)介绍
    理解smart pointer之一(auto_ptr介绍)
  • 原文地址:https://www.cnblogs.com/wxq1229/p/12245036.html
Copyright © 2011-2022 走看看