Code
#include<bits/stdc++.h>
#define endl '\n'
#define rep(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
using namespace std;
const int P=998244353,G=3,LIMIT=50;
union mi{
int w;
mi(){w=0;}
mi(int u){u%=P,w=u+(u<0?P:0);}
explicit operator int()const{return w;}
};
mi operator+(const mi&a,const mi&b){
return mi(a.w+b.w-(a.w+b.w>=P?P:0));
}
mi operator-(const mi&a,const mi&b){
return mi(a.w-b.w+(a.w<b.w?P:0));
}
mi operator*(const mi&a,const mi&b){
return mi(1ll*a.w*b.w%P);
}
bool operator==(const mi&a,const mi&b){
return a.w==b.w;
}
bool operator!=(const mi&a,const mi&b){
return a.w!=b.w;
}
bool operator<(const mi&a,const mi&b){
return a.w<b.w;
}
bool operator>(const mi&a,const mi&b){
return a.w>b.w;
}
bool operator>=(const mi&a,const mi&b){
return a.w>=b.w;
}
bool operator<=(const mi&a,const mi&b){
return a.w<=b.w;
}
mi&operator+=(mi&a,const mi&b){
a.w=a.w+b.w-(a.w+b.w>=P?P:0);
return a;
}
mi&operator-=(mi&a,const mi&b){
a.w=a.w-b.w+(a.w<b.w?P:0);
return a;
}
mi&operator*=(mi&a,const mi&b){
a.w=1ll*a.w*b.w%P;
return a;
}
mi operator-(const mi&a){
return mi(a.w?P-a.w:0);
}
mi&operator++(mi&a){
a.w=a.w+1-(a.w+1>=P?P:0);
return a;
}
mi&operator--(mi&a){
a.w=a.w-1+(a.w?0:P);
return a;
}
typedef vector<mi> vec;
struct IO_Tp {
const static int _I_Buffer_Size = 2 << 22;
char _I_Buffer[_I_Buffer_Size], *_I_pos = _I_Buffer;
const static int _O_Buffer_Size = 2 << 22;
char _O_Buffer[_O_Buffer_Size], *_O_pos = _O_Buffer;
IO_Tp() { fread(_I_Buffer, 1, _I_Buffer_Size, stdin); }
~IO_Tp() { fwrite(_O_Buffer, 1, _O_pos - _O_Buffer, stdout); }
IO_Tp &operator>>(int &res) {
int f=1;
while (!isdigit(*_I_pos)&&(*_I_pos)!='-') ++_I_pos;
if(*_I_pos=='-')f=-1,++_I_pos;
res = *_I_pos++ - '0';
while (isdigit(*_I_pos)) res = res * 10 + (*_I_pos++ - '0');
res*=f;
return *this;
}
IO_Tp &operator>>(mi &res) {
int f=1;
while (!isdigit(*_I_pos)&&(*_I_pos)!='-') ++_I_pos;
if(*_I_pos=='-')f=-1,++_I_pos;
res = *_I_pos++ - '0';
while (isdigit(*_I_pos)) res = res * 10 + (*_I_pos++ - '0');
res*=f;
return *this;
}
IO_Tp &operator<<(int n) {
if(n<0)*_O_pos++='-',n=-n;
static char _buf[10];
char *_pos = _buf;
do
*_pos++ = '0' + n % 10;
while (n /= 10);
while (_pos != _buf) *_O_pos++ = *--_pos;
return *this;
}
IO_Tp &operator<<(mi n) {
if(n<0)*_O_pos++='-',n=-n;
static char _buf[10];
char *_pos = _buf;
do
*_pos++ = '0' + n.w % 10;
while (n.w /= 10);
while (_pos != _buf) *_O_pos++ = *--_pos;
return *this;
}
IO_Tp &operator<<(char ch) {
*_O_pos++ = ch;
return *this;
}
} IO;
void chkmax(int &x,int y){if(x<y)x=y;}
void chkmin(int &x,int y){if(x>y)x=y;}
int qpow(int a,int k,int p=P){
int ret=1;
while(k){
if(k&1)ret=1ll*ret*a%p;
a=1ll*a*a%p;
k>>=1;
}
return ret;
}
int norm(int x){return x>=P?x-P:x;}
int reduce(int x){return x<0?x+P:x;}
void add(int&x,int y){if((x+=y)>=P)x-=P;}
struct Maths{
int n;
vec fac,invfac,inv;
void build(int n){
this->n=n;
fac.resize(n+1);
invfac.resize(n+1);
inv.resize(n+1);
fac[0]=1;
rep(k,1,n)fac[k]=fac[k-1]*k;
inv[1]=inv[0]=1;
rep(k,2,n)inv[k]=-(P/k)*inv[P%k];
invfac[0]=1;
rep(k,1,n)invfac[k]=invfac[k-1]*inv[k];
}
Maths(){build(1);}
void chk(int k){
int lmt=n;
if(k>lmt){while(k>lmt)lmt<<=1;build(lmt);}
}
mi cfac(int k){return chk(k),fac[k];}
mi cifac(int k){return chk(k),invfac[k];}
mi cinv(int k){return chk(k),inv[k];}
mi binom(int n,int m){
if(m<0||m>n)return 0;
return cfac(n)*cifac(m)*cifac(n-m);
}
}math;
struct NumberTheory{
mt19937 rng;
NumberTheory():rng(chrono::steady_clock::now().time_since_epoch().count()){}
void exgcd(int a,int b,int&x,int&y){
if(!b)return x=1,y=0,void();
exgcd(b,a%b,y,x);
y-=a/b*x;
}
int inv(int a,int p=P){
int x,y;
exgcd(a,p,x,y);
if(x<0)x+=p;
return x;
}
template<class Integer>
bool quadRes(Integer a,Integer b){
if(a<=1)return 1;
while(a%4==0)a/=4;
if(a%2==0)return (b%8==1||b%8==7)==quadRes(a/2,b);
return ((a-1)%4==0||(b-1)%4==0)==quadRes(b%a,a);
}
int sqrt(int x,int p=P){
if(p==2||x<=1)return x;
int w,v,k=(p+1)/2;
do w=rng()%p;while(quadRes((v=(1ll*w*w%p-x+p)%p),p));
pair<int,int>res(1,0),a(w,1);
while(k){
if(k&1)res=make_pair((1ll*res.first*a.first%p+1ll*res.second*a.second%p*v%p)%p,
(1ll*res.first*a.second%p+1ll*res.second*a.first%p)%p);
if(k>>=1)
a=make_pair((1ll*a.first*a.first%p+1ll*a.second*a.second%p*v%p)%p,
2ll*a.first*a.second%p);
}
return min(res.first,p-res.first);
}
}nt;
struct NTT{
int L,brev[1<<11];
vec root;
NTT():L(-1){
rep(i,1,(1<<11)-1)brev[i]=brev[i>>1]>>1|((i&1)<<10);
}
void preproot(int l){
L=l;
root.resize(2<<L);
rep(i,0,L){
mi*w=root.data()+(1<<i);
w[0]=1;
int omega=qpow(G,(P-1)>>i);
rep(j,1,(1<<i)-1)w[j]=w[j-1]*omega;
}
}
void dft(mi*a,int lgn,int d=1){
if(L<lgn)preproot(lgn);
int n=1<<lgn;
rep(i,0,n-1){
int rev=(brev[i>>11]|(brev[i&((1<<11)-1)]<<11))>>((11<<1)-lgn);
if(i<rev)swap(a[i],a[rev]);
}
for(int i=1;i<n;i<<=1){
mi*w=root.data()+(i<<1);
for(int j=0;j<n;j+=i<<1)rep(k,0,i-1){
mi aa=w[k]*a[i+j+k];
a[i+j+k]=a[j+k]-aa;
a[j+k]+=aa;
}
}
if(d==-1){
reverse(a+1,a+n);
int inv=nt.inv(n);
rep(i,0,n-1)a[i]*=inv;
}
}
}ntt;
struct poly{
vec a;
poly(mi v):a(1){a[0]=v;}
poly(int v=0):a(1){a[0]=v;}
poly(const vec&a):a(a){}
poly(initializer_list<mi>init):a(init){}
mi operator[](int k)const{return k<a.size()?a[k]:0;}
mi&operator[](int k){
if(k>=a.size())a.resize(k+1);
return a[k];
}
int deg()const{return a.size()-1;}
void redeg(int d){a.resize(d+1);}
poly slice(int d)const{
if(d<a.size())return vec(a.begin(),a.begin()+d+1);
vec res(a);
res.resize(d+1);
return res;
}
mi*base(){return a.data();}
const mi*base()const{return a.data();}
poly println(FILE* fp=stdout)const{
fprintf(fp,"%d",a[0]);
rep(i,1,a.size()-1)fprintf(fp," %d",a[i]);
fputc('\n',fp);
return *this;
}
poly operator+(const poly&rhs)const{
vec res(max(a.size(),rhs.a.size()));
rep(i,0,res.size()-1)res[i]=operator[](i)+rhs[i];
return res;
}
poly operator-()const{
poly ret(a);
rep(i,0,a.size()-1)ret[i]=-ret[i];
return ret;
}
poly operator-(const poly&rhs)const{return operator+(-rhs);}
poly operator*(const poly&rhs)const;//done
poly operator/(const poly&rhs)const;//done
poly operator%(const poly&rhs)const;//done
poly der()const;//done
poly integ()const;//done
poly inv()const;//done
poly sqrt()const;//done
poly ln()const;//done
poly exp()const;//done
poly sin()const;//done
poly cos()const;//done
poly arcsin()const;//done
poly arctan()const;//done
pair<poly,poly>sqrti()const;//done
pair<poly,poly>expi()const;//done
poly quo(const poly&rhs)const;//done
pair<poly,poly>iquo(const poly&rhs)const;
pair<poly,poly>div(const poly&rhs)const;//done
poly taylor(int k)const;
poly pow(int k)const;//done
poly exp(int k)const;//done
poly shift(int k)const;//done
mi LHRR(const vec&a,int n)const;//done
mi LNRR(const poly&P,const vec&a,int n)const;
};
poly zeroes(int deg){return vec(deg+1);}
struct Newton{
void inv(const poly&f,const poly&nttf,poly&g,const poly&nttg,int t){
int n=1<<t;
poly prod=nttf;
rep(i,0,(n<<1)-1)prod[i]=prod[i]*nttg[i];
ntt.dft(prod.base(),t+1,-1);
rep(i,0,n-1)prod[i]=0;
ntt.dft(prod.base(),t+1,1);
rep(i,0,(n<<1)-1)prod[i]=prod[i]*nttg[i];
ntt.dft(prod.base(),t+1,-1);
rep(i,0,n-1)prod[i]=0;
g=g-prod;
}
void inv(const poly&f,const poly&nttf,poly&g,int t){
poly nttg=g;
nttg.redeg((2<<t)-1);
ntt.dft(nttg.base(),t+1,1);
inv(f,nttf,g,nttg,t);
}
void inv(const poly&f,poly&g,int t){
poly nttg=g;
nttg.redeg((2<<t)-1);
ntt.dft(nttg.base(),t+1,1);
poly nttf=f;
nttf.redeg((2<<t)-1);
ntt.dft(nttf.base(),t+1,1);
inv(f,nttf,g,nttg,t);
}
void sqrt(const poly&f,poly&g,poly&nttg,poly&h,int t){
rep(i,0,(1<<t)-1)nttg[i]=qpow(int(nttg[i]),2);
ntt.dft(nttg.base(),t,-1);
nttg=nttg-f;
rep(i,0,(1<<t)-1)nttg[i+(1<<t)]+=nttg[i];
memset(nttg.base(),0,sizeof(mi)<<t);
ntt.dft(nttg.base(),t+1,1);
poly tmp=h;
tmp.redeg((2<<t)-1);
ntt.dft(tmp.base(),t+1,1);
rep(i,0,(2<<t)-1)tmp[i]*=nttg[i];
ntt.dft(tmp.base(),t+1,-1);
memset(tmp.base(),0,sizeof(mi)<<t);
g=g-tmp*mi(499122177);
}
void exp(const poly&f,poly&g,poly&nttg,poly&h,int t){
poly ntth(h);
ntt.dft(ntth.base(),t,1);
poly dg=g.der().slice((1<<t)-1);
ntt.dft(dg.base(),t,1);
poly tmp=zeroes((1<<t)-1);
rep(i,0,(1<<t)-1){
tmp[i]=nttg[i<<1]*ntth[i];
dg[i]=dg[i]*ntth[i];
}
ntt.dft(tmp.base(),t,-1);
ntt.dft(dg.base(),t,-1);
if(--tmp[0]<0)tmp[0]=P-1;
dg.redeg((2<<t)-1);
poly df0=f.der().slice((1<<t)-1);
df0[(1<<t)-1]=0;
rep(i,0,(1<<t)-1)if((dg[i|1<<t]=dg[i]-df0[i])<0)dg[i|1<<t]+=P;
memcpy(dg.base(),df0.base(),sizeof(mi)*((1<<t)-1));
tmp.redeg((2<<t)-1);
ntt.dft(tmp.base(),t+1,1);
df0.redeg((2<<t)-1);
ntt.dft(df0.base(),t+1,1);
rep(i,0,(2<<t)-1)df0[i]*=df0[i]*tmp[i];
ntt.dft(df0.base(),t+1,-1);
memcpy(df0.base()+(1<<t),df0.base(),sizeof(mi)<<t);
memset(df0.base(),0,sizeof(mi)<<t);
dg=(dg-df0).integ().slice((2<<t)-1)-f;
ntt.dft(dg.base(),t+1,1);
rep(i,0,(2<<t)-1)tmp[i]=dg[i]*nttg[i];
ntt.dft(tmp.base(),t+1,-1);
g.redeg((2<<t)-1);
rep(i,1<<t,(2<<t)-1)if(tmp[i]!=0)g[i]=P-tmp[i];
}
}nit;
struct SemiRelaxedConvolution{
template<class Function>
void run(const vec&a,vec&b,const Function&relax){
int n=a.size()-1;
function<void(int,int)>cdq = [&](int l,int r){
if(r-l<=LIMIT){
rep(i,l,r){
rep(j,l,i-1)b[i]+=b[j]*a[i-j];
relax(i);
}
return;
}
int lg=31-__builtin_clz(r-l),d=(r-l)/lg+1,lgd=0;
while((1<<lgd)<d)++lgd;++lgd;
vec top(lg<<(lgd+1));
rep(i,0,lg-1){
copy(a.begin()+i*d,a.begin()+min((i+2)*d,n+1),top.begin()+(i<<lgd));
ntt.dft(top.data()+(i<<lgd),lgd,1);
}
rep(i,0,lg-1){
if(i)ntt.dft(top.data()+((lg+i)<<lgd),lgd,-1);
rep(j,0,min(d,r-l-i*d+1)-1)b[l+i*d+j]+=top[((lg+i)<<lgd)+d+j];
cdq(l+i*d,min(l+(i+1)*d-1,r));
if(i+1<lg){
copy(b.begin()+l+i*d,b.begin()+min(l+(i+1)*d,n+1),top.begin()+((lg+i)<<lgd));
fill(top.data()+((lg+i)<<lgd)+d,top.data()+((lg+i+1)<<lgd),0);
ntt.dft(top.data()+((lg+i)<<lgd),lgd,1);
}
rep(j,i+1,lg-1)rep(k,0,(1<<lgd)-1)
top[((lg+j)<<lgd)+k]+=top[((j-i-1)<<lgd)+k]*top[((lg+i)<<lgd)+k];
}
};
cdq(0,n);
}
}src;
struct Transposition{
vec _mul(int l,vec res,const poly&b){
vec tmp(1<<l);
memcpy(tmp.data(),b.a.data(),sizeof(mi)*(b.deg()+1));
reverse(tmp.begin()+1,tmp.end());
ntt.dft(tmp.data(),l,1);
rep(i,0,(1<<l)-1)res[i]*=tmp[i];
ntt.dft(res.data(),l,-1);
return res;
}
poly bfmul(const poly&a,const poly&b){
int n=a.deg(),m=b.deg();
poly ret=zeroes(n-m);
rep(i,0,n-m)rep(j,0,m)ret[i]+=a[i+j]*b[j];
return ret;
}
poly mul(const poly&a,const poly&b){
if(a.deg()<b.deg())return 0;
if(a.deg()<=LIMIT)return bfmul(a,b);
int l=0;
while((1<<l)<=a.deg())++l;
vec res(1<<l);
memcpy(res.data(),a.a.data(),sizeof(mi)*(a.deg()+1));
ntt.dft(res.data(),l,1);
res=_mul(l,res,b);
res.resize(a.deg()-b.deg()+1);
return res;
}
pair<poly,poly>mul2(const poly&a,const poly&b1,const poly&b2){
if(a.deg()<=LIMIT)return make_pair(bfmul(a,b1),bfmul(a,b2));
int l=0;
while((1<<l)<=a.deg())++l;
vec fa(1<<l);
memcpy(fa.data(),a.a.data(),sizeof(mi)*(a.deg()+1));
ntt.dft(fa.data(),l,1);
vec res1=_mul(l,fa,b1),res2=_mul(l,fa,b2);
res1.resize(a.deg()-b1.deg()+1);
res2.resize(a.deg()-b2.deg()+1);
return make_pair(res1,res2);
}
vector<int>ls,rs,pos;
vector<poly>p,q;
void _build(int n){
ls.assign(n*2-1,0);
rs.assign(n*2-1,0);
p.assign(n*2-1,0);
q.assign(n*2-1,0);
pos.resize(n);
int cnt=0;
function<int(int,int)>dfs=[&](int l,int r){
if(l==r){
pos[l]=cnt;
return cnt++;
}
int ret=cnt++;
int mid=l+r>>1;
ls[ret]=dfs(l,mid);
rs[ret]=dfs(mid+1,r);
return ret;
};
dfs(0,n-1);
}
vec _eval(vec f,const vec&x){
int n=f.size();
_build(n);
rep(i,0,n-1)q[pos[i]]={1,-x[i]};
Rep(i,n*2-2,0)if(ls[i])q[i]=q[ls[i]]*q[rs[i]];
f.resize(n*2);
p[0]=mul(f,q[0].inv());
rep(i,0,n*2-2)if(ls[i])tie(p[ls[i]],p[rs[i]])=mul2(p[i],q[rs[i]],q[ls[i]]);
vec ret(n);
rep(i,0,n-1)ret[i]=p[pos[i]][0];
return ret;
}
vec eval(const poly&f,const vec&x){
int n=f.deg()+1,m=x.size();
vec tmpf=f.a,tmpx=x;
tmpf.resize(max(n,m));
tmpx.resize(max(n,m));
vec ret=_eval(tmpf,tmpx);
ret.resize(m);
return ret;
}
poly inter(const vec&x,const vec&y){
int n=x.size();
_build(n);
rep(i,0,n-1)q[pos[i]]={1,-x[i]};
Rep(i,n*2-2,0)if(ls[i])q[i]=q[ls[i]]*q[rs[i]];
poly tmp=q[0];
reverse(tmp.a.begin(),tmp.a.end());
vec f=tmp.der().a;
f.resize(n*2);
p[0]=mul(f,q[0].inv());
rep(i,0,n*2-2)if(ls[i])tie(p[ls[i]],p[rs[i]])=mul2(p[i],q[rs[i]],q[ls[i]]);
rep(i,0,n-1)p[pos[i]]=nt.inv((int)p[pos[i]][0])*y[i];
rep(i,0,n*2-2)reverse(q[i].a.begin(),q[i].a.end());
Rep(i,n*2-2,0)if(ls[i])p[i]=p[ls[i]]*q[rs[i]]+p[rs[i]]*q[ls[i]];
return p[0];
}
}tp;
poly operator "" _z(unsigned long long a){return {0,(int)a};}
poly operator+(int v,const poly&rhs){return poly(v)+rhs;}
poly operator*(int v,const poly&rhs){
poly ret=zeroes(rhs.deg());
rep(i,0,rhs.deg())ret[i]=rhs[i]*v;
return ret;
}
poly operator*(const poly&lhs,int v){return v*lhs;}
poly poly::operator*(const poly&rhs)const{
int n=deg(),m=rhs.deg();
if(n<=10||m<=10||n+m<=LIMIT){
poly ret=zeroes(n+m);
rep(i,0,n)rep(j,0,m)ret[i+j]+=a[i]*rhs[j];
return ret;
}
n+=m;
int l=0;
while((1<<l)<=n)++l;
vec res(1<<l),tmp(1<<l);
memcpy(res.data(),base(),a.size()*sizeof(mi));
ntt.dft(res.begin().base(),l,1);
memcpy(tmp.data(),rhs.base(),rhs.a.size()*sizeof(mi));
ntt.dft(tmp.begin().base(),l,1);
rep(i,0,(1<<l)-1)res[i]*=tmp[i];
ntt.dft(res.data(),l,-1);
res.resize(n+1);
return res;
}
poly poly::inv()const{
poly g=nt.inv((int)a[0]);
for(int t=0;(1<<t)<=deg();++t)nit.inv(slice((2<<t)-1),g,t);
g.redeg(deg());
return g;
}
poly poly::taylor(int k)const{
int n=deg();
poly t=zeroes(n);
math.chk(n);
rep(i,0,n)t[n-i]=a[i]*math.fac[i];
int pw=1;
poly help=vec(math.invfac.begin(),math.invfac.begin()+n+1);
rep(i,0,n){
help[i]*=pw;
pw*=k;
}
t=t*help;
rep(i,0,n)help[i]=t[n-i]*math.invfac[i];
return help;
}
poly poly::pow(int k)const{
if(k==0)return 1;
if(k==1)return *this;
int n=deg()*k;
int lgn=0;
while((1<<lgn)<=n)++lgn;
vec val=a;
val.resize(1<<lgn);
ntt.dft(val.data(),lgn,1);
rep(i,0,(1<<lgn)-1)val[i]=qpow((int)val[i],k);
ntt.dft(val.data(),lgn,-1);
return val;
}
poly poly::der()const{
if(deg()==0)return 0;
vec res(deg());
rep(i,0,deg()-1)res[i]=a[i+1]*mi(i+1);
return res;
}
poly poly::integ()const{
vec res(deg()+2);
rep(i,0,deg())res[i+1]=a[i]*math.cinv(i+1);
return res;
}
poly poly::quo(const poly&rhs)const{
if(rhs.deg()==0)return a[0]*(mi)nt.inv((int)rhs[0]);
poly g=(mi)nt.inv((int)rhs[0]);
int t=0,n;
for(n=1;(n<<1)<=rhs.deg();++t,n<<=1)nit.inv(rhs.slice((n<<1)-1),g,t);
poly nttg=g;
nttg.redeg((n<<1)-1);
ntt.dft(nttg.base(),t+1,1);
poly eps1=rhs.slice((n<<1)-1);
ntt.dft(eps1.base(),t+1,1);
rep(i,0,(n<<1)-1)eps1[i]*=nttg[i];
ntt.dft(eps1.base(),t+1,-1);
memcpy(eps1.base(),eps1.base()+n,sizeof(mi)<<t);
memset(eps1.base()+n,0,sizeof(mi)<<t);
ntt.dft(eps1.base(),t+1,1);
poly h0=slice(n-1);
h0.redeg((n<<1)-1);
ntt.dft(h0.base(),t+1,1);
poly h0g0=zeroes((n<<1)-1);
rep(i,0,(n<<1)-1)h0g0[i]=h0[i]*nttg[i];
ntt.dft(h0g0.base(),t+1,-1);
poly h0eps1=zeroes((n<<1)-1);
rep(i,0,(n<<1)-1)h0eps1[i]=h0[i]*eps1[i];
ntt.dft(h0eps1.base(),t+1,-1);
rep(i,0,n-1)h0eps1[i]=operator[](i+n)-h0eps1[i];
memset(h0eps1.base()+n,0,sizeof(mi)<<t);
ntt.dft(h0eps1.base(),t+1,1);
rep(i,0,(n<<1)-1)h0eps1[i]*=nttg[i];
ntt.dft(h0eps1.base(),t+1,-1);
memcpy(h0eps1.base()+n,h0eps1.base(),sizeof(mi)<<t);
memset(h0eps1.base(),0,sizeof(mi)<<t);
return (h0g0+h0eps1).slice(rhs.deg());
}
poly poly::ln()const{
if(deg()==0)return 0;
return der().quo(slice(deg()-1)).integ();
}
pair<poly,poly> poly::sqrti()const{
poly g=nt.sqrt((int)a[0]),h=nt.inv((int)g[0]),nttg=g;
for(int t=0;(1<<t)<=deg();++t){
nit.sqrt(slice((2<<t)-1),g,nttg,h,t);
if((2<<t)<=deg()){
nttg=g;
ntt.dft(nttg.base(),t+1,1);
nit.inv(g,nttg,h,t);
}
}
return make_pair(g.slice(deg()),h.slice(deg()));
}
poly poly::sqrt()const{
poly g=nt.sqrt((int)a[0]),h=nt.inv((int)g[0]),nttg=g;
for(int t=0;(1<<t)<=deg();++t){
nit.sqrt(slice((2<<t)-1),g,nttg,h,t);
if((2<<t)<=deg()){
nttg=g;
ntt.dft(nttg.base(),t+1,1);
nit.inv(g,nttg,h,t);
}
}
return g.slice(deg());
}
poly poly::exp()const{
vec der(a),ret(a.size());
rep(i,0,a.size()-1)der[i]*=i;
src.run(der,ret,[&](int i){
if(i==0)ret[0]=1;
else ret[i]*=math.cinv(i);
});
return ret;
}
pair<poly,poly>poly::expi()const{
poly g=1,h=1,nttg={1,1};
for(int t=0;(1<<t)<=deg();++t){
nit.exp(slice((2<<t)-1),g,nttg,h,t);
nttg=g;
nttg.redeg((4<<t)-1);
ntt.dft(nttg.base(),t+2);
poly f2n=zeroes((2<<t)-1);
rep(i,0,(2<<t)-1)f2n[i]=nttg[i<<1];
nit.inv(g,f2n,h,t);
}
return make_pair(g.slice(deg()),h.slice(deg()));
}
poly poly::exp(int k)const{
int lead,lz=0;
while(lz<deg()&&a[lz]==0)++lz;
if(lz==deg()&&a[lz]==0)return *this;
lead=(int)a[lz];
if(1ll*lz*k>deg())return zeroes(deg());
poly part=poly(vec(a.begin()+lz,a.begin()+deg()-lz*(k-1)+1))*nt.inv(lead);
part=(part.ln()*k).exp()*qpow(lead,k);
vec ret(deg()+1);
memcpy(ret.data()+lz*k,part.base(),sizeof(mi)*(deg()-lz*k+1));
return ret;
}
poly poly::sin()const{
int imag=qpow(G,(P-1)>>2);
poly g=operator*(imag),eif,ieif;
tie(eif,ieif)=g.expi();
return (eif-ieif)*nt.inv(imag*2%P);
}
poly poly::cos()const{
int imag=qpow(G,(P-1)>>2);
poly g=operator*(imag),eif,ieif;
tie(eif,ieif)=g.expi();
return (eif+ieif)*nt.inv(2);
}
poly poly::arcsin()const{
return der().quo((poly(1)-pow(2).slice(deg()-1)).sqrt()).integ();
}
poly poly::arctan()const{
return der().quo(poly(1)+pow(2).slice(deg()-1)).integ();
}
poly poly::operator/(const poly&rhs)const{
int n=deg(),m=rhs.deg();
if(n<m)return 0;
poly ta(vec(a.rbegin(),a.rend())),tb(vec(rhs.a.rbegin(),rhs.a.rend()));
ta.redeg(n-m),tb.redeg(n-m);
poly q=ta.quo(tb);
reverse(q.a.begin(),q.a.end());
return q;
}
pair<poly,poly>poly::div(const poly&rhs)const{
if(deg()<rhs.deg())return make_pair(0,*this);
int n=deg(),m=rhs.deg();
poly q=operator/(rhs),r;
int lgn=0;
while((1<<lgn)<rhs.deg())++lgn;
int t=(1<<lgn)-1;
r=zeroes(t);
poly tmp=zeroes(t);
rep(i,0,m)r[i&t]+=rhs[i];
rep(i,0,n-m)tmp[i&t]+=q[i];
ntt.dft(r.base(),lgn,1);
ntt.dft(tmp.base(),lgn,1);
rep(i,0,t)tmp[i]*=r[i];
ntt.dft(tmp.base(),lgn,-1);
memset(r.base(),0,sizeof(mi)<<lgn);
rep(i,0,n)r[i&t]+=a[i];
rep(i,0,m-1)r[i]-=tmp[i];
return make_pair(q,r.slice(m-1));
}
poly poly::shift(int k)const{
poly g=zeroes(deg()+k);
rep(i,0,k-1)g[i]=0;
rep(i,max(0,-k),deg())g[i+k]=a[i];
return g;
}
int quo2(int x){return ((x&1)?(x+P):x)>>1;}
mi poly::LHRR(const vec&a,int n)const{
int k=a.size();
poly p=a,q=zeroes(k);
q[0]=1;
rep(i,1,k)q[i]=-this->a[i-1];
p=(p*q).slice(k-1);
int logk=0;
while((1<<logk)<=k)++logk;++logk;
int h=1<<(logk-1);
vec sp,sq;
if(ntt.L<logk)ntt.preproot(logk);
while(n>=k){
vec pdft(1<<logk),qdft(1<<logk);
copy_n(p.base(),k,pdft.data());
copy_n(q.base(),k+1,qdft.data());
mi*w=ntt.root.data()+(1<<logk);
if(sp.empty()){
ntt.dft(pdft.data(),logk,1);
ntt.dft(qdft.data(),logk,1);
}else{
vec odd(h);
rep(i,0,k-1)odd[i]=w[i]*p[i];
ntt.dft(odd.data(),logk-1,1);
rep(i,0,h-1){
pdft[i<<1]=sp[i];
pdft[i<<1|1]=odd[i];
}
odd.assign(h,0);
rep(i,0,k)odd[i]=w[i]*q[i];
ntt.dft(odd.data(),logk-1,1);
rep(i,0,h-1){
qdft[i<<1]=sq[i];
qdft[i<<1|1]=odd[i];
}
}
sp.resize(h);
sq.resize(h);
if(n%2==0)rep(i,0,h-1)sp[i]=quo2(int(pdft[i]*qdft[i^h]+pdft[i^h]*qdft[i]));
else{
rep(i,0,h-1)sp[i]=quo2(int(pdft[i]*qdft[i^h]-pdft[i^h]*qdft[i]));
rep(i,1,h-1)sp[i]=sp[i]*w[(1<<logk)-i];
}
vec tmp=sp;
ntt.dft(tmp.data(),logk-1,-1);
copy_n(tmp.data(),k,p.base());
rep(i,0,h-1)sq[i]=qdft[i]*qdft[i^h];
tmp=sq;
ntt.dft(tmp.data(),logk-1,-1);
copy_n(tmp.data(),k+1,q.base());
n>>=1;
}
p=p.quo(q);
return p[n];
}
mi poly::LNRR(const poly&P,const vec&a,int n)const{
int m=P.deg(),k=a.size();
vec d(m+1);
rep(i,0,m)d[i]=i+k;
vec b=tp.eval(P,d);
b.resize(k+m+1);
Rep(i,k+m,k)b[i]=b[i-k];
memset(b.data(),0,sizeof(mi)*k);
auto res=operator*(poly(a)).shift(1);
rep(i,0,k-1)b[i]=a[i]-res[i];
poly B(b),F=shift(1),A=B.quo((poly(1)-F).slice(m+k));
F.redeg(m+1);
rep(i,0,m+1)
if((m+1-i)&1)F[i]=-math.binom(m+1,i);
else F[i]=math.binom(m+1,i);
int T=m+k;
poly f=((shift(1)-1)*F);
if(f[0]==1)f=-f;
return f.shift(-1).slice(T).LHRR(A.a,n);
}
poly poly::operator%(const poly&rhs)const{
if(deg()<rhs.deg())return *this;
return div(rhs).second;
}
template<class T>
IO_Tp& operator>>(IO_Tp& IO,vector<T>&v){
for(T&x:v)IO>>x;
return IO;
}
template<class T>
IO_Tp& operator<<(IO_Tp& IO,vector<T>&v){
for(T&x:v)IO<<x<<' ';
return IO;
}
//hodf 缝合版
