//Timestamp: 2024-08-11 17:27:26 #define DROP #ifdef ONLINE #undef LOCAL #endif #ifndef LOCAL #undef _GLIBCXX_DEBUG #undef _DEBUG #endif #include #include #include #include #include //#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include //#include "compiler_hint.cpp" template struct MDVecDef { using Type = std::vector::Type>; template static Type Make(int n, Args... args) { return Type(n, MDVecDef::Make(args...)); } }; template struct MDVecDef { using Type = T; static Type Make(T val = T()) { return val; } }; template using MDVec = typename MDVecDef::Type; #ifndef M_PI #define M_PI 3.14159265358979323851280895940618620443274267017841L #endif #ifndef M_E #define M_E 2.718281828459045235428168107993940338928950950503355L #endif #ifdef LOCAL #define Assert(x) assert(x) #define DebugRun(X) X #define DebugPoint int _x_ = 0; _x_++; #else #define Debug(...) 42 #define DebugFmtln(...) 42 #define Assert(x) 42 #define DebugRun(X) #define DebugPoint #endif #define Trace(x) DebugFmtln("Line %d: %s", __LINE__, #x) template inline T DebugRet(T x) { Debug(x); return x; } #define const_ref(T) const T & #define mut_ref(T) T & #define let auto #define var auto #define MEMSET0(X) std::memset(&X, 0, sizeof(X)) #define Size(T) int((T).size()) #define All(data) data.begin(), data.end() #define MakeUnique(data) data.resize(std::unique(All(data)) - data.begin()) #define MakeUniqueAndSort(data) Sort(All(data)); MakeUnique(data) #define MakeAttribute(struct_name, Type, attr_name) \ struct struct_name { \ using attr_name ## _type = Type; \ Type attr_name; \ mut_ref(Type) get_##attr_name() { return attr_name; } \ const_ref(Type) get_##attr_name() const { return attr_name; } \ }; #define MakeTemplateAttribute(struct_name, attr_name) \ template \ struct struct_name { \ using attr_name##_type = T; \ T attr_name; \ mut_ref(T) get_##attr_name() { return attr_name; } \ const_ref(T) get_##attr_name() const { return attr_name; } \ }; #define ImplDefaultEq(name) \ bool operator==(const name &a, const name &b) { \ return std::memcmp(&a, &b, sizeof(name)) == 0; \ } \ bool operator!=(const name &a, const name &b) { return !(a == b); } #define ImplDefaultComparision(name) \ bool operator>(const name &rhs) const { return rhs < *this; } \ bool operator<=(const name &rhs) const { return !(*this > rhs); } \ bool operator>=(const name &rhs) const { return !(*this < rhs); } #define ImplArithmeticAssignOperation(name) \ name &operator+=(const name &rhs) { return *this = (*this) + rhs; } \ name &operator-=(const name &rhs) { return *this = (*this) - rhs; } \ name &operator*=(const name &rhs) { return *this = (*this) * rhs; } \ name &operator/=(const name &rhs) { return *this = (*this) / rhs; } #define IsType(Type, param, ret_type) \ template \ enable_if_t && is_same_v, \ ret_type> #define IsBool(param, ret_type) \ template \ enable_if_t #define IsBoolStatic(param, ret_type) \ template \ static enable_if_t #define MakeAnnotation(name) \ template \ struct is_##name { \ static const bool value = false; \ }; \ template \ inline constexpr bool is_##name##_v = is_##name::value; #define AssignAnnotation(cls, annotation) \ template <> \ struct is_##annotation { \ static const bool value = true; \ }; #define AssignAnnotationTemplate(cls, annotation, type) \ template \ struct is_##annotation> { \ static const bool value = true; \ }; #define FunctionAlias(from, to) \ template \ inline auto to(Args &&...args) \ ->decltype(from(std::forward(args)...)) { \ return from(std::forward(args)...); \ } #define CastToScalar(field, type) \ operator type() const { return type(field); } #define CastToAllScalar(field) \ CastToScalar(field, i8); \ CastToScalar(field, u8); \ CastToScalar(field, i16); \ CastToScalar(field, u16); \ CastToScalar(field, i32); \ CastToScalar(field, u32); \ CastToScalar(field, i64); \ CastToScalar(field, u64); \ CastToScalar(field, f32); \ CastToScalar(field, f64); \ CastToScalar(field, f80); #define COMMA , #ifndef LOCAL std::mt19937 rng(std::chrono::steady_clock::now().time_since_epoch().count()); #else std::mt19937 rng(0); #endif template T random_choice(T l, T r, std::mt19937 &gen = rng) { std::uniform_int_distribution random(l, r); return random(gen); } namespace dalt { #ifndef LOCAL struct Timer {explicit Timer(const char* m) {}void stop() const {}}; #else #endif } using i8 = char; using i16 = short; using i32 = int; using i64 = long long; using u8 = unsigned char; using u16 = unsigned short; using u32 = unsigned int; using u64 = unsigned long long; using usize = size_t; using f32 = float; using f64 = double; // 16 exp, 64 precision using f80 = long double; FunctionAlias(std::lower_bound, LowerBound); FunctionAlias(std::upper_bound, UpperBound); FunctionAlias(std::unique, Unique); FunctionAlias(std::swap, Swap); FunctionAlias(std::min, Min); FunctionAlias(std::max, Max); FunctionAlias(std::abs, Abs); FunctionAlias(std::sin, Sin); FunctionAlias(std::asin, Asin); FunctionAlias(std::cos, Cos); FunctionAlias(std::acos, Acos); FunctionAlias(std::tan, Tan); FunctionAlias(std::atan, Atan); FunctionAlias(std::sort, Sort); FunctionAlias(std::fill, Fill); FunctionAlias(std::move, Move); FunctionAlias(std::reverse, Reverse); FunctionAlias(std::max_element, MaxElement); FunctionAlias(std::min_element, MinElement); FunctionAlias(std::make_tuple, MakeTuple); FunctionAlias(std::make_pair, MakePair); FunctionAlias(std::clamp, Clamp); FunctionAlias(std::shuffle, Shuffle); FunctionAlias(std::to_string, ToString); FunctionAlias(std::tie, Tie); FunctionAlias(std::get<0>, Get0); FunctionAlias(std::get<1>, Get1); FunctionAlias(std::get<2>, Get2); FunctionAlias(std::get<3>, Get3); FunctionAlias(std::get<4>, Get4); template using Function = std::function<_Signature>; template using Func = Function<_Signature>; using Str = std::string; using String = Str; using StringStream = std::stringstream; using IStream = std::istream; using OStream = std::ostream; using std::enable_if; using std::enable_if_t; using std::is_base_of; using std::is_base_of_v; using std::is_floating_point; using std::is_floating_point_v; using std::is_integral; using std::is_integral_v; using std::is_arithmetic; using std::is_arithmetic_v; using std::is_same; using std::is_same_v; using std::tie; auto &Stderr = std::cerr; auto &Stdin = std::cin; auto &Stdout = std::cout; template using Less = std::less; template using Greater = std::greater; template > using TreeMap = std::map<_Key, _Tp, _Compare>; template > using TreeSet = std::set<_Key, _Compare>; template , typename _Alloc = std::allocator<_Key>> using MultiTreeSet = std::multiset<_Key, _Compare, _Alloc>; template using Deque = std::deque; template using Queue = std::queue; template using Vec = std::vector; template using Reducer = Func; template using Comparator = Func; template using Indexer = Func; template using Indexer2 = Func; template using Adder = Func; template using Checker = Func; template using BiChecker = Func; template using Consumer = Func; template using Supplier = Func; template using BiConsumer = Func; template using Mapper = Func; template using MinHeap = std::priority_queue, Greater>; template using MaxHeap = std::priority_queue, Less>; template using Array = std::array; template using Tuple = std::tuple<_Elements...>; template >> using Complex = std::complex; template using Pair = std::pair; namespace dalt { template IStream& operator>>(IStream& is, Vec& val) { for (auto& v : val) { is >> v; } return is; } #define VEC_OP(op) \ template \ Vec& operator op(Vec& data, T x) { \ for (auto& v : data) { \ v op x; \ } \ return data; \ } VEC_OP(+=) VEC_OP(-=) VEC_OP(*=) VEC_OP(/=) VEC_OP(%=) VEC_OP(^=) VEC_OP(&=) VEC_OP(|=) VEC_OP(==) VEC_OP(!=) template int Compare(const Vec& lhs, const Vec& rhs) { for(int i = 0; i < Size(lhs) && i < Size(rhs); i++) { if(lhs[i] != rhs[i]) { return lhs[i] < rhs[i] ? -1 : 1; } } return Size(lhs) < Size(rhs) ? -1 : Size(lhs) > Size(rhs) ? 1 : 0; } template bool operator<(const Vec& lhs, const Vec& rhs) { return Compare(lhs, rhs) < 0; } template bool operator>(const Vec& lhs, const Vec& rhs) { return Compare(lhs, rhs) > 0; } template bool operator<=(const Vec& lhs, const Vec& rhs) { return Compare(lhs, rhs) <= 0; } template bool operator>=(const Vec& lhs, const Vec& rhs) { return Compare(lhs, rhs) >= 0; } } // namespace dalt //#include "array_adder.cpp" using namespace dalt; namespace dalt { template struct Optional { using Self = Optional; private: T val; bool show_up; public: Optional(const T &arg_val) : val(arg_val), show_up(true) {} Optional(const T &&arg_val) : val(arg_val), show_up(true) {} Optional() : show_up(false) {} const T &value() const { Assert(show_up); return val; } T &value() { Assert(show_up); return val; } T &operator*() { return value(); } const T &operator*() const { return value(); } bool is_some() const { return show_up; } bool is_none() const { return !show_up; } const T *operator->() const { return &value(); } T *operator->() { return &value(); } inline operator T() const { return value(); } T or_else(T def) const { if (is_some()) { return val; } else { return def; } } template Optional map(const Mapper &mapper) const { if (is_some()) { return mapper(value()); } else { return Optional(); } } bool operator==(const Self &b) const { return show_up == b.show_up && (!show_up || val == b.val); } }; template bool operator!=(const Optional &a, const Optional &b) { return !(a == b); } template OStream &operator<<(OStream &os, const Optional &v) { if (v.is_none()) { os << "{}"; } else { os << '{' << v.value() << '}'; } return os; } } // namespace dalt namespace dalt { template inline enable_if_t, T> Gcd(T a, T b) { while (b != 0) { a %= b; Swap(a, b); } return a; } // ret_value = [x, y, gcd(a,b)] that x * a + y * b = gcd(a, b) template inline enable_if_t, Array> ExtGcd(T a, T b) { if (b == 0) { return Array{1, 0, a}; } auto div = a / b; auto ans = ExtGcd(b, a - b * div); auto x = ans[0]; auto y = ans[1]; return Array{y, x - a / b * y, ans[2]}; } template inline enable_if_t, Optional> PossibleModInverse( T a, T modulus) { auto res = ExtGcd(a, modulus); if (res[2] == 1) { auto ans = res[0] % modulus; if (ans < 0) { ans += modulus; } return ans; } return {}; } } // namespace dalt namespace dalt { template T Modular(E val, T mod) { val = val % mod; if (val < 0) { val = val + mod; } return T(val); } template inline T MulMod(T a, T b, T modulus) { i64 res = i64(a) * i64(b) % modulus; return T(res); } template<> inline i64 MulMod(i64 a, i64 b, i64 modulus) { #ifdef __SIZEOF_INT128__ // do some fancy stuff here return __int128_t(a) * __int128_t(b) % modulus; #else // do some fallback stuff here i64 k = roundl((f80)a / modulus * b); i64 res = (a * b - k * modulus) % modulus; if (res < 0) { res += modulus; } return res; #endif } template inline enable_if_t, T> AddMod(T a, T b, T modulus) { T res = a + b; if (res >= modulus) { res -= modulus; } return res; } template inline enable_if_t && is_integral_v, T> PowMod(T x, E exp, T modulus) { Assert(exp >= E(0)); if (exp == E(0)) { return modulus == T(1) ? T(0) : T(1); } T ans = PowMod(x, exp >> 1, modulus); ans = MulMod(ans, ans, modulus); if (exp & T(1)) { ans = MulMod(ans, x, modulus); } return ans; } template inline enable_if_t, T> SubMod(T a, T b, T modulus) { T res = a - b; if (res < T(0)) { res += modulus; } return res; } } // namespace dalt namespace dalt { template inline A& Chmin(A& a, const B& b) { if (a > b) { a = b; } return a; } template inline T& Chmin(T& a, const T& b, const Comparator &comp) { if (comp(b, a)) { a = b; } return a; } template inline A& Chmax(A& a, const B& b) { if (a < b) { a = b; } return a; } template inline T& Chmax(T& a, const T& b, const Comparator& comp) { if (comp(a, b)) { a = b; } return a; } template inline T& AddTo(T& a, const T& b) { a = a + b; return a; } template inline T& MulTo(T& a, const T& b) { a = a * b; return a; } template inline T& SubFrom(T& a, const T& b) { a = a - b; return a; } template inline T& DivFrom(T& a, const T& b) { a = a / b; return a; } template constexpr enable_if_t, T> PowBinaryLift(T x, E n) { if (n == E(0)) { return T(1); } auto ans = PowBinaryLift(x, n >> 1); ans = ans * ans; if (n & 1) { ans = ans * x; } return ans; } template inline T MulLimit(T a, T b, T max, T def) { if (a == T(0) || b == T(0)) { return T(0); } // a * b <= max // a <= max / b // a <= floor(max / b) if (a <= max / b) { return a * b; } else { return def; } } template inline T MulLimit(T a, T b, T max) { return MulLimit(a, b, max, max); } template inline T AddLimit(T a, T b, T max, T def) { if (a <= max - b) { return a + b; } else { return def; } } template inline T AddLimit(T a, T b, T max) { return AddLimit(a, b, max, max); } i64 Round(f32 x) { return roundf(x); } i64 Round(f64 x) { return round(x); } i64 Round(f80 x) { return roundl(x); } //l + ... + r template T SumOfInterval(T l, T r) { if(l > r) { return T(0); } return (l + r) * (r - l + 1) / T(2); } template T Pow2(T x) { return x * x; } } // namespace dalt namespace dalt { // using Type = T; // static T modulus; // static T primitive_root; MakeAnnotation(modular); template struct StaticModular { static_assert(is_integral_v); const static T modulus = M; const static T primitive_root = PR; const static T phi = PHI; using Type = T; }; template struct is_modular> { static const bool value = true; }; using MOD998244353 = StaticModular; using MOD1004535809 = StaticModular; using MOD1000000007 = StaticModular; using MOD1000000009 = StaticModular; using MOD_BIG = StaticModular; // last used: -2 template struct DynamicModular { static_assert(is_integral_v); static T modulus; static T primitive_root; static T phi; static void Register(T _modulus, T _primitive_root = T(), T _phi = T()) { modulus = _modulus; primitive_root = _primitive_root; phi = _phi; } using Type = T; }; template T DynamicModular::modulus = T(); template T DynamicModular::primitive_root = T(); template T DynamicModular::phi = T(); template struct is_modular> { static const bool value = true; }; MakeAnnotation(modint); MakeAnnotation(modint_32); #define MOD MODULAR::modulus #define SELF ModInt #define TEMPLATE_ARGS template TEMPLATE_ARGS struct ModInt { using Modular = MODULAR; using Type = typename MODULAR::Type; static_assert(is_modular_v); static_assert(is_integral_v); Type value; using Self = SELF; ModInt() : value(0) {} ModInt(const Type &v) { value = v; if (value < 0 || value >= MOD) { value %= MOD; if (value < 0) { value += MOD; } } Assert(value >= 0); } static Self nil() { Self res; res.value = -1; return res; } bool is_nil() { return value == -1; } explicit operator Type() const { return value; } static Type modulus() { return MOD; } static Type primitive_root() { return Modular::primitive_root; } Self operator-() const { return Self(0) - *this; } template static enable_if_t, Self> of(F v) { v %= MOD; if (v < 0) { v += MOD; } return Self(v); } Optional possible_inv() const { auto raw_optional_inv = PossibleModInverse(value, MOD); if (raw_optional_inv.is_some()) { return Self(raw_optional_inv.value()); } else { return {}; } } Self operator+(const Self &rhs) const { auto res = value + rhs.value; if (res >= MOD || res < 0) { res -= MOD; } return res; } Self operator-(const Self &rhs) const { if(value >= rhs.value) { return value - rhs.value; } else { return value + MOD - rhs.value; } } Self operator/(const SELF &rhs) const { auto inv = Self(rhs.possible_inv().value()); return (*this) * inv; } bool operator==(const Self &b) const { return value == b.value; } bool operator!=(const Self &b) const { return value != b.value; } bool operator<(const Self &b) const { return value < b.value; } ImplDefaultComparision(Self); ImplArithmeticAssignOperation(Self); template enable_if_t, Self> pow(E n) { return PowBinaryLift(*this, n); } friend inline IStream &operator>>(IStream &is, Self &x) { Type val; is >> val; x = Self::of(val); return is; } friend inline OStream &operator<<(OStream &os, const Self &x) { os << x.value; return os; } }; TEMPLATE_ARGS inline enable_if_t, SELF> operator*( const SELF &a, const SELF &b) { return SELF(MulMod(a.value, b.value, MOD)); } TEMPLATE_ARGS inline enable_if_t, SELF> operator*( const SELF &x, const SELF &y) { static u64 mask = (u64(1) << 32) - 1; static u64 mod = MOD_BIG::modulus; u64 a = x.value; u64 b = y.value; u64 l1 = a & mask; u64 h1 = (a >> 32) & mask; u64 l2 = b & mask; u64 h2 = (b >> 32) & mask; u64 l = l1 * l2; u64 m = l1 * h2 + l2 * h1; u64 h = h1 * h2; u64 ret = (l & mod) + (l >> 61) + (h << 3) + (m >> 29) + ((m << 35) >> 3) + 1; ret = (ret & mod) + (ret >> 61); ret = (ret & mod) + (ret >> 61); return SELF(ret - 1); } TEMPLATE_ARGS struct is_modint> { static const bool value = true; }; TEMPLATE_ARGS struct is_modint_32> { static const bool value = is_same_v; }; #undef TEMPLATE_TYPE #undef MOD using ModInt998244353 = ModInt; using ModInt1000000007 = ModInt; using ModInt1000000009 = ModInt; using ModInt1004535809 = ModInt; } // namespace dalt #ifndef _builtin_clz inline i32 _builtin_clz(u32 i) { // HD, Count leading 0's if (i <= 0) return i == 0 ? 32 : 0; int n = 31; if (i >= 1 << 16) { n -= 16; i >>= 16; } if (i >= 1 << 8) { n -= 8; i >>= 8; } if (i >= 1 << 4) { n -= 4; i >>= 4; } if (i >= 1 << 2) { n -= 2; i >>= 2; } return n - (i >> 1); } #endif #ifndef _builtin_clzll inline i32 _builtin_clzll(u64 i) { u32 x = u32(i >> 32); return x == 0 ? 32 + _builtin_clz((int)i) : _builtin_clz(x); } #endif #ifndef _builtin_ctz inline i32 _builtin_ctz(u32 i) { // HD, Figure 5-14 int y; if (i == 0) return 32; int n = 31; y = i << 16; if (y != 0) { n = n - 16; i = y; } y = i << 8; if (y != 0) { n = n - 8; i = y; } y = i << 4; if (y != 0) { n = n - 4; i = y; } y = i << 2; if (y != 0) { n = n - 2; i = y; } return n - ((i << 1) >> 31); } #endif #ifndef _builtin_ctzll inline i32 _builtin_ctzll(u64 i) { // HD, Figure 5-14 int x, y; if (i == 0) return 64; int n = 63; y = (int)i; if (y != 0) { n = n - 32; x = y; } else x = (int)(i >> 32); y = x << 16; if (y != 0) { n = n - 16; x = y; } y = x << 8; if (y != 0) { n = n - 8; x = y; } y = x << 4; if (y != 0) { n = n - 4; x = y; } y = x << 2; if (y != 0) { n = n - 2; x = y; } return n - ((x << 1) >> 31); } #endif #ifndef _builtin_popcount inline i32 _builtin_popcount(u32 i) { // HD, Figure 5-2 i = i - ((i >> 1) & 0x55555555); i = (i & 0x33333333) + ((i >> 2) & 0x33333333); i = (i + (i >> 4)) & 0x0f0f0f0f; i = i + (i >> 8); i = i + (i >> 16); return i & 0x3f; } #endif #ifndef _builtin_popcountll inline i32 _builtin_popcountll(u64 i) { // HD, Figure 5-2 i = i - ((i >> 1) & 0x5555555555555555ll); i = (i & 0x3333333333333333ll) + ((i >> 2) & 0x3333333333333333ll); i = (i + (i >> 4)) & 0x0f0f0f0f0f0f0f0fll; i = i + (i >> 8); i = i + (i >> 16); i = i + (i >> 32); return (int)i & 0x7f; } #endif namespace dalt { inline i32 LeadingZeroNumber(u32 x) { if (x == 0) return 32; return _builtin_clz(x); } inline i32 LeadingZeroNumber(i32 x) { return LeadingZeroNumber(u32(x)); } inline i32 LeadingZeroNumber(u64 x) { if (x == 0) return 64; return _builtin_clzll(x); } inline i32 LeadingZeroNumber(i64 x) { return LeadingZeroNumber(u64(x)); } inline i32 TrailingZeroNumber(u32 x) { if (x == 0) return 32; return _builtin_ctz(x); } inline i32 TrailingZeroNumber(i32 x) { return TrailingZeroNumber(u32(x)); } inline i32 TrailingZeroNumber(u64 x) { if (x == 0) return 64; return _builtin_ctzll(x); } inline i32 TrailingZeroNumber(i64 x) { return TrailingZeroNumber(u64(x)); } inline i32 Log2Ceil(u32 x) { if (x == 0) { return 0; } return 32 - LeadingZeroNumber(x - 1); } inline i32 Log2Ceil(u64 x) { if (x == 0) { return 0; } return 64 - LeadingZeroNumber(x - 1); } inline i32 Log2Floor(u32 x) { if (x == 0) { return -1; } return 31 - LeadingZeroNumber(x); } inline i32 Log2Floor(u64 x) { if (x == 0) { return -1; } return 63 - LeadingZeroNumber(x); } inline i32 Log2Ceil(i32 x) { return Log2Ceil(u32(x)); } inline i32 Log2Ceil(i64 x) { return Log2Ceil(u64(x)); } inline i32 Log2Floor(i32 x) { return Log2Floor(u32(x)); } inline i32 Log2Floor(i64 x) { return Log2Floor(u64(x)); } inline i32 CountBit(u32 x) { return _builtin_popcount(x); } inline i32 CountBit(i32 x) { return CountBit(u32(x)); } inline i32 CountBit(u64 x) { return _builtin_popcountll(x); } inline i32 CountBit(i64 x) { return CountBit(u64(x)); } inline i32 HighestOneBitOffset(u32 x) { return Log2Floor(x); } inline i32 HighestOneBitOffset(i32 x) { return HighestOneBitOffset(u32(x)); } inline i32 HighestOneBitOffset(u64 x) { return Log2Floor(x); } inline i32 HighestOneBitOffset(i64 x) { return HighestOneBitOffset(u64(x)); } template inline enable_if_t, T> LowestOneBit(T x) { return x & -x; } template inline enable_if_t, T> HighestOneBit(T x) { if (x == 0) { return x; } return T(1) << HighestOneBitOffset(x); } template inline enable_if_t, i32> LowestOneBitOffset(T x) { if (x == 0) { return -1; } return HighestOneBitOffset(LowestOneBit(x)); } inline u32 HighestKOnes32(i32 k) { if (k == 0) { return 0; } return (~u32()) << (32 - k); } inline u64 HighestKOnes64(i32 k) { if (k == 0) { return 0; } return (~u64()) << (64 - k); } inline u32 LowestKOnes32(i32 k) { if (k == 0) { return 0; } return (~u32()) >> (32 - k); } inline u64 LowestKOnes64(i32 k) { if (k == 0) { return 0; } return (~u64()) >> (64 - k); } inline u64 IntervalOnes64(i32 l, i32 r) { if (l > r) { return 0; } u64 high = r < 63 ? (u64(-1) << r + 1) : 0; u64 low = u64(-1) << l; return high ^ low; } inline u32 IntervalOnes32(i32 l, i32 r) { if (l > r) { return 0; } u32 high = r < 31 ? (u32(-1) << r + 1) : 0; u32 low = u32(-1) << l; return high ^ low; } template inline enable_if_t, i32> KthBit(T x, i32 k) { return (x >> k) & 1; } template inline enable_if_t, T> SetBit(T x, i32 k) { return x | (T(1) << k); } template inline enable_if_t, T> ClearBit(T x, i32 k) { return x & ~(T(1) << k); } } // namespace dalt namespace dalt { static const u8 BitReverseTable[]{ 0, 128, 64, 192, 32, 160, 96, 224, 16, 144, 80, 208, 48, 176, 112, 240, 8, 136, 72, 200, 40, 168, 104, 232, 24, 152, 88, 216, 56, 184, 120, 248, 4, 132, 68, 196, 36, 164, 100, 228, 20, 148, 84, 212, 52, 180, 116, 244, 12, 140, 76, 204, 44, 172, 108, 236, 28, 156, 92, 220, 60, 188, 124, 252, 2, 130, 66, 194, 34, 162, 98, 226, 18, 146, 82, 210, 50, 178, 114, 242, 10, 138, 74, 202, 42, 170, 106, 234, 26, 154, 90, 218, 58, 186, 122, 250, 6, 134, 70, 198, 38, 166, 102, 230, 22, 150, 86, 214, 54, 182, 118, 246, 14, 142, 78, 206, 46, 174, 110, 238, 30, 158, 94, 222, 62, 190, 126, 254, 1, 129, 65, 193, 33, 161, 97, 225, 17, 145, 81, 209, 49, 177, 113, 241, 9, 137, 73, 201, 41, 169, 105, 233, 25, 153, 89, 217, 57, 185, 121, 249, 5, 133, 69, 197, 37, 165, 101, 229, 21, 149, 85, 213, 53, 181, 117, 245, 13, 141, 77, 205, 45, 173, 109, 237, 29, 157, 93, 221, 61, 189, 125, 253, 3, 131, 67, 195, 35, 163, 99, 227, 19, 147, 83, 211, 51, 179, 115, 243, 11, 139, 75, 203, 43, 171, 107, 235, 27, 155, 91, 219, 59, 187, 123, 251, 7, 135, 71, 199, 39, 167, 103, 231, 23, 151, 87, 215, 55, 183, 119, 247, 15, 143, 79, 207, 47, 175, 111, 239, 31, 159, 95, 223, 63, 191, 127, 255}; inline u8 ReverseBit(u8 x) { return BitReverseTable[x]; } inline i8 ReverseBit(i8 x) { return ReverseBit((u8)x); } #define DefineReverseBit(bit, low) \ inline u##bit ReverseBit(u##bit x) { \ static const u##bit mask = (u##bit(1) << low) - 1; \ return ReverseBit(u##low(x >> low)) | \ (u##bit(ReverseBit(u##low(x & mask))) << low); \ } \ inline i##bit ReverseBit(i##bit x) { return ReverseBit(u##bit(x)); } DefineReverseBit(16, 8); DefineReverseBit(32, 16); DefineReverseBit(64, 32); #undef DefineReverseBit } // namespace dalt namespace dalt { namespace poly { MakeAnnotation(convolution) } } // namespace dalt namespace dalt { struct uint128 { u64 high, low; uint128(u64 _low = 0, u64 _high = 0) : high(_high), low(_low) {} using Self = uint128; CastToAllScalar(low); friend Self operator+(const Self& a, const Self& b) { Self ans(a.low + b.low, a.high + b.high); if (ans.low < a.low) { ans.high++; } return ans; } friend Self operator-(const Self& a, const Self& b) { u64 h = a.high - b.high; u64 l = a.low - b.low; if(a.low < b.low) { l = -l; h--; } return Self(l, h); } friend OStream& operator<<(OStream &os, const Self& x) { if(x.high) { os << x.high; } os << x.low; return os; } friend IStream& operator>>(IStream &is, Self& x) { x.high = 0; is >> x.low; return is; } friend Self operator*(const Self& a, const Self& b) { static u64 mask = ((u64(1) << 32) - 1); u64 all = a.low & mask; u64 bll = b.low & mask; u64 alh = a.low >> 32; u64 blh = b.low >> 32; Self ans(0, a.high * b.low + a.low * b.high + alh * blh + (all * blh >> 32) + (alh * bll >> 32)); ans += all * bll; ans += ((all * blh) & mask) << 32; ans += ((bll * alh) & mask) << 32; return ans; } Self& operator+=(const u64& rhs) { low += rhs; if (low < rhs) { high++; } return *this; } bool operator<(const Self& rhs) const { return MakePair(high, low) < MakePair(rhs.high, rhs.low); } ImplDefaultComparision(Self); ImplArithmeticAssignOperation(Self); //ImplArithmeticAssignOperation(Self); bool operator==(const Self& rhs) const { return high == rhs.high && low == rhs.low; } bool operator!=(const Self& rhs) const { return !(*this == rhs); } // x < 2 ^ 31 u32 modular(u32 x) const { static u64 max = std::numeric_limits::max(); u64 ans = low < x ? low : low % x; if (high > 0) { ans = ans + (high >= x ? high % x : high) * (max % x + 1); ans %= x; } return ans; } u32 operator%(u32 x) const { return modular(x); } Self operator/(u32 x) const { static const u32 ALL_ONE = -1; u64 y = high; u64 top = y / x; y = y - top * x; y = (y << 32) | (low >> 32); u64 low_high = y / x; y = y - low_high * x; y = (y << 32) | (low & ALL_ONE); u64 low_low = y / x; return Self(low_low | (low_high << 32), top); } }; using u128 = uint128; } // namespace dalt namespace dalt { namespace poly { template struct BruteForceConv { using Type = T; template static enable_if_t && !(is_modint_v && is_same_v), Vec> conv(const Vec &a, const Vec &b) { int rank = Size(a) + Size(b) - 2; Vec c(rank + 1); for (int i = 0; i < Size(a); i++) { for (int j = 0; j < Size(b); j++) { c[i + j] += a[i] * b[j]; } } return c; } template static enable_if_t && is_modint_v && is_same_v, Vec> conv(const Vec &a, const Vec &b) { int rank = Size(a) + Size(b) - 2; Vec data(rank + 1); for(int i = 0; i < Size(a); i++) { for(int j = 0; j < Size(b); j++) { data[i + j] += u64(a[i].value) * b[j].value; } } Vec ans(rank + 1); i32 modulus = T::modulus(); for(int i = 0; i <= rank; i++) { ans[i] = data[i].modular(modulus); } return ans; } static Vec conv2(const Vec &a) { return conv(a, a); } static Vec inverse(Vec p, i32 n) { Extend(p, n); auto dfs = [&](auto &dfs, i32 m) -> Vec { if (m == 1) { return Vec{T(1) / p[0]}; } i32 prev_mod = (m + 1) / 2; auto C = dfs(dfs, prev_mod); auto AC = conv(p, m, C, prev_mod); AC.resize(m); for (int i = 0; i < m; i++) { AC[i] = T(0) - AC[i]; } AC[0] = AC[0] + T(2); auto res = conv(C, AC); res.resize(m); return res; }; auto ans = dfs(dfs, n); ans.resize(n); return ans; } static Array, 2> div_and_rem(Vec a, Vec b) { Reverse(All(b)); auto inv_first = T(1) / b[0]; i32 n = Size(a); i32 m = Size(b); if (n < m) { return {Vec{T(0)}, Vec{Move(a)}}; } Vec divisor(n - m + 1); for (int i = n - 1; i >= m - 1; i++) { if (a[i] == T(0)) { continue; } T factor = a[i] * inv_first; divisor[i - (m - 1)] = factor; for (int j = 0; j < m; j++) { a[i - j] = a[i - j] - b[j] * factor; } } return Array, 2>{Move(divisor), Move(a)}; } }; template struct is_convolution> { static const bool value = true; }; } // namespace poly } // namespace dalt namespace dalt { namespace poly { template Vec CopyAndExtend(const Vec &data, int n) { Vec res; res.reserve(n); if (n <= Size(data)) { res.insert(res.begin(), data.begin(), data.begin() + n); } else { res.insert(res.begin(), All(data)); res.resize(n); } return res; } template void Normalize(Vec &p) { while (!p.empty() && p.back() == T(0)) { p.pop_back(); } if (p.empty()) { p.push_back(T(0)); } } template void Extend(Vec &p, int cap) { while (Size(p) < cap) { p.push_back(T()); } } template Vec DotMul(const Vec &a, const Vec &b) { int n = Min(Size(a), Size(b)); Vec ans(n); for (int i = 0; i < n; i++) { ans[i] = a[i] * b[i]; } return ans; } template void DotMulInplace(Vec &a, const Vec &b) { int n = Size(b); a.resize(n); for (int i = 0; i < n; i++) { a[i] = a[i] * b[i]; } } template void DotMulPlus(Vec &a, const Vec &b, Vec &res) { int n = Min(Size(a), Size(b)); for (int i = 0; i < n; i++) { res[i] = res[i] + a[i] * b[i]; } } template T Apply(const Vec &a, T x) { T sum = 0; for (int i = Size(a) - 1; i >= 0; i--) { sum = sum * x + a[i]; } return sum; } } // namespace poly } // namespace dalt namespace dalt { namespace poly { template enable_if_t> fft(Vec> &p, bool inv) { static const F PI = std::asin((long double)1) * 2; using cpx = Complex; static Vec> multi_levels(30); int m = Log2Ceil(Size(p)); int n = 1 << m; Assert((1 << m) == Size(p)); int shift = 32 - TrailingZeroNumber(n); for (int i = 1; i < n; i++) { int j = ReverseBit(i << shift); if (i < j) { Swap(p[i], p[j]); } } for (int d = 0; d < m; d++) { int s = 1 << d; int s2 = s << 1; auto &level = multi_levels[d]; if (level.empty()) { // init level.resize(1 << d); for (int j = 0, s = 1 << d; j < s; j++) { level[j] = cpx(Cos(F(PI) / s * j), Sin(F(PI) / s * j)); } } for (int i = 0; i < n; i += s2) { for (int j = 0; j < s; j++) { auto a = i + j; auto b = a + s; auto t = level[j] * p[b]; p[b] = p[a] - t; p[a] = p[a] + t; } } } if (inv) { int i = 0; int j = 0; F tn = n; while (i <= j) { auto pj = p[j]; p[j] = p[i] / tn; if (i != j) { p[i] = pj / tn; } i += 1; j = n - i; } } } template struct FFTConv { static_assert(is_modint_v); using cpx = Complex; using mi = M; using Type = mi; static_assert(is_same_v); static Vec conv(const Vec &a, const Vec &b) { if (Size(a) <= 20 || Size(b) <= 20) { return BruteForceConv::conv(a, b); } if (&a == &b) { return conv2(a); } return conv(a, Size(a), b, Size(b)); } static Vec conv(const Vec &a, int na, const Vec &b, int nb) { let rank_a = na - 1; let rank_b = nb - 1; i32 step = 15; i32 mask = (1 << step) - 1; let n = 1 << Log2Ceil(rank_a + rank_b + 1); Vec a_cpx(n); Vec b_cpx(n); for (int i = 0; i < na; i++) { let x = a[i].value; a_cpx[i] = cpx(x & mask, x >> step); } for (int i = 0; i < nb; i++) { let x = b[i].value; b_cpx[i] = cpx(x & mask, x >> step); } fft(a_cpx, false); fft(b_cpx, false); i32 i = 0; i32 j = 0; while (i <= j) { let ari = a_cpx[i].real(); let aii = a_cpx[i].imag(); let bri = b_cpx[i].real(); let bii = b_cpx[i].imag(); let arj = a_cpx[j].real(); let aij = a_cpx[j].imag(); let brj = b_cpx[j].real(); let bij = b_cpx[j].imag(); let a1r = (ari + arj) / 2; let a1i = (aii - aij) / 2; let a2r = (aii + aij) / 2; let a2i = (arj - ari) / 2; let b1r = (bri + brj) / 2; let b1i = (bii - bij) / 2; let b2r = (bii + bij) / 2; let b2i = (brj - bri) / 2; a_cpx[i] = cpx(a1r * b1r - a1i * b1i - a2r * b2i - a2i * b2r, a1r * b1i + a1i * b1r + a2r * b2r - a2i * b2i); b_cpx[i] = cpx(a1r * b2r - a1i * b2i + a2r * b1r - a2i * b1i, a1r * b2i + a1i * b2r + a2r * b1i + a2i * b1r); if (i != j) { let a1r = (arj + ari) / 2; let a1i = (aij - aii) / 2; let a2r = (aij + aii) / 2; let a2i = (ari - arj) / 2; let b1r = (brj + bri) / 2; let b1i = (bij - bii) / 2; let b2r = (bij + bii) / 2; let b2i = (bri - brj) / 2; a_cpx[j] = cpx(a1r * b1r - a1i * b1i - a2r * b2i - a2i * b2r, a1r * b1i + a1i * b1r + a2r * b2r - a2i * b2i); b_cpx[j] = cpx(a1r * b2r - a1i * b2i + a2r * b1r - a2i * b1i, a1r * b2i + a1i * b2r + a2r * b1i + a2i * b1r); } i += 1; j = n - i; } fft(a_cpx, true); fft(b_cpx, true); i64 modulus = M::modulus(); Vec ans(n); for (int i = 0; i < n; i++) { i64 aa = Round(a_cpx[i].real()); i64 bb = Round(b_cpx[i].real()); i64 cc = Round(a_cpx[i].imag()); ans[i] = (aa + (bb % modulus << 15) + (cc % modulus << 30)) % modulus; } return ans; } static Vec conv2(const Vec &p) { int na, nb; na = nb = Size(p); let rank_a = na - 1; let rank_b = nb - 1; i32 step = 15; i32 mask = (1 << step) - 1; let n = 1 << Log2Ceil(rank_a + rank_b + 1); Vec a_cpx(n); for (int i = 0; i < na; i++) { let x = p[i].value; a_cpx[i] = cpx(x & mask, x >> step); } fft(a_cpx, false); auto b_cpx = a_cpx; i32 i = 0; i32 j = 0; while (i <= j) { let ari = a_cpx[i].real(); let aii = a_cpx[i].imag(); let bri = b_cpx[i].real(); let bii = b_cpx[i].imag(); let arj = a_cpx[j].real(); let aij = a_cpx[j].imag(); let brj = b_cpx[j].real(); let bij = b_cpx[j].imag(); let a1r = (ari + arj) / 2; let a1i = (aii - aij) / 2; let a2r = (aii + aij) / 2; let a2i = (arj - ari) / 2; let b1r = (bri + brj) / 2; let b1i = (bii - bij) / 2; let b2r = (bii + bij) / 2; let b2i = (brj - bri) / 2; a_cpx[i] = cpx(a1r * b1r - a1i * b1i - a2r * b2i - a2i * b2r, a1r * b1i + a1i * b1r + a2r * b2r - a2i * b2i); b_cpx[i] = cpx(a1r * b2r - a1i * b2i + a2r * b1r - a2i * b1i, a1r * b2i + a1i * b2r + a2r * b1i + a2i * b1r); if (i != j) { let a1r = (arj + ari) / 2; let a1i = (aij - aii) / 2; let a2r = (aij + aii) / 2; let a2i = (ari - arj) / 2; let b1r = (brj + bri) / 2; let b1i = (bij - bii) / 2; let b2r = (bij + bii) / 2; let b2i = (bri - brj) / 2; a_cpx[j] = cpx(a1r * b1r - a1i * b1i - a2r * b2i - a2i * b2r, a1r * b1i + a1i * b1r + a2r * b2r - a2i * b2i); b_cpx[j] = cpx(a1r * b2r - a1i * b2i + a2r * b1r - a2i * b1i, a1r * b2i + a1i * b2r + a2r * b1i + a2i * b1r); } i += 1; j = n - i; } fft(a_cpx, true); fft(b_cpx, true); i64 modulus = M::modulus(); Vec ans(n); for (int i = 0; i < n; i++) { i64 aa = Round(a_cpx[i].real()); i64 bb = Round(b_cpx[i].real()); i64 cc = Round(a_cpx[i].imag()); ans[i] = (aa + (bb % modulus << 15) + (cc % modulus << 30)) % modulus; } return ans; } static Vec inverse(Vec p, i32 n) { Extend(p, n); auto dfs = [&](auto &dfs, i32 m) -> Vec { if (m == 1) { return Vec{mi(1) / p[0]}; } i32 prev_mod = (m + 1) / 2; auto C = dfs(dfs, prev_mod); auto AC = conv(p, m, C, prev_mod); AC.resize(m); for (int i = 0; i < m; i++) { AC[i] = mi(0) - AC[i]; } AC[0] = AC[0] + mi(2); auto res = conv(C, AC); res.resize(m); return res; }; auto ans = dfs(dfs, n); ans.resize(n); return ans; } }; template struct is_convolution> { static const bool value = true; }; } // namespace poly } // namespace dalt namespace dalt { namespace math { template Vec InverseBatch(Vec batch) { int n = int(batch.size()); Vec origin = batch; T fact = 1; for (int i = 0; i < n; i++) { if (origin[i] == T(0)) { continue; } T val = batch[i]; batch[i] = fact; fact = fact * val; } T inv_fact = T(1) / fact; for (int i = n - 1; i >= 0; i--) { if (origin[i] == T(0)) { continue; } batch[i] = batch[i] * inv_fact; inv_fact = inv_fact * origin[i]; } return batch; } } // namespace math } // namespace dalt namespace dalt { namespace math { template struct Combination { static_assert(is_modint_v); using Modular = typename T::Modular; Vec fact; Vec inv_fact; Combination(int cap) { cap += 1; fact.resize(cap); inv_fact.resize(cap); fact[0] = T(1); for (int i = 1; i < cap; i++) { T v = T(i); if (v != T(0)) { fact[i] = fact[i - 1] * v; } else { fact[i] = fact[i]; } } T inv = T(1) / fact.back(); for (int i = cap - 1; i >= 0; i--) { T v = T(i); inv_fact[i] = inv; if (v != T(0)) { inv = inv * T(i); } } } T inverse(int x) { return inv_fact[x] * fact[x - 1]; } T combination(int a, int b) { if (a < b || b < 0) { return T(0); } return fact[a] * inv_fact[b] * inv_fact[a - b]; } template T combination_lucas(E a, E b) { if (a < b || b < 0) { return T(0); } if (a < Modular::modulus) { return fact[a] * inv_fact[b] * inv_fact[a - b]; } else { return combination(a % Modular::modulus, b % Modular::modulus) * combination_lucas(a / Modular::modulus, b / Modular::modulus); } } }; } // namespace math } // namespace dalt namespace dalt { namespace poly { const static int POLY_FAST_MAGIC_THRESHOLD = 64; MakeAnnotation(polynomial); template struct Polynomial { static_assert(is_convolution_v); using T = typename Conv::Type; using Type = T; using Seq = Vec; using Self = Polynomial; Seq data; Polynomial(T v = T(0)): Polynomial(Vec{v}) {} Polynomial(Vec _data) : data(Move(_data)) { Normalize(data); } T operator()(T x) const { return Apply(data, x); } T apply(T x) const { return (*this)(x); } Self integral() const { let rank = this->rank(); Vec range(rank + 1); for (int i = 0; i <= rank; i++) { range[i] = i + 1; } let inv = math::InverseBatch(move(range)); Vec ans(rank + 2); for (int i = 0; i <= rank; i++) { ans[i + 1] = inv[i] * data[i]; } return Self(ans); } Self differential() const { Vec ans(rank()); for (int i = 1; i <= Size(ans); i++) { ans[i - 1] = data[i] * T(i); } return Self(ans); } void self_modular(i32 n) { if(data.size() >= n) { data.resize(n); Normalize(data); } } Self modular(i32 n) const { return Self(CopyAndExtend(data, n)); } static Self of(T val) { return Self(Vec{val}); } Self ln(i32 n) const { Assert(data[0] == T(1)); let diff = differential().modular(n); let inv = inverse(n); return (diff * inv).integral().modular(n); } Self exp(i32 n) const { if (n == 0) { return Self::of(T(0)); } auto dfs = [&](auto &dfs, i32 n) -> Self { if (n == 1) { return Self::of(T(1)); } let ans = dfs(dfs, (n + 1) / 2); let ln = this->modular(n) - ans.ln(n); ln.data[0] = ln.data[0] + T(1); return (ans * ln).modular(n); }; return dfs(dfs, n); } int rank() const { return Size(data) - 1; } Self operator*(const Self &rhs) const { const Self &lhs = *this; // if (Min(lhs.rank(), rhs.rank()) < POLY_FAST_MAGIC_THRESHOLD) { // return Self(BruteForceConv::conv(lhs.data, rhs.data)); // } return Self(Conv::conv(lhs.data, rhs.data)); } Self operator*(const T &rhs) const { Vec res = data; for (int i = 0; i < Size(res); i++) { res[i] = res[i] * rhs; } return Self(res); } Self &operator*=(const T &rhs) { for (int i = 0; i < Size(data); i++) { data[i] = data[i] * rhs; } Normalize(data); return *this; } Self operator+(const T &rhs) const { Vec res = data; res[0] = res[0] + rhs; return Self(res); } Self operator+=(const T &rhs) const { data[0] = data[0] + rhs; Normalize(data); return data; } Self operator-(const T &rhs) const { Vec res = data; res[0] = res[0] - rhs; return Self(res); } Self operator-=(const T &rhs) const { data[0] = data[0] - rhs; Normalize(data); return data; } Self operator>>(i32 n) const { if (n < 0) { return *this << -n; } if (*this == Self::of(T(0))) { return Self::of(T(0)); } Vec res(Size(data) + n); for (int i = 0; i < Size(data); i++) { res[i + n] = data[i]; } return Self(res); } Self operator<<(i32 n) const { if (n < 0) { return *this >> -n; } if (Size(data) < n) { return Self::of(T(0)); } Vec res(Size(data) - n); for (int i = 0; i < Size(res); i++) { res[i] = data[i + n]; } return Self(res); } Self operator/(const Self &rhs) const { auto a = data; auto b = rhs.data; if (a.size() < b.size()) { return Self::of(T(0)); } if (b.size() <= POLY_FAST_MAGIC_THRESHOLD) { return BruteForceConv::div_and_rem(Move(a), Move(b))[0]; } Reverse(All(a)); Reverse(All(b)); i32 c_rank = Size(a) - Size(b); i32 proper_len = 1 << Log2Ceil(c_rank * 2 + 1); a.resize(proper_len); b.resize(proper_len); Vec c = Conv::inverse(move(b), c_rank + 1); Vec prod = Conv::conv(move(a), move(c)); prod.resize(c_rank + 1); Reverse(All(prod)); return Self(prod); } Self operator%(const Self &rhs) const { if (Min(rank(), rhs.rank()) < POLY_FAST_MAGIC_THRESHOLD) { return BruteForceConv::div_and_rem(data, rhs.data)[1]; } return *this - (*this / rhs) * rhs; } //return this(x + s) Self shift(T s) const { int r = rank(); var comb = math::Combination(r + 1); Vec A(r + 1), B(r + 1); T s_pow = 1; for(int i = 0; i <= r; i++, s_pow *= s) { A[i] = data[i] * comb.fact[i]; B[i] = s_pow * comb.inv_fact[i]; } var C = Self(Move(A)).delta_convolution(Self(Move(B))); for(int i = 0; i <= C.rank(); i++) { C.data[i] *= comb.inv_fact[i]; } Normalize(C.data); return C; } Self operator+(const Self &rhs) const { const Self &lhs = *this; int n = Size(lhs.data); int m = Size(rhs.data); Vec res(Max(n, m)); for (int i = 0; i < n; i++) { res[i] = lhs.data[i]; } for (int i = 0; i < m; i++) { res[i] = res[i] + rhs.data[i]; } return Self(move(res)); } Self operator-(const Self &rhs) const { const Self &lhs = *this; int n = Size(lhs.data); int m = Size(rhs.data); Vec res(Max(n, m)); for (int i = 0; i < n; i++) { res[i] = lhs.data[i]; } for (int i = 0; i < m; i++) { res[i] = res[i] - rhs.data[i]; } return Self(move(res)); } T operator[](int index) const { return get(index); } T get(int index) const { if (index < Size(data)) { return data[index]; } return T(0); } bool operator==(const_ref(Self) rhs) const { return data == rhs.data; } bool operator!=(const_ref(Self) rhs) const { return !(*this == rhs); } Vec to_vec() const { return data; } Self inverse(int n) const { return Self(Conv::inverse(data, n)); } static Self mulmod(const Self &a, const Self &b, int mod) { return (a * b).modular(mod); } //O(n ln n ln n) Self powmod_binary_lift(i64 n, i32 mod) const { if (n == 0) { return Self::of(T(1)).modular(mod); } Self res = powmod_binary_lift(n / 2, mod); res = (res * res).modular(mod); if (n % 2 == 1) { res = (res * *this).modular(mod); } return res; } //O(n ln n) Self powmod_fast(i32 n_mod_modulus, i32 n_mod_phi, i64 estimate, i32 mod) const { static_assert(is_modint_v); static_assert(is_same_v); if (estimate == 0) { return Self::of(T(1)).modular(mod); } if (*this == Self::of(T(0))) { return *this; } i32 k = 0; while (data[k] == T(0)) { k++; } if (MulLimit(k, estimate, mod) >= mod) { return Self::of(T(0)); } auto expln = [&](const Self &p, i32 n_mod_modulus, i32 mod) -> Self { return (p.ln(mod) *= T(n_mod_modulus)).exp(mod); }; auto expln_ext = [&](Self p, i32 n_mod_modulus, i32 n_mod_phi, i32 mod) -> Self { T val = p[0]; T inv = T(1) / p[0]; p *= inv; Self res = expln(p, n_mod_modulus, mod); res *= PowBinaryLift(val, n_mod_phi); return res; }; if (k == 0) { return expln_ext(*this, n_mod_modulus, n_mod_phi, mod); } Self trim = (*this) << k; Self res = expln_ext(Move(trim), n_mod_modulus, n_mod_phi, mod); return (res >> k * estimate).modular(mod); } static Self product(const Vec &data) { if (data.empty()) { return Self(Vec{Type(1)}); } auto dfs = [&](auto &dfs, int l, int r) { if (l == r) { return data[l]; } int m = (l + r) / 2; return dfs(dfs, l, m) * dfs(dfs, m + 1, r); }; return dfs(dfs, 0, Size(data) - 1); } //ret[i] = \sum_{j} this[i + j] * rhs[j] Self delta_convolution(const Self& rhs) const { Vec lhs = data; Reverse(All(lhs)); auto ans = Conv::conv(lhs, rhs.data); ans.resize(Size(lhs)); Reverse(All(ans)); return Self(Move(ans)); } }; AssignAnnotationTemplate(Polynomial, polynomial, class); } // namespace poly } // namespace dalt namespace dalt { template Indexer MakeIndexer(const C &data) { return [&](auto i) -> T { return data[i]; }; } template Indexer MakeReverseIndexer(const C &data) { return [&](auto i) -> T { return data[Size(data) - 1 - i]; }; } template Vec ExpandIndexer(int n, const Indexer &indexer) { Vec ans; ans.reserve(n); for (int i = 0; i < n; i++) { ans.push_back(indexer(i)); } return ans; } Indexer SelfIndexer() { return [](auto i) { return i; }; } template Indexer ConstantIndexer(const T &val) { return [=](auto i) { return val; }; } template Mapper ConstructorMapper() { return [&](auto a) { return B(a); }; } template Adder NaturalAdder() { return [](auto a, auto b) { return a + b; }; } template constexpr Adder EmptyAdder() { return [](auto a, auto b) { return C(); }; } template constexpr Adder ReturnLeftAdder() { return [](auto a, auto b) { return a; }; } template constexpr Adder ReturnRightAdder() { return [](auto a, auto b) { return b; }; } template Indexer BinaryIndexer(const T& val) { return [=](int i) {return int((val >> i) & 1);}; } template Indexer ReverseIndexer(int n, Indexer indexer) { return [=](int i) {return indexer(n - 1 - i);}; } Vec MakeIndexVec(int n) { Vec ans(n); for(int i = 0; i < n; i++) ans[i] = i; return ans; } } // namespace dalt namespace dalt { namespace poly { // Bostan-Mori Algorithm // [x^k] P/Q // time complexity: O(M(Poly) n) where M(Poly) is the time taken to perform polynomial mutiplication // consider from last bit to first bit (lower bit to higher bit), n - 1 means the lowest bit template typename Poly::Type KthTermOfInversePolynomial(int n, const Indexer &bit_indexer, const Poly &P, const Poly &Q) { static_assert(is_polynomial_v); using T = typename Poly::Type; if (n == 0) { return P[0] / Q[0]; } auto bit = bit_indexer(n - 1); //assert(P.rank() <= 1 && Q.rank() <= 1); Vec neg_Q_data = Q.data; for (int i = 1; i < Size(neg_Q_data); i += 2) { neg_Q_data[i] = T(0) - neg_Q_data[i]; } Poly neg_Q(Move(neg_Q_data)); let AB = P * neg_Q; let QQ = Q * neg_Q; Vec A; Vec C; A.reserve((AB.rank() + 1 + 1) / 2); C.reserve((QQ.rank() + 1 + 1) / 2); for (int i = bit; i <= AB.rank(); i += 2) { A.push_back(AB[i]); } for (int i = 0; i <= QQ.rank(); i += 2) { C.push_back(QQ[i]); } return KthTermOfInversePolynomial(n - 1, bit_indexer, Poly(Move(A)), Poly(Move(C))); } } // namespace poly } // namespace dalt using namespace std; using MOD = StaticModular; using Mi = ModInt; //using Mi = ModInt998244353; using Conv = poly::FFTConv; using Poly = poly::Polynomial; void SolveOne(int test_id, IStream &in, OStream &out) { i32 N; i64 M; in >> N >> M; Vec A(N); in >> A; Vec ps(N); for(int i = 0; i < N; i++) { Vec data(A[i] + 1); data[0] = 1; data[A[i]] = -1; ps[i] = Move(data); } var prod = Poly::product(ps); Vec bits; { i64 M2 = M; while(M2 > 0) { bits.push_back(M2 % 2); M2 /= 2; } Reverse(All(bits)); } var ans = poly::KthTermOfInversePolynomial(Size(bits), [&](int i) {return bits[i];}, Poly(1), prod); out << ans; } void SolveMulti(IStream &in, OStream &out) { //std::ifstream input("in"); int num_of_input = 1; //in >> num_of_input; for (int i = 0; i < num_of_input; i++) { //SolveOne(i + 1, input, out); SolveOne(i + 1, in, out); } } int main() { std::ios_base::sync_with_stdio(false); Stdin.tie(0); Stdout << std::setiosflags(std::ios::fixed); Stdout << std::setprecision(15); #ifdef STRESS stress::Stress(); #else SolveMulti(Stdin, Stdout); #endif return 0; }