結果
問題 | No.2747 Permutation Adjacent Sum |
ユーザー | shiomusubi496 |
提出日時 | 2024-04-21 02:10:45 |
言語 | C++17 (gcc 12.3.0 + boost 1.83.0) |
結果 |
MLE
|
実行時間 | - |
コード長 | 54,350 bytes |
コンパイル時間 | 5,220 ms |
コンパイル使用メモリ | 279,816 KB |
実行使用メモリ | 814,892 KB |
最終ジャッジ日時 | 2024-10-13 00:40:27 |
合計ジャッジ時間 | 6,901 ms |
ジャッジサーバーID (参考情報) |
judge5 / judge4 |
(要ログイン)
テストケース
テストケース表示入力 | 結果 | 実行時間 実行使用メモリ |
---|---|---|
testcase_00 | MLE | - |
testcase_01 | -- | - |
testcase_02 | -- | - |
testcase_03 | -- | - |
testcase_04 | -- | - |
testcase_05 | -- | - |
testcase_06 | -- | - |
testcase_07 | -- | - |
testcase_08 | -- | - |
testcase_09 | -- | - |
testcase_10 | -- | - |
testcase_11 | -- | - |
testcase_12 | -- | - |
testcase_13 | -- | - |
testcase_14 | -- | - |
testcase_15 | -- | - |
testcase_16 | -- | - |
testcase_17 | -- | - |
testcase_18 | -- | - |
testcase_19 | -- | - |
testcase_20 | -- | - |
testcase_21 | -- | - |
testcase_22 | -- | - |
testcase_23 | -- | - |
testcase_24 | -- | - |
testcase_25 | -- | - |
testcase_26 | -- | - |
testcase_27 | -- | - |
testcase_28 | -- | - |
testcase_29 | -- | - |
testcase_30 | -- | - |
testcase_31 | -- | - |
testcase_32 | -- | - |
testcase_33 | -- | - |
testcase_34 | -- | - |
testcase_35 | -- | - |
testcase_36 | -- | - |
testcase_37 | -- | - |
testcase_38 | -- | - |
testcase_39 | -- | - |
testcase_40 | -- | - |
testcase_41 | -- | - |
ソースコード
#line 2 "library/other/template.hpp" #include <bits/stdc++.h> #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<ll, ll>; template<class T> using prique = std::priority_queue<T, std::vector<T>, std::greater<T>>; template<class T> struct infinity { static constexpr T value = std::numeric_limits<T>::max() / 2; static constexpr T mvalue = std::numeric_limits<T>::lowest() / 2; static constexpr T max = std::numeric_limits<T>::max(); static constexpr T min = std::numeric_limits<T>::lowest(); }; #if __cplusplus <= 201402L template<class T> constexpr T infinity<T>::value; template<class T> constexpr T infinity<T>::mvalue; template<class T> constexpr T infinity<T>::max; template<class T> constexpr T infinity<T>::min; #endif #if __cpp_variable_templates >= 201304L template<class T> constexpr T INF = infinity<T>::value; #endif constexpr ll inf = infinity<ll>::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<class T, class... Args> struct function_traits_impl { using result_type = T; template<std::size_t idx> using argument_type = typename std::tuple_element<idx, std::tuple<Args...>>::type; using argument_tuple = std::tuple<Args...>; static constexpr std::size_t arg_size() { return sizeof...(Args); } }; template<class> struct function_traits_helper; template<class Res, class Tp, class... Args> struct function_traits_helper<Res (Tp::*)(Args...)> { using type = function_traits_impl<Res, Args...>; }; template<class Res, class Tp, class... Args> struct function_traits_helper<Res (Tp::*)(Args...)&> { using type = function_traits_impl<Res, Args...>; }; template<class Res, class Tp, class... Args> struct function_traits_helper<Res (Tp::*)(Args...) const> { using type = function_traits_impl<Res, Args...>; }; template<class Res, class Tp, class... Args> struct function_traits_helper<Res (Tp::*)(Args...) const&> { using type = function_traits_impl<Res, Args...>; }; #if __cpp_noexcept_function_type >= 201510L template<class Res, class Tp, class... Args> struct function_traits_helper<Res (Tp::*)(Args...) noexcept> { using type = function_traits_impl<Res, Args...>; }; template<class Res, class Tp, class... Args> struct function_traits_helper<Res (Tp::*)(Args...)& noexcept> { using type = function_traits_impl<Res, Args...>; }; template<class Res, class Tp, class... Args> struct function_traits_helper<Res (Tp::*)(Args...) const noexcept> { using type = function_traits_impl<Res, Args...>; }; template<class Res, class Tp, class... Args> struct function_traits_helper<Res (Tp::*)(Args...) const& noexcept> { using type = function_traits_impl<Res, Args...>; }; #endif template<class F> using function_traits = typename function_traits_helper< decltype(&std::remove_reference<F>::type::operator())>::type; template<class F> using function_result_type = typename function_traits<F>::result_type; template<class F, std::size_t idx> using function_argument_type = typename function_traits<F>::template argument_type<idx>; template<class F> using function_argument_tuple = typename function_traits<F>::argument_tuple; template<class T> using is_signed_int = std::integral_constant<bool, (std::is_integral<T>::value && std::is_signed<T>::value) || std::is_same<T, i128>::value>; template<class T> using is_unsigned_int = std::integral_constant<bool, (std::is_integral<T>::value && std::is_unsigned<T>::value) || std::is_same<T, u128>::value>; template<class T> using is_int = std::integral_constant<bool, is_signed_int<T>::value || is_unsigned_int<T>::value>; template<class T> using make_signed_int = typename std::conditional< std::is_same<T, i128>::value || std::is_same<T, u128>::value, std::common_type<i128>, std::make_signed<T>>::type; template<class T> using make_unsigned_int = typename std::conditional< std::is_same<T, i128>::value || std::is_same<T, u128>::value, std::common_type<u128>, std::make_unsigned<T>>::type; template<class T, class = void> struct is_range : std::false_type {}; template<class T> struct is_range< T, decltype(all(std::declval<typename std::add_lvalue_reference<T>::type>()), (void)0)> : std::true_type {}; template<class T, bool = is_range<T>::value> struct range_rank : std::integral_constant<std::size_t, 0> {}; template<class T> struct range_rank<T, true> : std::integral_constant<std::size_t, range_rank<typename T::value_type>::value + 1> {}; template<std::size_t size> 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<size <= 64, std::int_least64_t, i128>::type>::type>::type>::type; }; template<std::size_t size> using int_least_t = typename int_least<size>::type; template<std::size_t size> 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<size <= 64, std::uint_least64_t, u128>::type>::type>::type>::type; }; template<std::size_t size> using uint_least_t = typename uint_least<size>::type; template<class T> using double_size_int = int_least<std::numeric_limits<T>::digits * 2 + 1>; template<class T> using double_size_int_t = typename double_size_int<T>::type; template<class T> using double_size_uint = uint_least<std::numeric_limits<T>::digits * 2>; template<class T> using double_size_uint_t = typename double_size_uint<T>::type; template<class T> using double_size = typename std::conditional<is_signed_int<T>::value, double_size_int<T>, double_size_uint<T>>::type; template<class T> using double_size_t = typename double_size<T>::type; #line 7 "library/other/template.hpp" // #include "../template/in.hpp" // #include "../template/out.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<class T, class U, class Comp = std::less<>> 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<class T, class U, class Comp = std::less<>> 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 bool is_prime(ll N) { if (N <= 1) return false; for (ll i = 2; i * i <= N; ++i) { if (N % i == 0) return false; } return true; } inline std::vector<ll> prime_factor(ll N) { std::vector<ll> res; for (ll i = 2; i * i <= N; ++i) { while (N % i == 0) { res.push_back(i); N /= i; } } if (N != 1) res.push_back(N); return res; } inline CONSTEXPR ll my_pow(ll a, ll b) { ll res = 1; while (b) { if (b & 1) res *= a; b >>= 1; a *= a; } return res; } 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<class F> class RecLambda { private: F f; public: explicit constexpr RecLambda(F&& f_) : f(std::forward<F>(f_)) {} template<class... Args> constexpr auto operator()(Args&&... args) -> decltype(f(*this, std::forward<Args>(args)...)) { return f(*this, std::forward<Args>(args)...); } }; template<class F> inline constexpr RecLambda<F> rec_lambda(F&& f) { return RecLambda<F>(std::forward<F>(f)); } template<class Head, class... Tail> struct multi_dim_vector { using type = std::vector<typename multi_dim_vector<Tail...>::type>; }; template<class T> struct multi_dim_vector<T> { using type = T; }; template<class T, class Arg> constexpr std::vector<T> make_vec(int n, Arg&& arg) { return std::vector<T>(n, std::forward<Arg>(arg)); } template<class T, class... Args> constexpr typename multi_dim_vector<Args..., T>::type make_vec(int n, Args&&... args) { return typename multi_dim_vector<Args..., T>::type( n, make_vec<T>(std::forward<Args>(args)...)); } #line 2 "library/math/Combinatorics.hpp" #line 2 "library/math/ModInt.hpp" #line 4 "library/math/ModInt.hpp" template<class T, T mod> class StaticModInt { static_assert(std::is_integral<T>::value, "T must be integral"); static_assert(std::is_unsigned<T>::value, "T must be unsigned"); static_assert(mod > 0, "mod must be positive"); static_assert(mod <= std::numeric_limits<T>::max() / 2, "mod * 2 must be less than or equal to T::max()"); private: using large_t = typename double_size_uint<T>::type; using signed_t = typename std::make_signed<T>::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<class U, typename std::enable_if<std::is_integral<U>::value && std::is_signed<U>::value>::type* = nullptr> constexpr StaticModInt(U v) : val{} { v %= static_cast<signed_t>(mod); if (v < 0) v += static_cast<signed_t>(mod); val = static_cast<T>(v); } template<class U, typename std::enable_if< std::is_integral<U>::value && std::is_unsigned<U>::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<class Pr> void print(Pr& a) const { a.print(val); } template<class Pr> void debug(Pr& a) const { a.print(val); } template<class Sc> void scan(Sc& a) { ll v; a.scan(v); *this = v; } }; #if __cplusplus < 201703L template<class T, T mod> constexpr unsigned int StaticModInt<T, mod>::inv1000000007[]; template<class T, T mod> constexpr unsigned int StaticModInt<T, mod>::inv998244353[]; #endif template<unsigned int p> using static_modint = StaticModInt<unsigned int, p>; using modint1000000007 = static_modint<1000000007>; using modint998244353 = static_modint<998244353>; /** * @brief ModInt * @docs docs/math/ModInt.md */ #line 5 "library/math/Combinatorics.hpp" template<class T> class Combinatorics { private: static std::vector<T> factorial; static std::vector<T> 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<class T> std::vector<T> Combinatorics<T>::factorial = std::vector<T>(1, 1); template<class T> std::vector<T> Combinatorics<T>::factinv = std::vector<T>(1, 1); /** * @brief Combinatorics * @docs docs/math/Combinatorics.md */ #line 2 "library/math/poly/SamplingPointsShift.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 Engine> 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<class IntType = ll> IntType uniform(IntType l, IntType r) { static_assert(std::is_integral<IntType>::value, "template argument must be an integral type"); assert(l <= r); return std::uniform_int_distribution<IntType>{l, r}(rnd); } template<class RealType = double> RealType uniform_real(RealType l, RealType r) { static_assert(std::is_floating_point<RealType>::value, "template argument must be an floating point type"); assert(l <= r); return std::uniform_real_distribution<RealType>{l, r}(rnd); } bool uniform_bool() { return uniform<int>(0, 1) == 1; } template<class T = ll> std::pair<T, T> uniform_pair(T l, T r) { assert(l < r); T a, b; do { a = uniform<T>(l, r); b = uniform<T>(l, r); } while (a == b); if (a > b) swap(a, b); return {a, b}; } template<class T = ll> std::vector<T> choice(int n, T l, T r) { assert(l <= r); assert(T(n) <= (r - l + 1)); std::set<T> res; while ((int)res.size() < n) res.insert(uniform<T>(l, r)); return {res.begin(), res.end()}; } template<class Iter> void shuffle(const Iter& first, const Iter& last) { std::shuffle(first, last, rnd); } template<class T> std::vector<T> permutation(T n) { std::vector<T> res(n); rep (i, n) res[i] = i; shuffle(all(res)); return res; } template<class T = ll> std::vector<T> choice_shuffle(int n, T l, T r, bool sorted = true) { assert(l <= r); assert(T(n) <= (r - l + 1)); std::vector<T> 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<std::mt19937>; Random32 rand32; using Random64 = Random<std::mt19937_64>; Random64 rand64; /** * @brief Random * @docs docs/random/Random.md */ #line 2 "library/math/MontgomeryModInt.hpp" #line 4 "library/math/MontgomeryModInt.hpp" template<class T> class MontgomeryReduction { static_assert(std::is_integral<T>::value, "T must be integral"); static_assert(std::is_unsigned<T>::value, "T must be unsigned"); private: using large_t = typename double_size_uint<T>::type; static constexpr int lg = std::numeric_limits<T>::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<T>(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<T>::max() / 2); mod = v; r = (-static_cast<T>(mod)) % mod; r2 = (-static_cast<large_t>(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<large_t>(static_cast<T>(x) * minv) * mod) >> lg; return tmp >= mod ? tmp - mod : tmp; } T transform(large_t x) const { return reduce(x * r2); } }; template<class T, int id> class MontgomeryModInt { private: using large_t = typename double_size_uint<T>::type; using signed_t = typename std::make_signed<T>::type; T val; static MontgomeryReduction<T> mont; public: MontgomeryModInt() : val(0) {} template<class U, typename std::enable_if< std::is_integral<U>::value && std::is_unsigned<U>::value>::type* = nullptr> MontgomeryModInt(U x) : val(mont.transform( x < (static_cast<large_t>(mont.get_mod()) << mont.get_lg()) ? x : x % mont.get_mod())) {} template<class U, typename std::enable_if<std::is_integral<U>::value && std::is_signed<U>::value>::type* = nullptr> MontgomeryModInt(U x) : MontgomeryModInt(static_cast<typename std::make_unsigned<U>::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<large_t>(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<class Pr> void print(Pr& a) const { a.print(mont.reduce(val)); } template<class Pr> void debug(Pr& a) const { a.print(mont.reduce(val)); } template<class Sc> void scan(Sc& a) { ll v; a.scan(v); *this = v; } }; template<class T, int id> MontgomeryReduction<T> MontgomeryModInt<T, id>::mont = MontgomeryReduction<T>(998244353); using mmodint = MontgomeryModInt<unsigned int, -1>; /** * @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<class T> 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<MontgomeryModInt<unsigned int, -2>>( n, base_miller_rabin_int, 3); return miller_rabin<MontgomeryModInt<ull, -2>>(n, base_miller_rabin_ll, 7); } #if __cpp_variable_templates >= 201304L && __cpp_constexpr >= 201304L template<ull n> 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<class Cont, class Comp> std::vector<std::pair<typename Cont::value_type, int>> RunLength(const Cont& str, const Comp& cmp) { std::vector<std::pair<typename Cont::value_type, int>> 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<class Cont> std::vector<std::pair<typename Cont::value_type, int>> RunLength(const Cont& str) { return RunLength(str, std::equal_to<typename Cont::value_type>()); } /** * @brief RunLength(ランレングス圧縮) * @docs docs/string/RunLength.md */ #line 8 "library/math/PollardRho.hpp" template<class T, class Rnd> 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 T = MontgomeryModInt<ull, -3>, class Rnd = Random64> std::vector<ull> factorize(ull n, Rnd& rnd = rand64) { if (n == 1) return {}; std::vector<ull> res; std::vector<ull> 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>(t, rnd); st.push_back(f); st.push_back(t / f); } std::sort(all(res)); return res; } template<class T = MontgomeryModInt<ull, -3>, class Rnd = Random64> std::vector<std::pair<ull, int>> expfactorize(ull n, Rnd& rnd = rand64) { auto f = factorize<T, Rnd>(n, rnd); return RunLength(f); } std::vector<ll> divisors_pr(ll n) { std::vector<ll> res; auto r = expfactorize(n); int m = r.size(); rec_lambda([&](auto&& self, int k, ll d) -> void { if (k == m) { res.push_back(d); return; } ll t = d; rep (r[k].second) { self(k + 1, d); d *= r[k].first; } self(k + 1, d); d = t; })(0, 1); std::sort(all(res)); return res; } /** * @brief PollardRho(素因数分解) * @docs docs/math/PollardRho.md */ #line 9 "library/math/PrimitiveRoot.hpp" template<class T = MontgomeryModInt<ull, -4>> 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<unsigned int p> 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<unsigned int p> constexpr NthRoot<p> nth_root; template<class T> void number_theoretic_transform(std::vector<T>& 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<T::get_mod()>.get_rate(popcnt(j & ~(j + 1))); } } } } template<class T> void inverse_number_theoretic_transform(std::vector<T>& 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<T::get_mod()>.get_inv_rate(popcnt(j & ~(j + 1))); } } } T inv_n = T(1) / n; each_for (x : a) x *= inv_n; } template<class T> std::vector<T> convolution_naive(const std::vector<T>& a, const std::vector<T>& b) { int n = a.size(), m = b.size(); std::vector<T> c(n + m - 1); rep (i, n) rep (j, m) c[i + j] += a[i] * b[j]; return c; } template<class T> std::vector<T> convolution_pow2(std::vector<T> 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<T> a1(a.begin(), a.begin() + m), a2(a.begin() + m, a.end()); std::vector<T> c(n); std::vector<T> c1 = convolution_naive(a1, a1); std::vector<T> c2 = convolution_naive(a1, a2); std::vector<T> 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<class T> std::vector<T> convolution(std::vector<T> a, std::vector<T> b) { int n = a.size() + b.size() - 1; int lg = bitop::ceil_log2(n); 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<T> a1(a.begin(), a.begin() + m), a2(a.begin() + m, a.end()); std::vector<T> c(n); std::vector<T> c1 = convolution_naive(a1, b); std::vector<T> 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<unsigned int p> std::vector<static_modint<p>> convolution_for_any_mod(const std::vector<static_modint<p>>& a, const std::vector<static_modint<p>>& b); template<unsigned int p> std::vector<static_modint<p>> convolution(const std::vector<static_modint<p>>& a, const std::vector<static_modint<p>>& 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<unsigned int p> std::vector<ll> convolution(const std::vector<ll>& a, const std::vector<ll>& b) { int n = a.size(), m = b.size(); std::vector<static_modint<p>> 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<ll> c(n + m - 1); rep (i, n + m - 1) c[i] = c2[i].get(); return c; } template<unsigned int p> std::vector<static_modint<p>> convolution_for_any_mod(const std::vector<static_modint<p>>& a, const std::vector<static_modint<p>>& b) { int n = a.size(), m = b.size(); assert(n + m - 1 <= (1 << 26)); std::vector<ll> 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<MOD1>(a2, b2); auto c2 = convolution<MOD2>(a2, b2); auto c3 = convolution<MOD3>(a2, b2); std::vector<static_modint<p>> 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<p>(t1 + (t2 + t3 * MOD2) % p * MOD1); } return res; } template<class T> void ntt_doubling_(std::vector<T>& a) { int n = a.size(); auto b = a; inverse_number_theoretic_transform(b); const T z = internal::nth_root<T::get_mod()>.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<unsigned int p> 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/poly/FormalPowerSeries.hpp" #line 2 "library/math/SqrtMod.hpp" #line 5 "library/math/SqrtMod.hpp" template<class T> 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<ll>(x.get(), p - x.get()); } ll sqrt_mod(ll a, ll p) { if (p == 2) return a; using mint = MontgomeryModInt<unsigned int, 493174342>; mint::set_mod(p); return sqrt_mod<mint>(a); } /** * @brief SqrtMod(平方剰余) * @docs docs/math/SqrtMod.md * @see https://37zigen.com/tonelli-shanks-algorithm/ */ #line 7 "library/math/poly/FormalPowerSeries.hpp" template<class T> class FormalPowerSeries : public std::vector<T> { private: using Base = std::vector<T>; using Comb = Combinatorics<T>; 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<FormalPowerSeries, FormalPowerSeries> 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<bool AlwaysTrue = true, typename std::enable_if< AlwaysTrue && is_ntt_friendly<T::get_mod()>::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<bool AlwaysTrue = true, typename std::enable_if< AlwaysTrue && !is_ntt_friendly<T::get_mod()>::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<bool AlwaysTrue = true, typename std::enable_if< AlwaysTrue && is_ntt_friendly<T::get_mod()>::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<class T, class Comb = Combinatorics<T>> FormalPowerSeries<T> taylor_shift(FormalPowerSeries<T> f, T a) { const int n = f.size(); Comb::init(n); rep (i, n) f[i] *= Comb::fact(i); FormalPowerSeries<T> 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 8 "library/math/poly/SamplingPointsShift.hpp" template<class T, class Comb = Combinatorics<T>> std::vector<T> sampling_points_shift(std::vector<T> a, int m, T t) { const int n = a.size(); Comb::init(std::max(n, m)); FormalPowerSeries<T> f(n), g(n); rep (i, n) f[i] = i & 1 ? -Comb::finv(i) : Comb::finv(i); rep (i, n) g[i] = a[i] * Comb::finv(i); f = (f * g).prefix(n); rep (i, n) f[i] *= Comb::fact(i); T p = 1; rep (i, n) { g[n - 1 - i] = p * Comb::finv(i); p *= t--; } f = (f * g) >> (n - 1); rep (i, n) f[i] *= Comb::finv(i); g.resize(m); rep (i, m) g[i] = Comb::finv(i); f = (f * g).prefix(m); rep (i, m) f[i] *= Comb::fact(i); return std::vector<T>(f); } /** * @brief SamplingPointsShift(標本点シフト) * @docs docs/math/poly/SamplingPointsShift.md */ #line 4 "main.cpp" using namespace std; using mint = modint998244353; using comb = Combinatorics<mint>; int main() { ll N, K; cin >> N >> K; mint ans = 0; if (N <= K + 10) { rep (i, 1, N) ans += mint{i}.pow(K) * (N - i) * 2 * comb::fact(N - 1); cout << ans.get() << endl; return 0; } vector<mint> a(K + 3); rep (i, 1, K + 3) { a[i] = mint{i}.pow(K) * (N - i) * 2 * comb::fact(N - 1) + a[i - 1]; } cout << sampling_points_shift(a, 1, mint{N})[0].get() << endl; }