#line 1 "helper.cpp" // #define _GLIBCXX_DEBUG #include // #include // #include // #include // clang-format off std::ostream&operator<<(std::ostream&os,std::int8_t x){return os<<(int)x;} std::ostream&operator<<(std::ostream&os,std::uint8_t x){return os<<(int)x;} std::ostream&operator<<(std::ostream&os,const __int128_t &v){if(!v)os<<"0";__int128_t tmp=v<0?(os<<"-",-v):v;std::string s;while(tmp)s+='0'+(tmp%10),tmp/=10;return std::reverse(s.begin(),s.end()),os< #line 4 "/Users/hashimotoryoma/code/CompetitiveProgramming/Library/Library/src/DataStructure/CsrArray.hpp" template struct ListRange { using Iterator= typename std::vector::const_iterator; Iterator bg, ed; Iterator begin() const { return bg; } Iterator end() const { return ed; } size_t size() const { return std::distance(bg, ed); } const T &operator[](int i) const { return bg[i]; } }; template class CsrArray { std::vector csr; std::vector pos; public: CsrArray()= default; CsrArray(const std::vector &c, const std::vector &p): csr(c), pos(p) {} size_t size() const { return pos.size() - 1; } const ListRange operator[](int i) const { return {csr.cbegin() + pos[i], csr.cbegin() + pos[i + 1]}; } }; #line 10 "/Users/hashimotoryoma/code/CompetitiveProgramming/Library/Library/src/Graph/Tree.hpp" template class Tree { template struct Edge_B { int to; T cost; operator int() const { return to; } }; template struct Edge_B { int to; operator int() const { return to; } }; using Edge= Edge_B; std::vector, std::pair, std::tuple>> es; std::vector g; std::vector P, PP, D, I, L, R, pos; public: Tree(int n): P(n, -2) {} template , std::nullptr_t> = nullptr> void add_edge(int u, int v) { es.emplace_back(u, v), es.emplace_back(v, u); } template , std::nullptr_t> = nullptr> void add_edge(int u, int v, T c) { es.emplace_back(u, v, c), es.emplace_back(v, u, c); } template , std::is_convertible>, std::nullptr_t> = nullptr> void add_edge(int u, int v, T c, U d) /* c:u->v, d:v->u */ { es.emplace_back(u, v, c), es.emplace_back(v, u, d); } void build(int root= 0) { size_t n= P.size(); I.resize(n), PP.resize(n), std::iota(PP.begin(), PP.end(), 0), D.assign(n, 0), L.assign(n, 0), R.assign(n, 0), pos.resize(n + 1), g.resize(es.size()); for (const auto &e: es) ++pos[std::get<0>(e)]; std::partial_sum(pos.begin(), pos.end(), pos.begin()); if constexpr (std::is_same_v) for (const auto &[f, t]: es) g[--pos[f]]= {t}; else for (const auto &[f, t, c]: es) g[--pos[f]]= {t, c}; auto f= [&, i= 0, v= 0, t= 0](int r) mutable { for (P[r]= -1, I[t++]= r; i < t; ++i) for (int u: operator[](v= I[i])) if (P[v] != u) P[I[t++]= u]= v; }; f(root); for (size_t r= 0; r < n; ++r) if (P[r] == -2) f(r); std::vector Z(n, 1), nx(n, -1); for (int i= n, v; i--;) { if (P[v= I[i]] == -1) continue; if (Z[P[v]]+= Z[v]; nx[P[v]] == -1) nx[P[v]]= v; if (Z[nx[P[v]]] < Z[v]) nx[P[v]]= v; } for (int v: I) if (nx[v] != -1) PP[nx[v]]= v; for (int v: I) if (P[v] != -1) PP[v]= PP[PP[v]], D[v]= D[P[v]] + 1; for (int i= n; i--;) L[I[i]]= i; for (int v: I) { int ir= R[v]= L[v] + Z[v]; for (int u: operator[](v)) if (u != P[v] && u != nx[v]) L[u]= ir-= Z[u]; if (nx[v] != -1) L[nx[v]]= L[v] + 1; } for (int i= n; i--;) I[L[i]]= i; } size_t size() const { return P.size(); } const ListRange operator[](int v) const { return {g.cbegin() + pos[v], g.cbegin() + pos[v + 1]}; } int depth(int v) const { return D[v]; } int to_seq(int v) const { return L[v]; } int to_node(int i) const { return I[i]; } int parent(int v) const { return P[v]; } int root(int v) const { for (v= PP[v];; v= PP[P[v]]) if (P[v] == -1) return v; } bool connected(int u, int v) const { return root(u) == root(v); } int lca(int u, int v) const { for (;; v= P[PP[v]]) { if (L[u] > L[v]) std::swap(u, v); if (PP[u] == PP[v]) return u; } } int la(int v, int k) const { assert(k <= D[v]); for (int u;; k-= L[v] - L[u] + 1, v= P[u]) if (L[v] - k >= L[u= PP[v]]) return I[L[v] - k]; } int jump(int u, int v, int k) const { if (!k) return u; if (u == v) return -1; if (k == 1) return in_subtree(v, u) ? la(v, D[v] - D[u] - 1) : P[u]; int w= lca(u, v), d_uw= D[u] - D[w], d_vw= D[v] - D[w]; return k > d_uw + d_vw ? -1 : k <= d_uw ? la(u, k) : la(v, d_uw + d_vw - k); } int dist(int u, int v) const { return depth(u) + depth(v) - depth(lca(u, v)) * 2; } // u is in v bool in_subtree(int u, int v) const { return L[v] <= L[u] && L[u] < R[v]; } int subtree_size(int v) const { return R[v] - L[v]; } // half-open interval std::array subtree(int v) const { return std::array{L[v], R[v]}; } // sequence of closed intervals template std::vector> path(int u, int v) const { std::vector> up, down; while (PP[u] != PP[v]) { if (L[u] < L[v]) down.emplace_back(std::array{L[PP[v]], L[v]}), v= P[PP[v]]; else up.emplace_back(std::array{L[u], L[PP[u]]}), u= P[PP[u]]; } if (L[u] < L[v]) down.emplace_back(std::array{L[u] + edge, L[v]}); else if (L[v] + edge <= L[u]) up.emplace_back(std::array{L[u], L[v] + edge}); return up.insert(up.end(), down.rbegin(), down.rend()), up; } }; #line 3 "/Users/hashimotoryoma/code/CompetitiveProgramming/Library/Library/src/Graph/FunctionalGraph.hpp" class FunctionalGraph { std::vector to, rt; Tree<> tree; public: FunctionalGraph(int n): to(n, -1), rt(n, -1), tree(n + 1) {} void add_edge(int src, int dst) { assert(to[src] == -1), to[src]= dst; } void build() { const int n= to.size(); for (int u, w, v= n; v--;) if (rt[v] == -1) { for (rt[v]= -2, w= to[v];; rt[w]= -2, w= to[w]) if (assert(w != -1); rt[w] != -1) { if (rt[w] != -2) w= rt[w]; break; } for (u= v; rt[u] == -2; u= to[u]) rt[u]= w; } for (int v= n; v--;) if (rt[v] == v) tree.add_edge(v, n); else tree.add_edge(v, to[v]); tree.build(n); } template std::enable_if_t, int> jump(int v, Int k) const { int n= to.size(), d= tree.depth(v) - 1; if (k <= d) return tree.jump(v, n, (int)k); int b= to[v= rt[v]], l= (k-= d) % tree.depth(b); if (l == 0) return v; return tree.jump(b, n, l - 1); } // ((a_0,...,a_{i-1}) x 1, (a_i,...,a_{j-1}) x loop_num, (a_j,...,a_m) x 1) template std::enable_if_t, std::array, Int>, 3>> path(int v, Int k) const { std::array, Int>, 3> ret; int n= to.size(), d= tree.depth(v) - 1; if (ret[0].second= 1; k <= d) { for (int e= k; e--; v= to[v]) ret[0].first.push_back(v); return ret; } for (int e= d; e--; v= to[v]) ret[0].first.push_back(v); int b= to[v= rt[v]], c= tree.depth(b), l= (k-= d) % c; ret[1].second= k / c, ret[2].second= 1; for (int e= c; e--; v= to[v]) ret[1].first.push_back(v); for (int e= l; e--; v= to[v]) ret[2].first.push_back(v); return ret; } }; #line 2 "/Users/hashimotoryoma/code/CompetitiveProgramming/Library/Library/src/Internal/Remainder.hpp" namespace math_internal { using namespace std; using u8= uint8_t; using u32= uint32_t; using u64= uint64_t; using i64= int64_t; using u128= __uint128_t; #define CE constexpr #define IL inline #define NORM \ if (n >= mod) n-= mod; \ return n #define PLUS(U, M) \ CE IL U plus(U l, U r) const { \ if (l+= r; l >= M) l-= M; \ return l; \ } #define DIFF(U, C, M) \ CE IL U diff(U l, U r) const { \ if (l-= r; l >> C) l+= M; \ return l; \ } #define SGN(U) \ static CE IL U set(U n) { return n; } \ static CE IL U get(U n) { return n; } \ static CE IL U norm(U n) { return n; } template struct MP_Mo { u_t mod; CE MP_Mo(): mod(0), iv(0), r2(0) {} CE MP_Mo(u_t m): mod(m), iv(inv(m)), r2(-du_t(mod) % mod) {} CE IL u_t mul(u_t l, u_t r) const { return reduce(du_t(l) * r); } PLUS(u_t, mod << 1) DIFF(u_t, A, mod << 1) CE IL u_t set(u_t n) const { return mul(n, r2); } CE IL u_t get(u_t n) const { n= reduce(n); NORM; } CE IL u_t norm(u_t n) const { NORM; } private: u_t iv, r2; static CE u_t inv(u_t n, int e= 6, u_t x= 1) { return e ? inv(n, e - 1, x * (2 - x * n)) : x; } CE IL u_t reduce(const du_t &w) const { return u_t(w >> B) + mod - ((du_t(u_t(w) * iv) * mod) >> B); } }; struct MP_Na { u32 mod; CE MP_Na(): mod(0){}; CE MP_Na(u32 m): mod(m) {} CE IL u32 mul(u32 l, u32 r) const { return u64(l) * r % mod; } PLUS(u32, mod) DIFF(u32, 31, mod) SGN(u32) }; struct MP_Br { // mod < 2^31 u32 mod; CE MP_Br(): mod(0), s(0), x(0) {} CE MP_Br(u32 m): mod(m), s(95 - __builtin_clz(m - 1)), x(((u128(1) << s) + m - 1) / m) {} CE IL u32 mul(u32 l, u32 r) const { return rem(u64(l) * r); } PLUS(u32, mod) DIFF(u32, 31, mod) SGN(u32) private: u8 s; u64 x; CE IL u64 quo(u64 n) const { return (u128(x) * n) >> s; } CE IL u32 rem(u64 n) const { return n - quo(n) * mod; } }; struct MP_Br2 { // 2^20 < mod <= 2^41 u64 mod; CE MP_Br2(): mod(0), x(0) {} CE MP_Br2(u64 m): mod(m), x((u128(1) << 84) / m) {} CE IL u64 mul(u64 l, u64 r) const { return rem(u128(l) * r); } PLUS(u64, mod << 1) DIFF(u64, 63, mod << 1) static CE IL u64 set(u64 n) { return n; } CE IL u64 get(u64 n) const { NORM; } CE IL u64 norm(u64 n) const { NORM; } private: u64 x; CE IL u128 quo(const u128 &n) const { return (n * x) >> 84; } CE IL u64 rem(const u128 &n) const { return n - quo(n) * mod; } }; struct MP_D2B1 { u8 s; u64 mod, d, v; CE MP_D2B1(): s(0), mod(0), d(0), v(0) {} CE MP_D2B1(u64 m): s(__builtin_clzll(m)), mod(m), d(m << s), v(u128(-1) / d) {} CE IL u64 mul(u64 l, u64 r) const { return rem((u128(l) * r) << s) >> s; } PLUS(u64, mod) DIFF(u64, 63, mod) SGN(u64) private: CE IL u64 rem(const u128 &u) const { u128 q= (u >> 64) * v + u; u64 r= u64(u) - (q >> 64) * d - d; if (r > u64(q)) r+= d; if (r >= d) r-= d; return r; } }; template CE u_t pow(u_t x, u64 k, const MP &md) { for (u_t ret= md.set(1);; x= md.mul(x, x)) if (k & 1 ? ret= md.mul(ret, x) : 0; !(k>>= 1)) return ret; } #undef NORM #undef PLUS #undef DIFF #undef SGN #undef CE } #line 3 "/Users/hashimotoryoma/code/CompetitiveProgramming/Library/Library/src/Math/is_prime.hpp" namespace math_internal { template constexpr bool miller_rabin(Uint n) { const MP md(n); const Uint s= __builtin_ctzll(n - 1), d= n >> s, one= md.set(1), n1= md.norm(md.set(n - 1)); for (auto a: {args...}) if (Uint b= a % n; b) if (Uint p= md.norm(pow(md.set(b), d, md)); p != one) for (int i= s; p != n1; p= md.norm(md.mul(p, p))) if (!(--i)) return 0; return 1; } constexpr bool is_prime(u64 n) { if (n < 2 || n % 6 % 4 != 1) return (n | 1) == 3; if (n < (1 << 30)) return miller_rabin, 2, 7, 61>(n); if (n < (1ull << 62)) return miller_rabin, 2, 325, 9375, 28178, 450775, 9780504, 1795265022>(n); return miller_rabin(n); } } using math_internal::is_prime; #line 4 "/Users/hashimotoryoma/code/CompetitiveProgramming/Library/Library/src/Math/mod_inv.hpp" template constexpr inline Int mod_inv(Int a, Int mod) { static_assert(std::is_signed_v); Int x= 1, y= 0, b= mod; for (Int q= 0, z= 0, c= 0; b;) z= x, c= a, x= y, y= z - y * (q= a / b), a= b, b= c - b * q; return assert(a == 1), x < 0 ? mod - (-x) % mod : x % mod; } #line 4 "/Users/hashimotoryoma/code/CompetitiveProgramming/Library/Library/src/Math/ModInt.hpp" namespace math_internal { #define CE constexpr struct m_b {}; struct s_b: m_b {}; template CE bool is_modint_v= is_base_of_v; template CE bool is_staticmodint_v= is_base_of_v; template struct SB: s_b { protected: static CE MP md= MP(MOD); }; template struct MInt: public B { using Uint= U; static CE inline auto mod() { return B::md.mod; } CE MInt(): x(0) {} CE MInt(const MInt& r): x(r.x) {} template , nullptr_t> = nullptr> CE MInt(T v): x(B::md.set(v.val() % B::md.mod)) {} template , nullptr_t> = nullptr> CE MInt(T n): x(B::md.set((n < 0 ? ((n= (-n) % B::md.mod) ? B::md.mod - n : n) : n % B::md.mod))) {} CE MInt operator-() const { return MInt() - *this; } #define FUNC(name, op) \ CE MInt name const { \ MInt ret; \ ret.x= op; \ return ret; \ } FUNC(operator+(const MInt& r), B::md.plus(x, r.x)) FUNC(operator-(const MInt& r), B::md.diff(x, r.x)) FUNC(operator*(const MInt& r), B::md.mul(x, r.x)) FUNC(pow(u64 k), math_internal::pow(x, k, B::md)) #undef FUNC CE MInt operator/(const MInt& r) const { return *this * r.inv(); } CE MInt& operator+=(const MInt& r) { return *this= *this + r; } CE MInt& operator-=(const MInt& r) { return *this= *this - r; } CE MInt& operator*=(const MInt& r) { return *this= *this * r; } CE MInt& operator/=(const MInt& r) { return *this= *this / r; } CE bool operator==(const MInt& r) const { return B::md.norm(x) == B::md.norm(r.x); } CE bool operator!=(const MInt& r) const { return !(*this == r); } CE bool operator<(const MInt& r) const { return B::md.norm(x) < B::md.norm(r.x); } CE inline MInt inv() const { return mod_inv(val(), B::md.mod); } CE inline Uint val() const { return B::md.get(x); } friend ostream& operator<<(ostream& os, const MInt& r) { return os << r.val(); } friend istream& operator>>(istream& is, MInt& r) { i64 v; return is >> v, r= MInt(v), is; } private: Uint x; }; template using ModInt= conditional_t < (MOD < (1 << 30)) & MOD, MInt, MOD>>, conditional_t < (MOD < (1ull << 62)) & MOD, MInt, MOD>>, conditional_t>, conditional_t>, conditional_t>, MInt>>>>>>; #undef CE } using math_internal::ModInt, math_internal::is_modint_v, math_internal::is_staticmodint_v; template mod_t get_inv(int n) { static_assert(is_modint_v); static const auto m= mod_t::mod(); static mod_t dat[LM]; static int l= 1; if (l == 1) dat[l++]= 1; while (l <= n) dat[l++]= dat[m % l] * (m - m / l); return dat[n]; } #line 6 "/Users/hashimotoryoma/code/CompetitiveProgramming/Library/Library/src/FFT/NTT.hpp" namespace math_internal { #define CE constexpr #define ST static #define TP template #define BSF(_, n) __builtin_ctz##_(n) TP struct NTT { #define _DFT(a, b, c, ...) \ mod_t r, u, *x0, *x1; \ for (int a= n, b= 1, s, i; a>>= 1; b<<= 1) \ for (s= 0, r= I, x0= x;; r*= c[BSF(, s)], x0= x1 + p) { \ for (x1= x0 + (i= p); i--;) __VA_ARGS__; \ if (++s == e) break; \ } ST inline void dft(int n, mod_t x[]) { _DFT(p, e, r2, x1[i]= x0[i] - (u= r * x1[i]), x0[i]+= u); } ST inline void idft(int n, mod_t x[]) { _DFT(e, p, ir2, u= x0[i] - x1[i], x0[i]+= x1[i], x1[i]= r * u) for (const mod_t iv= I / n; n--;) x[n]*= iv; } #undef _DFT ST inline void even_dft(int n, mod_t x[]) { for (int i= 0, j= 0; i < n; i+= 2) x[j++]= iv2 * (x[i] + x[i + 1]); } ST inline void odd_dft(int n, mod_t x[], mod_t r= iv2) { for (int i= 0, j= 0;; r*= ir2[BSF(, ++j)]) if (x[j]= r * (x[i] - x[i + 1]); (i+= 2) == n) break; } ST inline void dft_doubling(int n, mod_t x[], int i= 0) { mod_t k= I, t= rt[BSF(, n << 1)]; for (copy_n(x, n, x + n), idft(n, x + n); i < n; ++i) x[n + i]*= k, k*= t; dft(n, x + n); } protected: ST CE u64 md= mod_t::mod(); static_assert(md & 1); static_assert(is_prime(md)); ST CE u8 E= BSF(ll, md - 1); ST CE mod_t w= [](u8 e) { for (mod_t r= 2;; r+= 1) if (auto s= r.pow((md - 1) / 2); s != 1 && s * s == 1) return r.pow((md - 1) >> e); return mod_t(); }(E); static_assert(w != mod_t()); ST CE mod_t I= 1, iv2= (md + 1) / 2, iw= w.pow((1ULL << E) - 1); ST CE auto roots(mod_t w) { array x= {}; for (u8 e= E; e; w*= w) x[e--]= w; return x[0]= w, x; } TP ST CE auto ras(const array& rt, const array& irt, int i= N) { array x= {}; for (mod_t ro= 1; i <= E; ro*= irt[i++]) x[i - N]= rt[i] * ro; return x; } ST CE auto rt= roots(w), irt= roots(iw); ST CE auto r2= ras<2>(rt, irt), ir2= ras<2>(irt, rt); }; TP struct NI: public B { using B::B; #define FUNC(op, name, HG, ...) \ inline void name(__VA_ARGS__) { \ HG(op, 1); \ if CE (t > 1) HG(op, 2); \ if CE (t > 2) HG(op, 3); \ if CE (t > 3) HG(op, 4); \ if CE (t > 4) HG(op, 5); \ } #define REP for (int i= b; i < e; ++i) #define DFT(fft, _) B::ntt##_::fft(e - b, this->dt##_ + b) #define ZEROS(op, _) fill_n(this->dt##_ + b, e - b, typename B::m##_()) #define SET(op, _) copy(x + b, x + e, this->dt##_ + b) #define SET_S(op, _) this->dt##_[i]= x; #define SUBST(op, _) copy(r.dt##_ + b, r.dt##_ + e, this->dt##_ + b) #define ASGN(op, _) REP this->dt##_[i] op##= r.dt##_[i] #define ASN(nm, op) TP FUNC(op, nm, ASGN, const NI& r, int b, int e) #define BOP(op, _) REP this->dt##_[i]= l.dt##_[i] op r.dt##_[i] #define OP(nm, op) TP FUNC(op, nm, BOP, const NI& l, const NI& r, int b, int e) OP(add, +) OP(dif, -) OP(mul, *) ASN(add, +) ASN(dif, -) ASN(mul, *) FUNC(dft, dft, DFT, int b, int e) FUNC(idft, idft, DFT, int b, int e) FUNC(__, zeros, ZEROS, int b, int e) FUNC(__, set, SET, const T x[], int b, int e) FUNC(__, set, SET_S, int i, T x) TP FUNC(__, subst, SUBST, const NI& r, int b, int e) inline void get(T x[], int b, int e) const { if CE (t == 1) copy(this->dt1 + b, this->dt1 + e, x + b); else REP x[i]= get(i); } #define TMP(_) B::iv##_##1 * (this->dt##_[i] - r1) inline T get(int i) const { if CE (t > 1) { u64 r1= this->dt1[i].val(), r2= (TMP(2)).val(); T a= 0; if CE (t > 2) { u64 r3= (TMP(3) - B::iv32 * r2).val(); if CE (t > 3) { u64 r4= (TMP(4) - B::iv42 * r2 - B::iv43 * r3).val(); if CE (t > 4) a= T(B::m4::mod()) * (TMP(5) - B::iv52 * r2 - B::iv53 * r3 - B::iv54 * r4).val(); a= (a + r4) * B::m3::mod(); } a= (a + r3) * B::m2::mod(); } return (a + r2) * B::m1::mod() + r1; } else return this->dt1[i]; } #undef TMP #undef DFT #undef ZEROS #undef SET #undef SET_S #undef SUBST #undef ASGN #undef ASN #undef BOP #undef OP #undef FUNC #undef REP }; #define ARR(_) \ using m##_= ModInt; \ using ntt##_= NTT; \ m##_ dt##_[LM]= {}; #define IV2 ST CE m2 iv21= m2(1) / m1::mod(); #define IV3 ST CE m3 iv32= m3(1) / m2::mod(), iv31= iv32 / m1::mod(); #define IV4 ST CE m4 iv43= m4(1) / m3::mod(), iv42= iv43 / m2::mod(), iv41= iv42 / m1::mod(); #define IV5 ST CE m5 iv54= m5(1) / m4::mod(), iv53= iv54 / m3::mod(), iv52= iv53 / m2::mod(), iv51= iv52 / m1::mod(); TP struct NB { ARR(1) }; TP struct NB<2, M1, M2, M3, M4, M5, LM, 0> { ARR(1) ARR(2) IV2 }; TP struct NB<3, M1, M2, M3, M4, M5, LM, 0> { ARR(1) ARR(2) ARR(3) IV2 IV3 }; TP struct NB<4, M1, M2, M3, M4, M5, LM, 0> { ARR(1) ARR(2) ARR(3) ARR(4) IV2 IV3 IV4 }; TP struct NB<5, M1, M2, M3, M4, M5, LM, 0> { ARR(1) ARR(2) ARR(3) ARR(4) ARR(5) IV2 IV3 IV4 IV5 }; #undef ARR #define VC(_) \ using m##_= ModInt; \ using ntt##_= NTT; \ vector bf##_; \ m##_* dt##_; #define RS resize TP struct NB<1, M1, M2, M3, M4, M5, LM, 1> { NB(): dt1(bf1.data()) {} void RS(int n) { bf1.RS(n), dt1= bf1.data(); } u32 size() const { return bf1.size(); } VC(1) }; TP struct NB<2, M1, M2, M3, M4, M5, LM, 1> { NB(): dt1(bf1.data()), dt2(bf2.data()) {} void RS(int n) { bf1.RS(n), dt1= bf1.data(), bf2.RS(n), dt2= bf2.data(); } u32 size() const { return bf1.size(); } VC(1) VC(2) IV2 }; TP struct NB<3, M1, M2, M3, M4, M5, LM, 1> { NB(): dt1(bf1.data()), dt2(bf2.data()), dt3(bf3.data()) {} void RS(int n) { bf1.RS(n), dt1= bf1.data(), bf2.RS(n), dt2= bf2.data(), bf3.RS(n), dt3= bf3.data(); } u32 size() const { return bf1.size(); } VC(1) VC(2) VC(3) IV2 IV3 }; TP struct NB<4, M1, M2, M3, M4, M5, LM, 1> { NB(): dt1(bf1.data()), dt2(bf2.data()), dt3(bf3.data()), dt4(bf4.data()) {} void RS(int n) { bf1.RS(n), dt1= bf1.data(), bf2.RS(n), dt2= bf2.data(), bf3.RS(n), dt3= bf3.data(), bf4.RS(n), dt4= bf4.data(); } u32 size() const { return bf1.size(); } VC(1) VC(2) VC(3) VC(4) IV2 IV3 IV4 }; TP struct NB<5, M1, M2, M3, M4, M5, LM, 1> { NB(): dt1(bf1.data()), dt2(bf2.data()), dt3(bf3.data()), dt4(bf4.data()), dt5(bf5.data()) {} void RS(int n) { bf1.RS(n), dt1= bf1.data(), bf2.RS(n), dt2= bf2.data(), bf3.RS(n), dt3= bf3.data(), bf4.RS(n), dt4= bf4.data(), bf5.RS(n), dt5= bf5.data(); } u32 size() const { return bf1.size(); } VC(1) VC(2) VC(3) VC(4) VC(5) IV2 IV3 IV4 IV5 }; #undef VC #undef IV2 #undef IV3 #undef IV4 #undef IV5 TP CE bool is_nttfriend() { if CE (!is_staticmodint_v) return 0; else return (T::mod() & is_prime(T::mod())) && LM <= (1ULL << BSF(ll, T::mod() - 1)); } TP, nullptr_t> = nullptr> CE u64 mv() { return numeric_limits::max(); } TP, nullptr_t> = nullptr> CE u64 mv() { return T::mod(); } TP CE u8 nt() { if CE (!is_nttfriend()) { CE u128 m= mv(), mm= m * m; if CE (mm <= M1 / LM) return 1; else if CE (mm <= u64(M1) * M2 / LM) return 2; else if CE (mm <= u128(M1) * M2 * M3 / LM) return 3; else if CE (mm <= u128(M1) * M2 * M3 * M4 / LM) return 4; else return 5; } else return 1; } #undef BSF #undef RS CE u32 MOD1= 998244353, MOD2= 897581057, MOD3= 880803841, MOD4= 754974721, MOD5= 645922817; TP CE u8 nttarr_type= nt(); TP CE u8 nttarr_cat= is_nttfriend() && (mv() > (1 << 30)) ? 0 : nttarr_type; TP using NTTArray= NI, conditional_t(), NB<1, mv(), 0, 0, 0, 0, LM, v>, NB, MOD1, MOD2, MOD3, MOD4, MOD5, LM, v>>>; #undef CE #undef ST #undef TP } using math_internal::is_nttfriend, math_internal::nttarr_type, math_internal::nttarr_cat, math_internal::NTT, math_internal::NTTArray; template struct GlobalNTTArray { static inline NTTArray bf; }; template struct GlobalNTTArray2D { static inline NTTArray bf[LM2]; }; template struct GlobalArray { static inline T bf[LM]; }; constexpr unsigned pw2(unsigned n) { return --n, n|= n >> 1, n|= n >> 2, n|= n >> 4, n|= n >> 8, n|= n >> 16, ++n; } #line 6 "/Users/hashimotoryoma/code/CompetitiveProgramming/Library/Library/src/FFT/fps_inv.hpp" namespace math_internal { template inline void inv_base(const mod_t p[], int n, mod_t r[], int i= 1, int l= -1) { static constexpr int t= nttarr_cat, TH= (int[]){64, 64, 128, 256, 512, 512}[t]; if (n <= i) return; if (l < 0) l= n; assert(((n & -n) == n)), assert(i && ((i & -i) == i)); const mod_t miv= -r[0]; for (int j, m= min(n, TH); i < m; r[i++]*= miv) for (r[i]= mod_t(), j= min(i + 1, l); --j;) r[i]+= r[i - j] * p[j]; static constexpr int lnR= 2 + (!t), R= (1 << lnR) - 1; using GNA1= GlobalNTTArray; using GNA2= GlobalNTTArray; for (auto gt1= GlobalNTTArray2D::bf, gt2= GlobalNTTArray2D::bf; i < n;) { mod_t* rr= r; const mod_t* pp= p; const int s= i, e= s << 1, ss= (l - 1) / s; for (int k= 0, j; i < n && k < R; ++k, i+= s, pp+= s) { if (j= min(e, l - k * s); j > 0) gt2[k].set(pp, 0, j), gt2[k].zeros(j, e), gt2[k].dft(0, e); for (gt1[k].set(rr, 0, s), gt1[k].zeros(s, e), gt1[k].dft(0, e), GNA2::bf.mul(gt1[k], gt2[0], 0, e), j= min(k, ss) + 1; --j;) GNA1::bf.mul(gt1[k - j], gt2[j], 0, e), GNA2::bf.add(GNA1::bf, 0, e); GNA2::bf.idft(0, e), GNA2::bf.zeros(0, s); if constexpr (!is_nttfriend()) GNA2::bf.get(rr, s, e), GNA2::bf.set(rr, s, e); for (GNA2::bf.dft(0, e), GNA2::bf.mul(gt1[0], 0, e), GNA2::bf.idft(0, e), GNA2::bf.get(rr, s, e), rr+= j= s; j--;) rr[j]= -rr[j]; } } } template void inv_(const mod_t p[], int n, mod_t r[]) { static constexpr u32 R= (1 << lnR) - 1, LM2= LM >> (lnR - 1); using GNA1= GlobalNTTArray; using GNA2= GlobalNTTArray; auto gt1= GlobalNTTArray2D::bf, gt2= GlobalNTTArray2D::bf; assert(n > 0), assert(p[0] != mod_t()); const int m= pw2(n) >> lnR, m2= m << 1, ed= (n - 1) / m; inv_base(p, m, r); for (int k= 0, l; k < ed; p+= m) { for (gt2[k].set(p, 0, l= min(m2, n - m * k)), gt2[k].zeros(l, m2), gt2[k].dft(0, m2), gt1[k].set(r, 0, m), gt1[k].zeros(m, m2), gt1[k].dft(0, m2), GNA2::bf.mul(gt1[k], gt2[0], 0, m2), l= k; l--;) GNA1::bf.mul(gt1[l], gt2[k - l], 0, m2), GNA2::bf.add(GNA1::bf, 0, m2); GNA2::bf.idft(0, m2), GNA2::bf.zeros(0, m); if constexpr (!is_nttfriend()) GNA2::bf.get(r, m, m2), GNA2::bf.set(r, m, m2); for (GNA2::bf.dft(0, m2), GNA2::bf.mul(gt1[0], 0, m2), GNA2::bf.idft(0, m2), GNA2::bf.get(r, m, m + (l= min(m, n - m * ++k))), r+= m; l--;) r[l]= -r[l]; } } template vector inv(const vector& p) { static constexpr int t= nttarr_cat, TH= (int[]){234, 106, 280, 458, 603, 861}[t]; mod_t *pp= GlobalArray::bf, *r= GlobalArray::bf; const int n= p.size(); copy_n(p.begin(), n, pp), assert(n > 0), assert(p[0] != mod_t()); if (const mod_t miv= -(r[0]= mod_t(1) / p[0]); n > TH) { const int l= pw2(n), l1= l >> 1, k= (n - l1 - 1) / (l1 >> 3), bl= __builtin_ctz(l1); int a= 4; if constexpr (!t) a= bl < 8 ? k > 5 ? 1 : 3 : bl < 9 ? k & 1 ? 3 : 4 : bl < 10 ? k & 1 && k > 4 ? 3 : 4 : bl < 11 ? k > 6 ? 3 : 4 : 4; else if constexpr (t < 2) a= bl < 7 ? 2 : bl < 9 ? k ? 3 : 4 : k & 1 ? 3 : 4; else if constexpr (t < 3) a= bl < 9 ? k > 5 ? 1 : k ? 3 : 4 : k & 1 ? 3 : 4; else if constexpr (t < 4) a= bl < 9 ? 1 : bl < 10 ? k > 5 ? 1 : !k ? 4 : k & 2 ? 2 : 3 : k & 1 ? 3 : 4; else if constexpr (t < 5) a= bl < 10 ? k & 2 ? 2 : 3 : k & 1 ? 3 : 4; else a= bl < 10 ? 1 : bl < 11 ? k > 5 ? 1 : !k ? 4 : k & 2 ? 2 : 3 : k & 1 ? 3 : 4; (a < 2 ? inv_<1, mod_t, LM> : a < 3 ? inv_<2, mod_t, LM> : a < 4 ? inv_<3, mod_t, LM> : inv_<4, mod_t, LM>)(pp, n, r); } else for (int j, i= 1; i < n; r[i++]*= miv) for (r[j= i]= mod_t(); j--;) r[i]+= r[j] * pp[i - j]; return vector(r, r + n); } } using math_internal::inv_base, math_internal::inv; #line 5 "/Users/hashimotoryoma/code/CompetitiveProgramming/Library/Library/src/FFT/convolve.hpp" template std::vector convolve(const std::vector& p, const std::vector& q) { mod_t *pp= GlobalArray::bf, *qq= GlobalArray::bf, *rr= GlobalArray::bf; static constexpr int t= nttarr_cat, TH= (int[]){70, 30, 70, 100, 135, 150}[t]; auto f= [](int l) -> int { static constexpr double B[]= {(double[]){8.288, 5.418, 7.070, 9.676, 11.713, 13.374}[t], (double[]){8.252, 6.578, 9.283, 12.810, 13.853, 15.501}[t]}; return std::round(std::pow(l, 0.535) * B[__builtin_ctz(l) & 1]); }; const int n= p.size(), m= q.size(), sz= n + m - 1; if (!n || !m) return std::vector(); if (std::min(n, m) < TH) { std::fill_n(rr, sz, mod_t()), std::copy(p.begin(), p.end(), pp), std::copy(q.begin(), q.end(), qq); for (int i= n; i--;) for (int j= m; j--;) rr[i + j]+= pp[i] * qq[j]; } else { const int rl= pw2(sz), l= pw2(std::max(n, m)), fl= f(l); static constexpr size_t LM2= LM >> 3; static constexpr bool b= nttarr_cat < t; if (b || (l + fl < sz && sz <= (rl >> 3) * 5)) { using GNA1= GlobalNTTArray; using GNA2= GlobalNTTArray; auto gt1= GlobalNTTArray2D::bf, gt2= GlobalNTTArray2D::bf; const int l= rl >> 4, l2= l << 1, nn= (n + l - 1) / l, mm= (m + l - 1) / l, ss= nn + mm - 1; for (int i= 0, k= 0, s; k < n; ++i, k+= l) gt1[i].set(p.data() + k, 0, s= std::min(l, n - k)), gt1[i].zeros(s, l2), gt1[i].dft(0, l2); if (&p != &q) for (int i= 0, k= 0, s; k < m; ++i, k+= l) gt2[i].set(q.data() + k, 0, s= std::min(l, m - k)), gt2[i].zeros(s, l2), gt2[i].dft(0, l2); else for (int i= nn; i--;) gt2[i].subst(gt1[i], 0, l2); GNA2::bf.mul(gt1[0], gt2[0], 0, l2), GNA2::bf.idft(0, l2), GNA2::bf.get(rr, 0, l2); for (int i= 1, k= l, j, ed; i < ss; ++i, k+= l) { for (j= std::max(0, i - nn + 1), ed= std::min(mm - 1, i), GNA2::bf.mul(gt1[i - ed], gt2[ed], 0, l2); j < ed; ++j) GNA1::bf.mul(gt1[i - j], gt2[j], 0, l2), GNA2::bf.add(GNA1::bf, 0, l2); for (GNA2::bf.idft(0, l2), GNA2::bf.get(pp, 0, j= std::min(l, sz - k)); j--;) rr[k + j]+= pp[j]; if (l < sz - k) GNA2::bf.get(rr + k, l, std::min(l2, sz - k)); } } else { using GNA1= GlobalNTTArray; using GNA2= GlobalNTTArray; const int len= sz <= l + fl ? l : rl; if (GNA1::bf.set(p.data(), 0, n), GNA1::bf.zeros(n, len), GNA1::bf.dft(0, len); &p != &q) GNA2::bf.set(q.data(), 0, m), GNA2::bf.zeros(m, len), GNA2::bf.dft(0, len), GNA1::bf.mul(GNA2::bf, 0, len); else GNA1::bf.mul(GNA1::bf, 0, len); if (GNA1::bf.idft(0, len), GNA1::bf.get(rr, 0, std::min(sz, len)); len < sz) { std::copy(p.begin() + len - m + 1, p.end(), pp + len - m + 1), std::copy(q.begin() + len - n + 1, q.end(), qq + len - n + 1); for (int i= len, j; i < sz; rr[i - len]-= rr[i], ++i) for (rr[i]= mod_t(), j= i - m + 1; j < n; ++j) rr[i]+= pp[j] * qq[i - j]; } } } return std::vector(rr, rr + sz); } #line 4 "/Users/hashimotoryoma/code/CompetitiveProgramming/Library/Library/src/FFT/bostan_mori.hpp" namespace div_at_internal { template int deg(const std::vector &p) { const K ZERO= 0; for (int n= p.size() - 1;; n--) if (n < 0 || p[n] != ZERO) return n; } template void div_at_na(std::vector &p, std::vector &q, std::uint64_t &k) { unsigned n= deg(p), nn, j; const unsigned m= deg(q), l= std::max(n, m) + 1; std::vector tmp(l); for (p.resize(l), q.resize(l); k > m; q.swap(p), p.swap(tmp)) { std::fill_n(tmp.begin(), (nn= (n + m - ((n ^ m ^ k) & 1)) >> 1) + 1, mod_t()); for (j= 0; j <= m; j+= 2) for (int i= k & 1; i <= n; i+= 2) tmp[(i + j) >> 1]+= p[i] * q[j]; for (j= 1; j <= m; j+= 2) for (int i= (~k) & 1; i <= n; i+= 2) tmp[(i + j) >> 1]-= p[i] * q[j]; for (std::fill_n(p.begin(), m + 1, mod_t()), j= 2; j <= m; j+= 2) for (int i= j; (i-= 2) >= 0;) p[(i + j) >> 1]+= q[i] * q[j]; for (k>>= 1, n= nn, j= 3; j <= m; j+= 2) for (int i= j; (i-= 2) >= 0;) p[(i + j) >> 1]-= q[i] * q[j]; for (int i= m; i >= 0; i--) p[i]+= p[i]; for (int i= 0; i <= m; i+= 2) p[i]+= q[i] * q[i]; for (int i= 1; i <= m; i+= 2) p[i]-= q[i] * q[i]; } p.resize(n + 1), q.resize(m + 1); } template void div_at_ntt(std::vector &p, std::vector &q, std::uint64_t &k) { static_assert(!is_nttfriend()); using GNA= GlobalNTTArray; using GNA1= GlobalNTTArray; using GNA2= GlobalNTTArray; const unsigned m= deg(q) + 1, offset= std::max(deg(p) + 1, m), len= pw2((offset + m) - 1); for (p.resize(len >> 1); k >= offset; k>>= 1) { GNA::bf.set(p.data(), 0, len >> 1), GNA::bf.zeros(len >> 1, len), GNA1::bf.set(q.data(), 0, m), GNA1::bf.zeros(m, len), GNA2::bf.zeros(m, len); for (int i= m; i--;) GNA2::bf.set(i, i & 1 ? -q[i] : q[i]); GNA::bf.dft(0, len), GNA1::bf.dft(0, len), GNA2::bf.dft(0, len), GNA::bf.mul(GNA2::bf, 0, len), GNA::bf.idft(0, len), GNA1::bf.mul(GNA2::bf, 0, len), GNA1::bf.idft(0, len); for (int i= k & 1; i < len; i+= 2) p[i >> 1]= GNA::bf.get(i); for (int i= m; i--;) q[i]= GNA1::bf.get(i << 1); } } template void div_at_ntt_fast(std::vector &p, std::vector &q, std::uint64_t &k) { static_assert(is_nttfriend()); using ntt= NTT; const unsigned m= deg(q) + 1, offset= std::max(deg(p) + 1, m), len= pw2((offset + m) - 1), len2= len >> 1; p.resize(len), q.resize(len), ntt::dft(len, p.data()), ntt::dft(len, q.data()); while (1) { for (int i= len; i--;) p[i]*= q[i ^ 1]; k & 1 ? ntt::odd_dft(len, p.data()) : ntt::even_dft(len, p.data()); for (int i= 0; i < len; i+= 2) q[i]= q[i + 1]= q[i] * q[i + 1]; ntt::even_dft(len, q.data()); if ((k>>= 1) < offset) break; ntt::dft_doubling(len2, p.data()), ntt::dft_doubling(len2, q.data()); } ntt::idft(len2, p.data()), ntt::idft(len2, q.data()); } } // namespace div_at_internal #define __FPS_DIVAT(Vec) \ template mod_t div_at(Vec p, Vec q, std::uint64_t k) { \ using namespace div_at_internal; \ const int n= deg(p) + 1, m= deg(q) + 1; \ assert(m != 0); \ mod_t ret= 0; \ if (n == 0) return ret; \ if (m == 1) return k <= (std::uint64_t)n ? p[k] / q[0] : ret; \ if (k >= m) { \ if constexpr (is_nttfriend()) m <= 44 ? div_at_na(p, q, k) : div_at_ntt_fast(p, q, k); \ else m <= 340 ? div_at_na(p, q, k) : div_at_ntt(p, q, k); \ } \ p.resize(k + 1, ret), q.resize(k + 1, ret), q= inv(q); \ for (int i= k; i >= 0; i--) ret+= q[i] * p[k - i]; \ return ret; \ } __FPS_DIVAT(std::vector) #ifdef __POLYNOMIAL __FPS_DIVAT(__POLYNOMIAL) #endif // a[n] = c[0] * a[n-1] + c[1] * a[n-2] + ... + c[d-1] * a[n-d] // return a[k] template mod_t linear_recurrence(std::vector c, std::vector a, std::uint64_t k) { const std::size_t d= c.size(); assert(d <= a.size()); for (auto &x: c) x= -x; c.insert(c.begin(), mod_t(1)), a.resize(d); auto p= convolve(c, a); return p.resize(d), div_at(p, c, k); } #line 18 "helper.cpp" using namespace std; signed main() { cin.tie(0); ios::sync_with_stdio(0); constexpr int N= 10000; using Mint= ModInt; string s; cin >> s; int a= stoi(s); long long M, L; cin >> M >> L; FunctionalGraph graph(N); for (int n= 0; n < N; ++n) { auto x= linear_recurrence({n, 1}, {0, 1}, M); if (M & 1) x-= 1; graph.add_edge(n, x.val()); } graph.build(); string ans= to_string(graph.jump(a, L)); ans= string(4 - ans.length(), '0') + ans; cout << ans << '\n'; return 0; }