#include using namespace std; using lint = long long; using pi = pair; #define sz(v) ((int)(v).size()) #define all(v) (v).begin(), (v).end() const int mod = 998244353; // 1e9 + 7;//993244853; // I don't remember the credit of modint, but it's not mine. // I don't remember the credit of FFT, but it's probably mine. // Polynomial library is due to adamant: // https://github.com/cp-algorithms/cp-algorithms-aux/blob/master/src/polynomial.cpp // To use polynomial sqrt, need to amend sqrt for modint. struct mint { int val; mint() { val = 0; } mint(const lint &v) { val = (-mod <= v && v < mod) ? v : v % mod; if (val < 0) val += mod; } friend ostream &operator<<(ostream &os, const mint &a) { return os << a.val; } friend bool operator==(const mint &a, const mint &b) { return a.val == b.val; } friend bool operator!=(const mint &a, const mint &b) { return !(a == b); } friend bool operator<(const mint &a, const mint &b) { return a.val < b.val; } mint operator-() const { return mint(-val); } mint &operator+=(const mint &m) { if ((val += m.val) >= mod) val -= mod; return *this; } mint &operator-=(const mint &m) { if ((val -= m.val) < 0) val += mod; return *this; } mint &operator*=(const mint &m) { val = (lint)val * m.val % mod; return *this; } friend mint ipow(mint a, lint p) { mint ans = 1; for (; p; p /= 2, a *= a) if (p & 1) ans *= a; return ans; } mint inv() const { return ipow(*this, mod - 2); } mint &operator/=(const mint &m) { return (*this) *= m.inv(); } friend mint operator+(mint a, const mint &b) { return a += b; } friend mint operator-(mint a, const mint &b) { return a -= b; } friend mint operator*(mint a, const mint &b) { return a *= b; } friend mint operator/(mint a, const mint &b) { return a /= b; } operator int64_t() const { return val; } }; namespace fft { using real_t = double; using base = complex; void fft(vector &a, bool inv) { int n = a.size(), j = 0; vector roots(n / 2); for (int i = 1; i < n; i++) { int bit = (n >> 1); while (j >= bit) { j -= bit; bit >>= 1; } j += bit; if (i < j) swap(a[i], a[j]); } real_t ang = 2 * acos(real_t(-1)) / n * (inv ? -1 : 1); for (int i = 0; i < n / 2; i++) { roots[i] = base(cos(ang * i), sin(ang * i)); } /* XOR Convolution : set roots[*] = 1. OR Convolution : It's easy, don't try to use FFT. https://codeforces.com/blog/entry/115438 */ for (int i = 2; i <= n; i <<= 1) { int step = n / i; for (int j = 0; j < n; j += i) { for (int k = 0; k < i / 2; k++) { base u = a[j + k], v = a[j + k + i / 2] * roots[step * k]; a[j + k] = u + v; a[j + k + i / 2] = u - v; } } } if (inv) for (int i = 0; i < n; i++) a[i] /= n; } template void ntt(vector &a, bool inv) { const int prr = 3; // primitive root int n = a.size(), j = 0; vector roots(n / 2); for (int i = 1; i < n; i++) { int bit = (n >> 1); while (j >= bit) { j -= bit; bit >>= 1; } j += bit; if (i < j) swap(a[i], a[j]); } T ang = ipow(T(prr), (mod - 1) / n); if (inv) ang = T(1) / ang; for (int i = 0; i < n / 2; i++) { roots[i] = (i ? (roots[i - 1] * ang) : T(1)); } for (int i = 2; i <= n; i <<= 1) { int step = n / i; for (int j = 0; j < n; j += i) { for (int k = 0; k < i / 2; k++) { T u = a[j + k], v = a[j + k + i / 2] * roots[step * k]; a[j + k] = u + v; a[j + k + i / 2] = u - v; } } } if (inv) { T rev = T(1) / T(n); for (int i = 0; i < n; i++) a[i] *= rev; } } template vector multiply_ntt(vector &v, const vector &w) { vector fv(all(v)), fw(all(w)); int n = 2; while (n < sz(v) + sz(w)) n <<= 1; fv.resize(n); fw.resize(n); ntt(fv, 0); ntt(fw, 0); for (int i = 0; i < n; i++) fv[i] *= fw[i]; ntt(fv, 1); vector ret(n); for (int i = 0; i < n; i++) ret[i] = fv[i]; return ret; } template vector multiply(vector &v, const vector &w) { vector fv(all(v)), fw(all(w)); int n = 2; while (n < sz(v) + sz(w)) n <<= 1; fv.resize(n); fw.resize(n); fft(fv, 0); fft(fw, 0); for (int i = 0; i < n; i++) fv[i] *= fw[i]; fft(fv, 1); vector ret(n); for (int i = 0; i < n; i++) ret[i] = (T)llround(fv[i].real()); return ret; } template vector multiply_mod(vector v, const vector &w) { int n = 2; while (n < sz(v) + sz(w)) n <<= 1; vector v1(n), v2(n), r1(n), r2(n); for (int i = 0; i < v.size(); i++) { v1[i] = base(v[i] >> 15, v[i] & 32767); } for (int i = 0; i < w.size(); i++) { v2[i] = base(w[i] >> 15, w[i] & 32767); } fft(v1, 0); fft(v2, 0); for (int i = 0; i < n; i++) { int j = (i ? (n - i) : i); base ans1 = (v1[i] + conj(v1[j])) * base(0.5, 0); base ans2 = (v1[i] - conj(v1[j])) * base(0, -0.5); base ans3 = (v2[i] + conj(v2[j])) * base(0.5, 0); base ans4 = (v2[i] - conj(v2[j])) * base(0, -0.5); r1[i] = (ans1 * ans3) + (ans1 * ans4) * base(0, 1); r2[i] = (ans2 * ans3) + (ans2 * ans4) * base(0, 1); } fft(r1, 1); fft(r2, 1); vector ret(n); for (int i = 0; i < n; i++) { T av = llround(r1[i].real()); T bv = llround(r1[i].imag()) + llround(r2[i].real()); T cv = llround(r2[i].imag()); av = av << 30; bv = bv << 15; ret[i] = av + bv + cv; } return ret; } template vector multiply_naive(vector v, const vector &w) { if (sz(v) == 0 || sz(w) == 0) return vector(); vector ret(sz(v) + sz(w) - 1); for (int i = 0; i < sz(v); i++) { for (int j = 0; j < sz(w); j++) { ret[i + j] += v[i] * w[j]; } } return ret; } } // namespace fft const int maxn = 1 << 17; const int magic = 250; // threshold for sizes to run the naive algo template T fact(int n) { static T F[maxn]; static bool init = false; if (!init) { F[0] = T(1); for (int i = 1; i < maxn; i++) { F[i] = F[i - 1] * T(i); } init = true; } return F[n]; } template T rfact(int n) { static T F[maxn]; static bool init = false; if (!init) { F[maxn - 1] = T(1) / fact(maxn - 1); for (int i = maxn - 2; i >= 0; i--) { F[i] = F[i + 1] * T(i + 1); } init = true; } return F[n]; } template struct poly { vector a; void normalize() { // get rid of leading zeroes while (!a.empty() && a.back() == T(0)) { a.pop_back(); } } poly() {} poly(T a0) : a{a0} { normalize(); } poly(const vector &t) : a(t) { normalize(); } poly operator-() const { auto t = *this; for (auto &it : t.a) { it = -it; } return t; } poly operator+=(const poly &t) { a.resize(max(a.size(), t.a.size())); for (size_t i = 0; i < t.a.size(); i++) { a[i] += t.a[i]; } normalize(); return *this; } poly operator-=(const poly &t) { a.resize(max(a.size(), t.a.size())); for (size_t i = 0; i < t.a.size(); i++) { a[i] -= t.a[i]; } normalize(); return *this; } poly operator+(const poly &t) const { return poly(*this) += t; } poly operator-(const poly &t) const { return poly(*this) -= t; } poly mod_xk(size_t k) const { // get first k coefficients return vector(begin(a), begin(a) + min(k, a.size())); } poly mul_xk(size_t k) const { // multiply by x^k auto res = a; res.insert(begin(res), k, 0); return res; } poly div_xk(size_t k) const { // drop first k coefficients return vector(begin(a) + min(k, a.size()), end(a)); } poly substr(size_t l, size_t r) const { // return mod_xk(r).div_xk(l) return vector(begin(a) + min(l, a.size()), begin(a) + min(r, a.size())); } poly operator*=(const poly &t) { a = fft::multiply_mod(a, t.a); normalize(); return *this; } poly operator*(const poly &t) const { return poly(*this) *= t; } poly reverse(size_t n) const { // computes x^n A(x^{-1}) auto res = a; res.resize(max(n, res.size())); return vector(res.rbegin(), res.rbegin() + n); } poly reverse() const { return reverse(deg() + 1); } pair divmod_slow(const poly &b) const { // when divisor or quotient is small vector A(a); vector res; while (A.size() >= b.a.size()) { res.push_back(A.back() / b.a.back()); if (res.back() != T(0)) { for (size_t i = 0; i < b.a.size(); i++) { A[A.size() - i - 1] -= res.back() * b.a[b.a.size() - i - 1]; } } A.pop_back(); } std::reverse(begin(res), end(res)); return {res, A}; } pair divmod_hint(poly const &b, poly const &binv) const { // when inverse is known assert(!b.is_zero()); if (deg() < b.deg()) { return {poly{0}, *this}; } int d = deg() - b.deg(); if (min(d, b.deg()) < magic) { return divmod_slow(b); } poly D = (reverse().mod_xk(d + 1) * binv.mod_xk(d + 1)).mod_xk(d + 1).reverse(d + 1); return {D, *this - D * b}; } pair divmod(const poly &b) const { // returns quotiend and remainder of a mod b assert(!b.is_zero()); if (deg() < b.deg()) { return {poly{0}, *this}; } int d = deg() - b.deg(); if (min(d, b.deg()) < magic) { return divmod_slow(b); } poly D = (reverse().mod_xk(d + 1) * b.reverse().inv(d + 1)).mod_xk(d + 1).reverse(d + 1); return {D, *this - D * b}; } // (ax+b) / (cx+d) struct transform { poly a, b, c, d; transform(poly a, poly b = T(1), poly c = T(1), poly d = T(0)) : a(a), b(b), c(c), d(d) {} transform operator*(transform const &t) { return {a * t.a + b * t.c, a * t.b + b * t.d, c * t.a + d * t.c, c * t.b + d * t.d}; } transform adj() { return transform(d, -b, -c, a); } auto apply(poly A, poly B) { return make_pair(a * A + b * B, c * A + d * B); } }; template static void concat(vector &a, vector const &b) { for (auto it : b) { a.push_back(it); } } // finds a transform that changes A/B to A'/B' such that // deg B' is at least 2 times less than deg A static pair, transform> half_gcd(poly A, poly B) { assert(A.deg() >= B.deg()); int m = (A.deg() + 1) / 2; if (B.deg() < m) { return {{}, {T(1), T(0), T(0), T(1)}}; } auto [ar, Tr] = half_gcd(A.div_xk(m), B.div_xk(m)); tie(A, B) = Tr.adj().apply(A, B); if (B.deg() < m) { return {ar, Tr}; } auto [ai, R] = A.divmod(B); tie(A, B) = make_pair(B, R); int k = 2 * m - B.deg(); auto [as, Ts] = half_gcd(A.div_xk(k), B.div_xk(k)); concat(ar, {ai}); concat(ar, as); return {ar, Tr * transform(ai) * Ts}; } // return a transform that reduces A / B to gcd(A, B) / 0 static pair, transform> full_gcd(poly A, poly B) { vector ak; vector trs; while (!B.is_zero()) { if (2 * B.deg() > A.deg()) { auto [a, Tr] = half_gcd(A, B); concat(ak, a); trs.push_back(Tr); tie(A, B) = trs.back().adj().apply(A, B); } else { auto [a, R] = A.divmod(B); ak.push_back(a); trs.emplace_back(a); tie(A, B) = make_pair(B, R); } } trs.emplace_back(T(1), T(0), T(0), T(1)); while (trs.size() >= 2) { trs[trs.size() - 2] = trs[trs.size() - 2] * trs[trs.size() - 1]; trs.pop_back(); } return {ak, trs.back()}; } static poly gcd(poly A, poly B) { if (A.deg() < B.deg()) { return full_gcd(B, A); } auto Tr = fraction(A, B); return Tr.d * A - Tr.b * B; } // Returns the characteristic polynomial // of the minimum linear recurrence for the sequence poly min_rec_slow(int d) const { auto R1 = mod_xk(d + 1).reverse(d + 1), R2 = xk(d + 1); auto Q1 = poly(T(1)), Q2 = poly(T(0)); while (!R2.is_zero()) { auto [a, nR] = R1.divmod(R2); // R1 = a*R2 + nR, deg nR < deg R2 tie(R1, R2) = make_tuple(R2, nR); tie(Q1, Q2) = make_tuple(Q2, Q1 + a * Q2); if (R2.deg() < Q2.deg()) { return Q2 / Q2.lead(); } } assert(0); } // calculate inv to *this modulo t // quadratic complexity optional inv_mod_slow(poly const &t) const { auto R1 = *this, R2 = t; auto Q1 = poly(T(1)), Q2 = poly(T(0)); int k = 0; while (!R2.is_zero()) { k ^= 1; auto [a, nR] = R1.divmod(R2); tie(R1, R2) = make_tuple(R2, nR); tie(Q1, Q2) = make_tuple(Q2, Q1 + a * Q2); } if (R1.deg() > 0) { return nullopt; } else { return (k ? -Q1 : Q1) / R1[0]; } } optional inv_mod(poly const &t) const { assert(!t.is_zero()); if (false && min(deg(), t.deg()) < magic) { return inv_mod_slow(t); } auto A = t, B = *this % t; auto [a, Tr] = full_gcd(A, B); auto g = Tr.d * A - Tr.b * B; if (g.deg() != 0) { return nullopt; } return -Tr.b / g[0]; }; poly operator/(const poly &t) const { return divmod(t).first; } poly operator%(const poly &t) const { return divmod(t).second; } poly operator/=(const poly &t) { return *this = divmod(t).first; } poly operator%=(const poly &t) { return *this = divmod(t).second; } poly operator*=(const T &x) { for (auto &it : a) { it *= x; } normalize(); return *this; } poly operator/=(const T &x) { for (auto &it : a) { it /= x; } normalize(); return *this; } poly operator*(const T &x) const { return poly(*this) *= x; } poly operator/(const T &x) const { return poly(*this) /= x; } poly conj() const { // A(x) -> A(-x) auto res = *this; for (int i = 1; i <= deg(); i += 2) { res[i] = -res[i]; } return res; } void print(int n) const { for (int i = 0; i < n; i++) { cout << (*this)[i] << ' '; } cout << "\n"; } void print() const { print(deg() + 1); } T eval(T x) const { // evaluates in single point x T res(0); for (int i = deg(); i >= 0; i--) { res *= x; res += a[i]; } return res; } T lead() const { // leading coefficient assert(!is_zero()); return a.back(); } int deg() const { // degree, -1 for P(x) = 0 return (int)a.size() - 1; } bool is_zero() const { return a.empty(); } T operator[](int idx) const { return idx < 0 || idx > deg() ? T(0) : a[idx]; } T &coef(size_t idx) { // mutable reference at coefficient return a[idx]; } bool operator==(const poly &t) const { return a == t.a; } bool operator!=(const poly &t) const { return a != t.a; } poly deriv(int k = 1) { // calculate derivative if (deg() + 1 < k) { return poly(T(0)); } vector res(deg() + 1 - k); for (int i = k; i <= deg(); i++) { res[i - k] = fact(i) * rfact(i - k) * a[i]; } return res; } poly integr() { // calculate integral with C = 0 vector res(deg() + 2); for (int i = 0; i <= deg(); i++) { res[i + 1] = a[i] / T(i + 1); } return res; } size_t trailing_xk() const { // Let p(x) = x^k * t(x), return k if (is_zero()) { return -1; } int res = 0; while (a[res] == T(0)) { res++; } return res; } poly log(size_t n) { // calculate log p(x) mod x^n assert(a[0] == T(1)); return (deriv().mod_xk(n) * inv(n)).integr().mod_xk(n); } poly exp(size_t n) { // calculate exp p(x) mod x^n if (is_zero()) { return T(1); } assert(a[0] == T(0)); poly ans = T(1); size_t a = 1; while (a < n) { poly C = ans.log(2 * a).div_xk(a) - substr(a, 2 * a); ans -= (ans * C).mod_xk(a).mul_xk(a); a *= 2; } return ans.mod_xk(n); } poly pow_bin(int64_t k, size_t n) { // O(n log n log k) if (k == 0) { return poly(1).mod_xk(n); } else { auto t = pow(k / 2, n); t *= t; return (k % 2 ? *this * t : t).mod_xk(n); } } // Do not compute inverse from scratch poly powmod_hint(int64_t k, poly const &md, poly const &mdinv) { if (k == 0) { return poly(1); } else { auto t = powmod_hint(k / 2, md, mdinv); t = (t * t).divmod_hint(md, mdinv).second; if (k % 2) { t = (t * *this).divmod_hint(md, mdinv).second; } return t; } } poly powmod(int64_t k, poly const &md) { auto mdinv = md.reverse().inv(md.deg() + 1); return powmod_hint(k, md, mdinv); } // O(d * n) with the derivative trick from // https://codeforces.com/blog/entry/73947?#comment-581173 poly pow_dn(int64_t k, size_t n) { if (n == 0) { return poly(T(0)); } assert((*this)[0] != T(0)); vector Q(n); Q[0] = ipow(a[0], k); for (int i = 1; i < (int)n; i++) { for (int j = 1; j <= min(deg(), i); j++) { Q[i] += a[j] * Q[i - j] * (T(k) * T(j) - T(i - j)); } Q[i] /= T(i) * a[0]; } return Q; } // calculate p^k(n) mod x^n in O(n log n) // might be quite slow due to high constant poly pow(int64_t k, size_t n) { if (is_zero()) { return k ? *this : poly(1); } int i = trailing_xk(); if (i > 0) { return k >= int64_t(n + i - 1) / i ? poly(T(0)) : div_xk(i).pow(k, n - i * k).mul_xk(i * k); } if (min(deg(), (int)n) <= magic) { return pow_dn(k, n); } if (k <= magic) { return pow_bin(k, n); } T j = a[i]; poly t = *this / j; return ((t.log(n) * T(k)).exp(n).mod_xk(n)) * ipow(j, k); } // returns nullopt if undefined optional sqrt(size_t n) const { if (is_zero()) { return *this; } int i = trailing_xk(); if (i % 2) { return nullopt; } else if (i > 0) { auto ans = div_xk(i).sqrt(n - i / 2); return ans ? ans->mul_xk(i / 2) : ans; } // do it or assume this to be 1 auto st = (*this)[0].sqrt(); if (st) { poly ans = *st; size_t a = 1; while (a < n) { a *= 2; ans -= (ans - mod_xk(a) * ans.inv(a)).mod_xk(a) / 2; } return ans.mod_xk(n); } return nullopt; } poly mulx(T a) const { // component-wise multiplication with a^k T cur = 1; poly res(*this); for (int i = 0; i <= deg(); i++) { res.coef(i) *= cur; cur *= a; } return res; } poly mulx_sq(T a) const { // component-wise multiplication with a^{k^2} T cur = a; T total = 1; T aa = a * a; poly res(*this); for (int i = 0; i <= deg(); i++) { res.coef(i) *= total; total *= cur; cur *= aa; } return res; } vector chirpz_even(T z, int n) const { // P(1), P(z^2), P(z^4), ..., P(z^2(n-1)) int m = deg(); if (is_zero()) { return vector(n, 0); } vector vv(m + n); T zi = T(1) / z; T zz = zi * zi; T cur = zi; T total = 1; for (int i = 0; i <= max(n - 1, m); i++) { if (i <= m) { vv[m - i] = total; } if (i < n) { vv[m + i] = total; } total *= cur; cur *= zz; } poly w = (mulx_sq(z) * vv).substr(m, m + n).mulx_sq(z); vector res(n); for (int i = 0; i < n; i++) { res[i] = w[i]; } return res; } vector chirpz(T z, int n) const { // P(1), P(z), P(z^2), ..., P(z^(n-1)) if (is_zero() || n == 0) { return vector(n); } if (z == T(0) || n == 1) { vector ans(n, a[0]); for (int i = 1; i <= deg(); i++) ans[0] += a[i]; return ans; } auto even = chirpz_even(z, (n + 1) / 2); auto odd = mulx(z).chirpz_even(z, n / 2); vector ans(n); for (int i = 0; i < n / 2; i++) { ans[2 * i] = even[i]; ans[2 * i + 1] = odd[i]; } if (n % 2 == 1) { ans[n - 1] = even.back(); } return ans; } static auto resultant(poly a, poly b) { // computes resultant of a and b if (b.is_zero()) { return 0; } else if (b.deg() == 0) { return bpow(b.lead(), a.deg()); } else { int pw = a.deg(); a %= b; pw -= a.deg(); auto mul = bpow(b.lead(), pw) * T((b.deg() & a.deg() & 1) ? -1 : 1); auto ans = resultant(b, a); return ans * mul; } } static poly xk(size_t n) { // P(x) = x^n return poly(T(1)).mul_xk(n); } static poly ones(size_t n) { // P(x) = 1 + x + ... + x^{n-1} return vector(n, 1); } static poly expx(size_t n) { // P(x) = e^x (mod x^n) return ones(n).borel(); } // [x^k] (a corr b) = sum_{i+j=k} ai*b{m-j} // = sum_{i-j=k-m} ai*bj static poly corr(poly a, poly b) { // cross-correlation return a * b.reverse(); } poly invborel() const { // ak *= k! auto res = *this; for (int i = 0; i <= deg(); i++) { res.coef(i) *= fact(i); } return res; } poly borel() const { // ak /= k! auto res = *this; for (int i = 0; i <= deg(); i++) { res.coef(i) *= rfact(i); } return res; } poly shift(T a) const { // P(x + a) return (expx(deg() + 1).mulx(a).reverse() * invborel()).div_xk(deg()).borel(); } poly x2() { // P(x) -> P(x^2) vector res(2 * a.size()); for (size_t i = 0; i < a.size(); i++) { res[2 * i] = a[i]; } return res; } // Return {P0, P1}, where P(x) = P0(x) + xP1(x) pair bisect() const { vector res[2]; res[0].reserve(deg() / 2 + 1); res[1].reserve(deg() / 2 + 1); for (int i = 0; i <= deg(); i++) { res[i % 2].push_back(a[i]); } return {res[0], res[1]}; } poly inv(int n) { poly q(a[0].inv()); for (int i = 1; i < n; i <<= 1) { poly p = poly(2) - q * mod_xk(i * 2); q = (p * q).mod_xk(i * 2); } return q.mod_xk(n); } // compute A(B(x)) mod x^n in O(n^2) static poly compose(poly A, poly B, int n) { int q = std::sqrt(n); vector Bk(q); auto Bq = B.pow(q, n); Bk[0] = poly(T(1)); for (int i = 1; i < q; i++) { Bk[i] = (Bk[i - 1] * B).mod_xk(n); } poly Bqk(1); poly ans; for (int i = 0; i <= n / q; i++) { poly cur; for (int j = 0; j < q; j++) { cur += Bk[j] * A[i * q + j]; } ans += (Bqk * cur).mod_xk(n); Bqk = (Bqk * Bq).mod_xk(n); } return ans; } vector eval(vector &tree, int v, int l, int r, vector &vec) { // auxiliary evaluation function if (r - l == 1) { return {eval(vec[l])}; } else { auto m = l + (r - l) / 2; auto A = (*this % tree[2 * v]).eval(tree, 2 * v, l, m, vec); auto B = (*this % tree[2 * v + 1]).eval(tree, 2 * v + 1, m, r, vec); A.insert(end(A), begin(B), end(B)); return A; } } vector eval(vector x) { // evaluate polynomial in (x1, ..., xn) int n = x.size(); if (is_zero()) { return vector(n, T(0)); } vector tree(4 * n); build(tree, 1, 0, sz(x), x); return eval(tree, 1, 0, sz(x), x); } poly inter(vector &tree, int v, int l, int r, int ly, int ry, vector &vecx, vector &vecy) { // auxiliary interpolation function if (r - l == 1) { return {vecy[ly] / a[0]}; } else { auto m = l + (r - l) / 2; auto my = ly + (ry - ly) / 2; auto A = (*this % tree[2 * v]).inter(tree, 2 * v, l, m, ly, my, vecx, vecy); auto B = (*this % tree[2 * v + 1]).inter(tree, 2 * v + 1, m, r, my, ry, vecx, vecy); return A * tree[2 * v + 1] + B * tree[2 * v]; } } static poly build(vector &res, int v, int L, int R, vector &vec) { // builds evaluation tree for (x-a1)(x-a2)...(x-an) if (R - L == 1) { return res[v] = vector{-vec[L], 1}; } else { int M = L + (R - L) / 2; return res[v] = build(res, 2 * v, L, M, vec) * build(res, 2 * v + 1, M, R, vec); } } static poly inter(vector x, vector y) { // interpolates minimum polynomial from (xi, yi) pairs int n = x.size(); vector tree(4 * n); return build(tree, 1, 0, sz(x), x).deriv().inter(tree, 1, 0, sz(x), 0, sz(y), x, y); } }; using polyn = poly; int main() { ios_base::sync_with_stdio(0); cin.tie(0); cout.tie(0); int n; cin >> n; mint ans = 0; for (int i = 2; i <= n - 2; i++) { mint p1 = fact(n) * rfact(i - 2) * rfact(n - 2 - i) / mint(1ll * (n - 1) * (n - i) * i); ans += p1 * p1; } cout << ans << "\n"; }