「题解」挑战多项式
题目
题解
跟着它说的算就是了...这里主要是保存模板
如果要看具体操作，看这里
#include<cstdio>
#include<vector>
#include<utility>
#include<algorithm>
using namespace std;
#define rep(i,__l,__r) for(int i=__l,i##_end_=__r;i<=i##_end_;++i)
#define fep(i,__l,__r) for(int i=__l,i##_end_=__r;i>=i##_end_;--i)
#define writc(a,b) fwrit(a),putchar(b)
#define mp(a,b) make_pair(a,b)
#define ft first
#define sd second
#define LL long long
#define ull unsigned long long
#define uint unsigned int
#define pii pair<int,int>
#define Endl putchar('\n')
// #define FILEOI
// #define int long long
#ifdef FILEOI
#define MAXBUFFERSIZE 500000
inline char fgetc(){
static char buf[MAXBUFFERSIZE+5],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,MAXBUFFERSIZE,stdin),p1==p2)?EOF:*p1++;
}
#undef MAXBUFFERSIZE
#define cg (c=fgetc())
#else
#define cg (c=getchar())
#endif
template<class T>inline void qread(T& x){
char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
if(f)x=-x;
}
inline int qread(){
int x=0;char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
return f?-x:x;
}
template<class T,class... Args>inline void qread(T& x,Args&... args){qread(x),qread(args...);}
template<class T>inline T Max(const T x,const T y){return x>y?x:y;}
template<class T>inline T Min(const T x,const T y){return x<y?x:y;}
template<class T>inline T fab(const T x){return x>0?x:-x;}
inline int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
inline void getInv(int inv[],const int lim,const int MOD){
inv[0]=inv[1]=1;for(int i=2;i<=lim;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
}
template<class T>void fwrit(const T x){
if(x<0)return (void)(putchar('-'),fwrit(-x));
if(x>9)fwrit(x/10);putchar(x%10^48);
}
inline LL mulMod(const LL a,const LL b,const LL mod){//long long multiplie_mod
return ((a*b-(LL)((long double)a/mod*b+1e-8)*mod)%mod+mod)%mod;
}
const int MAXN=1<<21;
namespace __poly{
const int MOD=998244353,__G=3;
struct Z{
uint x;
Z(const uint _x=0):x(_x){}
inline Z operator +(const Z &rhs)const{return x+rhs.x<MOD?x+rhs.x:x+rhs.x-MOD;}
inline Z operator -(const Z &rhs)const{return x<rhs.x?x-rhs.x+MOD:x-rhs.x;}
inline Z operator -()const{return x?MOD-x:0;}
inline Z operator *(const Z &rhs)const{return (ull)x*rhs.x%MOD;}
inline Z operator +=(const Z &rhs){return x=x+rhs.x<MOD?x+rhs.x:x+rhs.x-MOD,*this;}
inline Z operator -=(const Z &rhs){return x=x<rhs.x?x-rhs.x+MOD:x-rhs.x,*this;}
inline Z operator *=(const Z &rhs){return x=(ull)x*rhs.x%MOD,*this;}
inline bool operator ==(const Z &rhs){return x==rhs.x;}
};
Z w[MAXN+5],Inv[MAXN+5];
vector<Z>ans;
vector<vector<Z> >p;
inline Z qkpow(Z a,int n){
//快速幂
Z ans=1;
for(;n>0;n>>=1,a=a*a)if(n&1)ans=ans*a;
return ans;
}
inline void Init(){
//解决初始化, 包括原根的次方以及 MAXN 以内的逆元
for(uint i=1;i<MAXN;i<<=1){w[i]=1;
Z t=qkpow(__G,(MOD-1)/i/2);
for(uint j=1;j<i;++j)w[i+j]=w[i+j-1]*t;
}
Inv[0]=Inv[1]=1;
rep(i,2,MAXN)Inv[i]=Inv[MOD%i]*(MOD-MOD/i);
}
#define pzz pair<Z,Z>
inline pzz Mul(pzz x,pzz y,Z f){
//计算 a+bω 部分的 pzz 的乘法
return mp(x.ft*y.ft+x.sd*y.sd*f,x.ft*y.sd+x.sd*y.ft);
}
inline Z Quadratic_residue(Z a){
//常数在模意义下的开根
if(a.x<=1)return a.x;
//使用欧拉定律判断有无解
if(qkpow(a,(MOD-1)>>1).x!=1)return -1;
Z x,f;
//找到一个 x 使得 x^2 - a 不是二次剩余
do x=(((ull)rand()<<15)^rand())%(a.x-1)+1;
while(qkpow(f=x*x-a,(MOD-1)>>1)==Z(1));
//初始变量
pzz ans=mp(1,0),t=mp(x,1);
//类似于快速幂
for(uint i=(MOD+1)>>1;i>0;i>>=1,t=Mul(t,t,f))if(i&1)ans=Mul(ans,t,f);
//返回较小根
return Min(ans.ft.x,MOD-ans.ft.x);
}
inline int Get(const int x){int n=1;while(n<=x)n<<=1;return n;}
inline void ntt(vector<Z>& f,const int n){
//ntt 的标准写法
static ull F[MAXN+5];
if((int)f.size()!=n)f.resize(n);//避免 RE
for(int i=0,j=0;i<n;++i){
//复制数组, 同时更改系数
F[i]=f[j].x;
for(int k=n>>1;(j^=k)<k;k>>=1);
}
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,t;
//通过数组指针进行访问
for(int k=j;k<j+i;++k,++W,++F0,++F1){
t=(*F1)*(W->x)%MOD;
(*F1)=*F0+MOD-t,(*F0)+=t;
}
}
rep(i,0,n-1)f[i]=F[i]%MOD;//将数组复制回去
}
inline void Intt(vector<Z>& f,const int n){
//逆 ntt 的运算, 和一般的写法有点不一样, 稍微研究一下
f.resize(n),reverse(f.begin()+1,f.end());
ntt(f,n);
Z n_Inv=qkpow(n,MOD-2);
//要乘逆元并取模 (取模在 * 里面重载了)
rep(i,0,n-1)f[i]=f[i]*n_Inv;
}
inline vector<Z> operator +(const vector<Z>& f,const vector<Z>& g){
//重载多项式的加法
vector<Z>ans=f;
rep(i,0,(int)f.size()-1)ans[i]+=g[i];
return ans;
}
inline vector<Z> operator *(const vector<Z>& f,const vector<Z>& g){
if((ull)f.size()*g.size()<=1000){
//在这种情况下, 暴力算要更快一些...
vector<Z>ans;
ans.resize(f.size()+g.size()-1);
rep(i,0,(int)f.size()-1)rep(j,0,(int)g.size()-1)ans[i+j]+=f[i]*g[j];
return ans;
}
//保存数组
static vector<Z> F,G;
F=f,G=g;
int p=Get((int)f.size()+(int)g.size()-2);//为什么是 -2 ?
ntt(F,p),ntt(G,p);
rep(i,0,p-1)F[i]*=G[i];
Intt(F,p);
return F.resize((int)f.size()+(int)g.size()-1),F;
}
vector<Z>& polyinv(const vector<Z>& f,int n=-1){
//返回逆元多项式
//默认为数组的大小
if(n==-1)n=f.size();
if(n==1){
//如果到了最后一项, 则直接求逆元
static vector<Z> ans;
return ans.clear(),ans.push_back(qkpow(f[0],MOD-2)),ans;
}
//剩下的是标准操作
vector<Z> &ans=polyinv(f,(n+1)>>1);
vector<Z> tmp(&f[0],&f[0]+n);
int p=Get(n*2-2);//为什么 -2?
ntt(tmp,p),ntt(ans,p);
rep(i,0,p-1)ans[i]=((Z)2-ans[i]*tmp[i])*ans[i];
Intt(ans,p);
return ans.resize(n),ans;
}
inline void polydiv(const vector<Z>& a,const vector<Z>& b,vector<Z>& d,vector<Z>& r){
//多项式除法
if(b.size()>a.size())return d.clear(),(void)(r=a);
vector<Z>A=a,B=b,iB;
int n=a.size(),m=b.size();
reverse(A.begin(),A.end()),reverse(B.begin(),B.end());
B.resize(n-m+1),iB=polyinv(B);
d=A*iB;
d.resize(n-m+1),reverse(d.begin(),d.end());
r=b*d,r.resize(m-1);
rep(i,0,m-2)r[i]=a[i]-r[i];
}
inline vector<Z> Derivative(const vector<Z>& a){
//函数的导数
vector<Z>ans((int)a.size()-1);
rep(i,1,(int)a.size()-1)ans[i-1]=a[i]*i;
return ans;
}
inline vector<Z> Integral(const vector<Z>& a){
//计算微分
vector<Z>ans(a.size()+1);
rep(i,0,a.size()-1)ans[i+1]=a[i]*Inv[i+1];
return ans;
}
inline vector<Z>polyln(const vector<Z>& f){
//自然对数
vector<Z>ans=Derivative(f)*polyinv(f);
ans.resize((int)f.size()-1);
return Integral(ans);
}
vector<Z> polyexp(const vector<Z>& f,int n=-1){
if(n==-1)n=f.size();
if(n==1)return {1};
vector<Z> ans=polyexp(f,(n+1)>>1),tmp;
ans.resize(n),tmp=polyln(ans);
for(vector<Z>::iterator it=tmp.begin();it!=tmp.end();++it)
(*it)=-(*it);//此处 C++11 的写法在 luoguOJ 上的评测结果不同...
// for(Z &i:tmp)i=-i;
++tmp[0].x;
ans=ans*(tmp+f);
return ans.resize(n),ans;
}
vector<Z> polysqrt(const vector<Z>& f,int n=-1){
if(n==-1)n=f.size();
vector<Z>ans;
if(n==1)return ans.push_back(Quadratic_residue(f[0])),ans;
ans=polysqrt(f,(n+1)>>1);
vector<Z>tmp(&f[0],&f[0]+n);
ans.resize(n),ans=ans+tmp*polyinv(ans);
for(vector<Z>::iterator it=ans.begin();it!=ans.end();++it)
it->x=(it->x&1)?(it->x+MOD)>>1:it->x>>1;
return ans;
}
inline vector<Z> polyqkpow(const vector<Z>& f,const Z k){
vector<Z>ans=polyln(f);
for(vector<Z>::iterator it=ans.begin();it!=ans.end();++it)
(*it)=(*it)*k;//此处也有 C11 的问题
// for(Z &i:ans)i=i*k;
return polyexp(ans);
}
void Evaluate_Interpolate_Init(int l, int r, int t, const vector<Z> &a){
if(l==r)return p[t].clear(),p[t].push_back(-a[l]),p[t].push_back(1);
int mid=(l+r)/2,k=t<<1;
Evaluate_Interpolate_Init(l,mid,k,a),Evaluate_Interpolate_Init(mid+1,r,k|1,a);
p[t]=p[k]*p[k|1];
}
void Evaluate(int l,int r,int t,const vector<Z> &f,const vector<Z> &a){
if(r-l+1<=512){
for(int i=l;i<=r;++i){
Z x=0,a1=a[i],a2=a[i]*a[i], a3=a[i]*a2, a4=a[i]*a3, a5=a[i]*a4, a6=a[i]*a5, a7=a[i]*a6, a8=a[i]*a7;
int j=f.size();
while(j>=8)x=x*a8+f[j-1]*a7+f[j-2]*a6+f[j-3]*a5+f[j-4]*a4+f[j-5]*a3+f[j-6]*a2+f[j-7]*a1+f[j-8], j-=8;
while(j--)x=x*a[i]+f[j];
ans.push_back(x);
}
return;
}
vector<Z>tmp;
polydiv(f,p[t],tmp,tmp);
Evaluate(l,(l+r)/2,t<<1,tmp,a),Evaluate((l+r)/2+1,r,t<<1|1,tmp, a);
}
inline vector<Z> Evaluate(const vector<Z> &f,const vector<Z> &a,int flag=-1){
if(flag==-1)p.resize(a.size()<<2),Evaluate_Interpolate_Init(0,a.size()-1,1,a);
return ans.clear(),Evaluate(0,a.size()-1,1,f,a),ans;
}
vector<Z> Interpolate(int l,int r,int t,const vector<Z> &x,const vector<Z> &f){
if(l==r)return {f[l]};
int mid=(l+r)/2,k=t<<1;
return Interpolate(l,mid,k,x,f)*p[k|1]+Interpolate(mid+1,r,k|1,x,f)*p[k];
}
inline vector<Z> Interpolate(const vector<Z> &x,const vector<Z> &y){
int n=x.size();
p.resize(n<<2),Evaluate_Interpolate_Init(0,n-1,1,x);
vector<Z> f=Evaluate(Derivative(p[1]),x,0);
for(int i=0;i<n;++i)f[i]=y[i]*qkpow(f[i],MOD-2);
return Interpolate(0,n-1,1,x,f);
}
inline vector<Z> solve(const vector<Z>&f,const int k){
vector<Z>ans=polyexp(Integral(polyinv(polysqrt(f))));
rep(i,1,f.size()-1)ans[i]=f[i]-ans[i];
ans=polyln(ans);
++ans[0].x;
return Derivative(polyqkpow(ans,k));
}
};
int n,k;
vector<__poly::Z>f;
signed main(){
#ifdef FILEOI
freopen("file.in","r",stdin);
freopen("file.out","w",stdout);
#endif
__poly::Init();
qread(n,k);
rep(i,0,n)f.push_back(qread());
f=solve(f,k);f.resize(n);
while(!f[f.size()-1].x)f.erase(f.end()-1);
rep(i,0,n-1)writc(f[i].x,' ');
Endl;
return 0;
}
