【学习笔记】多项式 3：多项式运算

多项式求导与积分

\[(ax^n)'=a\times nx^{n-1} \]

\[\int ax^n\mathrm{d}x=\dfrac{a}{n+1}x^{n+1} \]

单独对每一项操作后加和。

预处理逆元，时间复杂度 \(O(n)\)。

inline Poly Differntial(int n){
    Poly F;
    F.set(n);
    for(int i=1;i<n;++i) F[i-1]=f[i]*i%mod;
    return F;
} 
inline Poly Integral(int n){
    Poly F;
    F.set(n+1);
    for(int i=0;i<n;++i) F[i+1]=f[i]*inv[i+1]%mod;
    return F;
}

多项式牛顿迭代

牛顿强啊！

设函数 \(G\)，要求一多项式 \(f(x)\)，满足：

\[G(f(x))\equiv 0\pmod {x^n} \]

假定已经求得 \(f_0(x)\)，满足：

\[G(f_0(x))\equiv 0\pmod {x^{\left\lceil n/2\right\rceil}} \]

泰勒展开：

\[G(f(x))=\sum_{i=0}^{\infty} \dfrac{G^{(i)}(f_0(x))}{i!}(f(x)-f_0(x))^i \]

由于右侧 \(f(x)-f_0(x)\) 的最低非零次项次数为 \(\left\lceil n/2\right\rceil\)，因此：

\[\forall i\ge 2,(f(x)-f_0(x))^i\equiv 0\pmod {x^n} \]

于是式子实际上是：

\[G(f_0(x))+G'(f_0(x))(f(x)-f_0(x))\equiv 0\pmod {x^n} \]

整理一下就是：

\[f(x)\equiv f_0(x)-\dfrac{G(f_0(x))}{G'(f_0(x))}\pmod {x^n} \]

于是我们只需要求出 \([x^0]f(x)\) 再倍增即可。

多项式求逆

给定 \(F(x)\)，求 \(G(x)\) 使得：

\[F(x)G(x)\equiv 1\pmod {x^n} \]

常规倍增法

假设已经求得 \(G_0\) 满足：

\[F(x)G_0(x)\equiv 1\pmod {x^{\left\lceil n/2\right\rceil}} \]

接下来推一下式子：

\[\begin{aligned} G(x)&\equiv G_0(x)&\pmod {x^{\left\lceil n/2\right\rceil}}\\ G(x)-G_0(x)&\equiv 0&\pmod {x^{\left\lceil n/2\right\rceil}}\\ (G(x)-G_0(x))^2&\equiv 0&\pmod {x^n}\\ F(x)(G^2(x)-2G(x)G_0(x)+G_0^2(x))&\equiv 0&\pmod {x^n}\\ G(x)-2G_0(x)+F(x)G_0^2(x)&\equiv 0&\pmod {x^n}\\ \end{aligned}\]

于是得到：

\[G(x)\equiv G_0(x)(2-F(x)G_0(x))\pmod {x^n} \]

牛顿迭代法

构造：

\[G(f(x))=\dfrac{1}{f(x)}-F(x) \]

代入牛顿迭代式：

\[f(x)=f_0(x)-\dfrac{\dfrac{1}{f_0(x)}-F(x)}{-\dfrac{1}{f_0^2(x)}} \]

整理一下就是：

\[f(x)\equiv f_0(x)(2-F(x)f_0(x))\pmod {x^n} \]

inline Poly Inv(int n){
    Poly F,G;
    G.set(1);
    G[0]=q_pow(f[0],mod-2,mod);
    for(int L=2;L<2*n;L<<=1){
        F.set(2*L),G.set(2*L);
        F.clear(0,2*L-1);
        for(int i=0;i<L;++i) F[i]=f[i];
        F.NTT(2*L,1),G.NTT(2*L,1);
        for(int i=0;i<2*L;++i) G[i]=1ll*(int)G[i]*(2ll-1ll*(int)G[i]*(int)F[i]%mod+mod)%mod;
        G.NTT(2*L,0);
        G.clear(min(n,L),2*L-1);
    }
    return G;
}

多项式开根

给定 \(F(x)\)，求 \(G(x)\) 使得：

\[F(x)\equiv G^2(x)\pmod {x^n} \]

常规倍增法

假设已经求得 \(G_0\) 满足：

\[F(x)\equiv G_0^2(x)\pmod {x^{\left\lceil n/2\right\rceil}} \]

接下来推一下式子：

\[\begin{aligned} G(x)&\equiv G_0(x)&\pmod {x^{\left\lceil n/2\right\rceil}}\\ G(x)-G_0(x)&\equiv 0&\pmod {x^{\left\lceil n/2\right\rceil}}\\ (G(x)-G_0(x))^2&\equiv 0&\pmod {x^n}\\ G^2(x)-2G(x)G_0(x)+G_0^2(x)&\equiv 0&\pmod {x^n}\\ F(x)-2G(x)G_0(x)+G_0^2(x)&\equiv 0&\pmod {x^n}\\ \end{aligned}\]

于是得到：

\[G(x)\equiv \dfrac{F(x)+G_0^2(x)}{2G_0(x)}\pmod {x^n} \]

即：

\[G(x)\equiv \dfrac{F(x)}{2G_0(x)}+\dfrac{G_0(x)}{2}\pmod {x^n} \]

牛顿迭代法

构造：

\[G(f(x))=f^2(x)-F(x) \]

代入牛顿迭代式：

\[f(x)=f_0(x)-\dfrac{f_0^2(x)-F(x)}{2f_0(x)} \]

整理一下就是：

\[f(x)\equiv \dfrac{F(x)+f_0^2(x)}{2f_0(x)}\pmod {x^n} \]

即：

\[f(x)\equiv \dfrac{F(x)}{2f_0(x)}+\dfrac{f_0(x)}{2}\pmod {x^n} \]

inline Poly Sqrt(int n){
    Poly F,G,invG,H;
    G.set(1);
    G[0]=1;
    for(int L=2;L<2*n;L<<=1){
        F.set(2*L),G.set(2*L),H.set(2*L);
        F.clear(0,2*L-1);
        for(int i=0;i<L;++i) F[i]=f[i];
        invG=G.Inv(L);
        F.NTT(2*L,1),invG.NTT(2*L,1);
        for(int i=0;i<2*L;++i) H[i]=F[i]*invG[i]%mod*inv[2]%mod;
        H.NTT(2*L,0);
        for(int i=0;i<2*L;++i) G[i]=(H[i]+G[i]*inv[2]%mod)%mod;
        G.clear(min(n,L),2*L-1);
    }
    return G;
}

多项式对数函数

给定 \(F(x)\)，求 \(G(x)\) 使得：

\[\ln F(x)\equiv G(x)\pmod {x^n} \]

左侧相当于一个复合函数，对左右求导：

\[(\ln F(x))'\times F'(x)\equiv G'(x)\pmod {x^n} \]

根据简单初等函数的导数得：

\[\dfrac{F'(x)}{F(x)}\equiv G'(x)\pmod {x^n} \]

求出 \(G'(x)\) 后积分得到 \(G(x)\)：

\[G(x)=\int G'(x) \mathrm{d}x \]

inline Poly Ln(int n){
    Poly dF,invF,G,H;
    dF=Differntial(n),invF=Inv(n);
    int L=1;
    while(L<2*n) L<<=1;
    dF.set(L),invF.set(L),G.set(L);
    dF.NTT(L,1),invF.NTT(L,1);
    for(int i=0;i<L;++i) G[i]=dF[i]*invF[i]%mod;
    G.NTT(L,0);
    H=G.Integral(L);
    H.clear(n,L);
    return H;
}

注意多项式牛顿迭代和多项式 \(\ln\) 在求导时的区别，前者的 \(G(f(x))\) 是关于 \(f(x)\) 的函数，也即 \(G'(f(x))=\dfrac{\mathrm{d} G(f(x))}{\mathrm{d}f(x)}\)；后者 \(G(x)\) 是关于 \(x\) 的函数，也即 \(G'(x)=\dfrac{\mathrm{d}G(x)}{\mathrm{d}x}\)，而左侧的 \(\ln F(x)\) 是关于 \(x\) 的复合函数，所以要应用复合函数求导 \((\ln F(x))'=\dfrac{\mathrm{d} \ln F(x)}{\mathrm{d} F(x)}\times\dfrac{\mathrm{d} F(x)}{\mathrm{d}x}\)。

多项式指数函数

给定 \(F(x)\)，求 \(G(x)\) 使得：

\[\exp F(x)\equiv G(x)\pmod {x^n} \]

使用牛顿迭代法，构造：

\[G(f(x))=\ln f(x)-F(x) \]

代入牛顿迭代式：

\[f(x)=f_0(x)-\dfrac{\ln f_0(x)-F(x)}{\dfrac{1}{f_0(x)}} \]

整理一下就是：

\[f(x)\equiv f_0(x)(1-\ln f_0(x)+F(x))\pmod {x^n} \]

inline Poly Exp(int n){
    Poly F,G,lnG;
    G.set(1);
    G[0]=1;
    for(int L=2;L<2*n;L<<=1){
        F.set(2*L),G.set(2*L);
        F.clear(0,2*L-1);
        for(int i=0;i<L;++i) F[i]=f[i];
        lnG=G.Ln(L);
        F.NTT(2*L,1),G.NTT(2*L,1),lnG.NTT(2*L,1);
        for(int i=0;i<2*L;++i) G[i]=G[i]*(1ull-lnG[i]+F[i]+mod)%mod;
        G.NTT(2*L,0);
        G.clear(min(n,L),2*L-1);
    }
    return G;
}

多项式半家桶封装模板

点击查看代码

inline int q_pow(int A,int B,int P){
    int res=1;
    while(B){
        if(B&1) res=1ll*res*A%P;
        A=1ll*A*A%P;
        B>>=1;
    }
    return res;
}

int inv[maxn];
int rev[maxn];
int base,w[maxn];
struct Poly{
    const static int g=3;
    int deg;
    vector<ull> f;
    ull &operator[](const int &i){return f[i];}
    ull operator[](const int &i)const{return f[i];}
    inline void set(int L){deg=L;f.resize(L);}
    inline void clear(int L,int R){for(int i=L;i<=R;++i)f[i]=0;}
    inline void output(int L){for(int i=0;i<L;++i)printf("%llu ",f[i]);printf("\n");}
    inline void NTT(int L,bool type){
        set(L);
        for(int i=1;i<L;++i){
            rev[i]=(rev[i>>1]>>1)+(i&1?L>>1:0);
            if(i<rev[i]) swap(f[i],f[rev[i]]);
        }
        for(int d=1;d<L;d<<=1){
            base=q_pow(type?g:q_pow(g,mod-2,mod),(mod-1)/(2*d),mod);
            w[0]=1;
            for(int i=1;i<d;++i) w[i]=1ll*w[i-1]*base%mod;
            for(int i=0;i<L;i+=d<<1){
                for(int j=0;j<d;++j){
                    ull x=f[i+j],y=f[i+d+j]*w[j]%mod;
                    f[i+j]=x+y,f[i+d+j]=x-y+mod;
                }
            }
        }
        for(int i=0;i<L;++i) f[i]%=mod;
        if(!type){
            int inv_L=q_pow(L,mod-2,mod);
            for(int i=0;i<L;++i) f[i]=f[i]*inv_L%mod;
        }
    }
    inline Poly Differntial(int n){
        Poly F;
        F.set(n);
        for(int i=1;i<n;++i) F[i-1]=f[i]*i%mod;
        return F;
    } 
    inline Poly Integral(int n){
        Poly F;
        F.set(n+1);
        for(int i=0;i<n;++i) F[i+1]=f[i]*inv[i+1]%mod;
        return F;
    }
    inline Poly Inv(int n){
        Poly F,G;
        G.set(1);
        G[0]=q_pow(f[0],mod-2,mod);
        for(int L=2;L<2*n;L<<=1){
            F.set(2*L),G.set(2*L);
            F.clear(0,2*L-1);
            for(int i=0;i<L;++i) F[i]=f[i];
            F.NTT(2*L,1),G.NTT(2*L,1);
            for(int i=0;i<2*L;++i) G[i]=1ll*(int)G[i]*(2ll-1ll*(int)G[i]*(int)F[i]%mod+mod)%mod;
            G.NTT(2*L,0);
            G.clear(min(n,L),2*L-1);
        }
        return G;
    }
    inline Poly Sqrt(int n){
        Poly F,G,invG,H;
        G.set(1);
        G[0]=1;
        for(int L=2;L<2*n;L<<=1){
            F.set(2*L),G.set(2*L),H.set(2*L);
            F.clear(0,2*L-1);
            for(int i=0;i<L;++i) F[i]=f[i];
            invG=G.Inv(L);
            F.NTT(2*L,1),invG.NTT(2*L,1);
            for(int i=0;i<2*L;++i) H[i]=F[i]*invG[i]%mod*inv[2]%mod;
            H.NTT(2*L,0);
            for(int i=0;i<2*L;++i) G[i]=(H[i]+G[i]*inv[2]%mod)%mod;
            G.clear(min(n,L),2*L-1);
        }
        return G;
    }
    inline Poly Ln(int n){
        Poly dF,invF,G,H;
        dF=Differntial(n),invF=Inv(n);
        int L=1;
        while(L<2*n) L<<=1;
        dF.set(L),invF.set(L),G.set(L);
        dF.NTT(L,1),invF.NTT(L,1);
        for(int i=0;i<L;++i) G[i]=dF[i]*invF[i]%mod;
        G.NTT(L,0);
        H=G.Integral(L);
        H.clear(n,L);
        return H;
    }
    inline Poly Exp(int n){
        Poly F,G,lnG;
        G.set(1);
        G[0]=1;
        for(int L=2;L<2*n;L<<=1){
            F.set(2*L),G.set(2*L);
            F.clear(0,2*L-1);
            for(int i=0;i<L;++i) F[i]=f[i];
            lnG=G.Ln(L);
            F.NTT(2*L,1),G.NTT(2*L,1),lnG.NTT(2*L,1);
            for(int i=0;i<2*L;++i) G[i]=G[i]*(1ull-lnG[i]+F[i]+mod)%mod;
            G.NTT(2*L,0);
            G.clear(min(n,L),2*L-1);
        }
        return G;
    }
};

参考资料

• OI Wiki