#pragma GCC target("avx2") #pragma GCC optimize("O3") #pragma GCC optimize("unroll-loops") #include #include using namespace std; #if __has_include() #include using namespace atcoder; #endif using ll = long long; using ld = long double; using ull = long long; #define REP3(i, m, n) for (int i = (m); (i) < int(n); ++ (i)) #define rep(i,a,b) for(int i=(int)(a);i<(int)(b);i++) #define ALL(x) begin(x), end(x) #define all(s) (s).begin(),(s).end() //#define rep2(i, m, n) for (int i = (m); i < (n); ++i) //#define rep(i, n) rep2(i, 0, n) #define drep2(i, m, n) for (int i = (m)-1; i >= (n); --i) #define drep(i, n) drep2(i, n, 0) #define rever(vec) reverse(vec.begin(), vec.end()) #define sor(vec) sort(vec.begin(), vec.end()) #define fi first #define se second #define pb push_back #define pll pair #define REP(i, n) for (int i = 0; i < (n); ++i) #define in scanner.read_int() const ll mod = 998244353; //const ll mod = 1000000007; const ll inf = 2000000000000000000ll; static const long double pi = 3.141592653589793; templatevoid vcin(vector &n){for(int i=0;i>n[i];} templatevoid vcout(vector &n){for(int i=0;ibool chmax(T &a, const T &b) { if (abool chmin(T &a, const T &b) { if (b void ifmin(T t,T u){if(t>u){cout<<-1< void ifmax(T t,T u){if(t>u){cout<<-1<auto make_vector(T x,int arg,Args ...args){if constexpr(sizeof...(args)==0)return vector(arg,x);else return vector(arg,make_vector(x,args...));} ll modPow(ll a, ll n, ll mod) { ll ret = 1; ll p = a % mod; while (n) { if (n & 1) ret = ret * p % mod; p = p * p % mod; n >>= 1; } return ret; } const ll t0=1,t1=10,t2=100,t3=1000,t4=10000,t5=100000,t6=1000000,t7=10000000,t8=t7*10,t9=t8*10,t10=t9*10,t11=t10*10,t12=t11*10,t13=t12*10,t14=t13*10,t15=t14*10,t16=t15*10,t17=t16*10,t18=t17*10; class Scanner { vector buffer; ssize_t n_written; ssize_t n_read; public: Scanner(): buffer(1024*1024) { do_read(); } int64_t read_int() { int64_t ret = 0, sgn = 1; int ch = current_char(); while (isspace(ch)) { ch = next_char(); } if (ch == '-') { sgn = -1; ch = next_char(); } for (; isdigit(ch); ch = next_char()) ret = (ret * 10) + (ch - '0'); return sgn * ret; } private: void do_read() { ssize_t r = read(0, &buffer[0], buffer.size()); if (r < 0) { throw runtime_error(strerror(errno)); } n_written = r; n_read = 0; } inline int next_char() { ++n_read; if (n_read == n_written) { do_read(); } return current_char(); } inline int current_char() { return (n_read == n_written) ? EOF : buffer[n_read]; } }; //Scanner scanner; //void vin(vector &n){for(int i=0;i struct FormalPowerSeries : std::vector { using std::vector::vector; using std::vector::size; using std::vector::resize; using F = FormalPowerSeries; F &operator+=(const F &g){ for(int i=0;i>=(const int d) { int n=(*this).size(); (*this).erase((*this).begin(),(*this).begin()+min(n, d)); (*this).resize(n); return *this; } F &operator=(const std::vector &v) { int n = (*this).size(); for(int i = 0; i < n; ++i) (*this)[i] = v[i]; return *this; } F operator-() const { F ret = *this; return ret * -1; } F &operator*=(const F &g) { if(mode==FAST) { int n=(*this).size(); auto tmp=atcoder::convolution(*this,g); for(int i=0;i= 0; --i) { (*this)[i] *= g[0]; for(int j = 1; j < std::min(i + 1, m); j++) (*this)[i] += (*this)[i - j] * g[j]; } return *this; } } F &operator/=(const F &g) { if(mode == FAST){ int n = (*this).size(); (*this) = atcoder::convolution(*this, g.inv()); return *this; } else{ assert(g[0] != T(0)); T ig0 = g[0].inv(); int n = (*this).size(), m = g.size(); for(int i = 0; i < n; ++i) { for(int j = 1; j < std::min(i + 1, m); ++j) (*this)[i] -= (*this)[i - j] * g[j]; (*this)[i] *= ig0; } return *this; } } F &operator%=(const F &g) { return *this-=*this/g*g; } F operator*(const T &g) const { return F(*this)*=g;} F operator-(const T &g) const { return F(*this)-=g;} F operator*(const F &g) const { return F(*this)*=g;} F operator-(const F &g) const { return F(*this)-=g;} F operator+(const F &g) const { return F(*this)+=g;} F operator/(const F &g) const { return F(*this)/=g;} F operator%(const F &g) const { return F(*this)%=g;} F operator<<(const int d) const { return F(*this)<<=d;} F operator>>(const int d) const { return F(*this)>>=d;} void onemul(const int d,const T c){ int n=(*this).size(); for(int i=n-d-1;i>=0;i--){ (*this)[i+d]+=(*this)[i]*c; } } void onediv(const int d,const T c){ int n=(*this).size(); for(int i=0;i 0); F res{(*this)[0].inv()}; while(int(res.size()) < deg) { int m = res.size(); F f((*this).begin(), (*this).begin() + std::min(n, m * 2)), r(res); f.resize(m * 2), atcoder::internal::butterfly(f); r.resize(m * 2), atcoder::internal::butterfly(r); for(int i = 0; i < m * 2; ++i) f[i] *= r[i]; atcoder::internal::butterfly_inv(f); f.erase(f.begin(), f.begin() + m); f.resize(m * 2), atcoder::internal::butterfly(f); for(int i = 0; i < m * 2; ++i) f[i] *= r[i]; atcoder::internal::butterfly_inv(f); T iz = T(m * 2).inv(); iz *= -iz; for(int i = 0; i < m; ++i) f[i] *= iz; res.insert(res.end(), f.begin(), f.begin() + m); } res.resize(deg); return res; } F &diff_inplace() { int n = (*this).size(); for(int i = 1; i < n; ++i) (*this)[i - 1] = (*this)[i] * i; (*this)[n - 1] = 0; return *this; } F diff() const { F(*this).diff_inplace();} F &integral_inplace() { int n = (*this).size(), mod = T::mod(); std::vector inv(n); { inv[1] = 1; for(int i = 2; i < n; ++i) inv[i] = T(mod) - inv[mod % i] * (mod / i); } for(int i = n - 2; i >= 0; --i) (*this)[i + 1] = (*this)[i] * inv[i + 1]; (*this)[0] = 0; return *this; } F integral() const { return F(*this).integral_inplace(); } F &log_inplace() { int n = (*this).size(); assert(n and (*this)[0] == 1); F f_inv = (*this).inv(); (*this).diff_inplace(); (*this) *= f_inv; (*this).integral_inplace(); return *this; } F log() const { return F(*this).log_inplace(); } F &deriv_inplace() { int n = (*this).size(); assert(n); for(int i = 2; i < n; ++i) (*this)[i] *= i; (*this).erase((*this).begin()); (*this).push_back(0); return *this; } F deriv() const { return F(*this).deriv_inplace(); } F &exp_inplace() { int n = (*this).size(); assert(n and (*this)[0] == 0); F g{1}; (*this)[0] = 1; F h_drv((*this).deriv()); for(int m = 1; m < n; m *= 2) { F f((*this).begin(), (*this).begin() + m); f.resize(2 * m), atcoder::internal::butterfly(f); auto mult_f = [&](F &p) { p.resize(2 * m); atcoder::internal::butterfly(p); for(int i = 0; i < 2 * m; ++i) p[i] *= f[i]; atcoder::internal::butterfly_inv(p); p /= 2 * m; }; if(m > 1) { F g_(g); g_.resize(2 * m), atcoder::internal::butterfly(g_); for(int i = 0; i < 2 * m; ++i) g_[i] *= g_[i] * f[i]; atcoder::internal::butterfly_inv(g_); T iz = T(-2 * m).inv(); g_ *= iz; g.insert(g.end(), g_.begin() + m / 2, g_.begin() + m); } F t((*this).begin(), (*this).begin() + m); t.deriv_inplace(); { F r{h_drv.begin(), h_drv.begin() + m - 1}; mult_f(r); for(int i = 0; i < m; ++i) t[i] -= r[i] + r[m + i]; } t.insert(t.begin(), t.back()); t.pop_back(); t *= g; F v((*this).begin() + m, (*this).begin() + std::min(n, 2 * m)); v.resize(m); t.insert(t.begin(), m - 1, 0); t.push_back(0); t.integral_inplace(); for(int i = 0; i < m; ++i) v[i] -= t[m + i]; mult_f(v); for(int i = 0; i < std::min(n - m, m); ++i) (*this)[m + i] = v[i]; } return *this; } F exp() const { return F(*this).exp_inplace(); } F &pow_inplace(long long k) { int n = (*this).size(), l = 0; assert(k >= 0); if(!k){ for(int i = 0; i < n; ++i) (*this)[i] = !i; return *this; } while(l < n and (*this)[l] == 0) ++l; if(l > (n - 1) / k or l == n) return *this = F(n); T c = (*this)[l]; (*this).erase((*this).begin(), (*this).begin() + l); (*this) /= c; (*this).log_inplace(); (*this).resize(n - l * k); (*this) *= k; (*this).exp_inplace(); (*this) *= c.pow(k); (*this).insert((*this).begin(), l * k, 0); return *this; } F pow(const long long k) const { return F(*this).pow_inplace(); } void spacemul(vector> g) { int n = (*this).size(); auto [d, c] = g.front(); if (d == 0) g.erase(g.begin()); else c = 0; for(int i=n-1;i>=0;i--){ (*this)[i] *= c; for (auto &[j, b] : g) { if (j > i) break; (*this)[i] += (*this)[i-j] * b; } } } void spacediv(vector> g) { int n = (*this).size(); auto [d, c] = g.front(); assert(d == 0 && c != T(0)); T ic = c.inv(); g.erase(g.begin()); for(int i=0;i i) break; (*this)[i] -= (*this)[i-j] * b; } (*this)[i] *= ic; } } }; using fps = FormalPowerSeries; template void GaussJordan(vector> &A,bool is_extended = false){ ll m=A.size(),n=A[0].size(); ll rank=0; for(int i=0;i void linear_equation(vector> a, vector b, vector &res) { ll m=a.size(),n=a[0].size(); vector> M(m,vector(n+1)); for(int i=0;i pair Characteristic_equation(const F &a) { using T=typename F::value_type; ll n=a.size(); ll p=n/2; ll u=p+(p+1); vector> f(u,vector(u)); f[0][0]=1; for(int i=1;i<=p;i++){ f[i][i-1]=-1; } for(int i=p;i b(u); b[0]=1; vector res(u); linear_equation(f,b,res); F X(p),Y(p+1); for(int i=0;i T getK(FormalPowerSeries p, FormalPowerSeries q,ll k){ if(k<0) return T(0); ll d=q.size(); while(k){ auto qn=q; for(int i=1;i>n; cout<