#include using namespace std; using ll = long long; using ld = long double; using ull = unsigned long long; using pll = pair; using tlll = tuple; constexpr ll INF = 1LL << 60; template bool chmin(T& a, T b) {if (a > b) {a = b; return true;} return false;} template bool chmax(T& a, T b) {if (a < b) {a = b; return true;} return false;} ll safemod(ll A, ll M) {ll res = A % M; if (res < 0) res += M; return res;} ll divfloor(ll A, ll B) {if (B < 0) A = -A, B = -B; return (A - safemod(A, B)) / B;} ll divceil(ll A, ll B) {if (B < 0) A = -A, B = -B; return divfloor(A + B - 1, B);} ll pow_ll(ll A, ll B) {if (A == 0 || A == 1) {return A;} if (A == -1) {return B & 1 ? -1 : 1;} ll res = 1; for (int i = 0; i < B; i++) {res *= A;} return res;} ll mul_limited(ll A, ll B, ll M = INF) { return B == 0 ? 0 : A > M / B ? M : A * B; } ll pow_limited(ll A, ll B, ll M = INF) { if (A == 0 || A == 1) {return A;} ll res = 1; for (int i = 0; i < B; i++) {if (res > M / A) return M; res *= A;} return res;} template void unique(vector &V) {V.erase(unique(V.begin(), V.end()), V.end());} template void sortunique(vector &V) {sort(V.begin(), V.end()); V.erase(unique(V.begin(), V.end()), V.end());} #define FINALANS(A) do {cout << (A) << '\n'; exit(0);} while (false) template void printvec(const vector &V) {int _n = V.size(); for (int i = 0; i < _n; i++) cout << V[i] << (i == _n - 1 ? "" : " ");cout << '\n';} template void printvect(const vector &V) {for (auto v : V) cout << v << '\n';} template void printvec2(const vector> &V) {for (auto &v : V) printvec(v);} //* #include #include #include #include using namespace atcoder; //*/ // http://drken1215.hatenablog.com/entry/2018/06/08/210000 template class binom { public: vector fac, finv, inv; binom(int M) { fac.resize(M + 1); finv.resize(M + 1); inv.resize(M + 1); //* fac[0] = T(1); for (int i = 1; i <= M; i++) fac[i] = fac[i - 1] * T::raw(i); finv[M] = fac[M].inv(); for (int i = M - 1; i >= 0; i--) finv[i] = finv[i + 1] * T::raw(i + 1); for (int i = 1; i <= M; i++) inv[i] = fac[i - 1] * finv[i]; //*/ /* fac[0] = T(1), finv[0] = T(1); fac[1] = T(1), finv[1] = T(1), inv[1] = T(1); for (int i = 2; i <= M; i++) { fac[i] = fac[i - 1] * i; inv[i] = -inv[T::mod() % i] * (T::mod() / i); finv[i] = finv[i - 1] * inv[i]; } //*/ } T P(int N, int K) { if (N < K) return 0; if (N < 0 || K < 0) return 0; return fac[N] * finv[N - K]; } T C(int N, int K) { if (N < K) return 0; if (N < 0 || K < 0) return 0; return fac[N] * finv[K] * finv[N - K]; } T H(int N, int K) { if (N == 0 && K == 0) return 1; return C(N + K - 1, K); } }; // http://drken1215.hatenablog.com/entry/2018/06/08/210000 template class binom_mut { private: vector fac, finv, inv; void calc(int n) { int i = fac.size(); if (n < i) return; fac.resize(n + 1), finv.resize(n + 1), inv.resize(n + 1); for (; i <= n; i++) { fac[i] = fac[i - 1] * i; inv[i] = -inv[T::mod() % i] * (T::mod() / i); finv[i] = finv[i - 1] * inv[i]; } } public: binom_mut() { fac = {1, 1}, finv = {1, 1}, inv = {0, 1}; } T get_fac(int n) { assert(n >= 0); calc(n); return fac[n]; } T get_finv(int n) { assert(n >= 0); calc(n); return finv[n]; } T get_inv(int n) { assert(n > 0); calc(n); return inv[n]; } T P(int N, int K) { if (N < K) return 0; if (N < 0 || K < 0) return 0; calc(N); return fac[N] * finv[N - K]; } T C(int N, int K) { if (N < K) return 0; if (N < 0 || K < 0) return 0; calc(N); return fac[N] * finv[K] * finv[N - K]; } T H(int N, int K) { if (N == 0 && K == 0) return 1; return C(N + K - 1, K); } }; // https://qiita.com/taiyaki8926/items/f62f534d43ff006129f7 ll sqrt_mod(ll n, int p) // p は素数 { n %= p; if (n == 0) return 0; if (p == 2) return n; if (pow_mod(n, (p - 1) / 2, p) == p - 1) // 平方非剰余 return -1; if (p % 4 == 3) return pow_mod(n, (p + 1) / 4, p); internal::barrett ba(p); int q = p - 1, s = 0; while (q % 2 == 0) q /= 2, s++; int z = 2; while (pow_mod(z, (p - 1) / 2, p) != p - 1) z++; int m = s; ll c = pow_mod(z, q, p); ll t = pow_mod(n, q, p); ll r = pow_mod(n, (q + 1) / 2, p); while (t != 1) { int m2 = 1; for (ll tmp = ba.mul(t, t); tmp != 1; tmp = ba.mul(tmp, tmp), m2++); ll b = pow_mod(c, 1 << (m - m2 - 1), p); m = m2, c = ba.mul(b, b), t = ba.mul(t, c), r = ba.mul(r, b); } return r; } template vector convolution_anymod(const vector &A, const vector &B) { int N = A.size(), M = B.size(); if (min(N, M) <= 300) { using mint = static_modint; vector A2(N), B2(M); for (int i = 0; i < N; i++) A2[i] = A[i]; for (int j = 0; j < M; j++) B2[j] = B[j]; vector C2(N + M - 1, 0); for (int i = 0; i < N; i++) for (int j = 0; j < M; j++) C2[i + j] += A2[i] * B2[j]; vector C(N + M - 1); for (int i = 0; i < N + M - 1; i++) C[i] = C2[i].val(); return C; } constexpr ll MOD1 = 167772161, MOD2 = 469762049, MOD3 = 1224736769; using mint2 = static_modint; using mint3 = static_modint; using mint4 = static_modint; constexpr int i1_2 = internal::inv_gcd(MOD1, MOD2).second; constexpr int i12_3 = internal::inv_gcd(MOD1 * MOD2, MOD3).second; constexpr int m12_4 = MOD1 * MOD2 % MOD; auto C1 = convolution(A, B); auto C2 = convolution(A, B); auto C3 = convolution(A, B); vector C(N + M - 1); for (ll i = 0; i < N + M - 1; i++) { int c1 = C1[i], c2 = C2[i], c3 = C3[i]; int t1 = (mint2(c2 - c1) * mint2::raw(i1_2)).val(); mint3 x2_m3 = mint3::raw(c1) + mint3::raw(t1) * mint3::raw(MOD1); mint4 x2_m = mint4::raw(c1) + mint4::raw(t1) * mint4::raw(MOD1); int t2 = ((mint3::raw(c3) - x2_m3) * mint3::raw(i12_3)).val(); C[i] = (x2_m + mint4::raw(t2) * mint4::raw(m12_4)).val(); } return C; } template vector convolution_anymod(const vector &A, const vector &B, const int MOD) { int N = A.size(), M = B.size(); if (min(N, M) <= 300) { using mint = dynamic_modint<100>; mint::set_mod(MOD); vector A2(N), B2(M); for (int i = 0; i < N; i++) A2[i] = A[i]; for (int j = 0; j < M; j++) B2[j] = B[j]; vector C2(N + M - 1, 0); for (int i = 0; i < N; i++) for (int j = 0; j < M; j++) C2[i + j] += A2[i] * B2[j]; vector C(N + M - 1); for (int i = 0; i < N + M - 1; i++) C[i] = C2[i].val(); return C; } constexpr ll MOD1 = 167772161, MOD2 = 469762049, MOD3 = 1224736769; using mint2 = static_modint; using mint3 = static_modint; using mint4 = dynamic_modint<100>; mint4::set_mod(MOD); constexpr int i1_2 = internal::inv_gcd(MOD1, MOD2).second; constexpr int i12_3 = internal::inv_gcd(MOD1 * MOD2, MOD3).second; auto C1 = convolution(A, B); auto C2 = convolution(A, B); auto C3 = convolution(A, B); vector C(N + M - 1); for (ll i = 0; i < N + M - 1; i++) { int c1 = C1[i], c2 = C2[i], c3 = C3[i]; int t1 = (mint2(c2 - c1) * mint2::raw(i1_2)).val(); mint3 x2_m3 = mint3::raw(c1) + mint3::raw(t1) * mint3::raw(MOD1); mint4 x2_m = mint4::raw(c1) + mint4::raw(t1) * mint4::raw(MOD1); int t2 = ((mint3::raw(c3) - x2_m3) * mint3::raw(i12_3)).val(); C[i] = (x2_m + mint4::raw(t2) * mint4::raw(MOD1) * mint4::raw(MOD2)).val(); } return C; } template vector> convolution_anymod(const vector> &A, const vector> &B) { int N = A.size(), M = B.size(); vector A2(N), B2(M); for (int i = 0; i < N; i++) A2[i] = A[i].val(); for (int i = 0; i < M; i++) B2[i] = B[i].val(); vector C2 = convolution_anymod(A2, B2); vector> C(N + M - 1); for (int i = 0; i < N + M - 1; i++) C[i] = static_modint::raw(C2[i]); return C; } template vector> convolution_anymod(const vector> &A, const vector> &B) { int N = A.size(), M = B.size(); vector A2(N), B2(M); for (int i = 0; i < N; i++) A2[i] = A[i].val(); for (int i = 0; i < M; i++) B2[i] = B[i].val(); vector C2 = convolution_anymod(A2, B2, dynamic_modint::mod()); vector> C(N + M - 1); for (int i = 0; i < N + M - 1; i++) C[i] = dynamic_modint::raw(C2[i]); return C; } template struct LagrangeInterpolation { int D; vector Y, fac, finv, prodl, prodr; template LagrangeInterpolation(const vector &y) { D = (int)y.size() - 1; Y.resize(D + 1); for (int i = 0; i <= D; i++) { Y[i] = y[i]; } fac.resize(D + 1), finv.resize(D + 1); fac[0] = 1; for (int i = 1; i <= D; i++) fac[i] = fac[i - 1] * i; finv[D] = fac[D].inv(); for (int i = D - 1; i >= 0; i--) finv[i] = finv[i + 1] * (i + 1); prodl.resize(D + 2), prodr.resize(D + 2); } T1 eval(T1 x) { prodl[0] = 1; for (int i = 0; i <= D; i++) { prodl[i + 1] = prodl[i] * (x - i); } prodr[D + 1] = 1; for (int i = D; i >= 0; i--) { prodr[i] = prodr[i + 1] * (x - i); } T1 res = 0; for (int i = 0; i <= D; i++) { T1 tmp = Y[i] * prodl[i] * prodr[i + 1] * finv[i] * finv[D - i]; if ((D - i) % 2 == 0) res += tmp; else res -= tmp; } return res; } }; // https://opt-cp.com/fps-implementation/ // https://qiita.com/hotman78/items/f0e6d2265badd84d429a // https://opt-cp.com/fps-fast-algorithms/ // https://maspypy.com/%E5%A4%9A%E9%A0%85%E5%BC%8F%E3%83%BB%E5%BD%A2%E5%BC%8F%E7%9A%84%E3%81%B9%E3%81%8D%E7%B4%9A%E6%95%B0-%E9%AB%98%E9%80%9F%E3%81%AB%E8%A8%88%E7%AE%97%E3%81%A7%E3%81%8D%E3%82%8B%E3%82%82%E3%81%AE template struct FormalPowerSeries : vector { private: static vector fac, finv, invmint; void calc(int n) { while ((int)fac.size() <= n) { int i = fac.size(); fac.emplace_back(fac[i - 1] * i); invmint.emplace_back(-invmint[T::mod() % i] * (T::mod() / i)); finv.emplace_back(finv[i - 1] * invmint[i]); } } public: T get_fac(int n) { calc(n); return fac[n]; } T get_finv(int n) { calc(n); return finv[n]; } T get_invmint(int n) { calc(n); return invmint[n]; } using vector::vector; using vector::operator=; using F = FormalPowerSeries; using S = vector>; FormalPowerSeries(const S &f, int n = -1) { if (n == -1) n = f.back().first + 1; (*this).assign(n, T(0)); for (auto [d, a] : f) (*this)[d] += a; } F operator-() const { F res(*this); for (auto &a : res) a = -a; return res; } F operator*=(const T &k) { for (auto &a : *this) a *= k; return *this; } F operator*(const T &k) const { return F(*this) *= k; } friend F operator*(const T k, const F &f) { return f * k; } F operator/=(const T &k) { *this *= k.inv(); return *this; } F operator/(const T &k) const { return F(*this) /= k; } F &operator+=(const F &g) { int n = (*this).size(), m = g.size(); (*this).resize(max(n, m), T(0)); for (int i = 0; i < m; i++) (*this)[i] += g[i]; return *this; } F operator+(const F &g) const { return F(*this) += g; } F &operator-=(const F &g) { int n = (*this).size(), m = g.size(); (*this).resize(max(n, m), T(0)); for (int i = 0; i < m; i++) (*this)[i] -= g[i]; return *this; } F operator-(const F &g) const { return F(*this) -= g; } F &operator<<=(const ll d) { int n = (*this).size(); (*this).insert((*this).begin(), min(ll(n), d), T(0)); (*this).resize(n); return *this; } F operator<<(const ll d) const { return F(*this) <<= d; } F &operator>>=(const ll d) { int n = (*this).size(); (*this).erase((*this).begin(), (*this).begin() + min(ll(n), d)); (*this).resize(n, T(0)); return *this; } F operator>>(const ll d) const { return F(*this) >>= d; } F &operator*=(const S &g) { int n = (*this).size(); auto [d, c] = g.front(); if (d != 0) c = 0; for (int i = n - 1; i >= 0; i--) { (*this)[i] *= c; for (auto &[j, b] : g) { if (j == 0) continue; if (j > i) break; (*this)[i] += (*this)[i - j] * b; } } return *this; } F operator*(const S &g) const { return F(*this) *= g; } F &operator/=(const S &g) { int n = (*this).size(); auto [d, c] = g.front(); assert(d == 0 && c != T(0)); T inv_c = c.inv(); for (int i = 0; i < n; i++) { for (auto &[j, b] : g) { if (j == 0) continue; if (j > i) break; (*this)[i] -= (*this)[i - j] * b; } (*this)[i] *= inv_c; } return *this; } F operator/(const S &g) const { return F(*this) /= g; } // (1 + cx^d) を掛ける F multiply(const int d, const T c) { int n = (*this).size(); if (c == T(1)) { for (int i = n - 1 - d; i >= 0; i--) (*this)[i + d] += (*this)[i]; } else if (c == T(-1)) { for (int i = n - 1 - d; i >= 0; i--) (*this)[i + d] -= (*this)[i]; } else { for (int i = n - 1 - d; i >= 0; i--) (*this)[i + d] += (*this)[i] * c; } return *this; } F multiplication(const int d, const T c) const { return multiply(F(*this)); } // (1 + cx^d) で割る F divide(const int d, const T c) { int n = (*this).size(); if (c == T(1)) { for (int i = 0; i < n - d; i++) (*this)[i + d] -= (*this)[i]; } else if (c == T(-1)) { for (int i = 0; i < n - d; i++) (*this)[i + d] += (*this)[i]; } else { for (int i = 0; i < n - d; i++) (*this)[i + d] -= (*this)[i] * c; } return *this; } F division(const int d, const T c) const { return divide(F(*this)); } template F convolution2(const vector> &A, const vector> &B, const int d = -1) const { F res; if (is_ntt_friendly) res = convolution(A, B); else res = convolution_anymod(A, B); if (d != -1 && (int)res.size() > d) res.resize(d); return res; } template F convolution2(const vector> &A, const vector> &B, const int d = -1) const { F res; res = convolution_anymod(A, B); if (d != -1 && (int)res.size() > d) res.resize(d); return res; } F &operator*=(const F &g) { int n = (*this).size(); if (n == 0) return *this; *this = convolution2(*this, g, n); return *this; } F operator*(const F &g) const { return F(*this) *= g; } template void butterfly2(FormalPowerSeries, true> &A) const { internal::butterfly(A); } template void butterfly2(FormalPowerSeries, false> &A) const { assert(false); } template void butterfly2(FormalPowerSeries, false> &A) const { assert(false); } template void butterfly_inv2(FormalPowerSeries, true> &A) const { internal::butterfly_inv(A); } template void butterfly_inv2(FormalPowerSeries, false> &A) const { assert(false); } template void butterfly_inv2(FormalPowerSeries, false> &A) const { assert(false); } // mod (x^n - 1) をとったものを返す F circular_mod(int n) const { F res(n, T(0)); for (int i = 0; i < (int)(*this).size(); i++) res[i % n] += (*this)[i]; return res; } F inv(int d = -1) const { int n = (*this).size(); assert(!(*this).empty() && (*this).at(0) != T(0)); if (d == -1) d = n; //assert(d > 0); F f, g2; F g{(*this).front().inv()}; while ((int)g.size() < d) { if (is_ntt_friendly) { int m = g.size(); f = F{(*this).begin(), (*this).begin() + min(n, 2 * m)}; g2 = F(g); f.resize(2 * m, T(0)), butterfly2(f); g2.resize(2 * m, T(0)), butterfly2(g2); for (int i = 0; i < 2 * m; i++) f[i] *= g2[i]; butterfly_inv2(f); f.erase(f.begin(), f.begin() + m); f.resize(2 * m, T(0)), butterfly2(f); for (int i = 0; i < 2 * m; i++) f[i] *= g2[i]; butterfly_inv2(f); T iz = T(2 * m).inv(); iz *= -iz; for (int i = 0; i < m; i++) f[i] *= iz; g.insert(g.end(), f.begin(), f.begin() + m); } else { g.resize(2 * g.size(), T(0)); g *= F{T(2)} - g * (*this); } } return {g.begin(), g.begin() + d}; } F &operator/=(const F &g) { *this *= g.inv((*this).size()); return *this; } F operator/(const F &g) const { return F(*this) *= g.inv((*this).size()); } F differentiate() { *this >>= 1; for (int i = 0; i < int((*this).size()) - 1; i++) (*this)[i] *= i + 1; return *this; } F differential() const { return F(*this).differentiate(); } F integrate() { int n = (*this).size(); vector minv(n); minv[1] = T(1); *this <<= 1; for (int i = 2; i < n; i++) { minv[i] = -minv[T::mod() % i] * (T::mod() / i); (*this)[i] *= minv[i]; } return *this; } F integral() const { return F(*this).integrate(); } F log() const { assert((*this).front() == T(1)); return ((*this).differential() / (*this)).integral(); } F exp() const // https://arxiv.org/pdf/1301.5804.pdf { int n = (*this).size(); assert(n != 0 && (*this).front() == T(0)); //* if (is_ntt_friendly) { F f{T(1)}, g{T(1)}; F dh = (*this).differential(); F f2, g2, f3, q, s, h, u; g2 = {T(0)}; while ((int)f.size() < n) { int m = f.size(); T im = T(m).inv(), i2m = T(2 * m).inv(); f2 = F(f); f2.resize(2 * m), butterfly2(f2); // a F f3(f); butterfly2(f3); for (int i = 0; i < m; i++) f3[i] *= g2[i]; butterfly_inv2(f3); f3.erase(f3.begin(), f3.begin() + m / 2); f3.resize(m, T(0)), butterfly2(f3); for (int i = 0; i < m; i++) f3[i] *= g2[i]; butterfly_inv2(f3); for (int i = 0; i < m / 2; i++) f3[i] *= -im * im; g.insert(g.end(), f3.begin(), f3.begin() + m / 2); g2 = F(g), g2.resize(2 * m), butterfly2(g2); // b, c q = F(dh); q.resize(2 * m); for (int i = m - 1; i < 2 * m; i++) q[i] = T(0); butterfly2(q); for (int i = 0; i < 2 * m; i++) q[i] *= f2[i]; butterfly_inv2(q); q = q.circular_mod(m); for (int i = 0; i < m; i++) q[i] *= i2m; // d, e q.resize(m + 1); s = ((f.differential() - q) << 1).circular_mod(m); s.resize(2 * m); butterfly2(s); for (int i = 0; i < 2 * m; i++) s[i] *= g2[i]; butterfly_inv2(s); for (int i = 0; i < m; i++) s[i] *= i2m; s.resize(m); // f, g h = (*this); h.resize(2 * m), s.resize(2 * m); u = (h - (s << (m - 1)).integral()) >> m; butterfly2(u); for (int i = 0; i < 2 * m; i++) u[i] *= f2[i]; butterfly_inv2(u); for (int i = 0; i < m; i++) u[i] *= i2m; u.resize(m); // h f.insert(f.end(), u.begin(), u.end()); } return {f.begin(), f.begin() + n}; } else //*/ { F f{T(1)}, g{T(1)}; while ((int)f.size() < n) { int m = f.size(); g = convolution2(g, F{T(2)} - f * g, m); F q = (*this).differential(); q.resize(m - 1); F r = f.convolution2(f, q).circular_mod(m); r.resize(m + 1); F s = ((f.differential() - r) << 1).circular_mod(m); F t = g * s; F h = (*this); h.resize(2 * m), t.resize(2 * m); F u = (h - (t << (m - 1)).integral()) >> m; F v = f * u; f.insert(f.end(), v.begin(), v.end()); } return {f.begin(), f.begin() + n}; /* F f{T(1)}; while ((int)f.size() < n) { int m = f.size(); f.resize(min(n, 2 * m), T(0)); f *= (*this) + F{T(1)} - f.log(); } return f; //*/ } } F pow(const ll k) const { if (k == 0) { F res((*this).size(), T(0)); res[0] = T(1); return res; } int n = (*this).size(), d; for (d = 0; d < n; d++) { if ((*this)[d] != T(0)) break; } if (d == n) return F(n, 0); F res = F(*this) >> d; T c = res[0]; res /= c; res = (res.log() * T(k)).exp(); res *= c.pow(k), res <<= (d != 0 && k > n ? n : d * k); return res; } F powmod(ll k, const F &g) const { F res(2 * g.size(), 0); res.front() = 1; F tmp = (*this) % g; tmp.resize(g.size()); while (k > 0) { if (k & 1) { res *= tmp; res %= g; res.resize(2 * g.size()); } tmp = tmp.convolution2(tmp, tmp); tmp %= g; tmp.resize(g.size()); k >>= 1; } return res; } // f(x)^k mod (x^n - 1) F powmod_circular(ll k, ll n) const { F res(n, 0); res.front() = 1; F tmp = (*this).circular_mod(n); while (k > 0) { if (k & 1) res = res.convolution2(res, tmp).circular_mod(n); tmp = tmp.convolution2(tmp, tmp).circular_mod(n); k >>= 1; } return res; } // 素数 mod を要求 // 存在しないなら空配列を返す F sqrt() const { int n = (*this).size(), d; for (d = 0; d < n; d += 2) { if ((*this)[d] != 0) break; if (d + 1 < n && (*this)[d + 1] != 0) return F(0); } if (d >= n) return F(n, 0); T a = (*this)[d]; int p = T::mod(); int r = sqrt_mod(a.val(), p); if (r == -1) return F(0); T inv_2 = T(2).inv(); F f = F(*this) >> d, res = F{r}; while (res.size() < f.size()) { res.resize(min(f.size(), 2 * res.size()), T(0)); res = (res + res.inv() * f) * inv_2; } res <<= d / 2; return res; } F div_poly(const F &g) const { F f2 = F(*this), g2 = F(g); while (!f2.empty() && f2.back() == T(0)) f2.pop_back(); while (!g2.empty() && g2.back() == T(0)) g2.pop_back(); int n = f2.size(), m = g2.size(); int k = n - m + 1; if (k <= 0) return F{}; reverse(f2.begin(), f2.end()); reverse(g2.begin(), g2.end()); f2.resize(k, T(0)), g2.resize(k, T(0)); F q = f2 / g2; reverse(q.begin(), q.end()); while (!q.empty() && q.back() == T(0)) q.pop_back(); return q; } pair divmod(const F &g) const { int m = g.size(); assert(m != 0); F q = (*this).div_poly(g); F f3 = F(*this), g3 = F(g), q3 = F(q); f3.resize(m - 1, T(0)), g3.resize(m - 1, T(0)), q3.resize(m - 1, T(0)); F r = f3 - q3 * g3; while (!r.empty() && r.back() == T(0)) r.pop_back(); return make_pair(q, r); } F operator%(const F &g) const { return (*this).divmod(g).second; } F &operator%=(const F &g) { return (*this) = (*this) % g; } F div_poly(const S &g) const { F f2 = F(*this); while (!f2.empty() && f2.back() == T(0)) f2.pop_back(); assert(!g.empty()); int n = f2.size(), m = g.back().first + 1; int k = n - m + 1; if (k <= 0) return F{}; reverse(f2.begin(), f2.end()); S g2(g.size()); for (int i = 0; i < (int)g.size(); i++) g2[(int)g.size() - 1 - i] = make_pair(m - 1 - g[i].first, g[i].second); f2.resize(k, T(0)); F q = f2 / g2; reverse(q.begin(), q.end()); while (!q.empty() && q.back() == T(0)) q.pop_back(); return q; } pair divmod(const S &g) const { assert(!g.empty()); int m = g.back().first + 1; F q = (*this).div_poly(g); F f3 = F(*this), q3 = F(q); f3.resize(m - 1, T(0)), q3.resize(m - 1, T(0)); F r = f3 - q3 * g; while (!r.empty() && r.back() == T(0)) r.pop_back(); return make_pair(q, r); } F operator%(const S &g) const { return (*this).divmod(g).second; } F &operator%=(const S &g) { return (*this) = (*this) % g; } T eval(const T &x) const { T res(0); for (int i = (int)(*this).size() - 1; i >= 0; i--) { res *= x; res += (*this)[i]; } return res; } // 各係数 a_n を n! で割る F to_egf() { for (int i = 0; i < (int)(*this).size(); i++) (*this)[i] *= get_finv(i); return (*this); } // 各係数 a_n を n! で割ったものを返す F get_egf() const { return F(*this).to_egf(); } // 各係数 a_n に n! をかける F to_ogf() { for (int i = 0; i < (int)(*this).size(); i++) (*this)[i] *= get_fac(i); return (*this); } // 各係数 a_n に n! をかけたものを返す F get_ogf() const { return F(*this).to_ogf(); } F taylor_shift(const T &c) const { int n = (*this).size(); F f = F(*this).get_ogf(); reverse(f.begin(), f.end()); F g = F(n); g[0] = 1; for (int i = 1; i < n; i++) g[i] = c * g[i - 1]; g.to_egf(); F h = f * g; reverse(h.begin(), h.end()); return h.to_egf(); } vector eval_multipoint(const vector &xs) const { int m0 = xs.size(), m = 1; while (m < m0) m <<= 1; vector node(2 * m, F{1}); for (int i = 0; i < m0; i++) node[m + i] = {-xs[i], T(1)}; for (int i = m - 1; i > 0; i--) node[i] = convolution2(node[i << 1], node[(i << 1) | 1]); node[1] = (*this).divmod(node[1]).second; for (int i = 2; i < m + m0; i++) node[i] = node[i >> 1].divmod(node[i]).second; vector res(m0); for (int i = 0; i < m0; i++) res[i] = node[m + i].empty() ? T(0) : node[m + i][0]; return res; } // i = 0..m-1 に対する f(ar^i) // https://noshi91.github.io/algorithm-encyclopedia/chirp-z-transform vector eval_multipoint_geo(int m, T a, T r) const { if (r == 0) { vector res(m, (*this).eval(0)); res.front() = (*this).eval(a); return res; } auto calc_pw = [&](T x, int k) -> vector { vector res(k); res.front() = 1; for (int i = 1; i < k; i++) res[i] = res[i - 1] * x; return res; }; auto get_pw_tri = [&](const vector &pw, int i) -> T { if (i == 0) return 1; return i % 2 == 0 ? pw[i - 1].pow(i / 2) : pw[i].pow((i - 1) / 2); }; int n = (*this).size(); T invr = r.inv(); vector pwa = calc_pw(a, n), pwr = calc_pw(r, n + m), pwir = calc_pw(invr, max(n, m)); vector s(n), t(n + m); for (int i = 0; i < n; i++) s[n - 1 - i] = (*this)[i] * pwa[i] * get_pw_tri(pwir, i); for (int i = 0; i < n + m; i++) t[i] = get_pw_tri(pwr, i); vector u = convolution2(s, t, n + m - 1); u.erase(u.begin(), u.begin() + n - 1); for (int i = 0; i < m; i++) u[i] *= get_pw_tri(pwir, i); return u; } }; // (次数, 係数) を昇順に並べたもの template struct SparseFormalPowerSeries : vector> { using vector>::vector; using vector>::operator=; using F = FormalPowerSeries; using S = SparseFormalPowerSeries; F to_fps(int n) const { F res(n, T(0)); for (auto [d, a] : (*this)) res[d] += a; return res; } SparseFormalPowerSeries(const F &f) { (*this).clear(); for (int i = 0; i < (int)f.size(); i++) { if (f[i] != T(0)) (*this).emplace_back(make_pair(i, f[i])); } } S operator-() const { S res(*this); for (auto &[d, a] : res) a = -a; return res; } S operator*=(const T &k) { for (auto &[d, a] : (*this)) a *= k; return (*this); } S operator/=(const T &k) { (*this) *= k.inv(); return (*this); } S operator*(const T &k) const { return S(*this) *= k; } S operator/(const T &k) const { return S(*this) /= k; } friend S operator*(const T k, const S &f) { return f * k; } S operator+(const S &g) const { S res; int n = (*this).size(), m = g.size(), i = 0, j = 0; while (i < n || j < m) { pair tmp; if (j == m || (i != n && (*this)[i].first <= g[j].first)) tmp = (*this)[i++]; else tmp = g[j++]; if (!res.empty() && res.back().first == tmp.first) res.back().second += tmp.second; else res.emplace_back(tmp); } return res; } S operator-(const S &g) const { S res; int n = (*this).size(), m = g.size(), i = 0, j = 0; while (i < n || j < m) { pair tmp; if (j == m || (i != n && (*this)[i].first <= g[j].first)) tmp = (*this)[i++]; else { tmp = g[j++]; tmp.second = -tmp.second; } if (!res.empty() && res.back().first == tmp.first) res.back().second += tmp.second; else res.emplace_back(tmp); } return res; } S operator*(const S &g) const { S res; for (auto [d, a] : (*this)) for (auto [e, b] : g) res.emplace_back(make_pair(d + e, a * b)); sort(res.begin(), res.end(), [&](pair p1, pair p2) { return p1.first < p2.first; }); S res2; for (auto da : res) { auto [d, a] = da; if (res2.empty() || res2.back().first != d) res2.emplace_back(da); else res2.back().second += a; } return res2; } S operator+=(const S &g) { return (*this) = (*this) + g; } S operator-=(const S &g) { return (*this) = (*this) - g; } S operator*=(const S &g) { return (*this) = (*this) * g; } S operator<<=(ll k) { for (auto &[d, a] : (*this)) d += k; return (*this); } S operator<<(ll k) const { return (*this) <<= k; } S operator>>(ll k) const { S res; for (auto [d, a] : (*this)) { d -= k; if (d >= 0) res.emplace_back(make_pair(d, a)); } return res; } S operator>>=(ll k) { return (*this) = (*this) >> k; } F inv(int n) const { F f(n, T(0)); f.front() = T(1); return f / (*this); } S differentiate() { for (auto &[d, a] : (*this)) a *= d--; if (!(*this).empty() && (*this).front().first == -1) (*this).erase((*this).begin()); return (*this); } S differential() const { return S(*this).differentiate(); } S integrate() { for (auto &[d, a] : (*this)) a /= T(++d); return (*this); } S integral() const { return S(*this).integrate(); } F log(int n) const { F f = (*this).to_fps(n); return (f.differential() / (*this)).integral(); } // 微分方程式 a(x)F'(x) + b(x)F(x) = 0, [x^0]F(x) = 1 を満たす F を n 項まで求める // [x^0]a(x) = 1 である必要がある F diffeq(const S &a, const S &b, int n) const { assert(a.front().first == 0 && a.front().second == 1); vector minv(n); minv[1] = T(1); for (int i = 2; i < n; i++) minv[i] = -minv[T::mod() % i] * (T::mod() / i); F f(n, T(0)); f[0] = T(1); for (int k = 0; k < n - 1; k++) { for (auto [i, ai] : a) { if (0 <= k - i + 1 && k - i + 1 < k + 1) f[k + 1] -= ai * (k - i + 1) * f[k - i + 1]; } for (auto [j, bj] : b) { if (0 <= k - j && k - j < k + 1) f[k + 1] -= bj * f[k - j]; } f[k + 1] *= minv[k + 1]; } return f; } F exp(int n) const { return diffeq(S{{0, 1}}, -((*this).differential()), n); } // m >= 0 のときは O(nk) (k: sparse の non-zero の個数) // m < 0 のときは O((n + d_0 m)k) F pow(ll m, int n) const { S f(*this); if (f.empty()) { F res(n, T(0)); if (m == 0) res.front() = T(1); return res; } auto [d0, a0] = f.front(); T a0_inv = a0.inv(); for (auto &[d, a] : f) d -= d0, a *= a0_inv; if (m >= 0) { F g = diffeq(f, -m * f.differential(), n); return (g * a0.pow(m)) << mul_limited(d0, m); } else { F g = diffeq(f, -m * f.differential(), n + (d0 * (-m))); F h = (g * a0_inv.pow(-m)) >> (d0 * (-m)); h.resize(n); return h; } } // 素数 mod を要求 // 存在しないなら空配列を返す F sqrt(int n) const { S f(*this); if (f.empty()) return F(n, T(0)); auto [d0, a0] = f.front(); if (d0 % 2 != 0) return F(0); if (d0 >= n) return F(n, T(0)); int p = T::mod(); int r = sqrt_mod(a0.val(), p); if (r == -1) return F(0); T a0_inv = a0.inv(); T inv_2 = T(2).inv(); for (auto &[d, a] : f) d -= d0, a *= a0_inv; F g = diffeq(f, -inv_2 * f.differential(), n); return ((g * r) << (d0 / 2)); } }; template vector FormalPowerSeries::fac{1, 1}; template vector FormalPowerSeries::finv{1, 1}; template vector FormalPowerSeries::invmint{0, 1}; template struct RationalFormalPowerSeries { using F = FormalPowerSeries; using R = RationalFormalPowerSeries; F num, den; R operator-() const { R res(*this); res.num = -res.num; return res; } R operator*=(const T &k) { (*this).num *= k; return *this; } R operator*(const T &k) const { return R(*this) *= k; } friend R operator*(const T k, const R &r) { return r * k; } R operator/=(const T &k) { (*this).den *= k; return k; } R operator/(const T &k) const { return R(*this) /= k; } R &operator+=(const R &r) { // ここうまくやると FFT の回数が節約できる気がする // うまくやらないと次数に偏りがある場合にかえって遅くなったりしそうで面倒 F f, g; f = f.convolution2((*this).num, r.den); g = g.convolution2((*this).den, r.num); (*this).num = f + g; (*this).den = (*this).den.convolution2((*this).den, r.den); return *this; } R operator+(const R &r) const { return R(*this) += r; } R &operator-=(const R &r) { F f, g; f = f.convolution2((*this).num, r.den); g = g.convolution2((*this).den, r.num); (*this).num = f - g; (*this).den = (*this).den.convolution2((*this).den, r.den); return *this; } R operator-(const R &r) const { return R(*this) -= r; } R operator*=(const R &r) { (*this).num = (*this).num.convolution2((*this).num, r.num); (*this).den = (*this).den.convolution2((*this).den, r.den); return *this; } R operator*(const R &r) const { return R(*this) *= r; } R operator/=(const R &r) { (*this).num = (*this).num.convolution2((*this).num, r.den); (*this).den = (*this).den.convolution2((*this).den, r.num); return *this; } R operator/(const R &r) const { return R(*this) /= r; } R inv() { R res(*this); swap(res.num, res.den); return res; } }; template FormalPowerSeries convolution_many(const vector> &fs, int d = -1) { using F = FormalPowerSeries; if (fs.empty()) { if (d == -1) d = 0; F res(d + 1, T(0)); res.front() = T(1); return res; } deque deq; for (auto f : fs) deq.push_back(f); while ((int)deq.size() > 1) { F f = deq.front(); deq.pop_front(); F g = deq.front(); deq.pop_front(); f = f.convolution2(f, g, d); deq.push_back(f); } if (d != -1) deq.front().resize(d); return deq.front(); } template RationalFormalPowerSeries rational_sum(const vector> &rs, int d = -1) { using R = RationalFormalPowerSeries; if (rs.empty()) return R{{1}, {1}}; deque deq; for (auto &r : rs) deq.emplace_back(r); while ((int)deq.size() > 1) { R r1 = deq.front(); deq.pop_front(); R r2 = deq.front(); deq.pop_front(); R r3 = r1 + r2; if (d != -1) { if ((int)r3.num.size() > d) r3.num.resize(d); if ((int)r3.den.size() > d) r3.den.resize(d); } deq.emplace_back(r3); } if (d != -1) deq.front().num.resize(d), deq.front().den.resize(d); return deq.front(); } template FormalPowerSeries interpolation(const vector &xs, const vector &ys) { using F = FormalPowerSeries; using R = RationalFormalPowerSeries; int n = xs.size(); assert(n == ys.size()); vector fs(n); for (int i = 0; i < n; i++) fs[i] = F{-xs[i], T(1)}; F g = convolution_many(fs); F h = g.differential(); vector a = h.eval_multipoint(xs); vector rs(n); for (int i = 0; i < n; i++) rs[i] = R{F{ys[i] / a[i]}, fs[i]}; R q = rational_sum(rs, n); return q.num; } // f(i) = ys[i] で定まる多項式 f(x) について f(c), …, f(c + M - 1) を求める template vector sample_points_shift(const vector &ys, int M, T c) { using F = FormalPowerSeries; F f; int N = ys.size(); vector a; { vector p(N), q(N); for (int i = 0; i < N; i++) { p[i] = ys[i] * f.get_finv(i); q[i] = i % 2 == 0 ? f.get_finv(i) : -f.get_finv(i); } a = f.convolution2(p, q); a.resize(N); } vector b; { vector p(N), q(N); T tmp = 1; for (int i = 0; i < N; i++) { p[i] = a[i] * f.get_fac(i); q[i] = tmp * f.get_finv(i); tmp *= c - i; } reverse(q.begin(), q.end()); b = f.convolution2(p, q); b.erase(b.begin(), b.begin() + N - 1); for (int i = 0; i < N; i++) b[i] *= f.get_finv(i); } vector res; { vector p(M); for (int i = 0; i < M; i++) p[i] = f.get_finv(i); res = f.convolution2(b, p); res.resize(M); for (int i = 0; i < M; i++) res[i] *= f.get_fac(i); } return res; } // https://suisen-kyopro.hatenablog.com/entry/2023/11/22/201600 // 前計算 O(K 2^K + (P/2^K) log K), クエリ O(2^K) template struct FactorialFast { private: const int P, K; vector Y, Z, fac; public: FactorialFast(const int K = 9) : P(T::mod()), K(K) { Y = {1}; for (int i = 0; i < K; i++) { Z = sample_points_shift(Y, (1 << (i + 2)) - (1 << i), 1 << i); Z.insert(Z.begin(), Y.begin(), Y.end()); Y.resize(1 << (i + 1)); for (int j = 0; j < (1 << (i + 1)); j++) Y[j] = Z[2 * j] * Z[2 * j + 1] * T::raw((1 << i) * (2 * j + 1)); } if ((1 << K) <= P / (1 << K)) { Z = sample_points_shift(Y, P / (1 << K), 1 << K); Y.insert(Y.end(), Z.begin(), Z.end()); } fac.resize(P / (1 << K) + 1); fac.at(0) = 1; for (int i = 0; i < P / (1 << K); i++) fac[i + 1] = fac[i] * Y[i] * T::raw((1 + i) * (1 << K)); } T query(ll n) { if (n >= T::mod()) return 0; T res = fac.at(n / (1 << K)); for (int j = n / (1 << K) * (1 << K) + 1; j <= n; j++) res *= T::raw(j); return res; } }; // f(x)/prod[i](1-a[i]x) = sum[i] c[i]/(1-a[i]x) なる c を求める template vector partial_fraction_decomposition(const FormalPowerSeries &f, const vector &as) { using F = FormalPowerSeries; int N = as.size(); assert((int)f.size() <= N); vector bs(N); for (int i = 0; i < N; i++) bs[i] = as[i].inv(); vector gs(N); for (int i = 0; i < N; i++) gs[i] = F{T(1), -as[i]}; F g = convolution_many(gs); F dg = g.differential(); vector ys = f.eval_multipoint(bs), zs = dg.eval_multipoint(bs); vector cs(N); for (int i = 0; i < N; i++) cs[i] = -as[i] * ys[i] / zs[i]; return cs; } // n = 0, 1, …, N-1 に対して n^k を列挙 template vector enum_pow(int N, int k) { vector minfactor(N, -1); for (int i = 2; i < N; i++) { if (minfactor[i] != -1) continue; for (int j = 2 * i; j < N; j += i) minfactor[j] = i; } vector power(N); for (int i = 0; i < N; i++) { if (minfactor[i] == -1) power[i] = T(i).pow(k); else power[i] = power[minfactor[i]] * power[i / minfactor[i]]; } return power; } // sum_[i = 0..infty] r^i i^d template T sum_of_exp_times_poly_limit(T r, int d) { using F = FormalPowerSeries; vector pws = enum_pow(d + 2, d); vector pwr(d + 2, 1); for (int i = 0; i < d + 1; i++) pwr[i + 1] = pwr[i] * r; F f(d + 2), g(d + 2); f.front() = pws.front() * pwr.front(); for (int i = 0; i < d + 1; i++) f[i + 1] = f[i] + pws[i + 1] * pwr[i + 1]; for (int i = 0; i <= d + 1; i++) g[i] = (i % 2 == 0 ? 1 : -1) * pwr[i] * g.get_fac(d + 1) * g.get_finv(i) * g.get_finv(d + 1 - i); T c = 0; for (int i = 0; i <= d + 1; i++) c += f[i] * g[d + 1 - i]; c /= accumulate(g.begin(), g.end(), T(0)); return c; } // sum_[i = 0..n-1] r^i i^d template T sum_of_exp_times_poly(T r, int d, ll n) { using F = FormalPowerSeries; if (n == 0) return 0; if (r == 0) return d == 0 ? 1 : 0; vector pws = enum_pow(d + 2, d); vector pwr(d + 2, 1); for (int i = 0; i < d + 1; i++) pwr[i + 1] = pwr[i] * r; F f(d + 2), g(d + 2); f.front() = pws.front() * pwr.front(); for (int i = 0; i < d + 1; i++) f[i + 1] = f[i] + pws[i + 1] * pwr[i + 1]; if (r == 1) { LagrangeInterpolation lag(f); return lag.eval(n - 1); } for (int i = 0; i <= d + 1; i++) g[i] = (i % 2 == 0 ? 1 : -1) * pwr[i] * g.get_fac(d + 1) * g.get_finv(i) * g.get_finv(d + 1 - i); T c = 0; for (int i = 0; i <= d + 1; i++) c += f[i] * g[d + 1 - i]; c /= accumulate(g.begin(), g.end(), T(0)); F h = f - F(d + 2, c); { T rinv = r.inv(); T pwrinv = 1; for (int i = 0; i <= d + 1; i++) { h[i] *= pwrinv; pwrinv *= rinv; } } LagrangeInterpolation lag(h); return c + r.pow(n - 1) * lag.eval(n - 1); } // prod[d in D](1 + cx^d) を M 次の項まで求める template FormalPowerSeries multiply_many(const int &M, const T &c, const vector &D) { using F = FormalPowerSeries; vector cnt(M + 1, 0); for (auto d : D) { if (d < 0 || M < d) continue; cnt[d]++; } F f(M + 1, 0); for (int k = 1; k <= M; k++) { T pw = 1; for (int i = 1; k * i <= M; i++) { pw *= c; if (i & 1) f[k * i] += T::raw(cnt[k]) * pw * f.get_invmint(i); else f[k * i] -= T::raw(cnt[k]) * pw * f.get_invmint(i); } } return f.exp(); } // 多重集合 S の要素から何個か選んで総和を 0, 1, …, M にする方法の数 template FormalPowerSeries subset_sum(const int &M, const vector &S) { return multiply_many(M, T(1), S); } // 集合 S の各要素が無限個ある集合 T から何個か選んで総和を 0, 1, …, M にする方法の数 template FormalPowerSeries partition(const int &M, const vector &S) { return multiply_many(M, T(-1), S).inv(); } template FormalPowerSeries stirling1(const int &N) { using F = FormalPowerSeries; using S = SparseFormalPowerSeries; if (N == 0) return {1}; if (N == 1) return {0, 1}; if (N & 1) { F f = stirling1(N - 1); f.resize(N + 1, T(0)); return f * S{{0, 1 - N}, {1, 1}}; } else { F f = stirling1(N / 2); f.resize(N + 1, T(0)); F g = f.taylor_shift(-(N / 2)); return f * g; } } template FormalPowerSeries stirling2(const int &N) { using F = FormalPowerSeries; vector power = enum_pow(N + 1, N); F A(N + 1), B(N + 1); for (int i = 0; i <= N; i++) { A[i] = power[i] * A.get_finv(i); B[i] = (i & 1) ? -A.get_finv(i) : A.get_finv(i); } return A * B; } template FormalPowerSeries bernoulli_number(const int &N) { using F = FormalPowerSeries; F f(N + 1, T(0)); for (int i = 0; i <= N; i++) f[i] = f.get_finv(i + 1); return f.inv().to_ogf(); } // [x^N] P(x)/Q(x) を求める(P の次数は Q の次数より小さい) template T bostan_mori(const FormalPowerSeries &P, const FormalPowerSeries &Q, ll N) { using F = FormalPowerSeries; int d = (int)Q.size() - 1; assert((int)P.size() <= d); if (is_ntt_friendly) { int z = 1; while (z < 2 * d + 1) z <<= 1; T iz = T(z).inv(); F U = F(P), V = F(Q); U.resize(z), V.resize(z); while (N > 0) { U.butterfly2(U), V.butterfly2(V); for (int i = 0; i < z; i += 2) { T x = V[i + 1], y = V[i]; U[i] *= x, V[i] *= x; U[i + 1] *= y, V[i + 1] *= y; } U.butterfly_inv2(U), V.butterfly_inv2(V); for (int i = 0; i < (z >> 1); i++) { U[i] = U[2 * i + (N & 1)] * iz; V[i] = V[2 * i] * iz; } for (int i = (z >> 1); i < z; i++) U[i] = 0, V[i] = 0; N >>= 1; } return U.front() / V.front(); } else { F U = F(P), V = F(Q); U.resize(d), V.resize(d + 1); while (N > 0) { F U2 = F(U), V2 = F(V), V3 = F(V); for (int i = 1; i <= d; i += 2) V3[i] = -V3[i]; U2 *= V3, V2 *= V3; for (int i = 0; i <= d; i++) { U[i] = U2[2 * i + (N & 1)]; V[i] = V2[2 * i]; } N >>= 1; } return U.front() / V.front(); } } // a_n = sum[i = 1..d] c_i a_{n-i}(n ≥ d)を満たすとき、a_N を求める(A は 0-indexed で C は 1-indexed) template T linear_recurrence(const vector &A, const vector &C, ll N) { using F = FormalPowerSeries; int d = C.size(); assert((int)A.size() >= d); F Ga(d), Q(d + 1); Q[0] = 1; for (int i = 0; i < d; i++) Ga[i] = A[i], Q[i + 1] = -C[i]; F P = Ga * Q; return bostan_mori(P, Q, N); } // (P の次数) < (Q の次数) とする // P/Q = R + x^N (P'/Q) を満たす P' (R は N 次未満、P' は d 次未満) // [x^{N+n}](P/Q) = [x^n](P'/Q) 線形漸化的数列のシフト // 高速化の余地あり template FormalPowerSeries fiduccia(const FormalPowerSeries &P, const FormalPowerSeries &Q, ll N) { using F = FormalPowerSeries; assert(P.size() < Q.size()); F xinv = -(Q >> 1); if (Q[0] != 1) xinv /= Q[0]; return xinv.powmod(N, Q) * P % Q; } // a_n = sum[i = 1..d] c_i a_{n-i}(n ≥ d)を満たすとき、a_N, …, a_{N+d-1} を求める(A は 0-indexed で C は 1-indexed) template vector linear_recurrence_many(const vector &A, const vector &C, ll N) { using F = FormalPowerSeries; int d = C.size(); assert((int)A.size() >= d); F Ga(d), Q(d + 1); Q[0] = 1; for (int i = 0; i < d; i++) Ga[i] = A[i], Q[i + 1] = -C[i]; F P = Ga * Q; F P2 = fiduccia(P, Q, N); P2.resize(d); F Gb = P2 / Q + (P.div_poly(Q) >> N); return Gb; } // Σ[i = 0..M-1] a_i exp(b_i x) を N 項まで求める template FormalPowerSeries exp_sum(int N, const vector &a, const vector &b) { using F = FormalPowerSeries; using R = RationalFormalPowerSeries; assert(a.size() == b.size()); int M = a.size(); vector gs(M); for (int i = 0; i < M; i++) gs[i] = R{F{a[i]}, F{1, -b[i]}}; R g = rational_sum(gs, N); return (g.num / g.den).to_egf(); } // f(exp(kx)) を N 項まで求める template FormalPowerSeries eval_exp(FormalPowerSeries &f, T k, int N = -1) { if (N == -1) N = (int)f.size(); vector b(f.size()); for (int i = 0; i < (int)f.size(); i++) b[i] = k * i; return exp_sum(N, f, b); } /* using mint = modint998244353; const bool ntt = true; //*/ //* using mint = modint1000000007; const bool ntt = false; //*/ /* using mint = modint; const bool ntt = false; //*/ using fps = FormalPowerSeries; using sfps = SparseFormalPowerSeries; using rfps = RationalFormalPowerSeries; int main() { ll N; cin >> N; FactorialFast fac(13); cout << fac.query(N).val() << endl; }