#include #include #include #include using namespace std; #include "atcoder/modint.hpp" #include "atcoder/convolution.hpp" template struct Polynomial : vector { using vector::vector; using poly = Polynomial; void shrink() { while(!this->empty() and this->back() == mint(0)) this->pop_back(); } poly rev() const { poly res(*this); reverse(res.begin(), res.end()); return res; } poly pre(int sz) const { poly res(this->begin(), this->begin() + min(this->size(), (size_t)sz)); if (res.size() < (size_t)sz) res.resize(sz); return res; } poly &operator>>=(int sz) { if (this->size() <= (size_t)sz) { this->clear(); return *this; } this->erase(this->begin(), this->begin() + (size_t)sz); return *this; } poly &operator<<=(int sz) { this->insert(this->begin(), sz, mint::raw(0)); return *this; } poly &operator+=(const mint &s) { if (this->empty()) this->resize(1); (*this)[0] += s; return *this; } poly &operator+=(const poly &p) { if (this->size() < p.size()) this->resize(p.size()); for (size_t i = 0; i < p.size(); i++) (*this)[i] += p[i]; return *this; } poly &operator-=(const mint &s) { if (this->empty()) this->resize(1); (*this)[0] -= s; return *this; } poly &operator-=(const poly &p) { if (this->size() < p.size()) this->resize(p.size()); for (size_t i = 0; i < p.size(); i++) (*this)[i] -= p[i]; return *this; } poly &operator*=(const mint &s) { for (size_t i = 0; i < this->size(); i++) (*this)[i] *= s; return *this; } poly operator*=(const poly &p) { if (this->empty() or p.empty()) { this->clear(); return *this; } vector prod = atcoder::convolution(*this, p); return *this = poly(prod.begin(), prod.end()); } poly &operator/=(const mint &s) { mint s_inv = s.inv(); for (size_t i = 0; i < this->size(); i++) (*this)[i] *= s_inv; return *this; } poly operator/=(const poly &p) { if (this->size() < p.size()) { this->clear(); return *this; } size_t q_sz = this->size() - p.size() + 1; return *this = (this->rev().pre(q_sz) * p.rev().inv(q_sz)).pre(q_sz).rev(); } poly operator%=(const poly &p) { *this -= (*this / p) * p; shrink(); assert(this->size() < p.size()); return *this; } poly operator+(const mint &s) const { return poly(*this) += s; } poly operator+(const poly &p) const { return poly(*this) += p; } poly operator-(const mint &s) const { return poly(*this) -= s; } poly operator-(const poly &p) const { return poly(*this) -= p; } poly operator*(const mint &s) const { return poly(*this) *= s; } poly operator*(const poly &p) const { return poly(*this) *= p; } poly operator/(const mint &s) const { return poly(*this) /= s; } poly operator/(const poly &p) const { return poly(*this) /= p; } poly operator%(const poly &p) const { return poly(*this) %= p; } poly operator>>(int s) const { return poly(*this) >>= s; } poly operator<<(int s) const { return poly(*this) <<= s; } poly operator-() const { poly res(this->size()); for (size_t i = 0; i < this->size(); i++) res[i] = -(*this)[i]; return res; } poly dot(const poly &p) const { poly res(min(this->size(), p.size())); for (size_t i = 0; i < res.size(); i++) res[i] = (*this)[i] * p[i]; return res; } mint eval(mint x) const { mint res = 0, pw = 1; for (auto &c : *this) res += c * pw, pw *= x; return res; } void ntt() { atcoder::internal::butterfly(*this); } void intt() { atcoder::internal::butterfly_inv(*this); mint sz_inv = mint::raw(this->size()).inv(); for (auto &c : *this) c *= sz_inv; } poly inv(int deg = -1) const { assert(!this->empty() and (*this)[0] != 0); if (deg == -1) deg = (int)this->size(); poly res(deg); res[0] = (*this)[0].inv(); for (int d = 1; d < deg; d <<= 1) { poly 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]; // atcoder::internal::butterfly(f); f.ntt(); // atcoder::internal::butterfly(g); g.ntt(); for (int j = 0; j < 2 * d; j++) f[j] *= g[j]; // atcoder::internal::butterfly_inv(f); f.intt(); for (int j = 0; j < d; j++) f[j] = 0; // atcoder::internal::butterfly(f); f.ntt(); for (int j = 0; j < 2 * d; j++) f[j] *= g[j]; // atcoder::internal::butterfly_inv(f); f.intt(); for (int j = d; j < min(2 * d, deg); j++) res[j] = -f[j]; } return res.pre(deg); } poly differentiate() const { const size_t n = this->size(); poly res(n == 0 ? 0 : n - 1); mint coef = 1; for (size_t i = 1; i < n; i++) { res[i - 1] = (*this)[i] * coef; coef += 1; } return res; } poly integrate() const { const size_t n = this->size(); poly res(n + 1); res[0] = mint::raw(0); if (n > 0) res[1] = mint(1); int mod = mint::mod(); for (size_t i = 2; i <= n; i++) res[i] = (-res[mod % i]) * (mod / i); for (size_t i = 0; i < n; i++) res[i + 1] *= (*this)[i]; return res; } poly log(int deg = -1) const { assert(!this->empty() and (*this)[0] == 1); if (deg == -1) deg = (int)this->size(); return (this->differentiate() * this->inv(deg)).pre(deg - 1).integrate(); } poly exp(int deg = -1) const { assert(this->empty() or (*this)[0] == mint(0)); if (deg == -1) deg = (int)this->size(); vector inv; inv.reserve(deg + 1); inv.emplace_back(mint::raw(0)); inv.emplace_back(mint::raw(1)); poly b{1, 1 < (int)this->size() ? (*this)[1] : 0}; poly c{1}, z1, z2{1, 1}; int mod = mint::mod(); for (int d = 2; d < deg ; d <<= 1) { poly y = b; y.resize(2 * d); // atcoder::internal::butterfly(y); y.ntt(); z1 = z2; poly z(d); for (int i = 0; i < d; i++) z[i] = y[i] * z1[i]; // atcoder::internal::butterfly_inv(z); z.intt(); fill(z.begin(), z.begin() + d / 2, mint::raw(0)); // atcoder::internal::butterfly(z); z.ntt(); for (int i = 0; i < d; i++) z[i] *= -z1[i]; // atcoder::internal::butterfly_inv(z); z.intt(); c.insert(c.end(), z.begin() + d / 2, z.end()); z2 = c; z2.resize(2 * d); // atcoder::internal::butterfly(z2); z2.ntt(); poly x(this->begin(), this->begin() + min(this->size(), d)); x.resize(d); { x.erase(x.begin()); mint coef = 1; for (int i = 0; i < d-1; i++) x[i] *= coef, coef += 1; } x.emplace_back(mint::raw(0)); // atcoder::internal::butterfly(x); x.ntt(); for (int i = 0; i < d; i++) x[i] *= y[i]; // atcoder::internal::butterfly_inv(x); x.intt(); x -= b.differentiate(); x.resize(2 * d); for (int i = 0; i < d - 1; i++) x[d + i] = x[i], x[i] = mint::raw(0); // atcoder::internal::butterfly(x); x.ntt(); for (int i = 0; i < 2 * d; i++) x[i] *= z2[i]; // atcoder::internal::butterfly_inv(x); x.intt(); x.pop_back(); { int sz; while((sz = (int)inv.size()) <= (int)x.size()) { inv.emplace_back((-inv[mod % sz]) * (mod / sz)); } x.insert(x.begin(), mint::raw(0)); for (size_t i = 1; i < x.size(); i++) x[i] *= inv[i]; } for (int i = d; i < min((int)this->size(), 2 * d); i++) x[i] += (*this)[i]; fill(x.begin(), x.begin() + d, mint::raw(0)); // atcoder::internal::butterfly(x); x.ntt(); for (int i = 0; i < 2 * d; i++) x[i] *= y[i]; // atcoder::internal::butterfly_inv(x); x.intt(); b.insert(b.end(), x.begin() + d, x.end()); } return b.pre(deg); } poly pow(long long k, int deg = -1) const { const int n = (int)this->size(); if (deg == -1) deg = n; if (k == 0) { poly res(deg); if (deg > 0) res[0] = 1; return res; } int l = 0; while (l < n and (*this)[l] == 0){ l++; if (l >= (deg + k - 1) / k) return poly(deg, mint::raw(0)); } if (l == n) return poly(deg, mint::raw(0)); mint rev = (*this)[l].inv(); poly res = (((*this * rev) >> l).log(deg) * k).exp(deg); res *= (*this)[l].pow(k); res = (res << (k * l)).pre(deg); if ((int)res.size() < deg) res.resize(deg, mint::raw(0)); return res; } }; template struct Binomial { vector fact_, inv_, finv_; constexpr Binomial() {} constexpr Binomial(int n) noexcept : fact_(n, 1), inv_(n, 1), finv_(n, 1) { init(n); } constexpr void init(int n) noexcept { constexpr int mod = T::mod(); fact_.assign(n, 1), inv_.assign(n, 1), finv_.assign(n, 1); for(int i = 2; i < n; i++){ fact_[i] = fact_[i-1] * i; inv_[i] = -inv_[mod%i] * (mod/i); finv_[i] = finv_[i-1] * inv_[i]; } } constexpr T com(int n, int k) const noexcept { if (n < k || n < 0 || k < 0) return 0; return fact_[n] * finv_[k] * finv_[n-k]; } constexpr T perm(int n, int k) const noexcept { if (n < k || n < 0 || k < 0) return 0; return fact_[n] * finv_[n-k]; } constexpr T fact(int n) const noexcept { if (n < 0) return 0; return fact_[n]; } constexpr T inv(int n) const noexcept { if (n < 0) return 0; return inv_[n]; } constexpr T finv(int n) const noexcept { if (n < 0) return 0; return finv_[n]; } constexpr T com_naive(int n, int k) const noexcept { if (n < 0 || k < 0 || n < k) return 0; T res = T(1); k = min(k, n-k); for (int i = 1; i <= k; i++)res *= (n--) * inv(i); return res; } template constexpr T multi(const vector &v) const noexcept { static_assert(is_integral::value); I n = 0; for (auto& x : v) { if (x < 0) return 0; n += x; } T res = fact(n); for (auto &x : v) res *= finv(x); return res; } // [x^k] (1-x)^{-n} = com(n+k-1, k) constexpr T neg(int n, int k) const noexcept { if (n < 0 || k < 0) return 0; return k == 0 ? 1 : com(n+k-1, k); } }; using mint = atcoder::modint998244353; using poly = Polynomial; Binomial bc(300000); // deg(f) = n-1 // [x^{n-1}] g(x)f(x)^i for i=0,1,...,m template Polynomial power_projection(Polynomial f, Polynomial g = {1}, int m = -1) { using poly = Polynomial; if (f.empty()) return poly(m+1,0); if (f[0] != 0) { mint c = f[0]; f[0] = 0; poly a = power_projection(f, g, m); for (int i = 0; i <= m; i++) a[i] *= bc.finv(i); poly b(m+1); mint pc = 1; for (int i = 0; i <= m; i++) b[i] = bc.finv(i) * pc, pc *= c; a *= b; a.resize(m+1); for (int i = 0; i <= m; i++) a[i] *= bc.fact(i); return a; } int fn = f.size(); if(m == -1) m = fn - 1; g.resize(fn); int n = 1; while (n < fn) n *= 2; f.resize(n); reverse(g.begin(), g.end()); g.resize(n); reverse(g.begin(), g.end()); int k = 1; poly p(n * 2), q(n * 2), r(n * 2); for (int i = 0; i < n; i++) p[i] = g[i], q[i] = -f[i]; while (n > 1) { for (int i = 0; i < 2 * n * k; i++) r[i] = (i % 2 == 0 ? q[i] : -q[i]); poly pq = p * r; poly qq = q * r; pq.resize(4 * n * k); qq.resize(4 * n * k); for (int i = 0; i < 2 * n * k; i++) { pq[2 * n * k + i] += p[i]; qq[2 * n * k + i] += q[i] + r[i]; } fill(p.begin(), p.end(), 0); fill(q.begin(), q.end(), 0); for (int i = 0; i < 2 * k; i++) { for (int j = 0; j < n / 2; j++) { p[n * i + j] = pq[2 * n * i + 2 * j + 1]; q[n * i + j] = qq[2 * n * i + 2 * j + 0]; } } n /= 2; k *= 2; } poly pk(k),qk(k+1); for (int i = 0; i < k; i++) pk[i] = p[2 * i]; for (int i = 0; i < k; i++) qk[i] = q[2 * i]; qk[k] = 1; reverse(pk.begin(), pk.end()); reverse(qk.begin(), qk.end()); return (pk * qk.inv(m+1)).pre(m+1); } // g(f(x)) template Polynomial composition(Polynomial f, Polynomial g, int deg = -1) { using poly = Polynomial; auto middle_product = [&] (poly x, poly y) -> poly { int nm_1 = x.size(), m = y.size(); int l = bit_ceil(x.size()); reverse(y.begin(), y.end()); x.resize(l); y.resize(l); x.ntt(); y.ntt(); for (int i = 0; i < l; i++) x[i] *= y[i]; x.intt(); return poly{x.begin() + m - 1, x.begin() + nm_1}; }; auto rec = [&] (this auto self, int n, int k, poly q) -> poly { if (n == 1) { reverse(g.begin(), g.end()); poly p(2 * k); for (int i = 0; i < k; i++) p[2 * i] = g[i]; return p; } poly r(2 * n * k); for (int i = 0; i < 2 * n * k; i++) r[i] = (i % 2 == 0 ? q[i] : -q[i]); poly qq = q * r; qq.resize(4 * n * k); for (int i = 0; i < 2 * n * k; i++) qq[2 * n * k + i] += q[i] + r[i]; poly nq(2 * n * k); for (int i = 0; i < k * 2; i++) { for (int j = 0; j < n / 2; j++) { nq[n * i + j] = qq[2 * n * i + 2 * j]; } } poly np = self(n / 2, k * 2, nq); poly pq(4 * n * k); for (int i = 0; i < k * 2; i++) { for (int j = 0; j < n / 2; j++) { pq[2 * n * i + 2 * j + 1] += np[n * i + j]; } } poly p(2 * n * k); for (int i = 0; i < 2 * n * k; i++) p[i] += pq[2 * n * k + i]; pq.pop_back(); poly mp = middle_product(pq, r); for (int i = 0; i < 2 * n * k; i++) p[i] += mp[i]; return p; }; if (deg == -1) deg = max(f.size(), g.size()); int n = 1; while (n < deg) n *= 2; f.resize(n); g.resize(n); poly q(n * 2); for (int i = 0; i < n; i++) q[i] = -f[i]; poly p = rec(n, 1, q); return p.pre(n).rev().pre(deg); } template Polynomial power_projection_ntt(Polynomial f, Polynomial g = {1}, int m = -1) { using poly = Polynomial; if (f.empty()) return poly(m+1,0); if (f[0] != 0) { mint c = f[0]; f[0] = 0; poly a = power_projection(f, g, m); for (int i = 0; i <= m; i++) a[i] *= bc.finv(i); poly b(m+1); mint pc = 1; for (int i = 0; i <= m; i++) b[i] = bc.finv(i) * pc, pc *= c; a *= b; a.resize(m+1); for (int i = 0; i <= m; i++) a[i] *= bc.fact(i); return a; } int fn = f.size(); if(m == -1) m = fn - 1; g.resize(fn); int n = 1; while (n < fn) n *= 2; f.resize(n); reverse(g.begin(), g.end()); g.resize(n); reverse(g.begin(), g.end()); vector br(2 * n); const int log = 31 - countl_zero(static_cast(2 * n)); for (int i = 0; i < 2 * n; i++) br[i] = (br[i >> 1] >> 1) + ((i & 1) << (log - 1)); constexpr int mod = mint::mod(); constexpr int inv2 = atcoder::internal::pow_mod_constexpr(2, mod - 2, mod); constexpr int pr = atcoder::internal::primitive_root_constexpr(mod); constexpr int trz = countr_zero(static_cast(mod - 1)); constexpr int rt = atcoder::internal::pow_mod_constexpr(pr, (mod - 1) >> trz, mod); constexpr int inv_rt = atcoder::internal::pow_mod_constexpr(rt, mod - 2, mod); const int w = atcoder::internal::pow_mod_constexpr(inv_rt, (1 << trz) / (4 * n), mod); vector wp(n * 2); mint wi = 1; for (int i : br) wp[i] = wi, wi *= w; int k = 1; poly p(n * 2), q(n * 2); for (int i = 0; i < n; i++) p[i] = g[i], q[i] = -f[i]; while (n > 1) { p.resize(4 * n * k); q.resize(4 * n * k); q[2 * n * k] = 1; p.ntt(); q.ntt(); for (int i = 0; i < 2 * n * k; i++) { p[i] = wp[i] * inv2 * (p[2*i] * q[2*i+1] - p[2*i+1] * q[2*i]); q[i] = q[2*i] * q[2*i+1]; } p.resize(2 * n * k); q.resize(2 * n * k); p.intt(); q.intt(); for (int i = 0; i < 2 * k; i++) { for (int j = n / 2; j < n; j++) { p[n * i + j] = q[n * i + j] = 0; } } q[0] = 0; n /= 2; k *= 2; } poly pk(k),qk(k+1); for (int i = 0; i < k; i++) pk[i] = p[2 * i]; for (int i = 0; i < k; i++) qk[i] = q[2 * i]; qk[k] = 1; reverse(pk.begin(), pk.end()); reverse(qk.begin(), qk.end()); return (pk * qk.inv(m+1)).pre(m+1); } // g(f(x)) template Polynomial composition_ntt(Polynomial f, Polynomial g, int deg = -1) { using poly = Polynomial; if (deg == -1) deg = max(f.size(), g.size()); int n = 1; while (n < deg) n *= 2; f.resize(n); g.resize(n); vector br(2 * n); const int log = countr_zero(static_cast(2 * n)); for (int i = 0; i < 2 * n; i++) br[i] = (br[i >> 1] >> 1) + ((i & 1) << (log - 1)); constexpr int mod = mint::mod(); constexpr int inv2 = atcoder::internal::pow_mod_constexpr(2, mod - 2, mod); constexpr int pr = atcoder::internal::primitive_root_constexpr(mod); constexpr int trz = countr_zero(static_cast(mod - 1)); constexpr int rt = atcoder::internal::pow_mod_constexpr(pr, (mod - 1) >> trz, mod); constexpr int inv_rt = atcoder::internal::pow_mod_constexpr(rt, mod - 2, mod); const int w = atcoder::internal::pow_mod_constexpr(inv_rt, (1 << trz) / (4 * n), mod); vector wp(2 * n); mint wi = 1; for (int i : br) wp[i] = wi, wi *= w; auto transposed_ntt = [&] (poly &p) -> void { int sz = p.size(); p.intt(); reverse(p.begin() + 1, p.end()); for (mint &x : p) x *= sz; }; auto transposed_intt = [&] (poly &p) -> void { mint sz_inv = mint::raw(p.size()).inv(); reverse(p.begin() + 1, p.end()); p.ntt(); for (mint &x : p) x *= sz_inv; }; auto rec = [&] (this auto self, int n, int k, poly q) -> poly { if (n == 1) { reverse(g.begin(), g.end()); poly p(2 * k); for (int i = 0; i < k; i++) p[2 * i] = g[i]; return p; } q.resize(4 * n * k); q[2 * n * k] = 1; q.ntt(); poly nq(2 * n * k); for (int i = 0; i < 2 * n * k; i++) nq[i] = q[2*i] * q[2*i+1]; nq.intt(); for (int i = 0; i < 2 * k; i++) { for (int j = n / 2; j < n; j++) { nq[i * n + j] = 0; } } nq[0] = 0; poly p = self(n / 2, k * 2, nq); for (int i = 0; i < 2 * k; i++) { for (int j = n / 2; j < n; j++) { p[i * n + j] = 0; } } transposed_intt(p); p.resize(4 * n * k); for (int i = 2 * n * k - 1; i >= 0; i--) { p[2 * i + 1] = wp[i] * -inv2 * q[2 * i] * p[i]; p[2 * i] = wp[i] * inv2 * q[2 * i + 1] * p[i]; } transposed_ntt(p); p.resize(2 * n * k); return p; }; poly q(n * 2); for (int i = 0; i < n; i++) q[i] = -f[i]; poly p = rec(n, 1, q); return p.pre(n).rev().pre(deg); } int main() { long n,m,k,x; cin >> n >> m >> k >> x; poly f(n); for (int i = 1; i < n; i++) f[i] = bc.finv(i); f = f.pow(k, n); for (int i = 0; i < n; i++) f[i] *= bc.fact(i) * bc.finv(k); mint ff = 1; for (int i = 0; i < k; i++) ff *= m-1-i; for (int i = 0; i < n; i++) f[i] *= ff; poly c(n); for (int i = 0; i < n; i++) c[i] = mint(i+1).pow(x) - mint(i).pow(x); // poly h = composition(f,c,n); poly h = composition_ntt(f,c,n); vector ans(n); ans[0] = h[0]; for (int i = 1; i < n; i++) ans[i] = ans[i-1] * m + h[i]; for (int i = 0; i < n; i++) ans[i] *= m; for (int i = 0; i < n; i++) cout << ans[i].val() << "\n"; }