「学习笔记」多项式 II

多项式复合

这东西 (O(n^2+nsqrt nlog n)) 过了 (20000)

思路是非常精巧的,令 (L=sqrt n) 考虑每个次幂的对系数的贡献那么得到一个式子

[sum_{i=0}^L G^{iL}sum_{j=0}^{L-1}G^jF[x^{iL+j}] ]

分块预处理 (G) 的大块和小块就行了

那么这个算是数据结构维护分段多项式的入门题

多项式复合逆

考虑复合的过程和复合逆的过程是类似的

[G(x)=sum_{i=0}^Lfrac{1}{iL}(frac{x}{F(x)})^{iL} sum_{j=0}^{L-1} frac1j (frac{x}{F(x)})^{j} ]

这个做法的正确性如下

Lagrange反演&扩展Lagrange反演

普通反演看懂了,扩展反演没人写证明

式子如下:

[G(F(x))=x o [x^n]G(x)=frac1{n} [x^{n-1}] (frac{x}{F(x)})^n ]

证明并不是 (trivial) 的,考虑将 (G(F(x))=x) 的式子展开并求导得到

[sum F^{i-1}(x) ia_iF'(x)=1 ]

考虑反演本身的式子可以消掉 (x^n) 那么保留 (x^{-1}) 和常数项来观察

由于多项式函数求导以后 (x^{-1}) 的系数总是 (0)(导数定义),那么考虑对展开的式子除掉 (F^n(x)) 并得到 (x^{-1}) 的系数有:

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

对于左边而言,接着展开两个函数的乘积就可以得到 (x^{-1}) 的系数必然为 (1)

带入得到原反演式子,证毕

关于扩展的证明貌似是非常复杂的,所以先咕掉了,式子如下

[H(F(x))[x^n]=frac1n [x^{-1}] H'(x)frac1 {G^n(x)} ]

多项式除法

(F_R(x)=x^nF(frac{1}x)),其实也就是系数翻转(没错,(R) 就是 (Reverse)

[F(x)=G(x)Q(x)+R(x) ]

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

配上 (x^n) 可以得到:

[x^nF(x)=x^mG(frac1 x)x^{n-m} Q(frac{1}x)+x^{n-m-1} x^{m-1} R(frac 1 x) ]

换式子:

[F_R(x)=G_R(x)Q_R(x)+x^{n-m-1} R_R(x) ]

[F_R(x)equiv G_R(x)Q_R(x)mod x^{n-m-1} ]

那么求逆得到 (Q_R(x)) 之后减法做出来 (R(x)) 就行了

写得时候记得清空一些数组

取模和除法是并行的

多项式多点求值

(F(x)=G(x)Q(x)+R(x)),那么对于任意一个 (x_0) 都会有

[F(x_0)=G(x_0)Q(x_0)+R(x_0) ]

(G(x_0)=0)(F(x_0)=R(x_0))

那么设 (P_0(x)=prod_{i=0}^{n/2} x-x_i,P_1(x)=prod_{i=n/2+1}^n x-x_i)

(F(x)) 不停地对去区间的乘积取模就可以最后剩下常数作为答案

分治乘法作预处理,再 (dfs) 一次得到所有答案

这里取模的时候千万注意清空

代码如下:

#include<bits/stdc++.h>
using namespace std;
#define reg register
#define int long long
namespace yspm{
    inline int read(){
        int res=0,f=1; char k;
        while(!isdigit(k=getchar())) if(k=='-') f=-1;
        while(isdigit(k)) res=res*10+k-'0',k=getchar();
        return res*f;
    }
    const int N=3e5+10,mod=998244353;
    inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
    inline int del(int x,int y){return x-y<0?x-y+mod:x-y;}
    inline int mul(int x,int y){return x*y-x*y/mod*mod;}
    inline int ksm(int x,int y){int res=1; for(;y;y>>=1,x=mul(x,x)) if(y&1) res=mul(res,x); return res;}
    int r[N],W[N];
    inline void NTT(int *f,int lim,int opt){
        for(reg int i=0;i<lim;++i) r[i]=r[i>>1]>>1|((i&1)?(lim>>1):0);
        for(reg int i=0;i<lim;++i) if(i<r[i]) swap(f[i],f[r[i]]); 
        for(reg int p=2;p<=lim;p<<=1){
            int len=p>>1,tmp=ksm(3,(mod-1)/p); if(opt==-1) tmp=ksm(tmp,mod-2); 
            W[0]=1; for(reg int i=1;i<len;++i) W[i]=mul(W[i-1],tmp);
            for(reg int k=0;k<lim;k+=p) 
                  for(reg int tt,l=k;l<k+len;++l) tt=mul(W[l-k],f[l+len]),f[l+len]=del(f[l],tt),f[l]=add(f[l],tt);
        }if(opt==-1){int tp=ksm(lim,mod-2); for(reg int i=0;i<lim;++i) f[i]=mul(f[i],tp);} return ;
    }
    int H[N],m,n,Ar[N],Br[N],I[N],C[N],x[N],y[N],z[N];
    inline void Juan(int *F,int *G,int *res,int lim){
        for(reg int i=0;i<lim;++i) x[i]=y[i]=z[i]=0;
        for(reg int i=0;i<(lim>>1);++i) x[i]=F[i],y[i]=G[i];
        NTT(x,lim,1); NTT(y,lim,1); for(reg int i=0;i<lim;++i) z[i]=mul(x[i],y[i]); NTT(z,lim,-1);
        for(reg int i=0;i<lim;++i) res[i]=z[i]; return ;
    }
    inline void Inv(int *F,int *G,int len){
        if(len==1) return G[0]=ksm(F[0],mod-2),void(); Inv(F,G,(len+1)>>1);
        int lim=1; while(lim<=(len<<1)) lim<<=1; 
        for(reg int i=0;i<len;++i) H[i]=F[i]; for(reg int i=len;i<lim;++i) H[i]=0; 
        NTT(H,lim,1); NTT(G,lim,1);
        for(reg int i=0;i<lim;++i) G[i]=mul(G[i],del(2,mul(G[i],H[i]))); NTT(G,lim,-1); 
        for(reg int i=len;i<lim;++i) G[i]=0;  return ;
    }
    inline void Mod(int *F,vector<int> G,int *res,int n,int m){
        int lim=1; while(lim<=n) lim<<=1;
        for(reg int i=0;i<lim;++i) I[i]=Ar[i]=Br[i]=C[i]=0;  
        for(reg int i=0;i<=n;++i) Ar[i]=F[n-i]; for(reg int i=0;i<=m;++i) Br[i]=G[m-i];
        for(reg int i=n-m+1;i<=m;++i) Br[i]=0; Inv(Br,I,n-m+1); 
        lim=1; while(lim<=(n-m+1)*2) lim<<=1; Juan(Ar,I,C,lim); 
        reverse(C,C+n-m+1); lim=1; while(lim<=((n+1)<<1)) lim<<=1; 
        for(reg int i=n-m+1;i<lim;++i) C[i]=0;
        for(reg int i=0;i<=m;++i) Br[i]=G[i]; 
        //这里br要清空
        for(reg int i=m+1;i<lim;++i) Br[i]=0;
        
        NTT(Br,lim,1); NTT(C,lim,1); for(reg int i=0;i<lim;++i) C[i]=mul(Br[i],C[i]);
        NTT(C,lim,-1);
        for(reg int i=0;i<m;++i) res[i]=del(F[i],C[i]); 
        for(reg int i=m;i<=n;++i) res[i]=0; return ;
    }
    int a[N],f[N]; vector<int> P[N];
    inline void prework(int p,int l,int r){
        if(l==r){P[p].resize(2); P[p][0]=mod-a[l],P[p][1]=1; return ;}
        
        int mid=(l+r)>>1; prework(p<<1,l,mid); prework(p<<1|1,mid+1,r);
        int len=P[p<<1].size()+P[p<<1|1].size(); P[p].resize(len-1); 
        int lim=1; while(lim<len) lim<<=1;    
        for(reg int i=0;i<lim;++i){
            x[i]=(i<P[p<<1].size())?P[p<<1][i]:0;
            y[i]=(i<P[p<<1|1].size())?P[p<<1|1][i]:0;
        }
        NTT(x,lim,1); NTT(y,lim,1); 
        for(reg int i=0;i<lim;++i) x[i]=mul(x[i],y[i]); NTT(x,lim,-1); 
        for(reg int i=0;i<P[p].size();++i) P[p][i]=x[i];
        return ;
    }
    int g[30][N],bin[31],top;
    inline void solve(int p,int l,int r,int *f){
        if(l==r) return printf("%lld
",f[0]),void(); int mid=(l+r)>>1;
        int ls=bin[top--]; 
        Mod(f,P[p<<1],g[ls],P[p].size()-2,P[p<<1].size()-1);
        //谨慎思考多项式的系数是size-2
        solve(p<<1,l,mid,g[ls]); bin[++top]=ls;
        int rs=bin[top--]; Mod(f,P[p<<1|1],g[rs],P[p].size()-2,P[p<<1|1].size()-1); 
        solve(p<<1|1,mid+1,r,g[rs]); bin[++top]=rs;
        return ;
    }
    signed main(){
        for(reg int i=0;i<30;++i) bin[++top]=i; 
        n=read(); m=read(); for(reg int i=0;i<=n;++i) f[i]=read(); 
        for(reg int i=1;i<=m;++i) a[i]=read(); if(!m) return 0; 
        prework(1,1,m); if(n>=m) Mod(f,P[1],f,n,m); solve(1,1,m,f);
        return 0;
    }
}
signed main(){return yspm::main();}

这样 (1s 6e4) 都跑不了,可行的优化手段,设

[Mul(F(x),G(x))=sum_i (sum_j f_{i+j}g_j)x^i ]

所以 (F(x_i)=[z^0]Mul(F(z),frac{1}{1-x_iz})),还是分治,那么每次传入的变成了:

[Mul(F(z),frac1{prod (1-x_iz)}) ]

那么每次递归下去的是时候利用 乘另外半边的乘积并且保留区间长度的项数就可以了

用多项式乘法代替了取模,快了很多

这个方法据说可以 (1s 5e5) ,代码留坑了

多项式快速插值

根据朴素的 (Lagrange) 插值法可以得到一个公式:

[F(x)=sum_{i=1} ^n y_i prod_{j eq i}frac{x-x_j}{x_i-x_j} ]

这个公式是 (trivial) 的,因为每次带入 (x_i) 的话其它的项都会消掉

但是显然它是 (O(n^2)) 的,那么考虑分治(为啥全是分治)

设多项式 (p(l,r)=prod_{i=l}^r x-x_i,v(l,r)=sum_{i=l}^r y_iprod_{j eq i} frac{x-x_j}{x_i-x_j})

求是 (trivial) 的,因为每个不等于只用判一下,(v(l,r)=v(l,mid)p(mid+1,r)+v(mid+1,r)p(l,mid))

但是瓶颈在于求这个 (frac{1}{prod_{j eq i} x-x_j})

根据洛必达法则:

[G(x)=frac{1}{prod x-x_i},G(x_i)=G'(x) ]

多点求值即可(代码中的多点求值参考了 (1s 5e5) 版的)

#include<bits/stdc++.h>
using namespace std;
#define reg register
namespace yspm{
    inline int read(){
        int res=0,f=1; char k;
        while(!isdigit(k=getchar())) res=res*10+k-'0',k=getchar();
        while(isdigit(k)) res=res*10+k-'0',k=getchar();
        return res*f;
    }
    const int N=3e5+10,mod=998244353;
    inline int add(reg int x,reg int y){return x+y>=mod?x+y-mod:x+y;}
    inline int del(reg int x,reg int y){return x-y<0?x-y+mod:x-y;}
    inline int mul(reg int x,reg int y){return 1ll*x*y-1ll*x*y/mod*mod;}
    inline int ksm(reg int x,reg int y){reg int res=1; for(;y;y>>=1,x=mul(x,x)) if(y&1) res=mul(res,x); return res;}
    int r[N],W[N];
    inline void NTT(vector<int> &f,reg int lim,reg int opt){ f.resize(lim);
        for(reg int i=0;i<lim;++i) r[i]=r[i>>1]>>1|((i&1)?(lim>>1):0);
        for(reg int i=0;i<lim;++i) if(i<r[i]) swap(f[i],f[r[i]]);
        for(reg int p=2;p<=lim;p<<=1){
            int len=p>>1; W[0]=1; W[1]=ksm(3,(mod-1)/p); 
            if(opt==-1) W[1]=ksm(W[1],mod-2); for(reg int i=2;i<len;++i) W[i]=mul(W[i-1],W[1]);
            for(reg int k=0;k<lim;k+=p) 
                  for(reg int l=k,tt;l<k+len;++l) tt=mul(W[l-k],f[l+len]),f[l+len]=del(f[l],tt),f[l]=add(f[l],tt);
        }if(opt==-1){int tp=ksm(lim,mod-2); for(reg int i=0;i<lim;++i) f[i]=mul(f[i],tp);} return ;
    }
    vector<int> H;
    inline void Inv(vector<int> A,vector<int> &B,reg int len){
        if(len==1){B.resize(1); B[0]=ksm(A[0],mod-2); return ;} Inv(A,B,(len+1)>>1); 
        int lim=1; while(lim<(len<<1)) lim<<=1; 
        H.clear(); H.resize(lim); for(reg int i=0;i<len;++i) H[i]=A[i]; for(reg int i=len;i<lim;++i) H[i]=0; 
        NTT(H,lim,1); B.resize(lim); NTT(B,lim,1); 
        for(reg int i=0;i<lim;++i) B[i]=mul(B[i],del(2,mul(B[i],H[i]))); NTT(B,lim,-1); 
        for(reg int i=len;i<lim;++i) B[i]=0; B.resize(len);
        return ;
    }
    struct _1s5e5{
        inline vector<int> Mult(vector<int> f,vector<int> g,reg int ned){
            reg int sz=f.size();
            for(reg int i=0;i<sz;++i) f[i]=mul(f[i],g[i]); NTT(f,sz,1); 
            for(reg int i=ned;i<sz;++i) f[i]=0; f.resize(ned); return f;
        }
        vector<int> P[N<<1],arr,ans;
        inline void prework(reg int p,reg int l,reg int r){
            if(l==r) return P[p].resize(2),P[p][0]=1,P[p][1]=mod-arr[l],void();
            reg int mid=(l+r)>>1,len=r-l+2,lim=1; while(lim<len) lim<<=1;   
            prework(p<<1,l,mid); prework(p<<1|1,mid+1,r); P[p].resize(lim);
            NTT(P[p<<1],lim,1); NTT(P[p<<1|1],lim,1); 
            for(reg int i=0;i<lim;++i) P[p][i]=mul(P[p<<1][i],P[p<<1|1][i]); NTT(P[p],lim,-1); 
            for(reg int i=len;i<lim;++i) P[p][i]=0; P[p].resize(len);
            return ;
        }
        inline void dfs(int p,int l,int r,vector<int> F){
            if(l==r) return ans[l]=F[0],void();
            reg int mid=(l+r)>>1,lim=1; while(lim<r-l+2) lim<<=1; NTT(F,lim,-1); 
            dfs(p<<1,l,mid,Mult(F,P[p<<1|1],mid-l+1)); 
            dfs(p<<1|1,mid+1,r,Mult(F,P[p<<1],r-mid)); 
            return ;
        }
        vector<int> tmp; 
        inline vector<int> query(vector<int> poly,vector<int> ask){
            reg int n=poly.size(),m=ask.size(); arr=ask; prework(1,1,m-1);
            ans.resize(ask.size()); Inv(P[1],tmp,n); 
            reg int lim=1; while(lim<(n<<1)) lim<<=1; 
            NTT(tmp,lim,1); NTT(poly,lim,-1);
            dfs(1,1,m-1,Mult(poly,tmp,n)); 
            return ans;
        } 
    }P;
    struct KSC{
        vector<int> X,Y,T[N<<1],ans[N<<1],res;
        inline void prework(int p,int l,int r){
            if(l==r){
                T[p].resize(2); T[p][0]=mod-X[l]; T[p][1]=1; 
                return ; 
            }int mid=(l+r)>>1,len=r-l+2,lim=1; prework(p<<1,l,mid); prework(p<<1|1,mid+1,r); 
            while(lim<len) lim<<=1; T[p].resize(lim); 
            NTT(T[p<<1],lim,1); NTT(T[p<<1|1],lim,1); 
            for(reg int i=0;i<lim;++i) T[p][i]=mul(T[p<<1][i],T[p<<1|1][i]); NTT(T[p],lim,-1);
            for(reg int i=len;i<lim;++i) T[p][i]=0; T[p].resize(len);  
            return ;
        }
        inline void dfs(int p,int l,int r){
            if(l==r){
                ans[p].resize(1); 
                ans[p][0]=mul(Y[l],ksm(res[l],mod-2));
                return ;
            }int mid=(l+r)>>1,len=r-l+2,lim=1; while(lim<len) lim<<=1; 
            dfs(p<<1,l,mid); dfs(p<<1|1,mid+1,r); ans[p].resize(lim); 
            NTT(ans[p<<1],lim,1); NTT(ans[p<<1|1],lim,1); 
            for(reg int i=0;i<lim;++i) ans[p][i]=add(mul(ans[p<<1][i],T[p<<1|1][i]),mul(ans[p<<1|1][i],T[p<<1][i]));
            NTT(ans[p],lim,-1);
            for(reg int i=len;i<lim;++i) ans[p][i]=0; ans[p].resize(len);
            return ;
        }
        inline vector<int> query(vector<int> x,vector<int> y){
            X=x; Y=y; int n=x.size()-1; prework(1,1,n);
            for(reg int i=1;i<=n;++i) T[1][i-1]=mul(T[1][i],i); T[1][n]=0; 
            res=P.query(T[1],x);
            dfs(1,1,n); 
            return ans[1];
        }
    }PO; int n;
    signed main(){
        n=read(); vector<int> a,b; a.resize(n+1); b.resize(n+1);
        for(reg int i=1;i<=n;++i) a[i]=read(),b[i]=read(); 
        vector<int> ans=PO.query(a,b); 
        for(reg int i=0;i<n;++i) printf("%d ",ans[i]); puts("");
        return 0;
    }
}
signed main(){return yspm::main();}

本文式子可以很快看懂,但是写出来不易,博主中途自闭了很久,写了一整天

所以多项式的使用到这里也就不会再有太多扩展了

剩下的就是熟练各种技巧了,慢慢来

原文地址:https://www.cnblogs.com/yspm/p/14369875.html