多项式复合
这东西 (O(n^2+nsqrt nlog n)) 过了 (20000)
思路是非常精巧的,令 (L=sqrt n) 考虑每个次幂的对系数的贡献那么得到一个式子
分块预处理 (G) 的大块和小块就行了
那么这个算是数据结构维护分段多项式的入门题
多项式复合逆
考虑复合的过程和复合逆的过程是类似的
这个做法的正确性如下
Lagrange反演&扩展Lagrange反演
普通反演看懂了,扩展反演没人写证明
式子如下:
证明并不是 (trivial) 的,考虑将 (G(F(x))=x) 的式子展开并求导得到
考虑反演本身的式子可以消掉 (x^n) 那么保留 (x^{-1}) 和常数项来观察
由于多项式函数求导以后 (x^{-1}) 的系数总是 (0)(导数定义),那么考虑对展开的式子除掉 (F^n(x)) 并得到 (x^{-1}) 的系数有:
对于左边而言,接着展开两个函数的乘积就可以得到 (x^{-1}) 的系数必然为 (1)
带入得到原反演式子,证毕
关于扩展的证明貌似是非常复杂的,所以先咕掉了,式子如下
多项式除法
设 (F_R(x)=x^nF(frac{1}x)),其实也就是系数翻转(没错,(R) 就是 (Reverse))
配上 (x^n) 可以得到:
换式子:
那么求逆得到 (Q_R(x)) 之后减法做出来 (R(x)) 就行了
写得时候记得清空一些数组
取模和除法是并行的
多项式多点求值
设 (F(x)=G(x)Q(x)+R(x)),那么对于任意一个 (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) 都跑不了,可行的优化手段,设
所以 (F(x_i)=[z^0]Mul(F(z),frac{1}{1-x_iz})),还是分治,那么每次传入的变成了:
那么每次递归下去的是时候利用 乘另外半边的乘积并且保留区间长度的项数就可以了
用多项式乘法代替了取模,快了很多
这个方法据说可以 (1s 5e5) ,代码留坑了
多项式快速插值
根据朴素的 (Lagrange) 插值法可以得到一个公式:
这个公式是 (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})
根据洛必达法则:
多点求值即可(代码中的多点求值参考了 (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();}
本文式子可以很快看懂,但是写出来不易,博主中途自闭了很久,写了一整天
所以多项式的使用到这里也就不会再有太多扩展了
剩下的就是熟练各种技巧了,慢慢来