当前位置: 代码迷 >> 综合 >> hdu6900 Residual Polynomial(分治NTT/启发式NTT)
  详细解决方案

hdu6900 Residual Polynomial(分治NTT/启发式NTT)

热度:5   发布时间:2024-02-24 01:25:48.0

题目

n(3<=n<=1e5)个多项式,

,现给出的第0次项到第n次项系数,第i项是ai满足0<=ai<998244353

对于,,其中bi,ci满足0<=bi,ci<998244353

现在要求输出的第0项到第n项系数,每个系数均对998244353取模

T(T<=100)组样例,但只有最多3组样例,满足n>1000

思路来源

jyt、tzy

题解

一道看似多项式求导,实则分治NTT/启发式NTT的题目

小范围暴力,大范围NTT,事实上,都NTT也都没事

 

考虑f1的第i项最后通过求导落到了fn的第j项里,

首先满足i>=j,其次经过了i-j次求导,

递推式可以看成求导或者不求导二选一,求导的话*bi,不求导的话*ci,

这个方案贡献就是(b1x+c1)*(b2x+c2)*...*(bnx+cn),中间x^i-j的系数就是求了i-j次导,

这个可以分治/启发式NTT求一求

然后x^i本身的系数是ai,x^i求导求成x^j,又乘了i*(i-1)*...*(j+1)即i!/j!,

所以,,

 

xs[]先分治NTT求出来,然后这是一个减法卷积,

常用套路就是把xs[]数组反过来,i下标变n-1-i,再做卷积即可

最后,res[j]即ans[n-1+j]除以j!即可

心得

orz我好菜啊惨遭学弟学妹吊打

push自己换了个板子,3s的题自己的板子跑了2800ms,学弟的板子900ms

代码

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned ll
const int N = 1<<18, P = 998244353, Primitive_root = 3;
struct Z{int x;Z(const int _x=0):x(_x){}Z operator +(const Z &r)const{ return x+r.x<P?x+r.x:x+r.x-P;}Z operator -(const Z &r)const{ return x<r.x?x-r.x+P:x-r.x;}Z operator -()const{ return x?P-x:0;}Z operator *(const Z &r)const{ return static_cast<ull>(x)*r.x%P;}Z operator +=(const Z &r){ return x=x+r.x<P?x+r.x:x+r.x-P, *this;}Z operator -=(const Z &r){ return x=x<r.x?x-r.x+P:x-r.x, *this;}Z operator *=(const Z &r){ return x=static_cast<ull>(x)*r.x%P, *this;}friend Z Pow(Z, int);pair<Z,Z> Mul(pair<Z,Z> x, pair<Z,Z> y, Z f)const{return make_pair(x.first*y.first+x.second*y.second*f,x.second*y.first+x.first*y.second);}Z Quadratic_residue()const{if(x<=1) return x;if(Pow((Z)x, (P-1)/2).x!=1) return -1;Z y, f;mt19937 rng(20030226);do y=rng()%(x-1)+1; while(Pow(f=y*y-x, (P-1)/2).x==1);pair<Z,Z> ans=make_pair(1, 0), t=make_pair(y, 1);for(int i=(P+1)/2; i; i>>=1, t=Mul(t, t, f)) if(i&1) ans=Mul(ans, t, f);return min(ans.first.x, P-ans.first.x);}
};
Z Pow(Z x, int y=P-2){Z ans=1;for(; y; y>>=1, x=x*x) if(y&1) ans=ans*x;return ans;
}
namespace Poly{Z w[N];Z Inv[N];vector<Z> ans;vector<vector<Z>> p;ull F[N];void Init(){for(int i=1; i<N; i<<=1){w[i]=1;Z t=Pow((Z)Primitive_root, (P-1)/i/2);for(int j=1; j<i; ++j) w[i+j]=w[i+j-1]*t;}Inv[1]=1;for(int i=2; i<N; ++i) Inv[i]=Inv[P%i]*(P-P/i);}int Get(int x){ int n=1; while(n<=x) n<<=1; return n;}int Mod(int x){ return x<P?x:x-P;}void DFT(vector<Z> &f, int n){if((int)f.size()!=n) f.resize(n);for(int i=0, j=0; i<n; ++i){F[i]=f[j].x;for(int k=n>>1; (j^=k)<k; k>>=1);}if(n<=4){for(int i=1; i<n; i<<=1) for(int j=0; j<n; j+=i<<1){Z *W=w+i;ull *F0=F+j, *F1=F+j+i;for(int k=j; k<j+i; ++k, ++W, ++F0, ++F1){ull t=(*F1)*(W->x)%P;(*F1)=*F0+P-t, (*F0)+=t;}}}else{for(int j=0; j<n; j+=2){int t=F[j+1];F[j+1]=Mod(F[j]+P-t), F[j]=Mod(F[j]+t);}for(int j=0; j<n; j+=4){int t0=F[j+2], t1=F[j+3]*w[3].x%P;F[j+2]=F[j]+P-t0, F[j]+=t0;F[j+3]=F[j+1]+P-t1, F[j+1]+=t1;}for(int i=4; i<n; i<<=1) for(int j=0; j<n; j+=i<<1){Z *W=w+i;ull *F0=F+j, *F1=F+j+i;for(int k=j; k<j+i; k+=4, W+=4, F0+=4, F1+=4){int t0=(W->x)**F1%P;int t1=(W+1)->x**(F1+1)%P;int t2=(W+2)->x**(F1+2)%P;int t3=(W+3)->x**(F1+3)%P;*F1=*F0+P-t0, *F0+=t0;*(F1+1)=*(F0+1)+P-t1, *(F0+1)+=t1;*(F1+2)=*(F0+2)+P-t2, *(F0+2)+=t2;*(F1+3)=*(F0+3)+P-t3, *(F0+3)+=t3;}}}for(int i=0; i<n; ++i) f[i]=F[i]%P;}void IDFT(vector<Z> &f, int n){f.resize(n), reverse(f.begin()+1, f.end()), DFT(f, n);Z I=1;for(int i=1; i<n; i<<=1) I*=(P+1)/2;for(int i=0; i<n; ++i) f[i]*=I;}vector<Z> operator +(const vector<Z> &f, const vector<Z> &g){vector<Z> ans=f;ans.resize(max(f.size(), g.size()));for(int i=0; i<(int)g.size(); ++i) ans[i]+=g[i];return ans;}vector<Z> operator *(const vector<Z> &f, const vector<Z> &g){static vector<Z> F, G;F=f, G=g;int p=Get(f.size()+g.size()-2);DFT(F, p), DFT(G, p);for(int i=0; i<p; ++i) F[i]*=G[i];IDFT(F, p);return F.resize(f.size()+g.size()-1), F;}vector<Z> operator *(const vector<Z> &f, Z g){vector<Z> ans=f;for(Z &i:ans) i*=g;return ans;}
}
using namespace Poly;
int T,n,b[N],c[N];
Z fac[N],ifac[N];
vector<Z>work(int l, int r){if(l==r){return{c[l],b[l]};} int mid=(l+r)>>1;return work(l,mid)*work(mid+1,r);
}
void init(int n){fac[0]=1;for(int i=1;i<=n;++i){fac[i]=fac[i-1]*i;} ifac[n]=Pow(fac[n]);for(int i=n;i;--i){ifac[i-1]=ifac[i]*i;}
}
vector<Z>a,d,res;
int main(){Init();init(N-1);scanf("%d",&T);while(T--){scanf("%d",&n);a.resize(n+1);for(int i=0;i<=n;++i){scanf("%d",&a[i].x);a[i]*=fac[i];} for(int i=2;i<=n;++i){scanf("%d",&b[i]);}for(int i=2;i<=n;++i){scanf("%d",&c[i]);}d=work(2,n);reverse(d.begin(),d.end());res=d*a;for(int i=0;i<=n;++i){printf("%d%c",(res[i+n-1]*ifac[i]).x," \n"[i==n]);} }return 0;
}

 

  相关解决方案