#line 2 "library/other/template.hpp" #include #line 2 "library/template/macros.hpp" #line 4 "library/template/macros.hpp" #ifndef __COUNTER__ #define __COUNTER__ __LINE__ #endif #define OVERLOAD5(a, b, c, d, e, ...) e #define REP1_0(b, c) REP1_1(b, c) #define REP1_1(b, c) \ for (ll REP_COUNTER_##c = 0; REP_COUNTER_##c < (ll)(b); ++REP_COUNTER_##c) #define REP1(b) REP1_0(b, __COUNTER__) #define REP2(i, b) for (ll i = 0; i < (ll)(b); ++i) #define REP3(i, a, b) for (ll i = (ll)(a); i < (ll)(b); ++i) #define REP4(i, a, b, c) for (ll i = (ll)(a); i < (ll)(b); i += (ll)(c)) #define rep(...) OVERLOAD5(__VA_ARGS__, REP4, REP3, REP2, REP1)(__VA_ARGS__) #define RREP2(i, a) for (ll i = (ll)(a)-1; i >= 0; --i) #define RREP3(i, a, b) for (ll i = (ll)(a)-1; i >= (ll)(b); --i) #define RREP4(i, a, b, c) for (ll i = (ll)(a)-1; i >= (ll)(b); i -= (ll)(c)) #define rrep(...) OVERLOAD5(__VA_ARGS__, RREP4, RREP3, RREP2)(__VA_ARGS__) #define REPS2(i, b) for (ll i = 1; i <= (ll)(b); ++i) #define REPS3(i, a, b) for (ll i = (ll)(a) + 1; i <= (ll)(b); ++i) #define REPS4(i, a, b, c) for (ll i = (ll)(a) + 1; i <= (ll)(b); i += (ll)(c)) #define reps(...) OVERLOAD5(__VA_ARGS__, REPS4, REPS3, REPS2)(__VA_ARGS__) #define RREPS2(i, a) for (ll i = (ll)(a); i > 0; --i) #define RREPS3(i, a, b) for (ll i = (ll)(a); i > (ll)(b); --i) #define RREPS4(i, a, b, c) for (ll i = (ll)(a); i > (ll)(b); i -= (ll)(c)) #define rreps(...) OVERLOAD5(__VA_ARGS__, RREPS4, RREPS3, RREPS2)(__VA_ARGS__) #define each_for(...) for (auto&& __VA_ARGS__) #define each_const(...) for (const auto& __VA_ARGS__) #define all(v) std::begin(v), std::end(v) #if __cplusplus >= 201402L #define rall(v) std::rbegin(v), std::rend(v) #else #define rall(v) v.rbegin(), v.rend() #endif #if __cpp_constexpr >= 201304L #define CONSTEXPR constexpr #else #define CONSTEXPR #endif #if __cpp_if_constexpr >= 201606L #define IF_CONSTEXPR constexpr #else #define IF_CONSTEXPR #endif #define IO_BUFFER_SIZE 2048 #line 2 "library/template/alias.hpp" #line 4 "library/template/alias.hpp" using ll = long long; using uint = unsigned int; using ull = unsigned long long; using i128 = __int128_t; using u128 = __uint128_t; using ld = long double; using PLL = std::pair; template using prique = std::priority_queue, std::greater>; template struct infinity { static constexpr T value = std::numeric_limits::max() / 2; static constexpr T mvalue = std::numeric_limits::lowest() / 2; static constexpr T max = std::numeric_limits::max(); static constexpr T min = std::numeric_limits::lowest(); }; #if __cplusplus <= 201402L template constexpr T infinity::value; template constexpr T infinity::mvalue; template constexpr T infinity::max; template constexpr T infinity::min; #endif #if __cpp_variable_templates >= 201304L template constexpr T INF = infinity::value; #endif constexpr ll inf = infinity::value; constexpr ld EPS = 1e-8; constexpr ld PI = 3.1415926535897932384626; #line 2 "library/template/type_traits.hpp" #line 5 "library/template/type_traits.hpp" template using is_signed_int = std::integral_constant::value && std::is_signed::value) || std::is_same::value>; template using is_unsigned_int = std::integral_constant::value && std::is_unsigned::value) || std::is_same::value>; template using is_int = std::integral_constant::value || is_unsigned_int::value>; template using make_signed_int = typename std::conditional< std::is_same::value || std::is_same::value, std::common_type, std::make_signed>::type; template using make_unsigned_int = typename std::conditional< std::is_same::value || std::is_same::value, std::common_type, std::make_unsigned>::type; template struct is_range : std::false_type {}; template struct is_range< T, decltype(all(std::declval::type>()), (void)0)> : std::true_type {}; template::value> struct range_rank : std::integral_constant {}; template struct range_rank : std::integral_constant::value + 1> {}; template struct int_least { static_assert(size <= 128, "size must be less than or equal to 128"); using type = typename std::conditional< size <= 8, std::int_least8_t, typename std::conditional< size <= 16, std::int_least16_t, typename std::conditional< size <= 32, std::int_least32_t, typename std::conditional::type>::type>::type>::type; }; template using int_least_t = typename int_least::type; template struct uint_least { static_assert(size <= 128, "size must be less than or equal to 128"); using type = typename std::conditional< size <= 8, std::uint_least8_t, typename std::conditional< size <= 16, std::uint_least16_t, typename std::conditional< size <= 32, std::uint_least32_t, typename std::conditional::type>::type>::type>::type; }; template using uint_least_t = typename uint_least::type; template using double_size_int = int_least::digits * 2 + 1>; template using double_size_int_t = typename double_size_int::type; template using double_size_uint = uint_least::digits * 2>; template using double_size_uint_t = typename double_size_uint::type; template using double_size = typename std::conditional::value, double_size_int, double_size_uint>::type; template using double_size_t = typename double_size::type; #line 7 "library/other/template.hpp" #line 2 "library/template/bitop.hpp" #line 6 "library/template/bitop.hpp" namespace bitop { #define KTH_BIT(b, k) (((b) >> (k)) & 1) #define POW2(k) (1ull << (k)) inline ull next_combination(int n, ull x) { if (n == 0) return 1; ull a = x & -x; ull b = x + a; return (x & ~b) / a >> 1 | b; } #define rep_comb(i, n, k) \ for (ull i = (1ull << (k)) - 1; i < (1ull << (n)); \ i = bitop::next_combination((n), i)) inline CONSTEXPR int msb(ull x) { int res = x ? 0 : -1; if (x & 0xFFFFFFFF00000000) x &= 0xFFFFFFFF00000000, res += 32; if (x & 0xFFFF0000FFFF0000) x &= 0xFFFF0000FFFF0000, res += 16; if (x & 0xFF00FF00FF00FF00) x &= 0xFF00FF00FF00FF00, res += 8; if (x & 0xF0F0F0F0F0F0F0F0) x &= 0xF0F0F0F0F0F0F0F0, res += 4; if (x & 0xCCCCCCCCCCCCCCCC) x &= 0xCCCCCCCCCCCCCCCC, res += 2; return res + ((x & 0xAAAAAAAAAAAAAAAA) ? 1 : 0); } inline CONSTEXPR int ceil_log2(ull x) { return x ? msb(x - 1) + 1 : 0; } inline CONSTEXPR ull reverse(ull x) { x = ((x & 0xAAAAAAAAAAAAAAAA) >> 1) | ((x & 0x5555555555555555) << 1); x = ((x & 0xCCCCCCCCCCCCCCCC) >> 2) | ((x & 0x3333333333333333) << 2); x = ((x & 0xF0F0F0F0F0F0F0F0) >> 4) | ((x & 0x0F0F0F0F0F0F0F0F) << 4); x = ((x & 0xFF00FF00FF00FF00) >> 8) | ((x & 0x00FF00FF00FF00FF) << 8); x = ((x & 0xFFFF0000FFFF0000) >> 16) | ((x & 0x0000FFFF0000FFFF) << 16); return (x >> 32) | (x << 32); } inline CONSTEXPR ull reverse(ull x, int n) { return reverse(x) >> (64 - n); } } // namespace bitop inline CONSTEXPR int popcnt(ull x) noexcept { #if __cplusplus >= 202002L return std::popcount(x); #endif x = (x & 0x5555555555555555) + ((x >> 1) & 0x5555555555555555); x = (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333); x = (x & 0x0f0f0f0f0f0f0f0f) + ((x >> 4) & 0x0f0f0f0f0f0f0f0f); x = (x & 0x00ff00ff00ff00ff) + ((x >> 8) & 0x00ff00ff00ff00ff); x = (x & 0x0000ffff0000ffff) + ((x >> 16) & 0x0000ffff0000ffff); return (x & 0x00000000ffffffff) + ((x >> 32) & 0x00000000ffffffff); } #line 2 "library/template/func.hpp" #line 6 "library/template/func.hpp" template> inline constexpr bool chmin(T& a, const U& b, Comp cmp = Comp()) noexcept(noexcept(cmp(b, a))) { return cmp(b, a) ? a = b, true : false; } template> inline constexpr bool chmax(T& a, const U& b, Comp cmp = Comp()) noexcept(noexcept(cmp(a, b))) { return cmp(a, b) ? a = b, true : false; } inline CONSTEXPR ll gcd(ll a, ll b) { if (a < 0) a = -a; if (b < 0) b = -b; while (b) { const ll c = a; a = b; b = c % b; } return a; } inline CONSTEXPR ll lcm(ll a, ll b) { return a / gcd(a, b) * b; } inline CONSTEXPR ll mod_pow(ll a, ll b, ll mod) { assert(mod > 0); if (mod == 1) return 0; a %= mod; ll res = 1; while (b) { if (b & 1) (res *= a) %= mod; b >>= 1; (a *= a) %= mod; } return res; } inline PLL extGCD(ll a, ll b) { const ll n = a, m = b; ll x = 1, y = 0, u = 0, v = 1; ll t; while (b) { t = a / b; std::swap(a -= t * b, b); std::swap(x -= t * u, u); std::swap(y -= t * v, v); } if (x < 0) { x += m; y -= n; } return {x, y}; } inline ll mod_inv(ll a, ll mod) { ll b = mod; ll x = 1, u = 0; ll t; while (b) { t = a / b; std::swap(a -= t * b, b); std::swap(x -= t * u, u); } if (x < 0) x += mod; assert(a == 1); return x; } #line 2 "library/template/util.hpp" #line 6 "library/template/util.hpp" template struct multi_dim_vector { using type = std::vector::type>; }; template struct multi_dim_vector { using type = T; }; template constexpr std::vector make_vec(int n, Arg&& arg) { return std::vector(n, std::forward(arg)); } template constexpr typename multi_dim_vector::type make_vec(int n, Args&&... args) { return typename multi_dim_vector::type( n, make_vec(std::forward(args)...)); } #line 2 "library/math/Combinatorics.hpp" #line 2 "library/math/ModInt.hpp" #line 4 "library/math/ModInt.hpp" template class StaticModInt { static_assert(std::is_integral::value, "T must be integral"); static_assert(std::is_unsigned::value, "T must be unsigned"); static_assert(mod > 0, "mod must be positive"); static_assert(mod <= std::numeric_limits::max() / 2, "mod * 2 must be less than or equal to T::max()"); private: using large_t = typename double_size_uint::type; using signed_t = typename std::make_signed::type; T val; static constexpr unsigned int inv1000000007[] = { 0, 1, 500000004, 333333336, 250000002, 400000003, 166666668, 142857144, 125000001, 111111112, 700000005}; static constexpr unsigned int inv998244353[] = { 0, 1, 499122177, 332748118, 748683265, 598946612, 166374059, 855638017, 873463809, 443664157, 299473306}; public: constexpr StaticModInt() : val(0) {} template::value && std::is_signed::value>::type* = nullptr> constexpr StaticModInt(U v) : val{} { v %= static_cast(mod); if (v < 0) v += static_cast(mod); val = static_cast(v); } template::value && std::is_unsigned::value>::type* = nullptr> constexpr StaticModInt(U v) : val(v % mod) {} T get() const { return val; } static constexpr T get_mod() { return mod; } static StaticModInt raw(T v) { StaticModInt res; res.val = v; return res; } StaticModInt inv() const { if IF_CONSTEXPR (mod == 1000000007) { if (val <= 10) return inv1000000007[val]; } else if IF_CONSTEXPR (mod == 998244353) { if (val <= 10) return inv998244353[val]; } return mod_inv(val, mod); } StaticModInt& operator++() { ++val; if (val == mod) val = 0; return *this; } StaticModInt operator++(int) { StaticModInt res = *this; ++*this; return res; } StaticModInt& operator--() { if (val == 0) val = mod; --val; return *this; } StaticModInt operator--(int) { StaticModInt res = *this; --*this; return res; } StaticModInt& operator+=(const StaticModInt& other) { val += other.val; if (val >= mod) val -= mod; return *this; } StaticModInt& operator-=(const StaticModInt& other) { if (val < other.val) val += mod; val -= other.val; return *this; } StaticModInt& operator*=(const StaticModInt& other) { large_t a = val; a *= other.val; a %= mod; val = a; return *this; } StaticModInt& operator/=(const StaticModInt& other) { *this *= other.inv(); return *this; } friend StaticModInt operator+(const StaticModInt& lhs, const StaticModInt& rhs) { return StaticModInt(lhs) += rhs; } friend StaticModInt operator-(const StaticModInt& lhs, const StaticModInt& rhs) { return StaticModInt(lhs) -= rhs; } friend StaticModInt operator*(const StaticModInt& lhs, const StaticModInt& rhs) { return StaticModInt(lhs) *= rhs; } friend StaticModInt operator/(const StaticModInt& lhs, const StaticModInt& rhs) { return StaticModInt(lhs) /= rhs; } StaticModInt operator+() const { return StaticModInt(*this); } StaticModInt operator-() const { return StaticModInt() - *this; } friend bool operator==(const StaticModInt& lhs, const StaticModInt& rhs) { return lhs.val == rhs.val; } friend bool operator!=(const StaticModInt& lhs, const StaticModInt& rhs) { return lhs.val != rhs.val; } StaticModInt pow(ll a) const { StaticModInt v = *this, res = 1; while (a) { if (a & 1) res *= v; a >>= 1; v *= v; } return res; } template void print(Pr& a) const { a.print(val); } template void debug(Pr& a) const { a.print(val); } template void scan(Sc& a) { ll v; a.scan(v); *this = v; } }; #if __cplusplus < 201703L template constexpr unsigned int StaticModInt::inv1000000007[]; template constexpr unsigned int StaticModInt::inv998244353[]; #endif template using static_modint = StaticModInt; using modint1000000007 = static_modint<1000000007>; using modint998244353 = static_modint<998244353>; template class DynamicModInt { static_assert(std::is_integral::value, "T must be integral"); static_assert(std::is_unsigned::value, "T must be unsigned"); private: using large_t = typename double_size_uint::type; using signed_t = typename std::make_signed::type; T val; static T mod; public: constexpr DynamicModInt() : val(0) {} template::value && std::is_signed::value>::type* = nullptr> constexpr DynamicModInt(U v) : val{} { v %= static_cast(mod); if (v < 0) v += static_cast(mod); val = static_cast(v); } template::value && std::is_unsigned::value>::type* = nullptr> constexpr DynamicModInt(U v) : val(v % mod) {} T get() const { return val; } static T get_mod() { return mod; } static void set_mod(T v) { assert(v > 0); assert(v <= std::numeric_limits::max() / 2); mod = v; } static DynamicModInt raw(T v) { DynamicModInt res; res.val = v; return res; } DynamicModInt inv() const { return mod_inv(val, mod); } DynamicModInt& operator++() { ++val; if (val == mod) val = 0; return *this; } DynamicModInt operator++(int) { DynamicModInt res = *this; ++*this; return res; } DynamicModInt& operator--() { if (val == 0) val = mod; --val; return *this; } DynamicModInt operator--(int) { DynamicModInt res = *this; --*this; return res; } DynamicModInt& operator+=(const DynamicModInt& other) { val += other.val; if (val >= mod) val -= mod; return *this; } DynamicModInt& operator-=(const DynamicModInt& other) { if (val < other.val) val += mod; val -= other.val; return *this; } DynamicModInt& operator*=(const DynamicModInt& other) { large_t a = val; a *= other.val; a %= mod; val = a; return *this; } DynamicModInt& operator/=(const DynamicModInt& other) { *this *= other.inv(); return *this; } friend DynamicModInt operator+(const DynamicModInt& lhs, const DynamicModInt& rhs) { return DynamicModInt(lhs) += rhs; } friend DynamicModInt operator-(const DynamicModInt& lhs, const DynamicModInt& rhs) { return DynamicModInt(lhs) -= rhs; } friend DynamicModInt operator*(const DynamicModInt& lhs, const DynamicModInt& rhs) { return DynamicModInt(lhs) *= rhs; } friend DynamicModInt operator/(const DynamicModInt& lhs, const DynamicModInt& rhs) { return DynamicModInt(lhs) /= rhs; } DynamicModInt operator+() const { return DynamicModInt(*this); } DynamicModInt operator-() const { return DynamicModInt() - *this; } friend bool operator==(const DynamicModInt& lhs, const DynamicModInt& rhs) { return lhs.val == rhs.val; } friend bool operator!=(const DynamicModInt& lhs, const DynamicModInt& rhs) { return lhs.val != rhs.val; } DynamicModInt pow(ll a) const { DynamicModInt v = *this, res = 1; while (a) { if (a & 1) res *= v; a >>= 1; v *= v; } return res; } template void print(Pr& a) const { a.print(val); } template void debug(Pr& a) const { a.print(val); } template void scan(Sc& a) { ll v; a.scan(v); *this = v; } }; template T DynamicModInt::mod = 998244353; template using dynamic_modint = DynamicModInt; using modint = dynamic_modint<-1>; /** * @brief ModInt * @docs docs/math/ModInt.md */ #line 5 "library/math/Combinatorics.hpp" template class Combinatorics { private: static std::vector factorial; static std::vector factinv; public: static void init(ll n) { const int b = factorial.size(); if (n < b) return; factorial.resize(n + 1); rep (i, b, n + 1) factorial[i] = factorial[i - 1] * i; factinv.resize(n + 1); factinv[n] = T(1) / factorial[n]; rreps (i, n, b) factinv[i - 1] = factinv[i] * i; } static T fact(ll x) { if (x < 0) return 0; init(x); return factorial[x]; } static T finv(ll x) { if (x < 0) return 0; init(x); return factinv[x]; } static T inv(ll x) { if (x <= 0) return 0; init(x); return factorial[x - 1] * factinv[x]; } static T perm(ll n, ll r) { if (r < 0 || r > n) return 0; init(n); return factorial[n] * factinv[n - r]; } static T comb(ll n, ll r) { if (n < 0) return 0; if (r < 0 || r > n) return 0; init(n); return factorial[n] * factinv[n - r] * factinv[r]; } static T homo(ll n, ll r) { return comb(n + r - 1, r); } static T small_perm(ll n, ll r) { if (r < 0 || r > n) return 0; T res = 1; reps (i, r) res *= n - r + i; return res; } static T small_comb(ll n, ll r) { if (r < 0 || r > n) return 0; chmin(r, n - r); init(r); T res = factinv[r]; reps (i, r) res *= n - r + i; return res; } static T small_homo(ll n, ll r) { return small_comb(n + r - 1, r); } }; template std::vector Combinatorics::factorial = std::vector(1, 1); template std::vector Combinatorics::factinv = std::vector(1, 1); /** * @brief Combinatorics * @docs docs/math/Combinatorics.md */ #line 2 "library/math/poly/FormalPowerSeries.hpp" #line 2 "library/math/convolution/Convolution.hpp" #line 2 "library/math/PrimitiveRoot.hpp" #line 2 "library/random/Random.hpp" #line 4 "library/random/Random.hpp" template class Random { private: Engine rnd; public: using result_type = typename Engine::result_type; Random() : Random(std::random_device{}()) {} Random(result_type seed) : rnd(seed) {} result_type operator()() { return rnd(); } template IntType uniform(IntType l, IntType r) { static_assert(std::is_integral::value, "template argument must be an integral type"); assert(l <= r); return std::uniform_int_distribution{l, r}(rnd); } template RealType uniform_real(RealType l, RealType r) { static_assert(std::is_floating_point::value, "template argument must be an floating point type"); assert(l <= r); return std::uniform_real_distribution{l, r}(rnd); } bool uniform_bool() { return uniform(0, 1) == 1; } template std::pair uniform_pair(T l, T r) { assert(l < r); T a, b; do { a = uniform(l, r); b = uniform(l, r); } while (a == b); if (a > b) swap(a, b); return {a, b}; } template std::vector choice(int n, T l, T r) { assert(l <= r); assert(T(n) <= (r - l + 1)); std::set res; while ((int)res.size() < n) res.insert(uniform(l, r)); return {res.begin(), res.end()}; } template void shuffle(const Iter& first, const Iter& last) { std::shuffle(first, last, rnd); } template std::vector permutation(T n) { std::vector res(n); rep (i, n) res[i] = i; shuffle(all(res)); return res; } template std::vector choice_shuffle(int n, T l, T r, bool sorted = true) { assert(l <= r); assert(T(n) <= (r - l + 1)); std::vector res(r - l + 1); rep (i, l, r + 1) res[i - l] = i; shuffle(all(res)); res.erase(res.begin() + n, res.end()); if (sorted) sort(all(res)); return res; } }; using Random32 = Random; Random32 rand32; using Random64 = Random; Random64 rand64; /** * @brief Random * @docs docs/random/Random.md */ #line 2 "library/math/MontgomeryModInt.hpp" #line 4 "library/math/MontgomeryModInt.hpp" template class MontgomeryReduction { static_assert(std::is_integral::value, "T must be integral"); static_assert(std::is_unsigned::value, "T must be unsigned"); private: using large_t = typename double_size_uint::type; static constexpr int lg = std::numeric_limits::digits; T mod; T r; T r2; // r^2 mod m T calc_minv() { T t = 0, res = 0; rep (i, lg) { if (~t & 1) { t += mod; res += static_cast(1) << i; } t >>= 1; } return res; } T minv; public: MontgomeryReduction(T v) { set_mod(v); } static constexpr int get_lg() { return lg; } void set_mod(T v) { assert(v > 0); assert(v & 1); assert(v <= std::numeric_limits::max() / 2); mod = v; r = (-static_cast(mod)) % mod; r2 = (-static_cast(mod)) % mod; minv = calc_minv(); } inline T get_mod() const { return mod; } inline T get_r() const { return r; } T reduce(large_t x) const { large_t tmp = (x + static_cast(static_cast(x) * minv) * mod) >> lg; return tmp >= mod ? tmp - mod : tmp; } T transform(large_t x) const { return reduce(x * r2); } }; template class MontgomeryModInt { private: using large_t = typename double_size_uint::type; using signed_t = typename std::make_signed::type; T val; static MontgomeryReduction mont; public: MontgomeryModInt() : val(0) {} template::value && std::is_unsigned::value>::type* = nullptr> MontgomeryModInt(U x) : val(mont.transform( x < (static_cast(mont.get_mod()) << mont.get_lg()) ? x : x % mont.get_mod())) {} template::value && std::is_signed::value>::type* = nullptr> MontgomeryModInt(U x) : MontgomeryModInt(static_cast::type>( x < 0 ? -x : x)) { if (x < 0 && val) val = mont.get_mod() - val; } T get() const { return mont.reduce(val); } static T get_mod() { return mont.get_mod(); } static void set_mod(T v) { mont.set_mod(v); } MontgomeryModInt operator+() const { return *this; } MontgomeryModInt operator-() const { MontgomeryModInt res; if (val) res.val = mont.get_mod() - val; return res; } MontgomeryModInt& operator++() { val += mont.get_r(); if (val >= mont.get_mod()) val -= mont.get_mod(); return *this; } MontgomeryModInt& operator--() { if (val < mont.get_r()) val += mont.get_mod(); val -= mont.get_r(); return *this; } MontgomeryModInt operator++(int) { MontgomeryModInt res = *this; ++*this; return res; } MontgomeryModInt operator--(int) { MontgomeryModInt res = *this; --*this; return res; } MontgomeryModInt& operator+=(const MontgomeryModInt& rhs) { val += rhs.val; if (val >= mont.get_mod()) val -= mont.get_mod(); return *this; } MontgomeryModInt& operator-=(const MontgomeryModInt& rhs) { if (val < rhs.val) val += mont.get_mod(); val -= rhs.val; return *this; } MontgomeryModInt& operator*=(const MontgomeryModInt& rhs) { val = mont.reduce(static_cast(val) * rhs.val); return *this; } MontgomeryModInt pow(ull n) const { MontgomeryModInt res = 1, x = *this; while (n) { if (n & 1) res *= x; x *= x; n >>= 1; } return res; } MontgomeryModInt inv() const { return pow(mont.get_mod() - 2); } MontgomeryModInt& operator/=(const MontgomeryModInt& rhs) { return *this *= rhs.inv(); } friend MontgomeryModInt operator+(const MontgomeryModInt& lhs, const MontgomeryModInt& rhs) { return MontgomeryModInt(lhs) += rhs; } friend MontgomeryModInt operator-(const MontgomeryModInt& lhs, const MontgomeryModInt& rhs) { return MontgomeryModInt(lhs) -= rhs; } friend MontgomeryModInt operator*(const MontgomeryModInt& lhs, const MontgomeryModInt& rhs) { return MontgomeryModInt(lhs) *= rhs; } friend MontgomeryModInt operator/(const MontgomeryModInt& lhs, const MontgomeryModInt& rhs) { return MontgomeryModInt(lhs) /= rhs; } friend bool operator==(const MontgomeryModInt& lhs, const MontgomeryModInt& rhs) { return lhs.val == rhs.val; } friend bool operator!=(const MontgomeryModInt& lhs, const MontgomeryModInt& rhs) { return lhs.val != rhs.val; } template void print(Pr& a) const { a.print(mont.reduce(val)); } template void debug(Pr& a) const { a.print(mont.reduce(val)); } template void scan(Sc& a) { ll v; a.scan(v); *this = v; } }; template MontgomeryReduction MontgomeryModInt::mont = MontgomeryReduction(998244353); using mmodint = MontgomeryModInt; /** * @brief MontgomeryModInt(モンゴメリ乗算) * @docs docs/math/MontgomeryModInt.md */ #line 2 "library/math/MillerRabin.hpp" #line 5 "library/math/MillerRabin.hpp" constexpr ull base_miller_rabin_int[3] = {2, 7, 61}; constexpr ull base_miller_rabin_ll[7] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022}; template CONSTEXPR bool miller_rabin(ull n, const ull base[], int s) { if (T::get_mod() != n) T::set_mod(n); ull d = n - 1; while (~d & 1) d >>= 1; T e{1}, re{n - 1}; rep (i, s) { ull a = base[i]; if (a >= n) return true; ull t = d; T y = T(a).pow(t); while (t != n - 1 && y != e && y != re) { y *= y; t <<= 1; } if (y != re && !(t & 1)) return false; } return true; } CONSTEXPR bool is_prime_mr(ll n) { if (n == 2) return true; if (n < 2 || n % 2 == 0) return false; if (n < (1u << 31)) return miller_rabin>( n, base_miller_rabin_int, 3); return miller_rabin>(n, base_miller_rabin_ll, 7); } #if __cpp_variable_templates >= 201304L && __cpp_constexpr >= 201304L template constexpr bool is_prime_v = is_prime_mr(n); #endif /** * @brief MillerRabin(ミラーラビン素数判定) * @docs docs/math/MillerRabin.md */ #line 2 "library/math/PollardRho.hpp" #line 2 "library/string/RunLength.hpp" #line 4 "library/string/RunLength.hpp" template std::vector> RunLength(const Cont& str, const Comp& cmp) { std::vector> res; if (str.size() == 0) return res; res.emplace_back(str[0], 1); rep (i, 1, str.size()) { if (cmp(res.back().first, str[i])) ++res.back().second; else res.emplace_back(str[i], 1); } return res; } template std::vector> RunLength(const Cont& str) { return RunLength(str, std::equal_to()); } /** * @brief RunLength(ランレングス圧縮) * @docs docs/string/RunLength.md */ #line 8 "library/math/PollardRho.hpp" template ull pollard_rho(ull n, Rnd& rnd) { if (~n & 1) return 2; if (T::get_mod() != n) T::set_mod(n); T c, one = 1; auto f = [&](T x) -> T { return x * x + c; }; constexpr int M = 128; while (1) { c = rnd.uniform(1ull, n - 1); T x = rnd.uniform(2ull, n - 1), y = x; ull g = 1; while (g == 1) { T p = one, tx = x, ty = y; rep (M) { x = f(x); y = f(f(y)); p *= x - y; } g = gcd(p.get(), n); if (g == 1) continue; rep (M) { tx = f(tx); ty = f(f(ty)); g = gcd((tx - ty).get(), n); if (g != 1) { if (g != n) return g; break; } } } } return -1; } template, class Rnd = Random64> std::vector factorize(ull n, Rnd& rnd = rand64) { if (n == 1) return {}; std::vector res; std::vector st = {n}; while (!st.empty()) { ull t = st.back(); st.pop_back(); if (t == 1) continue; if (is_prime_mr(t)) { res.push_back(t); continue; } ull f = pollard_rho(t, rnd); st.push_back(f); st.push_back(t / f); } std::sort(all(res)); return res; } template, class Rnd = Random64> std::vector> expfactorize(ull n, Rnd& rnd = rand64) { auto f = factorize(n, rnd); return RunLength(f); } /** * @brief PollardRho(素因数分解) * @docs docs/math/PollardRho.md */ #line 9 "library/math/PrimitiveRoot.hpp" template> ull primitive_root(ull p) { assert(is_prime_mr(p)); if (p == 2) return 1; if (T::get_mod() != p) T::set_mod(p); auto pf = factorize(p - 1); pf.erase(std::unique(all(pf)), pf.end()); each_for (x : pf) x = (p - 1) / x; T one = 1; while (1) { ull g = rand64.uniform(2ull, p - 1); bool ok = true; each_const (x : pf) { if (T(g).pow(x) == one) { ok = false; break; } } if (ok) return g; } } CONSTEXPR ull primitive_root_for_convolution(ull p) { if (p == 2) return 1; if (p == 998244353) return 3; if (p == 469762049) return 3; if (p == 1811939329) return 11; if (p == 2013265921) return 11; rep (g, 2, p) { if (mod_pow(g, (p - 1) >> 1, p) != 1) return g; } return -1; } /** * @brief PrimitiveRoot(原始根) * @docs docs/math/PrimitiveRoot.md */ #line 6 "library/math/convolution/Convolution.hpp" namespace internal { template class NthRoot { private: static constexpr unsigned int lg = bitop::msb((p - 1) & (1 - p)); unsigned int root[lg + 1]; unsigned int inv_root[lg + 1]; unsigned int rate[lg + 1]; unsigned int inv_rate[lg + 1]; public: constexpr NthRoot() : root{}, inv_root{}, rate{}, inv_rate{} { root[lg] = mod_pow(primitive_root_for_convolution(p), (p - 1) >> lg, p); inv_root[lg] = mod_pow(root[lg], p - 2, p); rrep (i, lg) { root[i] = (ull)root[i + 1] * root[i + 1] % p; inv_root[i] = (ull)inv_root[i + 1] * inv_root[i + 1] % p; } ull r = 1; rep (i, 2, lg + 1) { rate[i - 2] = r * root[i] % p; r = r * inv_root[i] % p; } r = 1; rep (i, 2, lg + 1) { inv_rate[i - 2] = r * inv_root[i] % p; r = r * root[i] % p; } } static constexpr unsigned int get_lg() { return lg; } constexpr unsigned int get(int n) const { return root[n]; } constexpr unsigned int inv(int n) const { return inv_root[n]; } constexpr unsigned int get_rate(int n) const { return rate[n]; } constexpr unsigned int get_inv_rate(int n) const { return inv_rate[n]; } }; template constexpr NthRoot

nth_root; template void number_theoretic_transform(std::vector& a) { int n = a.size(); int lg = bitop::msb(n - 1) + 1; rrep (i, lg) { T z = T(1); rep (j, 1 << (lg - i - 1)) { int offset = j << (i + 1); rep (k, 1 << i) { T x = a[offset + k]; T y = a[offset + k + (1 << i)] * z; a[offset + k] = x + y; a[offset + k + (1 << i)] = x - y; } if (j != (1 << (lg - i - 1)) - 1) { z *= nth_root.get_rate(popcnt(j & ~(j + 1))); } } } } template void inverse_number_theoretic_transform(std::vector& a) { int n = a.size(); int lg = bitop::msb(n - 1) + 1; rep (i, lg) { T z = T(1); rep (j, 1 << (lg - i - 1)) { int offset = j << (i + 1); rep (k, 1 << i) { T x = a[offset + k]; T y = a[offset + k + (1 << i)]; a[offset + k] = x + y; a[offset + k + (1 << i)] = (x - y) * z; } if (j != (1 << (lg - i - 1)) - 1) { z *= nth_root.get_inv_rate(popcnt(j & ~(j + 1))); } } } T inv_n = T(1) / n; each_for (x : a) x *= inv_n; } template std::vector convolution_naive(const std::vector& a, const std::vector& b) { int n = a.size(), m = b.size(); std::vector c(n + m - 1); rep (i, n) rep (j, m) c[i + j] += a[i] * b[j]; return c; } template std::vector convolution_pow2(std::vector a) { int n = a.size() * 2 - 1; int lg = bitop::msb(n - 1) + 1; if (n - (1 << (lg - 1)) <= 5) { --lg; int m = a.size() - (1 << (lg - 1)); std::vector a1(a.begin(), a.begin() + m), a2(a.begin() + m, a.end()); std::vector c(n); std::vector c1 = convolution_naive(a1, a1); std::vector c2 = convolution_naive(a1, a2); std::vector c3 = convolution_pow2(a2); rep (i, c1.size()) c[i] += c1[i]; rep (i, c2.size()) c[i + m] += c2[i] * 2; rep (i, c3.size()) c[i + m * 2] += c3[i]; return c; } int m = 1 << lg; a.resize(m); number_theoretic_transform(a); rep (i, m) a[i] *= a[i]; inverse_number_theoretic_transform(a); a.resize(n); return a; } template std::vector convolution(std::vector a, std::vector b) { int n = a.size() + b.size() - 1; int lg = bitop::msb(n - 1) + 1; int m = 1 << lg; if (n - (1 << (lg - 1)) <= 5) { --lg; if (a.size() < b.size()) std::swap(a, b); int m = n - (1 << lg); std::vector a1(a.begin(), a.begin() + m), a2(a.begin() + m, a.end()); std::vector c(n); std::vector c1 = convolution_naive(a1, b); std::vector c2 = convolution(a2, b); rep (i, c1.size()) c[i] += c1[i]; rep (i, c2.size()) c[i + m] += c2[i]; return c; } a.resize(m); b.resize(m); number_theoretic_transform(a); number_theoretic_transform(b); rep (i, m) a[i] *= b[i]; inverse_number_theoretic_transform(a); a.resize(n); return a; } } // namespace internal using internal::inverse_number_theoretic_transform; using internal::number_theoretic_transform; template std::vector> convolution_for_any_mod(const std::vector>& a, const std::vector>& b); template std::vector> convolution(const std::vector>& a, const std::vector>& b) { unsigned int n = a.size(), m = b.size(); if (n == 0 || m == 0) return {}; if (n <= 60 || m <= 60) return internal::convolution_naive(a, b); if (n + m - 1 > ((1 - p) & (p - 1))) return convolution_for_any_mod(a, b); if (n == m && a == b) return internal::convolution_pow2(a); return internal::convolution(a, b); } template std::vector convolution(const std::vector& a, const std::vector& b) { int n = a.size(), m = b.size(); std::vector> a2(n), b2(m); rep (i, n) a2[i] = a[i]; rep (i, m) b2[i] = b[i]; auto c2 = convolution(a2, b2); std::vector c(n + m - 1); rep (i, n + m - 1) c[i] = c2[i].get(); return c; } template std::vector> convolution_for_any_mod(const std::vector>& a, const std::vector>& b) { int n = a.size(), m = b.size(); assert(n + m - 1 <= (1 << 26)); std::vector a2(n), b2(m); rep (i, n) a2[i] = a[i].get(); rep (i, m) b2[i] = b[i].get(); static constexpr ll MOD1 = 469762049; static constexpr ll MOD2 = 1811939329; static constexpr ll MOD3 = 2013265921; static constexpr ll INV1_2 = mod_pow(MOD1, MOD2 - 2, MOD2); static constexpr ll INV1_3 = mod_pow(MOD1, MOD3 - 2, MOD3); static constexpr ll INV2_3 = mod_pow(MOD2, MOD3 - 2, MOD3); auto c1 = convolution(a2, b2); auto c2 = convolution(a2, b2); auto c3 = convolution(a2, b2); std::vector> res(n + m - 1); rep (i, n + m - 1) { ll t1 = c1[i]; ll t2 = (c2[i] - t1 + MOD2) * INV1_2 % MOD2; if (t2 < 0) t2 += MOD2; ll t3 = ((c3[i] - t1 + MOD3) * INV1_3 % MOD3 - t2 + MOD3) * INV2_3 % MOD3; if (t3 < 0) t3 += MOD3; res[i] = static_modint

(t1 + (t2 + t3 * MOD2) % p * MOD1); } return res; } template void ntt_doubling_(std::vector& a) { int n = a.size(); auto b = a; inverse_number_theoretic_transform(b); const T z = internal::nth_root.get(bitop::msb(n) + 1); T r = 1; rep (i, n) { b[i] *= r; r *= z; } number_theoretic_transform(b); std::copy(all(b), std::back_inserter(a)); } template struct is_ntt_friendly : std::false_type {}; template<> struct is_ntt_friendly<998244353> : std::true_type {}; /** * @brief Convolution(畳み込み) * @docs docs/math/convolution/Convolution.md */ #line 2 "library/math/SqrtMod.hpp" #line 5 "library/math/SqrtMod.hpp" template ll sqrt_mod(ll a) { const ll p = T::get_mod(); if (p == 2) return a; if (a == 0) return 0; if (T{a}.pow((p - 1) >> 1) != 1) return -1; T b = 2; while (T{b}.pow((p - 1) >> 1) == 1) ++b; ll s = 0, t = p - 1; while ((t & 1) == 0) t >>= 1, ++s; T x = T{a}.pow((t + 1) >> 1); T w = T{a}.pow(t); T v = T{b}.pow(t); while (w != 1) { ll k = 0; T y = w; while (y != 1) { y *= y; ++k; } T z = v; rep (s - k - 1) z *= z; x *= z; w *= z * z; } return std::min(x.get(), p - x.get()); } ll sqrt_mod(ll a, ll p) { if (p == 2) return a; using mint = MontgomeryModInt; mint::set_mod(p); return sqrt_mod(a); } /** * @brief SqrtMod(平方剰余) * @docs docs/math/SqrtMod.md * @see https://37zigen.com/tonelli-shanks-algorithm/ */ #line 7 "library/math/poly/FormalPowerSeries.hpp" template class FormalPowerSeries : public std::vector { private: using Base = std::vector; using Comb = Combinatorics; public: using Base::Base; FormalPowerSeries(const Base& v) : Base(v) {} FormalPowerSeries(Base&& v) : Base(std::move(v)) {} FormalPowerSeries& shrink() { while (!this->empty() && this->back() == T{0}) this->pop_back(); return *this; } T eval(T x) const { T res = 0; rrep (i, this->size()) { res *= x; res += (*this)[i]; } return res; } FormalPowerSeries prefix(int deg) const { assert(0 <= deg); if (deg < (int)this->size()) { return FormalPowerSeries(this->begin(), this->begin() + deg); } FormalPowerSeries res(*this); res.resize(deg); return res; } FormalPowerSeries operator+() const { return *this; } FormalPowerSeries operator-() const { FormalPowerSeries res(this->size()); rep (i, this->size()) res[i] = -(*this)[i]; return res; } FormalPowerSeries& operator<<=(int n) { this->insert(this->begin(), n, T{0}); return *this; } FormalPowerSeries& operator>>=(int n) { this->erase(this->begin(), this->begin() + std::min(n, (int)this->size())); return *this; } friend FormalPowerSeries operator<<(const FormalPowerSeries& lhs, int rhs) { return FormalPowerSeries(lhs) <<= rhs; } friend FormalPowerSeries operator>>(const FormalPowerSeries& lhs, int rhs) { return FormalPowerSeries(lhs) >>= rhs; } FormalPowerSeries& operator+=(const FormalPowerSeries& rhs) { if (this->size() < rhs.size()) this->resize(rhs.size()); rep (i, rhs.size()) (*this)[i] += rhs[i]; return *this; } FormalPowerSeries& operator-=(const FormalPowerSeries& rhs) { if (this->size() < rhs.size()) this->resize(rhs.size()); rep (i, rhs.size()) (*this)[i] -= rhs[i]; return *this; } friend FormalPowerSeries operator+(const FormalPowerSeries& lhs, const FormalPowerSeries& rhs) { return FormalPowerSeries(lhs) += rhs; } friend FormalPowerSeries operator-(const FormalPowerSeries& lhs, const FormalPowerSeries& rhs) { return FormalPowerSeries(lhs) -= rhs; } friend FormalPowerSeries operator*(const FormalPowerSeries& lhs, const FormalPowerSeries& rhs) { return FormalPowerSeries(convolution(lhs, rhs)); } FormalPowerSeries& operator*=(const FormalPowerSeries& rhs) { return *this = *this * rhs; } FormalPowerSeries& operator*=(const T& rhs) { rep (i, this->size()) (*this)[i] *= rhs; return *this; } friend FormalPowerSeries operator*(const FormalPowerSeries& lhs, const T& rhs) { return FormalPowerSeries(lhs) *= rhs; } friend FormalPowerSeries operator*(const T& lhs, const FormalPowerSeries& rhs) { return FormalPowerSeries(rhs) *= lhs; } FormalPowerSeries& operator/=(const T& rhs) { rep (i, this->size()) (*this)[i] /= rhs; return *this; } friend FormalPowerSeries operator/(const FormalPowerSeries& lhs, const T& rhs) { return FormalPowerSeries(lhs) /= rhs; } FormalPowerSeries rev() const { FormalPowerSeries res(*this); std::reverse(all(res)); return res; } friend FormalPowerSeries div(FormalPowerSeries lhs, FormalPowerSeries rhs) { lhs.shrink(); rhs.shrink(); if (lhs.size() < rhs.size()) { return FormalPowerSeries{}; } int n = lhs.size() - rhs.size() + 1; if (rhs.size() <= 32) { FormalPowerSeries res(n); T iv = rhs.back().inv(); rrep (i, n) { T d = lhs[i + rhs.size() - 1] * iv; res[i] = d; rep (j, rhs.size()) lhs[i + j] -= d * rhs[j]; } return res; } return (lhs.rev().prefix(n) * rhs.rev().inv(n)).prefix(n).rev(); } friend FormalPowerSeries operator%(FormalPowerSeries lhs, FormalPowerSeries rhs) { lhs.shrink(); rhs.shrink(); if (lhs.size() < rhs.size()) { return lhs; } int n = lhs.size() - rhs.size() + 1; if (rhs.size() <= 32) { T iv = rhs.back().inv(); rrep (i, n) { T d = lhs[i + rhs.size() - 1] * iv; rep (j, rhs.size()) lhs[i + j] -= d * rhs[j]; } return lhs.shrink(); } return (lhs - div(lhs, rhs) * rhs).shrink(); } friend std::pair divmod(FormalPowerSeries lhs, FormalPowerSeries rhs) { lhs.shrink(); rhs.shrink(); if (lhs.size() < rhs.size()) { return {FormalPowerSeries{}, lhs}; } int n = lhs.size() - rhs.size() + 1; if (rhs.size() <= 32) { FormalPowerSeries res(n); T iv = rhs.back().inv(); rrep (i, n) { T d = lhs[i + rhs.size() - 1] * iv; res[i] = d; rep (j, rhs.size()) lhs[i + j] -= d * rhs[j]; } return {res, lhs.shrink()}; } FormalPowerSeries q = div(lhs, rhs); return {q, (lhs - q * rhs).shrink()}; } FormalPowerSeries& operator%=(const FormalPowerSeries& rhs) { return *this = *this % rhs; } FormalPowerSeries diff() const { if (this->empty()) return {}; FormalPowerSeries res(this->size() - 1); rep (i, res.size()) res[i] = (*this)[i + 1] * (i + 1); return res; } FormalPowerSeries integral() const { FormalPowerSeries res(this->size() + 1); res[0] = 0; Comb::init(this->size()); rep (i, this->size()) res[i + 1] = (*this)[i] * Comb::inv(i + 1); return res; } template::value>::type* = nullptr> FormalPowerSeries inv(int deg = -1) const { assert(this->size() > 0 && (*this)[0] != 0); if (deg == -1) deg = this->size(); FormalPowerSeries res(1, (*this)[0].inv()); for (int m = 1; m < deg; m <<= 1) { FormalPowerSeries f(2 * m); for (int i = 0; i < std::min(2 * m, (int)this->size()); i++) f[i] = (*this)[i]; res.resize(2 * m); FormalPowerSeries dft = res; number_theoretic_transform(f); number_theoretic_transform(dft); rep (i, 2 * m) f[i] *= dft[i]; inverse_number_theoretic_transform(f); std::fill(f.begin(), f.begin() + m, T{0}); number_theoretic_transform(f); rep (i, 2 * m) dft[i] *= f[i]; inverse_number_theoretic_transform(dft); rep (i, m, 2 * m) res[i] = -dft[i]; } return res.prefix(deg); } template::value>::type* = nullptr> FormalPowerSeries inv(int deg = -1) const { assert(this->size() > 0 && (*this)[0] != 0); if (deg == -1) deg = this->size(); FormalPowerSeries res(1, (*this)[0].inv()); for (int m = 1; m < deg; m <<= 1) { res = (res * 2 - (res * res * this->prefix(2 * m)).prefix(2 * m)) .prefix(2 * m); } return res.prefix(deg); } FormalPowerSeries log(int deg = -1) const { assert(this->size() > 0 && (*this)[0] == 1); if (deg == -1) deg = this->size(); return (diff() * inv(deg)).prefix(deg - 1).integral(); } template::value>::type* = nullptr> FormalPowerSeries exp(int deg = -1) const { assert(this->size() > 0 && (*this)[0] == 0); if (deg == -1) deg = this->size(); FormalPowerSeries f(1, 1); FormalPowerSeries g(1, 1); FormalPowerSeries dft_f(1, 1); for (int m = 1; m < deg; m <<= 1) { FormalPowerSeries q = prefix(m).diff(); q.resize(m); number_theoretic_transform(q); rep (i, m) q[i] *= dft_f[i]; inverse_number_theoretic_transform(q); FormalPowerSeries s = f.diff(); s.resize(m); rep (i, m) s[i] -= q[i]; s.insert(s.begin(), (T)s.back()); s.pop_back(); FormalPowerSeries dft_g = g; s.resize(2 * m); dft_g.resize(2 * m); number_theoretic_transform(s); number_theoretic_transform(dft_g); rep (i, 2 * m) s[i] *= dft_g[i]; inverse_number_theoretic_transform(s); FormalPowerSeries u = (prefix(2 * m) - (s.prefix(m) << (m - 1)).integral()) >> m; u.resize(2 * m); FormalPowerSeries dft_f_2 = f; dft_f_2.resize(2 * m); number_theoretic_transform(u); number_theoretic_transform(dft_f_2); rep (i, 2 * m) u[i] *= dft_f_2[i]; inverse_number_theoretic_transform(u); f = f + (u.prefix(m) << m); if (2 * m < deg) { g.resize(2 * m); FormalPowerSeries dft_g_2 = g; FormalPowerSeries dft_f_2 = f; number_theoretic_transform(dft_g_2); number_theoretic_transform(dft_f_2); dft_f = dft_f_2; rep (i, 2 * m) dft_f_2[i] *= dft_g_2[i]; inverse_number_theoretic_transform(dft_f_2); std::fill(dft_f_2.begin(), dft_f_2.begin() + m, T{0}); number_theoretic_transform(dft_f_2); rep (i, 2 * m) dft_f_2[i] *= dft_g_2[i]; inverse_number_theoretic_transform(dft_f_2); rep (i, m, 2 * m) g[i] = -dft_f_2[i]; } } return f.prefix(deg); } template::value>::type* = nullptr> FormalPowerSeries exp(int deg = -1) const { assert(this->size() > 0 && (*this)[0] == 0); if (deg == -1) deg = this->size(); FormalPowerSeries res(1, 1); for (int m = 1; m < deg; m <<= 1) { res = (res * (prefix(2 * m) - res.log(2 * m)) + res).prefix(2 * m); } return res.prefix(deg); } FormalPowerSeries pow(ll k, int deg = -1) const { if (deg == -1) deg = this->size(); if (deg == 0) return {}; if (k == 0) { FormalPowerSeries res(deg); res[0] = 1; return res; } if (k == 1) return prefix(deg); if (k == 2) return (*this * *this).prefix(deg); T a; int d = -1; rep (i, this->size()) { if ((*this)[i] != 0) { a = (*this)[i]; d = i; break; } } if (d == -1) { FormalPowerSeries res(deg); return res; } if ((i128)(d)*k >= deg) { FormalPowerSeries res(deg); return res; } deg -= d * k; FormalPowerSeries res = (((*this >> d) / a).log(deg) * k).exp(deg); res *= a.pow(k); res <<= d * k; return res; } FormalPowerSeries sqrt(int deg = -1) const { if (deg == -1) deg = this->size(); T a; int d = -1; rep (i, this->size()) { if ((*this)[i] != 0) { a = (*this)[i]; d = i; break; } } if (d == -1) { FormalPowerSeries res(deg); return res; } if (d & 1) return {}; deg -= (d >> 1); if (deg <= 0) { FormalPowerSeries res(deg); return res; } FormalPowerSeries f = (*this >> d); T sq = sqrt_mod(a.get()); if (sq == -1) return {}; FormalPowerSeries g(1, sq); for (int m = 1; m < deg; m <<= 1) { g = (g + (f.prefix(2 * m) * g.inv(2 * m)).prefix(2 * m)) / 2; } g.resize(deg); return g << (d >> 1); } FormalPowerSeries compose(FormalPowerSeries g, int deg = -1) const { if (this->empty()) return {}; if (g.empty()) return {(*this)[0]}; assert(g[0] == 0); int n = deg == -1 ? this->size() : deg; int m = 1 << (bitop::ceil_log2(std::max(1, std::sqrt(n / std::log2(n)))) + 1); FormalPowerSeries p = g.prefix(m), q = g >> m; p.shrink(); q.shrink(); int l = (n + m - 1) / m; std::vector fs(this->size()); rep (i, this->size()) fs[i] = FormalPowerSeries{(*this)[i]}; FormalPowerSeries pd = p.diff(); int z = 0; while (z < (int)pd.size() && pd[z] == T{0}) z++; if (z == (int)pd.size()) { FormalPowerSeries ans; rrep (i, l) { ans = ((ans * q) << m).prefix(n - i * m) + FormalPowerSeries{(*this)[i]}; } return ans; } pd = (pd >> z).inv(n); FormalPowerSeries t = p; for (int k = 1; fs.size() > 1; k <<= 1) { std::vector nfs((fs.size() + 1) / 2); t.resize(1 << (bitop::ceil_log2(t.size()) + 1)); number_theoretic_transform(t); rep (i, fs.size() / 2) { nfs[i] = std::move(fs[2 * i]); fs[2 * i + 1].resize(t.size()); number_theoretic_transform(fs[2 * i + 1]); rep (j, t.size()) fs[2 * i + 1][j] *= t[j]; inverse_number_theoretic_transform(fs[2 * i + 1]); if ((int)fs[2 * i + 1].size() > n) fs[2 * i + 1].resize(n); nfs[i] += fs[2 * i + 1]; } if (fs.size() & 1) nfs.back() = std::move(fs.back()); fs = std::move(nfs); if (fs.size() > 1) { rep (i, t.size()) t[i] *= t[i]; inverse_number_theoretic_transform(t); if ((int)t.size() > n) t.resize(n); } } FormalPowerSeries fp = fs[0].prefix(n); FormalPowerSeries res = fp; int n2 = 1 << (bitop::ceil_log2(n) + 1); FormalPowerSeries qpow(n2); qpow[0] = 1; q.resize(n2); number_theoretic_transform(q); pd.resize(n2); number_theoretic_transform(pd); rep (i, 1, l) { if ((n - i * m) * 4 <= n2) { while ((n - i * m) * 4 <= n2) { n2 /= 2; } inverse_number_theoretic_transform(q); q.resize(n - i * m); q.resize(n2); number_theoretic_transform(q); inverse_number_theoretic_transform(pd); pd.resize(n - i * m); pd.resize(n2); number_theoretic_transform(pd); } qpow.resize(n - i * m); qpow.resize(n2); number_theoretic_transform(qpow); rep (j, n2) qpow[j] *= q[j]; inverse_number_theoretic_transform(qpow); qpow.resize(n - i * m); fp = fp.diff() >> z; fp.resize(n - i * m); fp.resize(n2); number_theoretic_transform(fp); rep (j, n2) fp[j] *= pd[j]; inverse_number_theoretic_transform(fp); fp.resize(n - i * m); res += ((qpow * fp).prefix(n - i * m) * Comb::finv(i)) << (i * m); } return res; } FormalPowerSeries compinv(int deg = -1) const { assert(this->size() >= 2 && (*this)[0] == 0 && (*this)[1] != 0); if (deg == -1) deg = this->size(); FormalPowerSeries fd = diff(); FormalPowerSeries x{0, 1}; FormalPowerSeries res{0, (*this)[1].inv()}; for (int m = 2; m < deg; m <<= 1) { auto tmp = prefix(2 * m).compose(res); auto d = tmp.diff(); auto gd = res.diff(); res -= ((tmp - x) * (d.inv(2 * m) * gd).prefix(2 * m)).prefix(2 * m); } return res.prefix(deg); } template::value>::type* = nullptr> FormalPowerSeries& ntt_doubling() { ntt_doubling_(*this); return *this; } }; /** * @brief FormalPowerSeries(形式的冪級数) * @docs docs/math/poly/FormalPowerSeries.md * @see https://nyaannyaan.github.io/library/fps/formal-power-series.hpp */ #line 2 "library/math/poly/TaylorShift.hpp" #line 7 "library/math/poly/TaylorShift.hpp" template> FormalPowerSeries taylor_shift(FormalPowerSeries f, T a) { const int n = f.size(); Comb::init(n); rep (i, n) f[i] *= Comb::fact(i); FormalPowerSeries g(n); T p = 1; rep (i, n) { g[n - 1 - i] = p * Comb::finv(i); p *= a; } f *= g; f >>= n - 1; rep (i, n) f[i] *= Comb::finv(i); return f; } /** * @brief TaylorShift * @docs docs/math/poly/TaylorShift.md */ #line 5 "main.cpp" using namespace std; using mint = modint998244353; using comb = Combinatorics; using fps = FormalPowerSeries; int main() { int N; cin >> N; vector A(N); rep (i, N) cin >> A[i]; map> mp; for (auto x : A) ++mp[x]; int M = A[0] + 1; mint sm = 0; fps B(M); fps f(M); vector> C; int i = 0; for (auto [a, b] : mp) { if ((++i) % 1500 == 0) { f += taylor_shift(B, mint{1}); B.assign(a + 1, 0); C.clear(); } mint t = f[a] + 1; for (auto [c, d] : C) t += comb::comb(c, a) * d; t *= mint{2}.pow(b) - 1; sm += t; B[A[i]] += t; C.emplace_back(a, t); } cout << sm.get() << endl; }