点击查看代码
namespace Poly{
int *r[23];
int rr[MAX<<1];
struct polyinit{
polyinit(){
int nw = 1;
int len = 1;
while(nw+len<=(MAX<<1)){
int pos = lg(len);
r[pos] = rr+nw;
rep(i,1,len)r[pos][i] = (r[pos][i>>1]>>1)|((i&1)<<(pos-1));
nw += len;
len <<= 1;
}
}
}polyinit;
template <class T = mint>
struct poly{
vT a;
T G=3,Gi=G.inv();
int len = 0;
poly(){}
explicit poly(int Len){len = Len;a.resize(Len);}
poly(vT v){len=v.size();a = v;}
T & operator [] (int x){return a[x];}
void set(int x){
len = 1;
a.resize(1);
a[0] = x;
return ;
}
const poly<T> operator % (const int Len)const {
poly<T> rlt;
rlt.len = Len;
if(Len>=len){rlt.a=a;rlt.a.resize(Len);}
else rlt.a.assign(a.begin(),a.begin()+Len);
return rlt;
}
poly<T> & operator %= (const int Len)
{a.resize(Len);len = Len;return *this;}
const poly<T> operator << (int x)const {
poly<T> rlt(len+x);
rep(i,x,len+x)rlt[i] = a[i-x];
return rlt;
}
poly<T> & operator <<= (int x){return *this=*this<<x;}
const poly<T> operator >> (int x)const {
poly<T> rlt(len-x);
rep(i,0,x)assert(a[i]==T(0));
rep(i,x,len)rlt[i-x] = a[i];
return rlt;
}
poly<T> & operator >>= (int x){return *this=*this>>x;}
const bool operator < (const poly b)const {assert(false);return 0;}
const bool operator > (const poly b)const {assert(false);return 0;}
const bool empty()const {return a[0].empty();}
const poly<T> operator + (poly<T> b)const {
poly<T> rlt = *this;
rlt %= max(len,b.len);
rep(i,0,b.len)rlt[i] += b[i];
return rlt;
}
poly<T> & operator += (poly<T> b){
*this %= max(len,b.len);
rep(i,0,b.len)a[i] += b[i];
return *this;
}
const poly<T> operator - (poly<T> b)const {
poly<T> rlt = *this;
rlt %= max(len,b.len);
rep(i,0,b.len)rlt[i] -= b[i];
return rlt;
}
poly<T> & operator -= (poly<T> b){
*this %= max(len,b.len);
rep(i,0,b.len)a[i] -= b[i];
return *this;
}
T operator () (T x)const {
T rlt; rlt.set(0);
T w; w.set(0);
rep(i,0,len){rlt += w*a[i];w *= x;}
return rlt;
}
void NTT(int N=1){
while(N<len)N <<= 1;
*this %= N;
int pos = lg(N);
rep(i,0,N)if(r[pos][i]<i)swap(a[r[pos][i]],a[i]);
a.resize(N);
T x,y,w,omg;
for(int mid=1;mid<N;mid<<=1){
omg = pw(G,(mod-1)/(mid<<1));
for(int i=0;i<N;i+=(mid<<1)){
w = 1;
rep(j,i,i+mid){
x = a[j];y = a[j+mid]*w;
a[j] = x+y;a[j+mid] = x-y;
w *= omg;
}
}
}
return ;
}
void INTT(int N=1){
while(N<len)N <<= 1;
*this %= N;
int pos = lg(N);
rep(i,0,N)if(r[pos][i]<i)swap(a[r[pos][i]],a[i]);
T x,y,w,omg;
for(int mid=1;mid<N;mid<<=1){
omg = pw(Gi,(mod-1)/(mid<<1));
for(int i=0;i<N;i+=(mid<<1)){
w = 1;
rep(j,i,i+mid){
x = a[j];y = a[j+mid]*w;
a[j] = x+y;a[j+mid] = x-y;
w *= omg;
}
}
}
return ;
}
const poly<T> operator * (poly<T> b)const {
poly<T> c = *this;
int N = 1;
int Len = len+b.len-1;
while(N<Len)N <<= 1;
c.NTT(N);b.NTT(N);
T Ni = T(N).inv();
rep(i,0,N)c[i] *= b[i]*Ni;
c.INTT(N);
c %= Len;
return c;
}
poly<T> & operator *= (poly<T> b){
int Len = len+b.len-1;
int N = 1;
while(N<Len)N <<= 1;
this->NTT(N);b.NTT(N);
T Ni = T(N).inv();
rep(i,0,N)a[i] *= b[i]*Ni;
this->INTT(N);
*this %= Len;
return *this;
}
const poly<T> operator * (T b)const {
poly rlt(len);
rep(i,0,len)rlt[i] = a[i]*b;
return rlt;
}
poly<T> & operator *= (T b){
rep(i,0,len)a[i] *= b;
return *this;
}
const poly<T> inv_iter(const poly<T> &f)const
{return (*this)*(poly(vT{2})-(*this)*(f%(len*2)))%(len*2);}
const poly<T> inv(int n)const {
int N = 1;
poly<T> rlt(1);
rlt[0] = a[0].inv();
while(N<n){
rlt = rlt.inv_iter(*this);
N <<= 1;
}
rlt %= n;
return rlt;
}
poly<T> & deSelf(){
rep(i,0,len-1)a[i] = a[i+1]*T(i+1);
*this %= len-1;
return *this;
}
const poly<T> de()const {
poly<T> rlt(len-1);
rep(i,0,len-1)rlt[i] = a[i+1]*T(i+1);
return rlt;
}
poly<T> &interSelf(){
*this %= len+1;
per(i,len,1)a[i] = a[i-1]*idig[i];
a[0] = 0;
return *this;
}
const poly<T> inter()const {
poly<T> rlt(len+1);
per(i,len+1,1)rlt[i] = a[i-1]*idig[i];
return rlt;
}
const poly<T> ln(int n)const {
assert(a[0]==T(1));
poly<T> rlt;
rlt = (this->de()*this->inv(n)%n).inter()%n;
return rlt;
}
const poly<T> ln(int n,poly<T> invf)const
{return (this->de()*invf.inv_iter(*this)%n).inter()%n;}
const poly<T> exp_iter(const poly<T> &f,const poly<T> &invg)const
{return ((*this)*((f%(len*2))-(*this).ln(len*2,invg)+poly(vT{1})))%(len*2);}
const poly<T> exp(int n)const {
assert(a[0]==T(0));
poly<T> g(vT{1});
poly<T> invg(vT{1});
int N = 1;
while(N<n){
g = g.exp_iter(*(this),invg);
invg = invg.inv_iter(g);
N <<= 1;
}
g %= n;
return g;
}
pair<poly<T>,poly<T> > exp_with_inv(int n){
// e^f(x) e^-f(x)
assert(a[0]==T(0));
poly<T> g(vT{1});
poly<T> invg(vT{1});
int N = 1;
while(N<n){
g = g.exp_iter(*(this),invg);
invg = invg.inv_iter(g);
N <<= 1;
}
g %= n;invg %= n;
return m_p(g,invg);
}
poly<T> sqrt_iter(const poly<T> &f,const poly<T> &invg)const
{return ((*this)+(f%(len*2)*invg.inv_iter(*this))%(len*2))*(T(2).inv());}
poly<T> sqrt(int n)const {
assert(a[0]==T(1));
poly<T> g(vT{1});
poly<T> invg(vT{1});
int N = 1;
while(N<n){
g = g.sqrt_iter(*this,invg);
invg = invg.inv_iter(g);
N <<= 1;
}
g %= n;
return g;
}
const poly<T> ksm(int n,ll d)const {
assert(a[0]==T(1));
poly<T> f = this->ln(n);
rep(i,0,n)f[i] *= d;
poly<T> rlt = f.exp(n);
return rlt;
}
};
template<class T>
ostream & operator << (ostream & cout,poly<T> f){
cout<<"(";
rep(i,0,f.len)cout<<f[i]<<",)"[i==f.len-1];
return cout;
}
}
using Poly::poly;
