#include using namespace std; using ll = long long; static const int MOD = 998244353; static const int G = 3; ll modpow(ll a, ll e) { ll r = 1; while (e > 0) { if (e & 1) r = r * a % MOD; a = a * a % MOD; e >>= 1; } return r; } void ntt(vector &a, bool inv) { int n = (int)a.size(); for (int i = 1, j = 0; i < n; i++) { int bit = n >> 1; while (j & bit) { j ^= bit; bit >>= 1; } j ^= bit; if (i < j) swap(a[i], a[j]); } for (int len = 2; len <= n; len <<= 1) { int wlen = (int)modpow(G, (MOD - 1) / len); if (inv) wlen = (int)modpow(wlen, MOD - 2); for (int i = 0; i < n; i += len) { ll w = 1; int half = len >> 1; for (int j = 0; j < half; j++) { int u = a[i + j]; int v = (int)(w * a[i + j + half] % MOD); int x = u + v; if (x >= MOD) x -= MOD; a[i + j] = x; x = u - v; if (x < 0) x += MOD; a[i + j + half] = x; w = w * wlen % MOD; } } } if (inv) { int inv_n = (int)modpow(n, MOD - 2); for (int &x : a) x = (ll)x * inv_n % MOD; } } vector convolution_self(vector a, vector b) { int need = (int)a.size() + (int)b.size() - 1; int n = 1; while (n < need) n <<= 1; a.resize(n); b.resize(n); ntt(a, false); ntt(b, false); for (int i = 0; i < n; i++) { a[i] = (ll)a[i] * b[i] % MOD; } ntt(a, true); a.resize(need); return a; } vector multiply_poly_full(const vector &a, const vector &b, int E) { vector c = convolution_self(a, b); c.resize(E + 1, 0); return c; } vector shift_poly_full(const vector &p, int shift, int E) { vector res(E + 1, 0); if (shift > E) return res; for (int i = 0; i + shift <= E; i++) { res[i + shift] = p[i]; } return res; } void add_scaled_full(vector &dst, const vector &src, int scale, int E) { if (scale == 0) return; for (int i = 0; i <= E; i++) { if (src[i] == 0) continue; dst[i] = (dst[i] + (ll)src[i] * scale) % MOD; } } struct Comb { vector fact, ifact; void init(int n) { fact.assign(n + 1, 1); ifact.assign(n + 1, 1); for (int i = 1; i <= n; i++) fact[i] = (ll)fact[i - 1] * i % MOD; ifact[n] = (int)modpow(fact[n], MOD - 2); for (int i = n; i >= 1; i--) ifact[i - 1] = (ll)ifact[i] * i % MOD; } int C(int n, int r) const { if (r < 0 || r > n) return 0; return (ll)fact[n] * ifact[r] % MOD * ifact[n - r] % MOD; } }; int main() { ios::sync_with_stdio(false); cin.tie(nullptr); int N, M, K; cin >> N >> M >> K; int E = N * (N - 1) / 2; int R = E - M; int S = N - 2; Comb comb; comb.init(E); // base_poly[s] = q^s - 1 vector> base_poly(S + 2, vector(E + 1, 0)); for (int s = 0; s <= S + 1; s++) { base_poly[s][0] = MOD - 1; if (s <= E) { base_poly[s][s]++; if (base_poly[s][s] >= MOD) base_poly[s][s] -= MOD; } } // pow_poly[s][t] = (q^s - 1)^t vector>> pow_poly( S + 2, vector>(S + 1, vector(E + 1, 0)) ); for (int s = 0; s <= S + 1; s++) { pow_poly[s][0][0] = 1; for (int t = 1; t <= S; t++) { pow_poly[s][t] = multiply_poly_full(pow_poly[s][t - 1], base_poly[s], E); } } // kernel[s][t] = (q^s - 1)^t * q^{t(t-1)/2} vector>> kernel( S + 2, vector>(S + 1, vector(E + 1, 0)) ); for (int s = 0; s <= S + 1; s++) { for (int t = 0; t <= S; t++) { int shift = t * (t - 1) / 2; kernel[s][t] = shift_poly_full(pow_poly[s][t], shift, E); } } // dp[r][s] is a polynomial in q. vector>> dp( S + 1, vector>(S + 2, vector(E + 1, 0)) ); vector> active(S + 1, vector(S + 2, 0)); dp[S][1][0] = 1; active[S][1] = 1; // Build layers 1,2,...,K-1 without reaching vertex N. for (int step = 0; step < K - 1; step++) { vector>> ndp( S + 1, vector>(S + 2, vector(E + 1, 0)) ); vector> nactive(S + 1, vector(S + 2, 0)); for (int r = 0; r <= S; r++) { for (int s = 0; s <= S + 1; s++) { if (!active[r][s]) continue; for (int t = 0; t <= r; t++) { int ways = comb.C(r, t); vector prod = multiply_poly_full(dp[r][s], kernel[s][t], E); add_scaled_full(ndp[r - t][t], prod, ways, E); nactive[r - t][t] = 1; } } } dp.swap(ndp); active.swap(nactive); } // H(q): generating polynomial for dist(1,N)=K in q=1+x. vector H(E + 1, 0); for (int r = 0; r <= S; r++) { for (int s = 0; s <= S + 1; s++) { if (!active[r][s]) continue; int base = s * r + r + r * (r - 1) / 2; // final factor: (q^s - 1) q^base vector final_factor = shift_poly_full(base_poly[s], base, E); vector prod = multiply_poly_full(dp[r][s], final_factor, E); add_scaled_full(H, prod, 1, E); } } // We need [x^R] H(1+x). // If H(q)=sum_c h_c q^c, then [x^R]H(1+x)=sum_c h_c*C(c,R). ll ans = 0; for (int c = R; c <= E; c++) { ans += (ll)H[c] * comb.C(c, R) % MOD; ans %= MOD; } cout << ans << '\n'; return 0; }