#include #include #define rep(i,n) for(int i=0;i vi; typedef vector vl; typedef vector> vvi; typedef vector> vvl; typedef long double ld; typedef pair P; ostream& operator<<(ostream& os, const modint& a) {os << a.val(); return os;} template ostream& operator<<(ostream& os, const static_modint& a) {os << a.val(); return os;} template ostream& operator<<(ostream& os, const dynamic_modint& a) {os << a.val(); return os;} template istream& operator>>(istream& is, vector& v){int n = v.size(); assert(n > 0); rep(i, n) is >> v[i]; return is;} template ostream& operator<<(ostream& os, const pair& p){os << p.first << ' ' << p.second; return os;} template ostream& operator<<(ostream& os, const vector& v){int n = v.size(); rep(i, n) os << v[i] << (i == n - 1 ? "\n" : " "); return os;} template ostream& operator<<(ostream& os, const vector>& v){int n = v.size(); rep(i, n) os << v[i] << (i == n - 1 ? "\n" : ""); return os;} template void chmin(T& a, T b){a = min(a, b);} template void chmax(T& a, T b){a = max(a, b);} using mint = modint998244353; // thanks for Nyaan-san's library // https://nyaannyaan.github.io/library/fps/formal-power-series.hpp template struct NTT { static constexpr uint32_t get_pr() { uint32_t _mod = mint::mod(); using u64 = uint64_t; u64 ds[32] = {}; int idx = 0; u64 m = _mod - 1; for (u64 i = 2; i * i <= m; ++i) { if (m % i == 0) { ds[idx++] = i; while (m % i == 0) m /= i; } } if (m != 1) ds[idx++] = m; uint32_t _pr = 2; while (1) { int flg = 1; for (int i = 0; i < idx; ++i) { u64 a = _pr, b = (_mod - 1) / ds[i], r = 1; while (b) { if (b & 1) r = r * a % _mod; a = a * a % _mod; b >>= 1; } if (r == 1) { flg = 0; break; } } if (flg == 1) break; ++_pr; } return _pr; }; static constexpr uint32_t mod = mint::mod(); static constexpr uint32_t pr = get_pr(); static constexpr int level = __builtin_ctzll(mod - 1); mint dw[level], dy[level]; void setwy(int k) { mint w[level], y[level]; w[k - 1] = mint(pr).pow((mod - 1) / (1 << k)); y[k - 1] = w[k - 1].inv(); for (int i = k - 2; i > 0; --i) w[i] = w[i + 1] * w[i + 1], y[i] = y[i + 1] * y[i + 1]; dw[1] = w[1], dy[1] = y[1], dw[2] = w[2], dy[2] = y[2]; for (int i = 3; i < k; ++i) { dw[i] = dw[i - 1] * y[i - 2] * w[i]; dy[i] = dy[i - 1] * w[i - 2] * y[i]; } } NTT() { setwy(level); } void fft4(vector &a, int k) { if ((int)a.size() <= 1) return; if (k == 1) { mint a1 = a[1]; a[1] = a[0] - a[1]; a[0] = a[0] + a1; return; } if (k & 1) { int v = 1 << (k - 1); for (int j = 0; j < v; ++j) { mint ajv = a[j + v]; a[j + v] = a[j] - ajv; a[j] += ajv; } } int u = 1 << (2 + (k & 1)); int v = 1 << (k - 2 - (k & 1)); mint one = mint(1); mint imag = dw[1]; while (v) { // jh = 0 { int j0 = 0; int j1 = v; int j2 = j1 + v; int j3 = j2 + v; for (; j0 < v; ++j0, ++j1, ++j2, ++j3) { mint t0 = a[j0], t1 = a[j1], t2 = a[j2], t3 = a[j3]; mint t0p2 = t0 + t2, t1p3 = t1 + t3; mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag; a[j0] = t0p2 + t1p3, a[j1] = t0p2 - t1p3; a[j2] = t0m2 + t1m3, a[j3] = t0m2 - t1m3; } } // jh >= 1 mint ww = one, xx = one * dw[2], wx = one; for (int jh = 4; jh < u;) { ww = xx * xx, wx = ww * xx; int j0 = jh * v; int je = j0 + v; int j2 = je + v; for (; j0 < je; ++j0, ++j2) { mint t0 = a[j0], t1 = a[j0 + v] * xx, t2 = a[j2] * ww, t3 = a[j2 + v] * wx; mint t0p2 = t0 + t2, t1p3 = t1 + t3; mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag; a[j0] = t0p2 + t1p3, a[j0 + v] = t0p2 - t1p3; a[j2] = t0m2 + t1m3, a[j2 + v] = t0m2 - t1m3; } xx *= dw[__builtin_ctzll((jh += 4))]; } u <<= 2; v >>= 2; } } void ifft4(vector &a, int k) { if ((int)a.size() <= 1) return; if (k == 1) { mint a1 = a[1]; a[1] = a[0] - a[1]; a[0] = a[0] + a1; return; } int u = 1 << (k - 2); int v = 1; mint one = mint(1); mint imag = dy[1]; while (u) { // jh = 0 { int j0 = 0; int j1 = v; int j2 = v + v; int j3 = j2 + v; for (; j0 < v; ++j0, ++j1, ++j2, ++j3) { mint t0 = a[j0], t1 = a[j1], t2 = a[j2], t3 = a[j3]; mint t0p1 = t0 + t1, t2p3 = t2 + t3; mint t0m1 = t0 - t1, t2m3 = (t2 - t3) * imag; a[j0] = t0p1 + t2p3, a[j2] = t0p1 - t2p3; a[j1] = t0m1 + t2m3, a[j3] = t0m1 - t2m3; } } // jh >= 1 mint ww = one, xx = one * dy[2], yy = one; u <<= 2; for (int jh = 4; jh < u;) { ww = xx * xx, yy = xx * imag; int j0 = jh * v; int je = j0 + v; int j2 = je + v; for (; j0 < je; ++j0, ++j2) { mint t0 = a[j0], t1 = a[j0 + v], t2 = a[j2], t3 = a[j2 + v]; mint t0p1 = t0 + t1, t2p3 = t2 + t3; mint t0m1 = (t0 - t1) * xx, t2m3 = (t2 - t3) * yy; a[j0] = t0p1 + t2p3, a[j2] = (t0p1 - t2p3) * ww; a[j0 + v] = t0m1 + t2m3, a[j2 + v] = (t0m1 - t2m3) * ww; } xx *= dy[__builtin_ctzll(jh += 4)]; } u >>= 4; v <<= 2; } if (k & 1) { u = 1 << (k - 1); for (int j = 0; j < u; ++j) { mint ajv = a[j] - a[j + u]; a[j] += a[j + u]; a[j + u] = ajv; } } } void ntt(vector &a) { if ((int)a.size() <= 1) return; fft4(a, __builtin_ctz(a.size())); } void intt(vector &a) { if ((int)a.size() <= 1) return; ifft4(a, __builtin_ctz(a.size())); mint iv = mint(a.size()).inv(); for (auto &x : a) x *= iv; } vector multiply(const vector &a, const vector &b) { int l = a.size() + b.size() - 1; if (min(a.size(), b.size()) <= 40) { vector s(l); for (int i = 0; i < (int)a.size(); ++i) for (int j = 0; j < (int)b.size(); ++j) s[i + j] += a[i] * b[j]; return s; } int k = 2, M = 4; while (M < l) M <<= 1, ++k; setwy(k); vector s(M); for (int i = 0; i < (int)a.size(); ++i) s[i] = a[i]; fft4(s, k); if (a.size() == b.size() && a == b) { for (int i = 0; i < M; ++i) s[i] *= s[i]; } else { vector t(M); for (int i = 0; i < (int)b.size(); ++i) t[i] = b[i]; fft4(t, k); for (int i = 0; i < M; ++i) s[i] *= t[i]; } ifft4(s, k); s.resize(l); mint invm = mint(M).inv(); for (int i = 0; i < l; ++i) s[i] *= invm; return s; } void ntt_doubling(vector &a) { int M = (int)a.size(); auto b = a; intt(b); mint r = 1, zeta = mint(pr).pow((mint::mod() - 1) / (M << 1)); for (int i = 0; i < M; i++) b[i] *= r, r *= zeta; ntt(b); copy(begin(b), end(b), back_inserter(a)); } }; template struct FormalPowerSeries : vector { using vector::vector; using FPS = FormalPowerSeries; FPS &operator+=(const FPS &r) { if (r.size() > this->size()) this->resize(r.size()); for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i]; return *this; } FPS &operator+=(const mint &r) { if (this->empty()) this->resize(1); (*this)[0] += r; return *this; } FPS &operator-=(const FPS &r) { if (r.size() > this->size()) this->resize(r.size()); for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i]; return *this; } FPS &operator-=(const mint &r) { if (this->empty()) this->resize(1); (*this)[0] -= r; return *this; } FPS &operator*=(const mint &v) { for (int k = 0; k < (int)this->size(); k++) (*this)[k] *= v; return *this; } FPS &operator/=(const FPS &r) { if (this->size() < r.size()) { this->clear(); return *this; } int n = this->size() - r.size() + 1; if ((int)r.size() <= 64) { FPS f(*this), g(r); g.shrink(); mint coeff = g.back().inv(); for (auto &x : g) x *= coeff; int deg = (int)f.size() - (int)g.size() + 1; int gs = g.size(); FPS quo(deg); for (int i = deg - 1; i >= 0; i--) { quo[i] = f[i + gs - 1]; for (int j = 0; j < gs; j++) f[i + j] -= quo[i] * g[j]; } *this = quo * coeff; this->resize(n, mint(0)); return *this; } return *this = ((*this).rev().pre(n) * r.rev().inv(n)).pre(n).rev(); } FPS &operator%=(const FPS &r) { *this -= *this / r * r; shrink(); return *this; } FPS operator+(const FPS &r) const { return FPS(*this) += r; } FPS operator+(const mint &v) const { return FPS(*this) += v; } FPS operator-(const FPS &r) const { return FPS(*this) -= r; } FPS operator-(const mint &v) const { return FPS(*this) -= v; } FPS operator*(const FPS &r) const { return FPS(*this) *= r; } FPS operator*(const mint &v) const { return FPS(*this) *= v; } FPS operator/(const FPS &r) const { return FPS(*this) /= r; } FPS operator%(const FPS &r) const { return FPS(*this) %= r; } FPS operator-() const { FPS ret(this->size()); for (int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i]; return ret; } void shrink() { while (this->size() && this->back() == mint(0)) this->pop_back(); } FPS rev() const { FPS ret(*this); reverse(begin(ret), end(ret)); return ret; } FPS dot(FPS r) const { FPS ret(min(this->size(), r.size())); for (int i = 0; i < (int)ret.size(); i++) ret[i] = (*this)[i] * r[i]; return ret; } // 前 sz 項を取ってくる。sz に足りない項は 0 埋めする FPS pre(int sz) const { FPS ret(begin(*this), begin(*this) + min((int)this->size(), sz)); if ((int)ret.size() < sz) ret.resize(sz); return ret; } FPS operator>>(int sz) const { if ((int)this->size() <= sz) return {}; FPS ret(*this); ret.erase(ret.begin(), ret.begin() + sz); return ret; } FPS operator<<(int sz) const { FPS ret(*this); ret.insert(ret.begin(), sz, mint(0)); return ret; } FPS diff() const { const int n = (int)this->size(); FPS ret(max(0, n - 1)); mint one(1), coeff(1); for (int i = 1; i < n; i++) { ret[i - 1] = (*this)[i] * coeff; coeff += one; } return ret; } FPS integral() const { const int n = (int)this->size(); FPS ret(n + 1); ret[0] = mint(0); if (n > 0) ret[1] = mint(1); auto mod = mint::mod(); for (int i = 2; i <= n; i++) ret[i] = (-ret[mod % i]) * (mod / i); for (int i = 0; i < n; i++) ret[i + 1] *= (*this)[i]; return ret; } mint eval(mint x) const { mint r = 0, w = 1; for (auto &v : *this) r += w * v, w *= x; return r; } FPS log(int deg = -1) const { assert(!(*this).empty() && (*this)[0] == mint(1)); if (deg == -1) deg = (int)this->size(); return (this->diff() * this->inv(deg)).pre(deg - 1).integral(); } FPS pow(int64_t k, int deg = -1) const { const int n = (int)this->size(); if (deg == -1) deg = n; if (k == 0) { FPS ret(deg); if (deg) ret[0] = 1; return ret; } for (int i = 0; i < n; i++) { if ((*this)[i] != mint(0)) { mint rev = mint(1) / (*this)[i]; FPS ret = (((*this * rev) >> i).log(deg) * k).exp(deg); ret *= (*this)[i].pow(k); ret = (ret << (i * k)).pre(deg); if ((int)ret.size() < deg) ret.resize(deg, mint(0)); return ret; } if (__int128_t(i + 1) * k >= deg) return FPS(deg, mint(0)); } return FPS(deg, mint(0)); } static void *ntt_ptr; static void set_fft(); FPS &operator*=(const FPS &r); void ntt(); void intt(); void ntt_doubling(); static int ntt_pr(); FPS inv(int deg = -1) const; FPS exp(int deg = -1) const; }; template void *FormalPowerSeries::ntt_ptr = nullptr; template void FormalPowerSeries::set_fft() { if (!ntt_ptr) ntt_ptr = new NTT; } template FormalPowerSeries& FormalPowerSeries::operator*=( const FormalPowerSeries& r) { if (this->empty() || r.empty()) { this->clear(); return *this; } set_fft(); auto ret = static_cast*>(ntt_ptr)->multiply(*this, r); return *this = FormalPowerSeries(ret.begin(), ret.end()); } template void FormalPowerSeries::ntt() { set_fft(); static_cast*>(ntt_ptr)->ntt(*this); } template void FormalPowerSeries::intt() { set_fft(); static_cast*>(ntt_ptr)->intt(*this); } template void FormalPowerSeries::ntt_doubling() { set_fft(); static_cast*>(ntt_ptr)->ntt_doubling(*this); } template int FormalPowerSeries::ntt_pr() { set_fft(); return static_cast*>(ntt_ptr)->pr; } template FormalPowerSeries FormalPowerSeries::inv(int deg) const { assert((*this)[0] != mint(0)); if (deg == -1) deg = (int)this->size(); FormalPowerSeries res(deg); res[0] = {mint(1) / (*this)[0]}; for (int d = 1; d < deg; d <<= 1) { FormalPowerSeries f(2 * d), g(2 * d); for (int j = 0; j < min((int)this->size(), 2 * d); j++) f[j] = (*this)[j]; for (int j = 0; j < d; j++) g[j] = res[j]; f.ntt(); g.ntt(); for (int j = 0; j < 2 * d; j++) f[j] *= g[j]; f.intt(); for (int j = 0; j < d; j++) f[j] = 0; f.ntt(); for (int j = 0; j < 2 * d; j++) f[j] *= g[j]; f.intt(); for (int j = d; j < min(2 * d, deg); j++) res[j] = -f[j]; } return res.pre(deg); } template FormalPowerSeries FormalPowerSeries::exp(int deg) const { using fps = FormalPowerSeries; assert((*this).size() == 0 || (*this)[0] == mint(0)); if (deg == -1) deg = this->size(); fps inv; inv.reserve(deg + 1); inv.push_back(mint(0)); inv.push_back(mint(1)); auto inplace_integral = [&](fps& F) -> void { const int n = (int)F.size(); auto mod = mint::mod(); while ((int)inv.size() <= n) { int i = inv.size(); inv.push_back((-inv[mod % i]) * (mod / i)); } F.insert(begin(F), mint(0)); for (int i = 1; i <= n; i++) F[i] *= inv[i]; }; auto inplace_diff = [](fps& F) -> void { if (F.empty()) return; F.erase(begin(F)); mint coeff = 1, one = 1; for (int i = 0; i < (int)F.size(); i++) { F[i] *= coeff; coeff += one; } }; fps b{1, 1 < (int)this->size() ? (*this)[1] : 0}, c{1}, z1, z2{1, 1}; for (int m = 2; m < deg; m *= 2) { auto y = b; y.resize(2 * m); y.ntt(); z1 = z2; fps z(m); for (int i = 0; i < m; ++i) z[i] = y[i] * z1[i]; z.intt(); fill(begin(z), begin(z) + m / 2, mint(0)); z.ntt(); for (int i = 0; i < m; ++i) z[i] *= -z1[i]; z.intt(); c.insert(end(c), begin(z) + m / 2, end(z)); z2 = c; z2.resize(2 * m); z2.ntt(); fps x(begin(*this), begin(*this) + min(this->size(), m)); x.resize(m); inplace_diff(x); x.push_back(mint(0)); x.ntt(); for (int i = 0; i < m; ++i) x[i] *= y[i]; x.intt(); x -= b.diff(); x.resize(2 * m); for (int i = 0; i < m - 1; ++i) x[m + i] = x[i], x[i] = mint(0); x.ntt(); for (int i = 0; i < 2 * m; ++i) x[i] *= z2[i]; x.intt(); x.pop_back(); inplace_integral(x); for (int i = m; i < min(this->size(), 2 * m); ++i) x[i] += (*this)[i]; fill(begin(x), begin(x) + m, mint(0)); x.ntt(); for (int i = 0; i < 2 * m; ++i) x[i] *= y[i]; x.intt(); b.insert(end(b), begin(x) + m, end(x)); } return fps{begin(b), begin(b) + deg}; } // [x^n] f(x)^i g(x) を i=0,1,...,m で列挙 // n = (f の次数) - 1 template FormalPowerSeries pow_enumerate(FormalPowerSeries f, FormalPowerSeries g = {1}, int m = -1) { using fps = FormalPowerSeries; int n = f.size() - 1, k = 1; g.resize(n + 1); if (m == -1) m = n; int h = 1; while (h < n + 1) h *= 2; fps P((n + 1) * k), Q((n + 1) * k), nP, nQ, buf, buf2; for (int i = 0; i <= n; i++) P[i * k + 0] = g[i]; for (int i = 0; i <= n; i++) Q[i * k + 0] = -f[i]; Q[0] += 1; while (n) { mint inv2 = mint{2}.inv(); mint w = mint{fps::ntt_pr()}.pow((mint::mod() - 1) / (2 * k)); mint iw = w.inv(); buf2.resize(k); auto ntt_doubling = [&]() { copy(begin(buf), end(buf), begin(buf2)); buf2.intt(); mint c = 1; for (int i = 0; i < k; i++) buf2[i] *= c, c *= w; buf2.ntt(); copy(begin(buf2), end(buf2), back_inserter(buf)); }; nP.clear(), nQ.clear(); for (int i = 0; i <= n; i++) { buf.resize(k); copy(begin(P) + i * k, begin(P) + (i + 1) * k, begin(buf)); ntt_doubling(); copy(begin(buf), end(buf), back_inserter(nP)); buf.resize(k); copy(begin(Q) + i * k, begin(Q) + (i + 1) * k, begin(buf)); if (i == 0) { for (int j = 0; j < k; j++) buf[j] -= 1; ntt_doubling(); for (int j = 0; j < k; j++) buf[j] += 1; for (int j = 0; j < k; j++) buf[k + j] -= 1; } else { ntt_doubling(); } copy(begin(buf), end(buf), back_inserter(nQ)); } nP.resize(2 * h * 2 * k); nQ.resize(2 * h * 2 * k); fps p(2 * h), q(2 * h); w = mint{fps::ntt_pr()}.pow((mint::mod() - 1) / (2 * h)); iw = w.inv(); vector btr; if (n % 2) { btr.resize(h); for (int i = 0, lg = __builtin_ctz(h); i < h; i++) { btr[i] = (btr[i >> 1] >> 1) + ((i & 1) << (lg - 1)); } } for (int j = 0; j < 2 * k; j++) { p.assign(2 * h, 0); q.assign(2 * h, 0); for (int i = 0; i < h; i++) { p[i] = nP[i * 2 * k + j], q[i] = nQ[i * 2 * k + j]; } p.ntt(), q.ntt(); for (int i = 0; i < 2 * h; i += 2) swap(q[i], q[i + 1]); for (int i = 0; i < 2 * h; i++) p[i] *= q[i]; for (int i = 0; i < h; i++) q[i] = q[i * 2] * q[i * 2 + 1]; if (n % 2 == 0) { for (int i = 0; i < h; i++) p[i] = (p[i * 2] + p[i * 2 + 1]) * inv2; } else { mint c = inv2; buf.resize(h); for (int i : btr) buf[i] = (p[i * 2] - p[i * 2 + 1]) * c, c *= iw; swap(p, buf); } p.resize(h), q.resize(h); p.intt(), q.intt(); for (int i = 0; i < h; i++) nP[i * 2 * k + j] = p[i]; for (int i = 0; i < h; i++) nQ[i * 2 * k + j] = q[i]; } nP.resize((n / 2 + 1) * 2 * k); nQ.resize((n / 2 + 1) * 2 * k); swap(P, nP), swap(Q, nQ); n /= 2, h /= 2, k *= 2; } fps S{begin(P), begin(P) + k}; fps T{begin(Q), begin(Q) + k}; S.intt(), T.intt(), T[0] -= 1; if (f[0] == 0) return S.rev().pre(m + 1); return (S.rev() * (T + (fps{1} << k)).rev().inv(m + 1)).pre(m + 1); } // g(f(x)) を計算 template FormalPowerSeries composition(FormalPowerSeries f, FormalPowerSeries g, int deg = -1) { using fps = FormalPowerSeries; auto dfs = [&](auto rc, fps Q, int n, int h, int k) -> fps { if (n == 0) { fps T{begin(Q), begin(Q) + k}; T.push_back(1); fps u = g * T.rev().inv().rev(); fps P(h * k); for (int i = 0; i < (int)g.size(); i++) P[k - 1 - i] = u[i + k]; return P; } fps nQ(4 * h * k), nR(2 * h * k); for (int i = 0; i < k; i++) { copy(begin(Q) + i * h, begin(Q) + i * h + n + 1, begin(nQ) + i * 2 * h); } nQ[k * 2 * h] += 1; nQ.ntt(); for (int i = 0; i < 4 * h * k; i += 2) swap(nQ[i], nQ[i + 1]); for (int i = 0; i < 2 * h * k; i++) nR[i] = nQ[i * 2] * nQ[i * 2 + 1]; nR.intt(); nR[0] -= 1; Q.assign(h * k, 0); for (int i = 0; i < 2 * k; i++) { for (int j = 0; j <= n / 2; j++) { Q[i * h / 2 + j] = nR[i * h + j]; } } auto P = rc(rc, Q, n / 2, h / 2, k * 2); fps nP(4 * h * k); for (int i = 0; i < 2 * k; i++) { for (int j = 0; j <= n / 2; j++) { nP[i * 2 * h + j * 2 + n % 2] = P[i * h / 2 + j]; } } nP.ntt(); for (int i = 1; i < 4 * h * k; i *= 2) { reverse(begin(nQ) + i, begin(nQ) + i * 2); } for (int i = 0; i < 4 * h * k; i++) nP[i] *= nQ[i]; nP.intt(); P.assign(h * k, 0); for (int i = 0; i < k; i++) { copy(begin(nP) + i * 2 * h, begin(nP) + i * 2 * h + n + 1, begin(P) + i * h); } return P; }; if (deg == -1) deg = max(f.size(), g.size()); f.resize(deg), g.resize(deg); int n = f.size() - 1, k = 1; int h = 1; while (h < n + 1) h *= 2; fps Q(h * k); for (int i = 0; i <= n; i++) Q[i] = -f[i]; fps P = dfs(dfs, Q, n, h, k); return P.pre(n + 1).rev(); } // f を入力として, f(g(x)) = x を満たす g(x) mod x^{deg} を返す template FormalPowerSeries compositional_inverse(FormalPowerSeries f, int deg = -1) { assert(int(f.size()) == deg); using fps = FormalPowerSeries; assert((int)f.size() >= 2 and f[1] != 0); if (deg == -1) deg = f.size(); if (deg < 2) return fps{0, f[1].inv()}.pre(deg); int n = deg - 1; fps h = pow_enumerate(f) * n; for (int k = 1; k <= n; k++) h[k] /= k; h = h.rev(); h *= h[0].inv(); fps g = (h.log() * mint{-n}.inv()).exp(); g *= f[1].inv(); return (g << 1).pre(deg); } // f(g(x)) = x を満たす g(x) mod x^{deg} を返す // calc_f(g, d) は f(g(x)) mod x^d を計算する関数 template fps compositional_inverse(function calc_f, int deg) { if (deg <= 2) { fps g = calc_f(fps{0, 1}, 2); assert(g[0] == 0 && g[1] != 0); g[1] = g[1].inv(); return g.pre(deg); } fps g = compositional_inverse(calc_f, (deg + 1) / 2); fps fg = calc_f(g, deg + 1); fps fdg = (fg.diff() * g.diff().inv(deg)).pre(deg); return (g - (fg - fps{0, 1}) * fdg.inv()).pre(deg); } // divide and conquer template struct Merge{ int n; using P = fps; using Comp = std::function; Comp comp = [](const P& a, const P& b){return a.size() > b.size();}; priority_queue, Comp> pq; Merge(int n = -1) : n(n), pq(comp){ pq.push(P{1}); } void push(P r){ pq.push(r); } P get(){ while(pq.size() > 1){ auto f = pq.top(); pq.pop(); auto g = pq.top(); pq.pop(); f *= g; if(n != -1) if(int(f.size()) > n) f.resize(n); pq.push(f); } P res = pq.top(); res.resize(n); return res; } }; using fps = FormalPowerSeries; int main(){ int n, k; cin >> n >> k; assert(1 <= k && k <= n && n <= 200000); auto pentagonal = [&](int x){ return x * (3 * x - 1) / 2; }; fps f(n + 1); { int i = 0; mint c = 1; while(pentagonal(i) <= n){ f[pentagonal(i)] += c; c *= -1; i++; } } { int i = -1; mint c = -1; while(pentagonal(i) <= n){ f[pentagonal(i)] += c; c *= -1; i--; } } fps g(n + 1); for(int i = 0; i * (k + 1) <= n; i++) g[i * (k + 1)] += f[i]; fps h = g *= f.inv(); for(int i = 1; i <= n; i++){ cout << h[i].val(); if(i == n) cout << "\n"; else cout << " "; } return 0; }