#line 2 "library/template/template.hpp" #include using namespace std; #line 2 "library/template/macro.hpp" #define rep(i, a, b) for (int i = (a); i < (int)(b); i++) #define rrep(i, a, b) for (int i = (int)(b) - 1; i >= (a); i--) #define ALL(v) (v).begin(), (v).end() #define UNIQUE(v) sort(ALL(v)), (v).erase(unique(ALL(v)), (v).end()) #define SZ(v) (int)v.size() #define MIN(v) *min_element(ALL(v)) #define MAX(v) *max_element(ALL(v)) #define LB(v, x) int(lower_bound(ALL(v), (x)) - (v).begin()) #define UB(v, x) int(upper_bound(ALL(v), (x)) - (v).begin()) #define YN(b) cout << ((b) ? "YES" : "NO") << "\n"; #define Yn(b) cout << ((b) ? "Yes" : "No") << "\n"; #define yn(b) cout << ((b) ? "yes" : "no") << "\n"; #line 6 "library/template/template.hpp" #line 2 "library/template/util.hpp" using uint = unsigned int; using ll = long long int; using ull = unsigned long long; using i128 = __int128_t; using u128 = __uint128_t; template S SUM(const vector& a) { return accumulate(ALL(a), S(0)); } template inline bool chmin(T& a, T b) { if (a > b) { a = b; return true; } return false; } template inline bool chmax(T& a, T b) { if (a < b) { a = b; return true; } return false; } template int popcnt(T x) { return __builtin_popcountll(x); } template int topbit(T x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); } template int lowbit(T x) { return (x == 0 ? -1 : __builtin_ctzll(x)); } #line 8 "library/template/template.hpp" #line 2 "library/template/inout.hpp" struct Fast { Fast() { cin.tie(nullptr); ios_base::sync_with_stdio(false); cout << fixed << setprecision(15); } } fast; template istream& operator>>(istream& is, pair& p) { return is >> p.first >> p.second; } template ostream& operator<<(ostream& os, const pair& p) { return os << p.first << " " << p.second; } template istream& operator>>(istream& is, vector& a) { for (auto& v : a) is >> v; return is; } template ostream& operator<<(ostream& os, const vector& a) { for (auto it = a.begin(); it != a.end();) { os << *it; if (++it != a.end()) os << " "; } return os; } template ostream& operator<<(ostream& os, const set& st) { os << "{"; for (auto it = st.begin(); it != st.end();) { os << *it; if (++it != st.end()) os << ","; } os << "}"; return os; } template ostream& operator<<(ostream& os, const map& mp) { os << "{"; for (auto it = mp.begin(); it != mp.end();) { os << it->first << ":" << it->second; if (++it != mp.end()) os << ","; } os << "}"; return os; } ostream& operator<<(ostream& os, __uint128_t x) { char buf[40]; size_t k = 0; while (x > 0) buf[k++] = (char)(x % 10 + '0'), x /= 10; if (k == 0) buf[k++] = '0'; while (k) cout << buf[--k]; return cout; } ostream& operator<<(ostream& os, __int128_t x) { return x < 0 ? (cout << '-' << (__uint128_t)(-x)) : (cout << (__uint128_t)x); } void in() {} template void in(T& t, U&... u) { cin >> t; in(u...); } void out() { cout << "\n"; } template void out(const T& t, const U&... u) { cout << t; if (sizeof...(u)) cout << sep; out(u...); } #line 10 "library/template/template.hpp" #line 2 "library/template/debug.hpp" #ifdef LOCAL #define debug 1 #define show(...) _show(0, #__VA_ARGS__, __VA_ARGS__) #else #define debug 0 #define show(...) true #endif template void _show(int i, T name) { cerr << '\n'; } template void _show(int i, const T1& a, const T2& b, const T3&... c) { for (; a[i] != ',' && a[i] != '\0'; i++) cerr << a[i]; cerr << ":" << b << " "; _show(i + 1, a, c...); } #line 2 "main.cpp" #line 2 "library/math/util.hpp" namespace Math { template T safe_mod(T a, T b) { assert(b != 0); if (b < 0) a = -a, b = -b; a %= b; return a >= 0 ? a : a + b; } template T floor(T a, T b) { assert(b != 0); if (b < 0) a = -a, b = -b; return a >= 0 ? a / b : (a + 1) / b - 1; } template T ceil(T a, T b) { assert(b != 0); if (b < 0) a = -a, b = -b; return a > 0 ? (a - 1) / b + 1 : a / b; } long long isqrt(long long n) { if (n <= 0) return 0; long long x = sqrt(n); while ((x + 1) * (x + 1) <= n) x++; while (x * x > n) x--; return x; } // return g=gcd(a,b) // a*x+b*y=g // - b!=0 -> 0<=x<|b|/g // - b=0 -> ax=g template T ext_gcd(T a, T b, T& x, T& y) { T a0 = a, b0 = b; bool sgn_a = a < 0, sgn_b = b < 0; if (sgn_a) a = -a; if (sgn_b) b = -b; if (b == 0) { x = sgn_a ? -1 : 1; y = 0; return a; } T x00 = 1, x01 = 0, x10 = 0, x11 = 1; while (b != 0) { T q = a / b, r = a - b * q; x00 -= q * x01; x10 -= q * x11; swap(x00, x01); swap(x10, x11); a = b, b = r; } x = x00, y = x10; if (sgn_a) x = -x; if (sgn_b) y = -y; if (b0 != 0) { a0 /= a, b0 /= a; if (b0 < 0) a0 = -a0, b0 = -b0; T q = x >= 0 ? x / b0 : (x + 1) / b0 - 1; x -= b0 * q; y += a0 * q; } return a; } constexpr long long inv_mod(long long x, long long m) { x %= m; if (x < 0) x += m; long long a = m, b = x; long long y0 = 0, y1 = 1; while (b > 0) { long long q = a / b; swap(a -= q * b, b); swap(y0 -= q * y1, y1); } if (y0 < 0) y0 += m / a; return y0; } long long pow_mod(long long x, long long n, long long m) { x = (x % m + m) % m; long long y = 1; while (n) { if (n & 1) y = y * x % m; x = x * x % m; n >>= 1; } return y; } constexpr long long pow_mod_constexpr(long long x, long long n, int m) { if (m == 1) return 0; unsigned int _m = (unsigned int)(m); unsigned long long r = 1; unsigned long long y = x % m; if (y >= m) y += m; while (n) { if (n & 1) r = (r * y) % _m; y = (y * y) % _m; n >>= 1; } return r; } constexpr bool is_prime_constexpr(int n) { if (n <= 1) return false; if (n == 2 || n == 7 || n == 61) return true; if (n % 2 == 0) return false; long long d = n - 1; while (d % 2 == 0) d /= 2; constexpr long long bases[3] = {2, 7, 61}; for (long long a : bases) { long long t = d; long long y = pow_mod_constexpr(a, t, n); while (t != n - 1 && y != 1 && y != n - 1) { y = y * y % n; t <<= 1; } if (y != n - 1 && t % 2 == 0) { return false; } } return true; } template constexpr bool is_prime = is_prime_constexpr(n); }; // namespace Math #line 3 "library/modint/modint.hpp" template struct ModInt { using mint = ModInt; static constexpr unsigned int get_mod() { return m; } static mint raw(int v) { mint x; x._v = v; return x; } ModInt() : _v(0) {} ModInt(int64_t v) { long long x = (long long)(v % (long long)(umod())); if (x < 0) x += umod(); _v = (unsigned int)(x); } unsigned int val() const { return _v; } mint& operator++() { _v++; if (_v == umod()) _v = 0; return *this; } mint& operator--() { if (_v == 0) _v = umod(); _v--; return *this; } mint operator++(int) { mint result = *this; ++*this; return result; } mint operator--(int) { mint result = *this; --*this; return result; } mint& operator+=(const mint& rhs) { _v += rhs._v; if (_v >= umod()) _v -= umod(); return *this; } mint& operator-=(const mint& rhs) { _v -= rhs._v; if (_v >= umod()) _v += umod(); return *this; } mint& operator*=(const mint& rhs) { unsigned long long z = _v; z *= rhs._v; _v = (unsigned int)(z % umod()); return *this; } mint& operator/=(const mint& rhs) { return *this *= rhs.inv(); } mint operator+() const { return *this; } mint operator-() const { return mint() - *this; } mint pow(long long n) const { assert(0 <= n); mint x = *this, r = 1; while (n) { if (n & 1) r *= x; x *= x; n >>= 1; } return r; } mint inv() const { if (is_prime) { assert(_v); return pow(umod() - 2); } else { auto inv = Math::inv_mod(_v, umod()); return raw(inv); } } friend mint operator+(const mint& lhs, const mint& rhs) { return mint(lhs) += rhs; } friend mint operator-(const mint& lhs, const mint& rhs) { return mint(lhs) -= rhs; } friend mint operator*(const mint& lhs, const mint& rhs) { return mint(lhs) *= rhs; } friend mint operator/(const mint& lhs, const mint& rhs) { return mint(lhs) /= rhs; } friend bool operator==(const mint& lhs, const mint& rhs) { return lhs._v == rhs._v; } friend bool operator!=(const mint& lhs, const mint& rhs) { return lhs._v != rhs._v; } friend istream& operator>>(istream& is, mint& x) { int64_t v; is >> v; x = mint(v); return is; } friend ostream& operator<<(ostream& os, const mint& x) { return os << x.val(); } private: unsigned int _v; static constexpr unsigned int umod() { return m; } static constexpr bool is_prime = Math::is_prime; }; #line 4 "main.cpp" using mint = ModInt<120586241>; #line 2 "library/fps/fps-ntt-friendly.hpp" #line 2 "library/fft/ntt.hpp" template struct NTT { static constexpr unsigned int mod = mint::get_mod(); static constexpr unsigned long long pow_constexpr(unsigned long long x, unsigned long long n, unsigned long long m) { unsigned long long y = 1; while (n) { if (n & 1) y = y * x % m; x = x * x % m; n >>= 1; } return y; } static constexpr unsigned int get_g() { unsigned long long x = 2; while (pow_constexpr(x, (mod - 1) >> 1, mod) == 1) x += 1; return x; } static constexpr unsigned int g = get_g(); static constexpr int rank2 = __builtin_ctzll(mod - 1); array root; array iroot; array rate2; array irate2; array rate3; array irate3; NTT() { root[rank2] = mint(g).pow((mod - 1) >> rank2); iroot[rank2] = root[rank2].inv(); for (int i = rank2 - 1; i >= 0; i--) { root[i] = root[i + 1] * root[i + 1]; iroot[i] = iroot[i + 1] * iroot[i + 1]; } { mint prod = 1, iprod = 1; for (int i = 0; i <= rank2 - 2; i++) { rate2[i] = root[i + 2] * prod; irate2[i] = iroot[i + 2] * iprod; prod *= iroot[i + 2]; iprod *= root[i + 2]; } } { mint prod = 1, iprod = 1; for (int i = 0; i <= rank2 - 3; i++) { rate3[i] = root[i + 3] * prod; irate3[i] = iroot[i + 3] * iprod; prod *= iroot[i + 3]; iprod *= root[i + 3]; } } } void ntt(vector& a) { int n = int(a.size()); int h = __builtin_ctzll((unsigned int)n); a.resize(1 << h); int len = 0; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed while (len < h) { if (h - len == 1) { int p = 1 << (h - len - 1); mint rot = 1; for (int s = 0; s < (1 << len); s++) { int offset = s << (h - len); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p] * rot; a[i + offset] = l + r; a[i + offset + p] = l - r; } if (s + 1 != (1 << len)) rot *= rate2[__builtin_ctzll(~(unsigned int)(s))]; } len++; } else { // 4-base int p = 1 << (h - len - 2); mint rot = 1, imag = root[2]; for (int s = 0; s < (1 << len); s++) { mint rot2 = rot * rot; mint rot3 = rot2 * rot; int offset = s << (h - len); for (int i = 0; i < p; i++) { auto mod2 = 1ULL * mint::get_mod() * mint::get_mod(); auto a0 = 1ULL * a[i + offset].val(); auto a1 = 1ULL * a[i + offset + p].val() * rot.val(); auto a2 = 1ULL * a[i + offset + 2 * p].val() * rot2.val(); auto a3 = 1ULL * a[i + offset + 3 * p].val() * rot3.val(); auto a1na3imag = 1ULL * mint(a1 + mod2 - a3).val() * imag.val(); auto na2 = mod2 - a2; a[i + offset] = a0 + a2 + a1 + a3; a[i + offset + 1 * p] = a0 + a2 + (2 * mod2 - (a1 + a3)); a[i + offset + 2 * p] = a0 + na2 + a1na3imag; a[i + offset + 3 * p] = a0 + na2 + (mod2 - a1na3imag); } if (s + 1 != (1 << len)) rot *= rate3[__builtin_ctzll(~(unsigned int)(s))]; } len += 2; } } } void intt(vector& a) { int n = int(a.size()); int h = __builtin_ctzll((unsigned int)n); a.resize(1 << h); int len = h; // a[i, i+(n>>len), i+2*(n>>len), ..] is transformed while (len) { if (len == 1) { int p = 1 << (h - len); mint irot = 1; for (int s = 0; s < (1 << (len - 1)); s++) { int offset = s << (h - len + 1); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p]; a[i + offset] = l + r; a[i + offset + p] = (unsigned long long)(mint::get_mod() + l.val() - r.val()) * irot.val(); } if (s + 1 != (1 << (len - 1))) irot *= irate2[__builtin_ctzll(~(unsigned int)(s))]; } len--; } else { // 4-base int p = 1 << (h - len); mint irot = 1, iimag = iroot[2]; for (int s = 0; s < (1 << (len - 2)); s++) { mint irot2 = irot * irot; mint irot3 = irot2 * irot; int offset = s << (h - len + 2); for (int i = 0; i < p; i++) { auto a0 = 1ULL * a[i + offset + 0 * p].val(); auto a1 = 1ULL * a[i + offset + 1 * p].val(); auto a2 = 1ULL * a[i + offset + 2 * p].val(); auto a3 = 1ULL * a[i + offset + 3 * p].val(); auto a2na3iimag = 1ULL * mint((mint::get_mod() + a2 - a3) * iimag.val()).val(); a[i + offset] = a0 + a1 + a2 + a3; a[i + offset + 1 * p] = (a0 + (mint::get_mod() - a1) + a2na3iimag) * irot.val(); a[i + offset + 2 * p] = (a0 + a1 + (mint::get_mod() - a2) + (mint::get_mod() - a3)) * irot2.val(); a[i + offset + 3 * p] = (a0 + (mint::get_mod() - a1) + (mint::get_mod() - a2na3iimag)) * irot3.val(); } if (s + 1 != (1 << (len - 2))) irot *= irate3[__builtin_ctzll(~(unsigned int)(s))]; } len -= 2; } } mint e = mint(n).inv(); for (auto& x : a) x *= e; } vector multiply(const vector& a, const vector& b) { if (a.empty() || b.empty()) return vector(); int n = a.size(), m = b.size(); int sz = n + m - 1; if (n <= 30 || m <= 30) { if (n > 30) return multiply(b, a); vector res(sz); for (int i = 0; i < n; i++) for (int j = 0; j < m; j++) res[i + j] += a[i] * b[j]; return res; } int sz1 = 1; while (sz1 < sz) sz1 <<= 1; vector res(sz1); for (int i = 0; i < n; i++) res[i] = a[i]; ntt(res); if (a == b) for (int i = 0; i < sz1; i++) res[i] *= res[i]; else { vector c(sz1); for (int i = 0; i < m; i++) c[i] = b[i]; ntt(c); for (int i = 0; i < sz1; i++) res[i] *= c[i]; } intt(res); res.resize(sz); return res; } // c[i]=sum[j]a[j]b[i+j] vector middle_product(const vector& a, const vector& b) { if (b.empty() || a.size() > b.size()) return {}; int n = a.size(), m = b.size(); int sz = m - n + 1; if (n <= 30 || sz <= 30) { vector res(sz); for (int i = 0; i < sz; i++) for (int j = 0; j < n; j++) res[i] += a[j] * b[i + j]; return res; } int sz1 = 1; while (sz1 < m) sz1 <<= 1; vector res(sz1), b2(sz1); reverse_copy(a.begin(), a.end(), res.begin()); copy(b.begin(), b.end(), b2.begin()); ntt(res); ntt(b2); for (int i = 0; i < res.size(); i++) res[i] *= b2[i]; intt(res); res.resize(m); res.erase(res.begin(), res.begin() + n - 1); return res; } void ntt_doubling(vector& a) { int n = (int)a.size(); auto b = a; intt(b); mint r = 1, zeta = mint(g).pow((mint::get_mod() - 1) / (n << 1)); for (int i = 0; i < n; i++) b[i] *= r, r *= zeta; ntt(b); copy(b.begin(), b.end(), back_inserter(a)); } }; /** * @brief NTT (数論変換) * @docs docs/fft/ntt.md */ #line 2 "library/fps/formal-power-series.hpp" template struct FormalPowerSeries : vector { using vector::vector; using FPS = FormalPowerSeries; FormalPowerSeries(const vector& r) : vector(r) {} FormalPowerSeries(vector&& r) : vector(std::move(r)) {} FPS& operator=(const vector& r) { vector::operator=(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 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.at(g.size() - 1).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; } FPS pre(int sz) const { return FPS(begin(*this), begin(*this) + min((int)this->size(), sz)); } FPS operator>>=(int sz) { assert(sz >= 0); if ((int)this->size() <= sz) this->clear(); else this->erase(this->begin(), this->begin() + sz); return *this; } 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) { assert(sz >= 0); this->insert(this->begin(), sz, mint(0)); return *this; } 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::get_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)[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_ntt(); FPS& operator*=(const FPS& r); FPS middle_product(const FPS& r) const; void ntt(); void intt(); void ntt_doubling(); static int ntt_root(); FPS inv(int deg = -1) const; FPS exp(int deg = -1) const; }; template void* FormalPowerSeries::ntt_ptr = nullptr; #line 5 "library/fps/fps-ntt-friendly.hpp" template void FormalPowerSeries::set_ntt() { if (!ntt_ptr) ntt_ptr = new NTT; } template FormalPowerSeries& FormalPowerSeries::operator*=(const FormalPowerSeries& r) { if (this->empty() || r.empty()) { this->clear(); return *this; } set_ntt(); auto ret = static_cast*>(ntt_ptr)->multiply(*this, r); return *this = FormalPowerSeries(ret.begin(), ret.end()); } template FormalPowerSeries FormalPowerSeries::middle_product(const FormalPowerSeries& r) const { set_ntt(); auto ret = static_cast*>(ntt_ptr)->middle_product(*this, r); return FormalPowerSeries(ret.begin(), ret.end()); } template void FormalPowerSeries::ntt() { set_ntt(); static_cast*>(ntt_ptr)->ntt(*this); } template void FormalPowerSeries::intt() { set_ntt(); static_cast*>(ntt_ptr)->intt(*this); } template void FormalPowerSeries::ntt_doubling() { set_ntt(); static_cast*>(ntt_ptr)->ntt_doubling(*this); } template int FormalPowerSeries::ntt_root() { set_ntt(); return static_cast*>(ntt_ptr)->g; } template FormalPowerSeries FormalPowerSeries::inv(int deg) const { assert((*this)[0] != mint(0)); if (deg == -1) deg = (*this).size(); FPS ret{mint(1) / (*this)[0]}; for (int i = 1; i < deg; i <<= 1) ret = (ret + ret - ret * ret * (*this).pre(i << 1)).pre(i << 1); return ret.pre(deg); } template FormalPowerSeries FormalPowerSeries::exp(int deg) const { assert((*this)[0] == mint(0)); if (deg == -1) deg = (*this).size(); FPS ret{mint(1)}; for (int i = 1; i < deg; i <<= 1) ret = (ret * ((*this).pre(i << 1) - ret.log(i << 1) + 1)).pre(i << 1); return ret.pre(deg); } #line 6 "main.cpp" using fps = FormalPowerSeries; #line 3 "library/fps/fps-rational.hpp" template struct FPSRational { using F = FormalPowerSeries; using R = FPSRational; F num, den; FPSRational() : num(F{}), den(F{1}) {} FPSRational(F f) : num(f), den(F{1}) {} FPSRational(F f, F g) : num(f), den(g) {} R& operator+=(const R& r) { num *= r.den; num += den * r.num; den *= r.den; return *this; } R& operator-=(const R& r) { num *= r.den; num -= den * r.num; den *= r.den; return *this; } R& operator*=(const R& r) { num *= r.num; den *= r.den; return *this; } R& operator/=(const R& r) { num *= r.den; den *= r.num; return *this; } R operator+(const R& r) const { return R(*this) += r; } R operator-(const R& r) const { return R(*this) -= r; } R operator*(const R& r) const { return R(*this) *= r; } R operator/(const R& r) const { return R(*this) /= r; } R inv() const { return {den, num}; } F approx(int deg) const { return (den * num.inv(deg)).pre(deg); } }; #line 8 "main.cpp" using rat = FPSRational; #line 3 "library/fps/bostan-mori.hpp" // [x^n]f(x)/g(x) template mint BostanMori(FormalPowerSeries f, FormalPowerSeries g, long long n) { g.shrink(); mint ret = 0; if (f.size() >= g.size()) { auto q = f / g; if (n < q.size()) ret += q[n]; f -= q * g; f.shrink(); } if (f.empty()) return ret; FormalPowerSeries::set_ntt(); if (!FormalPowerSeries::ntt_ptr) { f.resize(g.size() - 1); for (; n > 0; n >>= 1) { auto g1 = g; for (int i = 1; i < g1.size(); i += 2) g1[i] = -g1[i]; auto p = f * g1, q = g * g1; if (n & 1) { for (int i = 0; i < f.size(); i++) f[i] = p[(i << 1) | 1]; } else { for (int i = 0; i < f.size(); i++) f[i] = p[i << 1]; } for (int i = 0; i < g.size(); i++) g[i] = q[i << 1]; } return ret + f[0] / g[0]; } else { int m = 1, log = 0; while (m < g.size()) m <<= 1, log++; mint wi = mint(FormalPowerSeries::ntt_root()).inv().pow((mint::get_mod() - 1) >> (log + 1)); vector rev(m); for (int i = 0; i < rev.size(); i++) rev[i] = (rev[i / 2] / 2) | ((i & 1) << (log - 1)); vector pow(m, 1); for (int i = 1; i < m; i++) pow[rev[i]] = pow[rev[i - 1]] * wi; f.resize(m), g.resize(m); f.ntt(), g.ntt(); mint inv2 = mint(2).inv(); for (; n >= m; n >>= 1) { f.ntt_doubling(), g.ntt_doubling(); if (n & 1) { for (int i = 0; i < m; i++) f[i] = (f[i << 1] * g[(i << 1) | 1] - f[(i << 1) | 1] * g[i << 1]) * inv2 * pow[i]; } else { for (int i = 0; i < m; i++) f[i] = (f[i << 1] * g[(i << 1) | 1] + f[(i << 1) | 1] * g[i << 1]) * inv2; } for (int i = 0; i < m; i++) g[i] = g[i << 1] * g[(i << 1) | 1]; f.resize(m), g.resize(m); } f.intt(), g.intt(); return ret + (f * g.inv())[n]; } } /** * @brief Bostan-Mori * @docs docs/fps/bostan-mori.md */ #line 10 "main.cpp" rat RatSum(vector rs) { for (int i = rs.size() - 1; i >= 1; i--) rs[i - (i & -i)] += rs[i]; return rs[0]; } const mint w = mint(6).pow((mint::get_mod() - 1) / 10); void solve() { int n, k, t; in(n, k, t); vector a(n); in(a); vector r(k, vector(n)); in(r); mint invk = mint(k).inv(); vector pow10(n + 1, 1); rep(i, 0, n) pow10[i + 1] = pow10[i] * 10; int m = pow10[n]; mint invm = mint(m).inv(); vector poww(100, 1); rep(i, 1, poww.size()) poww[i] = poww[i - 1] * w; vector z(m), za(m); auto pre_calc = [&](auto rec, int k, int l, int r, vector f, mint va) -> void { if (k < 0) { z[l] = f[0]; za[l] = va; return; } else { rep(t, 0, 10) { int l1 = l + (r - l) / 10 * t; int r1 = l1 + (r - l) / 10; vector g(f.size() / 10); rep(i, 0, g.size()) rep(j, 0, 10) { g[i] += f[i + pow10[k] * j] * poww[j * t]; } rec(rec, k - 1, l1, r1, g, va * poww[a[k] * t]); } } }; { vector f0(m); rep(i, 0, k) { int x = 0; rep(j, 0, n) x += pow10[j] * r[i][j]; f0[x] = invk; } pre_calc(pre_calc, n - 1, 0, m, f0, 1); } rat f; { vector rs; rep(i, 0, m) rs.push_back(rat(fps{invm}, fps{za[i], -za[i] * z[i]})); f = RatSum(rs); } { vector rs; rep(i, 0, m) rs.push_back(rat(fps{invm}, fps{1, -z[i]})); f /= RatSum(rs); } f /= rat(fps{1, -1}); mint ans = BostanMori(f.num, f.den, t); out(ans); } int main() { int t = 1; // in(t); while (t--) solve(); }