#include using namespace std; struct Montgomery{ //2^62未満&奇数modのみ. //初めにsetmodする. using u64 = uint64_t; using u128 = __uint128_t; private: static u64 mod,N2,Rsq; //N*N2≡1(mod N); //Rsq = R^2modN; R=2^64. u64 v = 0; public: long long val(){return reduce(v);} static u64 getmod(){return mod;} static void setmod(u64 m){ assert(m<(1LL<<62)&&(m&1)); mod = m; N2 = mod; for(int i=0; i<5; i++) N2 *= 2-N2*mod; Rsq = (-u128(mod))%mod; } //reduce = T*R^-1modNを求める. u64 reduce(const u128 &T){ //T*R^-1≡(T+(T*(-N2))modR*N)/R 2N未満なので-N必要かだけで良い. u64 ret = (T+u128(((u64)T)*(-N2))*mod)>>64; if(ret >= mod) ret -= mod; return ret; } //初期値reduce(R^2)でok. Montgomery(){v = 0;} Montgomery(long long w):v(reduce(u128(w)*Rsq)){} Montgomery& operator=(const Montgomery &b) = default; Montgomery operator-()const{return Montgomery()-Montgomery(*this);} Montgomery operator+(const Montgomery &b)const{return Montgomery(*this)+=b;} Montgomery operator-(const Montgomery &b)const{return Montgomery(*this)-=b;} Montgomery operator*(const Montgomery &b)const{return Montgomery(*this)*=b;} Montgomery operator/(const Montgomery &b)const{return Montgomery(*this)/=b;} Montgomery& operator+=(const Montgomery &b){ v += b.v; if(v >= mod) v -= mod; return (*this); } Montgomery& operator-=(const Montgomery &b){ v += mod-b.v; if(v >= mod) v -= mod; return (*this); } Montgomery& operator*=(const Montgomery &b){ v = reduce(u128(v)*b.v); return (*this); } Montgomery& operator/=(const Montgomery &b){ (*this) *= b.inv(); return (*this); } Montgomery pow(u64 b)const{ Montgomery ret = 1,p = (*this); while(b){ if(b&1) ret *= p; p *= p; b >>= 1; } return ret; } Montgomery inv()const{return pow(mod-2);} bool operator!=(const Montgomery &b)const{return v!=b.v;} bool operator==(const Montgomery &b)const{return v==b.v;} }; typename Montgomery::u64 Montgomery::mod,Montgomery::N2,Montgomery::Rsq; using mint = Montgomery; namespace to_fold{ //纏めて折り畳むためのやつ 苦労したが結局ほぼACLの写経. __int128_t safemod(__int128_t a,long long m){a %= m; if(a < 0) a += m; return a;} pair invgcd(long long a,long long b){ //return {gcd(a,b),x} (xa≡g(mod b)) a = safemod(a,b); if(a == 0) return {b,0}; long long x = 0,y = 1,memob = b; while(a){ long long q = b/a; b -= a*q; swap(x,y); y -= q*x; swap(a,b); } if(x < 0) x += memob/b; return {b,x}; } long long Garner(vector &A,vector &M){ __int128_t mulM = 1,x = A.at(0)%M.at(0); //Mの要素のペア互いに素必須. for(int i=1; i root,iroot,rate2,irate2,rate3,irate3; fftinfo(mint g){ checkmod = mint::getmod(); int rank2 = countzero(mint::getmod()-1); root.resize(rank2+1),iroot = root; rate2.resize(max(0,rank2-2+1)); irate2 = rate2; rate3.resize(max(0,rank2-3+1)); irate3 = rate3; root[rank2] = g.pow((mint::getmod()-1)>>rank2); iroot[rank2] = root[rank2].inv(); for(int i=rank2-1; i>=0; i--){ root[i] = root[i+1]*root[i+1]; iroot[i] = iroot[i+1]*iroot[i+1]; } mint mul = 1,imul = 1; for(int i=0; i<=rank2-2; i++){ rate2[i] = root[i+2]*mul; irate2[i] = iroot[i+2]*imul; mul *= iroot[i+2],imul *= root[i+2]; } mul = 1,imul = 1; for(int i=0; i<=rank2-3; i++){ rate3[i] = root[i+3]*mul; irate3[i] = iroot[i+3]*imul; mul *= iroot[i+3],imul *= root[i+3]; } } }; mint findroot(){ if(mint::getmod() == 998244353) return mint(3); else if(mint::getmod() == 754974721) return mint(11); else if(mint::getmod() == 167772161) return mint(3); else if(mint::getmod() == 469762049) return mint(3); assert(false); //バグあったので. } fftinfo info(3); void NTT(vector &A){ long long Mod = mint::getmod(); if(info.getmod() != Mod) info = fftinfo(findroot()); int N = A.size(),ln = countzero(N); int dep = 0; while(dep < ln){ if(ln-dep == 1){ int p = 1<<(ln-dep-1); mint rot = 1; for(int o=0; o<(1< &A){ long long Mod = mint::getmod(); if(info.getmod() != Mod) info = fftinfo(findroot()); int N = A.size(),ln = countzero(N),dep = ln; while(dep){ if(dep == 1){ int p = 1<<(ln-dep); mint irot = 1; for(int o=0; o<(1<<(dep-1)); o++){ int offset = o<<(ln-dep+1); for(int i=0; i convolution(vector &A,vector &B){ int siza = A.size(),sizb = B.size(),sizc = siza+sizb-1,N = 1; if(siza == 0 || sizb == 0) return {}; if(min(siza,sizb) <= 50){ //naive. vector ret(sizc); for(int i=0; i ca = A,cb = B; ca.resize(N),cb.resize(N); assert((mint::getmod()-1)%N == 0); mint root = findroot(); NTT(ca); NTT(cb); for(int i=0; i convolution_ll(vector &A,vector &B){ //long longに収まる範囲. int siza = A.size(),sizb = B.size(),sizc = siza+sizb-1; if(siza == 0 || sizb == 0) return {}; vector ret(sizc); if(min(siza,sizb) <= 200){ //naive 200はやばい?. for(int i=0; i a(siza),b(sizb); for(int i=0; i C1 = convolution(a,b); mint::setmod(mod2); for(int i=0; i C2 = convolution(a,b); mint::setmod(mod3); for(int i=0; i C3 = convolution(a,b); vector offset = {0,0,m1m2m3,2*m1m2m3,3*m1m2m3}; for(int i=0; i convolution_llmod(vector &A,vector &B,long long memomod){ int siza = A.size(),sizb = B.size(),sizc = siza+sizb-1; if(siza == 0 || sizb == 0) return {}; vector ret(sizc); if(min(siza,sizb) <= 200){ for(int i=0; i a(siza),b(sizb); for(int i=0; i C1 = convolution(a,b); mint::setmod(mod2); for(int i=0; i C2 = convolution(a,b); mint::setmod(mod3); for(int i=0; i C3 = convolution(a,b); mint::setmod(memomod); for(int i=0; i A = {C1.at(i).val(),C2.at(i).val(),C3.at(i).val()}; vector M = {mod1,mod2,mod3}; ret.at(i) = Garner(A,M); } return ret; } vector convolution_int(vector &A,vector &B){ //intに収まる範囲. if(A.size() == 0 || B.size() == 0) return {}; vector ret; if(min(A.size(),B.size()) <= 50){ ret.resize(A.size()+B.size()-1); for(int i=0; i X,Y,Z; for(auto &a : A) X.push_back(a); for(auto &b : B) Y.push_back(b); Z = convolution(X,Y); for(auto &z : Z) ret.push_back(z.val()); } return ret; } } using namespace to_fold; template //実質mintだけ?. struct FormalPowerSeries:vector{ //NTT-friendly素数だけ じゃなくてもいいけど全部書き直せ!. using vector::vector; using fps = FormalPowerSeries; //重要なところは某のほぼパクリ. fps operator++(){*this += 1; return *this;} fps operator--(){*this -= 1; return *this;} fps operator++(int){*this += 1; return *this;} fps operator--(int){*this -= 1; return *this;} fps operator+(const fps &b) const {return fps(*this)+=b;} fps operator+(const T &b) const {return fps(*this)+=b;} fps operator-(const fps &b) const {return fps(*this)-=b;} fps operator-(const T &b) const {return fps(*this)-=b;} fps operator*(const fps &b){return fps(*this)*=b;} fps operator*(const T &b) const {return fps(*this)*=b;} fps operator/(const fps &b) const {return fps(*this)/=b;} fps operator%(const fps &b) const {return fps(*this)%=b;} fps operator>>(const unsigned int b) const {return fps(*this)>>=b;} fps operator<<(const unsigned int b) const {return fps(*this)<<=b;} fps operator-() const { fps ret = (*this); for(auto &v : ret) v = -v; return ret; } bool operator==(const fps &b)const{ if((*this).size() != b.size()) return false; for(int i=0; i<(*this).size(); i++) if((*this).at(i) != b.at(i)) return false; return true; } bool operator!=(const fps &b)const{return !((*this)==b);} fps &operator+=(const fps &b){ //Cix^i = (Ai+Bi)x^i. O(n). if((*this).size() < b.size()) (*this).resize(b.size(),0); for(int i=0; i f,g; for(auto v : (*this)) f.emplace_back(v); for(auto v : b) g.emplace_back(v); vector ret; if(mint::getmod() == 998244353) ret = convolution(f,g); //型を仕方なく合わせる C++わからん. else ret = convolution_llmod(f,g,mint::getmod()); //modがNTT-friendlyじゃない時用 直すのだるい 998以外はとりあえずこっち. (*this).resize(ret.size()); for(int i=0; i>=(const unsigned int &b){//b<0は対象外. 先頭b項を削除. O(n) if((*this).size() <= b) (*this).clear(); else (*this).erase((*this).begin(),(*this).begin()+b); return *this; } fps &operator<<=(const unsigned int &b){//b<0は対象外. 先頭b項に0を挿入. O(n) (*this).insert((*this).begin(),b,0); return *this; } fps &operator%=(const fps &b){ //多項式の余り. O(nlogn) (*this) -= (*this)/b*b; del0(); return (*this); } fps &operator/=(const fps &b){ //多項式としての除算 O(nlogn). assert(b.size() > 0); //分母の末尾0は駄目. T check = b.back(); assert(check != 0); long long mod = mint::getmod(); del0(); //分子の末尾0は消して許容. if((*this).size() < b.size()){ (*this).clear(); return *this; } int n = (*this).size()-b.size()+1; if(b.size() <= 64){ //愚直. fps G(b); assert(G.size() > 0); T div = G.back().inv(); for (auto &v : G) v *= div; int deg = (*this).size()-G.size()+1; fps Q(deg); for(int i=deg-1; i>=0; i--) { Q[i] = (*this).at(i+G.size()-1); for(int k=0; k> P; for(int i=0; i> P; if((*this).at(0) != 1) div = T((*this).at(0)).inv(); for(int i=1; i= k) ret.at(i) -= v*ret.at(i-k); //-xが+ret[i-1]に対応. if(div != 1) for(auto &v : ret) v *= div; return ret; } fps inv_sparse(const fps &b,int deg = -1)const{ //f/gを返す 1/fでは分母だがこれは分子に注意. int n = (*this).size(),m = b.size(); if(deg == -1) deg = n; assert(b.at(0) != 0); T div = 1; vector> P; if(b.at(0) != 1) div = T(b.at(0)).inv(); for(int i=1; i= k) ret.at(i) -= v*ret.at(i-k); if(div != 1) for(auto &v : ret) v *= div; return ret; } fps log_sparse(int deg = -1){ //log(f)を返す O(N*非0). //logf = ∫(f'/f) inv,1/f*(f')がO(N*非0) 他はO(N). assert((*this).size()&&(*this).at(0)==1); if(deg == -1) deg = (*this).size(); fps ret = (*this).diff(); ret = ((*this).inv_sparse(deg)).multi_sparse(ret); return ret.inte().prefix(deg); } fps exp_sparse(int deg = -1)const{ //exp(f)を返す O(N*非0). //(expf)'=(f')*exp(f)より低次から決まる. //[x^0]expf=1より左辺のx^0の係数が求まる->expfのx^1の係数が求まる->... if(deg == -1) deg = (*this).size(); fps ret(deg); if((*this).size() == 0){ if(deg > 0) ret.at(0) = 1; return ret; } assert((*this).at(0) == 0); long long mod = mint::getmod(); if(deg == 1) return fps{1}; ret.at(0) = 1; ret.at(1) = 1; for(int i=2; i> P; for(int i=1; i<(*this).size(); i++) if((*this).at(i) != 0) P.emplace_back(pair{i-1,(*this).at(i)*i}); for(int i=0; i 0) ret.at(0) = 1; return ret; } for(int t=0; t> P; for(int i=t+1; i 1) ret.at(1) = 1; for(int i=2; i n+1) break; if(i > 0 && i<=n+1) now += mulK*ret.at(n+1-i)*i*v; if(i > 0 && i <= n) now -= v*(n+1-i)*ret.at(n+1-i); } ret.at(n+1) *= now; } ret *= T((*this).at(t)).pow(K); return (ret<<(t*K)).prefix(deg); } if(K >= deg || (t+1)*K >= deg) break; } return fps(deg,0); } fps diff()const{ //微分 nx^(n-1) = (x^n)' O(n). int n = (*this).size(); if(n <= 1) return {}; fps ret(n-1); T multi = 1; for(int i=1; i divi; //invと衝突回避用 mintでiの逆元. divi.resize(deg*2); divi.at(1) = 1; for(int i=2; i void { //inplaceで積分. int n = f.size(); f.insert(f.begin(),0); for(int i=1; i<=n; i++) f.at(i) *= divi.at(i); }; auto differential = [&](fps &f) -> void { //inplaceで微分. if(f.size() == 0) return; f.erase(f.begin()); T multi = 0; for(int i=0; i 1) f.push_back((*this).at(1)); else f.push_back(0); for(int m=2; m 0) ret.at(0) = 1; return ret; } for(int i=0; i>i).log(deg)*(K%mod)).exp(deg); ret *= T((*this).at(i)).pow(K); //[x^i]f^Kの分. ret = (ret<<(i*K)).prefix(deg); //*x^(i*k)の分. return ret; } else{ T div = T(1)/(*this).at(i); //*([x^i]f)^Kと *x^(i*k)の分は後回し. fps ret = ((((*this)*div)>>i).log(deg,false)*(K%mod)).exp2(deg); ret *= T((*this).at(i)).pow(K); //[x^i]f^Kの分. ret = (ret<<(i*K)).prefix(deg); //*x^(i*k)の分. return ret; } } if(K >= deg || (i+1)*K >= deg) break; //((i+1)*K)乗未満は0確定 int128回避用にK>=deg(degがllはやばい). } return fps(deg,0); //fの係数全て0なら係数全て0. } fps prefix(int siz)const{ //先頭siz項を返す なかったら0埋め. O(siz). fps ret((*this).begin(),(*this).begin()+min((int)(*this).size(),siz)); if(ret.size() < siz) ret.resize(siz,0); return ret; } void del0(){ //末尾の0を消す O(n). while((*this).size() && (*this).back() == 0) (*this).pop_back(); } fps rev()const{ //ひっくり返す O(n). fps ret(*this); reverse(ret.begin(),ret.end()); return ret; } pair getQR(const fps &b)const{ //多項式の商と余りを同時に得る O(nlogn). fps Q = (*this)/b,R = (*this)-Q*b; R.del0(); return {Q,R}; } fps cumulativeNtimes(int N,T b,int deg=-1){ //1/(1-bx)^Nをdeg次まで返す 指定なしはN次まで. //負の二項定理を使う. 1/(1-x)^N=Σ[i=0~∞]((n+i-1) choose i)(bx^i); //fps{}.cumulativeNtime()で無理やり関数を呼び出す. assert(N <= 0); //N=0も駄目? {1}を返すべき所{0}になる. if(deg == -1) deg = N+1; long long mod = mint::getmod(); int Limit = N+deg; //Limit 必要なサイズ fac->x! facinv->1/x! inv->1/x. vector FAC(Limit+1,1); for(int i=1; i<=Limit; i++) FAC.at(i) = FAC.at(i-1)*i; vector FACinv(Limit+1); FACinv.at(min((int)mod-1,Limit)) = FAC.at(min((int)mod-1,Limit)).inv(); for(int i=min((int)mod-2,Limit-1); i>=0; i--) FACinv.at(i) = FACinv.at(i+1)*(i+1); auto nCr = [&](int n, int r) -> T { if(n < r || r < 0 || n < 0) return 0; return FAC.at(n)*FACinv.at(r)*FACinv.at(n-r); }; fps ret(deg); T value = 1; for(int i=0; ix! facinv->1/x! inv->1/x. vector fac(deg+1,1); for(int i=1; i<=deg; i++) fac.at(i) = fac.at(i-1)*i; vector facinv(deg+1); facinv.at(min((int)mod-1,deg)) = fac.at(min((int)mod-1,deg)).inv(); for(int i=min((int)mod-2,deg-1); i>=0; i--) facinv.at(i) = facinv.at(i+1)*(i+1); fps X(deg),Y(deg); for(int i=0; i; //mod998244353以外のNTT-friendlyの時は色々気を付ける 後で書き直すかも?. pair SumfpsFraction(vector A,vector B){ //次数の和=n O(nlog^2n). //ΣAi/Biを求める 分割統治でO(nlog^2n). //a/b+c/d = (ad+bc)/bd 面倒なのでvectorは用意せずfpsだけ; assert(A.size() == B.size() && A.size()); int n = A.size(); auto comp = [&](int a,int b) -> bool {return A[a].size()+B[a].size()>A[b].size()+B[b].size();}; priority_queue,decltype(comp)> Q(comp); for(int i=0; i 1){ int p = Q.top(); Q.pop(); int q = Q.top(); Q.pop(); fps num = A.at(p)*B.at(q)+B.at(p)*A.at(q); A.at(p) = num; B.at(p) *= B.at(q); Q.push(p); } int leader = Q.top(); return {A.at(leader),B.at(leader)}; } pair,vector> SummintFraction(vector> A,vector> B){ //次数の和=n O(nlog^2n). //ΣAi/Biを求める 分割統治でO(nlog^2n). //a/b+c/d = (ad+bc)/bd 面倒なのでvectorは用意せずfpsだけ; assert(A.size() == B.size() && A.size()); int n = A.size(); auto comp = [&](int a,int b) -> bool {return A[a].size()+B[a].size()>A[b].size()+B[b].size();}; priority_queue,decltype(comp)> Q(comp); for(int i=0; i 1){ int p = Q.top(); Q.pop(); int q = Q.top(); Q.pop(); vector num = convolution(A.at(p),B.at(q)); vector num2 = convolution(B.at(p),A.at(q)); if(num.size() < num2.size()) num.resize(num2.size()); for(int i=0; i> N >> S; vector> A(N),B(N); mint divS = mint(1)/S; for(int i=0; i> p; mint one = divS*p,two = one*one; A.at(i) = {two}; B.at(i) = {1,one}; } auto [n,d] = SummintFraction(A,B); mint answer = 0,fac = 1; for(int i=0; i