代码

#ifndef SOUPO_POLYNOMIAL
#define SOUPO_POLYNOMIAL
#ifndef SOUPO_MODINT
#define SOUPO_MODINT
#include<chrono>
#include<utility>
#include<cstdint>
#include<type_traits>
#include<iostream>
template<uint64_t P> class ModInt;
template<uint64_t P> constexpr ModInt<P> qpw(ModInt<P>,uint64_t);
template<uint64_t P> constexpr ModInt<P> inv(ModInt<P>);
template<uint64_t P> constexpr ModInt<P> sqrt(ModInt<P>);
template<uint64_t P>
class ModInt{
    static_assert(!(P>>62));
    static_assert(P&1);
private:
    static constexpr uint64_t PP=P<<1;
    uint64_t x;
    static constexpr uint64_t calc_pinv(){
        uint64_t inv=P;
        return inv*=2-P*inv,inv*=2-P*inv,inv*=2-P*inv,inv*=2-P*inv,inv*=2-P*inv,-inv;
    }
    static constexpr uint64_t calc_r2(){
        __uint128_t r=-1;
        return r%P+1;
    }
    static constexpr uint64_t PINV=calc_pinv(),R2=calc_r2();
    static constexpr uint64_t reduce(__uint128_t T){
        uint64_t m=uint64_t(T)*PINV;
        uint64_t t=(T+__uint128_t(m)*P)>>64;
        return t;
    }
    struct RawConstructTag{};
    constexpr ModInt(uint64_t raw_val,RawConstructTag):x(raw_val){}
public:
    constexpr ModInt():x(0){}
    constexpr ModInt(uint64_t _x):x(reduce((__uint128_t)(_x)*R2)){}
    constexpr uint64_t val()const{uint64_t t=reduce(x);return t-(P&-(t>=P));}
    constexpr ModInt operator+(const ModInt&rhs)const{
        uint64_t t=x+rhs.x;
        return ModInt(t-=(PP&-(t>=PP)),RawConstructTag{});
    }
    constexpr ModInt operator-(const ModInt&rhs)const{
        uint64_t t=x-rhs.x;
        return ModInt(t+=(PP&(int64_t(t)>>63)),RawConstructTag{});
    }
    constexpr ModInt operator*(const ModInt&rhs)const{
        return ModInt(reduce(__uint128_t(x)*rhs.x),RawConstructTag{});
    }
    constexpr ModInt operator-()const{
        return ModInt((PP-x)&((int64_t)(-x)>>63),RawConstructTag{});
    }
    constexpr ModInt&operator+=(const ModInt&rhs){
        return x+=rhs.x,x-=(PP&-(x>=PP)),*this;
    }
    constexpr ModInt&operator-=(const ModInt&rhs){
        return x-=rhs.x,x+=(PP&(int64_t(x)>>63)),*this;
    }
    constexpr ModInt&operator*=(const ModInt&rhs){
        return x=reduce(__uint128_t(x)*rhs.x),*this;
    }
    constexpr bool operator==(const ModInt&rhs)const{
        uint64_t a=x-(P&(-(x>=P))),b=rhs.x-(P&(-(rhs.x>=P)));
        return a==b;
    }
    constexpr bool operator!=(const ModInt&rhs)const{
        uint64_t a=x-(P&(-(x>=P))),b=rhs.x-(P&(-(rhs.x>=P)));
        return a!=b;
    }
    friend constexpr ModInt<P> qpw<>(ModInt<P>,uint64_t);
    friend constexpr ModInt<P> inv<>(ModInt<P>);
    friend constexpr ModInt<P> sqrt<>(ModInt<P>);
};
template<uint64_t P>
constexpr ModInt<P> qpw(ModInt<P> x,uint64_t y){
    ModInt<P> res{1};
    for(;y;y>>=1,x*=x)
        if(y&1)
            res*=x;
    return res;
}
template<uint64_t P>
constexpr ModInt<P> inv(ModInt<P> x){
    ModInt<P> res{1};
    for(uint64_t y=P-2;y;y>>=1,x*=x)
        if(y&1)
            res*=x;
    return res;
}
template<uint64_t P>
constexpr ModInt<P> sqrt(ModInt<P>x){
    if(x==ModInt<P>(0))return x;
    ModInt<P> a{},w2;
    for(;qpw(w2=a*a-x,(P-1)>>1)==ModInt<P>(1);a+=ModInt<P>(1));
    if(w2==ModInt<P>(0))return a;
    uint64_t y=(P+1)>>1;
    std::pair<ModInt<P>,ModInt<P>>res{ModInt<P>(1),ModInt<P>()},X{a,ModInt<P>(1)};
    for(;y;y>>=1,X={X.first*X.first+X.second*X.second*w2,X.first*X.second*2})if(y&1)res={res.first*X.first+res.second*X.second*w2,res.first*X.second+res.second*X.first};
    return res.first;
}
template<> class ModInt<2>;
template<> constexpr ModInt<2> qpw(ModInt<2>,uint64_t);
template<> constexpr ModInt<2> inv(ModInt<2>);
template<> constexpr ModInt<2> sqrt(ModInt<2>);
template<>
class ModInt<2>{
private:
    uint8_t x;
public:
    constexpr ModInt<2>():x(0){}
    constexpr ModInt<2>(uint64_t _x):x(_x&1){}
    constexpr uint64_t val(){return x;}
    constexpr ModInt<2> operator+(const ModInt<2>&rhs)const{
        return ModInt<2>(x^rhs.x);
    }
    constexpr ModInt<2> operator-(const ModInt<2>&rhs)const{
        return ModInt<2>(x^rhs.x);
    }
    constexpr ModInt<2> operator*(const ModInt<2>&rhs)const{
        return ModInt<2>(x&rhs.x);
    }
    constexpr ModInt<2> operator-()const{
        return ModInt<2>(x);
    }
    constexpr ModInt<2>&operator+=(const ModInt<2>&rhs){
        return x^=rhs.x,*this;
    }
    constexpr ModInt<2>&operator-=(const ModInt<2>&rhs){
        return x^=rhs.x,*this;
    }
    constexpr ModInt<2>&operator*=(const ModInt<2>&rhs){
        return x&=rhs.x,*this;
    }
    constexpr bool operator==(const ModInt<2>&rhs)const{return x==rhs.x;}
    constexpr bool operator!=(const ModInt<2>&rhs)const{return x!=rhs.x;}
    constexpr friend ModInt<2> qpw<>(ModInt<2>,uint64_t);
    constexpr friend ModInt<2> inv<>(ModInt<2>);
    constexpr friend ModInt<2> sqrt<>(ModInt<2>);
};
template<>
constexpr ModInt<2> qpw(ModInt<2>x,uint64_t y){
    return ModInt<2>(x.x|(!y));
}
template<>
constexpr ModInt<2> inv(ModInt<2>x){
    return x;
}
template<>
constexpr ModInt<2> sqrt(ModInt<2>x){
    return x;
}
#endif
#include<stdint.h>
#include<vector>
#include<algorithm>
constexpr uint64_t P=998244353,primitive_root=3;
using mint=ModInt<P>;
using STDVECTORMINT=std::vector<mint>;
class poly: public STDVECTORMINT{
public:
    using STDVECTORMINT::STDVECTORMINT;
protected:
    static vector<mint>GW[22];
    static void init(int t){
        if(GW[t].size())return;
        auto&gw=GW[t];
        gw.resize((1<<t)+1);
        gw[0]=1,gw[1<<t]=qpw(mint(31),1<<(21-t));
        for(int i=t;i;--i)gw[1<<(i-1)]=gw[1<<i]*gw[1<<i];
        for(int i=1;i<(1<<t);++i)gw[i]=gw[i&(i-1)]*gw[i&-i];
    }
    friend void DIF(poly&a,int t);
    friend void DIT(poly&a,int t);
public:
    friend poly operator*(poly a,poly b);
    friend poly&operator*=(poly&a,poly b);
};
poly operator+(poly a,const poly&b){
    if(b.size()>a.size())a.resize(b.size());
    for(int i=0;i<b.size();++i)a[i]+=b[i];
    return a;
}
poly&operator+=(poly&a,const poly&b){
    if(b.size()>a.size())a.resize(b.size());
    for(int i=0;i<b.size();++i)a[i]+=b[i];
    return a;
}
poly operator-(poly a,const poly&b){
    if(b.size()>a.size())a.resize(b.size());
    for(int i=0;i<b.size();++i)a[i]-=b[i];
    return a;
}
poly&operator-=(poly&a,const poly&b){
    if(b.size()>a.size())a.resize(b.size());
    for(int i=0;i<b.size();++i)a[i]-=b[i];
    return a;
}
poly operator-(poly a){
    for(auto&x:a)x=-x;
    return a;
}
void DIT(poly&a,int t){
    poly::init(t);
    auto&gw=poly::GW[t];
    int lim=1<<t;
    a.resize(lim);
    for(int l=1;l<lim;l<<=1){
        auto k=a.begin();
        for(auto g=gw.begin();k<a.end();k+=(l<<1),++g)
            for(auto x=k;x<k+l;++x){
                auto _=*(x+l);
                (*(x+l))=((*x)-_)*(*g),(*x)+=_;
            }
    }
    uint64_t iv=P-((P-1)>>t);
    for(auto&x:a)x*=iv;
    reverse(a.begin()+1,a.end());
}
void DIF(poly&a,int t){
    poly::init(t);
    auto&gw=poly::GW[t];
    int lim=1<<t;
    a.resize(lim);
    for(int l=(lim>>1);l;l>>=1){
        auto k=a.begin();
        for(auto g=gw.begin();k<a.end();k+=(l<<1),++g)
            for(auto x=k;x<k+l;++x){
                auto _=(*(x+l))*(*g);
                (*(x+l))=(*x)-_,(*x)+=_;
            }
    }
}
poly operator*(poly a,poly b){
    int _=a.size()+b.size()-1;
    int t=0;
    for(;(1<<t)<_;++t);
    --t;
    DIF(a,t),DIF(b,t);
    for(int i=0;i<(1<<t);++i)a[i]*=b[i];
    DIT(a,t);
    a.resize(_);
    return a;
}
poly&operator*=(poly&a,poly b){
    int _=a.size()+b.size()-1;
    int t=0;
    for(;(1<<t)<_;++t);
    --t;
    DIF(a,t),DIF(b,t);
    for(int i=0;i<(1<<t);++i)a[i]*=b[i];
    DIT(a,t);
    a.resize(_);
    return a;
}
poly operator*(poly a,const mint&b){
    for(auto&x:a)x*=b;
    return a;
}
poly&operator*=(poly&a,const mint&b){
    for(auto&x:a)x*=b;
    return a;
}
poly Product(std::vector<poly>&a){
    if(a.size()==1)return a[0];
    auto mid=a.size()/2;
    std::vector<poly>a0(a.begin(),a.begin()+mid),a1(a.begin()+mid,a.end());
    a.clear();
    return Product(a0)*Product(a1);
}
poly Inv(const poly&f,int N){
    poly g({inv(f[0])});
    for(int n=1;n<N;n<<=1){
        poly h(f.begin(),f.begin()+std::min({n<<1,N,(int)f.size()}));
        h*=-g,h[0]+=2,(g*=h).resize(std::min(n<<1,N));
    }return g;
}
poly Div(poly f,poly g){
    if(f.size()<g.size())return poly({mint(0)});
    auto fR=f,gR=g;
    std::reverse(fR.begin(),fR.end()),std::reverse(gR.begin(),gR.end());
    auto q=fR*Inv(gR,f.size()-g.size()+1);
    q.resize(f.size()-g.size()+1);
    std::reverse(q.begin()+1,q.end());
    return q;
}
poly Mod(poly f,poly g){
    if(f.size()<g.size())return f;
    if(g.empty())return (poly){mint(0)};
    auto fR=f,gR=g;
    std::reverse(fR.begin(),fR.end()),std::reverse(gR.begin(),gR.end());
    auto q=fR*Inv(gR,f.size()-g.size()+1);
    q.resize(f.size()-g.size()+1);
    std::reverse(q.begin(),q.end());
    auto r=f-q*g;r.resize(g.size()-1);
    return r;
}
poly Ln(const poly&f,int N){
    auto g=f;
    for(int i=1;i<g.size();++i)g[i-1]=g[i]*i;
    g.pop_back();
    g*=Inv(f,N),g.resize(N);
    for(int i=N;--i;g[i]=g[i-1]*inv(mint(i)));
    g[0]=0;
    return g;
}
poly Exp(const poly&f,int N){
    poly g{1};
    for(int n=1;n<N;n<<=1){
        poly h(f.begin(),f.begin()+std::min({n<<1,N,(int)f.size()}));
        h-=Ln(g,std::min(n<<1,N)),h[0]+=1;
        (g*=h).resize(std::min(n<<1,N));
    }
    return g;
}
poly Pow(const poly&f,long long k,int N){
    if(k==0){
        poly g(N);g[0]=1;
        return g;
    }
    int cnt=0;
    for(;cnt<f.size()&&f[cnt]==0;++cnt);
    if(cnt==f.size()||(cnt&&(k>=N||k*cnt>=N)))return poly(N);
    N-=(cnt*k);
    poly g(f.begin()+cnt,f.begin()+std::min(cnt+N,(int)f.size()));
    cnt*=k;
    mint c=g[0],ic=inv(c);
    for(auto&x:g)x*=ic;
    mint _k=k;
    g=Ln(g,N);
    for(auto&x:g)x*=_k;
    g=Exp(g,N);
    int __k=k%(P-1);if(__k<0)__k+=P-1;
    c=qpw(c,__k);
    for(auto&x:g)x*=c;
    g.insert(g.begin(),cnt,mint{});
    return g;
}
#endif
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c){
    poly A(n+1),B(m+1);
    for(int i=0;i<=n;++i)A[i]=a[i];
    for(int i=0;i<=m;++i)B[i]=b[i];
    poly C=A*B;
    for(int i=0;i<=n+m;++i)c[i]=C[i].val();
}

评测结果