結果
| 問題 | No.3799 1 + 6 = 7 (Another version) |
| ユーザー |
 miscalc
|
| 提出日時 | 2026-11-14 00:45:58 |
| 言語 | C++23 (gcc 15.3.0 + boost 1.92.0 + ACL) |
| 結果 |
AC
不安定
|
| 実行時間 | 1 ms / 2,000 ms |
| + 729µs | |
| コード長 | 109,978 bytes |
| 記録 | |
| コンパイル時間 | 6,618 ms |
| コンパイル使用メモリ | 442,832 KB |
| 実行使用メモリ | 9,924 KB |
| 最終ジャッジ日時 | 2026-11-14 00:46:07 |
| 合計ジャッジ時間 | 8,979 ms |
|
ジャッジサーバーID （参考情報） |
judge3_0 / judge5_0 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 3 |
| other | AC * 36 |
ソースコード
// https://github.com/miscalculation53/library
#define SINGLE_TESTCASE
// #define MULTI_TESTCASE
// #define AOJ_TESTCASE
#define FAST_IO
// #define FAST_CIO
// #define INTERACTIVE
#define INF 4'000'000'000'000'000'037LL
#define EPS 1e-11
#line 30
#include <bits/stdc++.h>
using namespace std;
#if __cplusplus > 202002L && defined(__STDCPP_FLOAT128_T__)
#include <stdfloat>
using f128 = std::float128_t;
#elif defined(__SIZEOF_FLOAT128__)
using f128 = __float128;
#else
#error "f128 requires a compiler with IEEE binary128 support, such as GCC on x86_64."
#endif
#define LIBRARY_HAS_FLOAT128 1
#include <cassert>
using namespace std;
using ld = decltype(EPS);
using ll = long long;
using uint = unsigned int;
using ull = unsigned long long;
using pll = pair<ll, ll>;
using tlll = tuple<ll, ll, ll>;
using tllll = tuple<ll, ll, ll, ll>;
#define vc vector
template <class T>
using vvc = vc<vc<T>>;
template <class T>
using vvvc = vc<vc<vc<T>>>;
using vb = vc<bool>;
using vl = vc<ll>;
using vpll = vc<pll>;
using vtlll = vc<tlll>;
using vtllll = vc<tllll>;
using vstr = vc<string>;
using vvb = vvc<bool>;
using vvl = vvc<ll>;
template <class T>
using pql = priority_queue<T, vc<T>, greater<T>>;
template <class T>
using pqg = priority_queue<T>;
using i128 = __int128_t;
using u128 = __uint128_t;
#define cauto const auto
#line 103
//! https://gcc.gnu.org/onlinedocs/libquadmath/Math-Library-Routines.html
#line 289
i128 stoi128(const string &s)
{
const bool neg = s.front() == '-';
u128 res = 0;
for (size_t i = size_t(neg); i < s.size(); i++)
res = 10 * res + unsigned(s[i] - '0');
if (neg)
return -i128(res - 1) - 1;
return i128(res);
}
inline string u128tos(u128 x)
{
string res;
do
{
res += char('0' + x % 10);
x /= 10;
} while (x > 0);
reverse(res.begin(), res.end());
return res;
}
string i128tos(i128 x)
{
return x < 0 ? "-" + u128tos(u128(-(x + 1)) + 1) : u128tos(u128(x));
}
istream &operator>>(istream &is, i128 &a)
{
string s;
if (!(is >> s)) return is;
a = stoi128(s);
return is;
}
ostream &operator<<(ostream &os, const i128 &a)
{
os << i128tos(a);
return os;
}
#line 337
#define overload4(_1, _2, _3, _4, name, ...) name
#define LIBRARY_REPEAT_LOOP(T, i, l, r, d) \
for (auto [i, library_repeat_end, library_repeat_step] = array<T, 3>{T(l), T(r), T(d)}; \
library_repeat_step > 0 ? i < library_repeat_end : i > library_repeat_end; \
i += library_repeat_step)
#define rep1(i, n) LIBRARY_REPEAT_LOOP(ll, i, 0, n, 1)
#define rep2(i, l, r) LIBRARY_REPEAT_LOOP(ll, i, l, r, 1)
#define rep3(i, l, r, d) LIBRARY_REPEAT_LOOP(ll, i, l, r, d)
#define rep(...) overload4(__VA_ARGS__, rep3, rep2, rep1)(__VA_ARGS__)
#define repi1(i, n) LIBRARY_REPEAT_LOOP(int, i, 0, n, 1)
#define repi2(i, l, r) LIBRARY_REPEAT_LOOP(int, i, l, r, 1)
#define repi3(i, l, r, d) LIBRARY_REPEAT_LOOP(int, i, l, r, d)
#define repi(...) overload4(__VA_ARGS__, repi3, repi2, repi1)(__VA_ARGS__)
#define fe(...) for (auto __VA_ARGS__)
#define fec(...) for (cauto &__VA_ARGS__)
#define fem(...) for (auto &__VA_ARGS__)
#line 408
template <class T, class U, class P>
constexpr pair<T, U> &operator+=(pair<T, U> &a, const P &b)
{
a.first += b.first;
a.second += b.second;
return a;
}
template <class T, class U, class P>
constexpr pair<T, U> operator+(pair<T, U> a, const P &b) { return a += b; }
template <class T, class U, class P>
constexpr pair<T, U> &operator-=(pair<T, U> &a, const P &b)
{
a.first -= b.first;
a.second -= b.second;
return a;
}
template <class T, class U, class P>
constexpr pair<T, U> operator-(pair<T, U> a, const P &b) { return a -= b; }
template <class T, class U>
constexpr pair<T, U> operator-(pair<T, U> a)
{
a.first = -a.first;
a.second = -a.second;
return a;
}
template <class T, size_t n, class A>
constexpr array<T, n> &operator+=(array<T, n> &a, const A &b)
{
for (size_t i = 0; i < n; i++)
a[i] += b[i];
return a;
}
template <class T, size_t n, class A>
constexpr array<T, n> operator+(array<T, n> a, const A &b) { return a += b; }
template <class T, size_t n, class A>
constexpr array<T, n> &operator-=(array<T, n> &a, const A &b)
{
for (size_t i = 0; i < n; i++)
a[i] -= b[i];
return a;
}
template <class T, size_t n, class A>
constexpr array<T, n> operator-(array<T, n> a, const A &b) { return a -= b; }
template <class T, size_t n>
constexpr array<T, n> operator-(array<T, n> a)
{
for (auto &ai : a)
ai = -ai;
return a;
}
namespace internal
{
template <size_t... I, class A, class B>
constexpr auto &tuple_add_impl(A &a, const B &b, const index_sequence<I...>)
;
template <size_t... I, class A, class B>
constexpr auto &tuple_sub_impl(A &a, const B &b, const index_sequence<I...>)
;
template <size_t... I, class A>
constexpr auto &tuple_neg_impl(A &a, const index_sequence<I...>)
;
};
template <class... Ts, class Tp>
constexpr tuple<Ts...> &operator+=(tuple<Ts...> &a, const Tp &b)
{ return internal::tuple_add_impl(a, b, make_index_sequence<tuple_size_v<tuple<Ts...>>>{}); }
template <class... Ts, class Tp>
constexpr tuple<Ts...> operator+(tuple<Ts...> a, const Tp &b) { return a += b; }
template <class... Ts, class Tp>
constexpr tuple<Ts...> &operator-=(tuple<Ts...> &a, const Tp &b)
{ return internal::tuple_sub_impl(a, b, make_index_sequence<tuple_size_v<tuple<Ts...>>>{}); }
template <class... Ts, class Tp>
constexpr tuple<Ts...> operator-(tuple<Ts...> a, const Tp &b) { return a -= b; }
template <class... Ts>
constexpr tuple<Ts...> operator-(tuple<Ts...> a)
{
internal::tuple_neg_impl(a, make_index_sequence<sizeof...(Ts)>{});
return a;
}
#line 512
template <class T>
constexpr bool is_integral_ext = is_integral_v<T> || is_same_v<T, i128> || is_same_v<T, u128>;
template <class T>
constexpr bool is_signed_ext = is_signed_v<T> || is_same_v<T, i128>;
template <class T>
constexpr bool is_unsigned_ext = is_unsigned_v<T> || is_same_v<T, u128>;
#line 531
template <class T>
struct make_unsigned_ext
{
using type = make_unsigned_t<T>;
};
template <>
struct make_unsigned_ext<i128>
{
using type = u128;
};
template <>
struct make_unsigned_ext<u128>
{
using type = u128;
};
template <class T>
using make_unsigned_ext_t = typename make_unsigned_ext<T>::type;
template <class T>
struct make_signed_ext
{
using type = make_signed_t<T>;
};
template <>
struct make_signed_ext<u128>
{
using type = i128;
};
template <class T>
using make_signed_ext_t = typename make_signed_ext<T>::type;
template <class T>
constexpr bool is_negative(T x)
{
static_assert(is_integral_ext<T>);
if constexpr (is_signed_ext<T>) return x < 0;
return false;
}
template <class T>
constexpr make_unsigned_ext_t<T> unsigned_abs(T x)
;
template <class T, bool = is_integral_ext<T>>
struct unsigned_arithmetic { using type = T; };
template <class T>
struct unsigned_arithmetic<T, true>
{
static_assert(!is_same_v<T, bool>);
using type = common_type_t<unsigned int, make_unsigned_ext_t<T>>;
};
template <class T>
using unsigned_arithmetic_t = typename unsigned_arithmetic<T>::type;
template <class R, class T>
constexpr R wrapping_cast(T x)
;
#line 625
template <class T>
constexpr decltype(auto) default_infty();
namespace default_infty_detail
{
template <class T>
struct is_product : false_type {};
template <class T, class U>
struct is_product<pair<T, U>> : true_type {};
template <class... Ts>
struct is_product<tuple<Ts...>> : true_type {};
template <class T, size_t n>
struct is_product<array<T, n>> : true_type {};
template <class T, size_t... I>
constexpr T product_value(index_sequence<I...>)
;
template <class T, class = void>
struct has_infty : false_type {};
template <class T>
struct has_infty<T, void_t<decltype(T::infty())>> : true_type {};
template <class T>
const T &custom_value()
;
template <class T>
inline constexpr bool unsupported = false;
}
template <class T>
constexpr decltype(auto) default_infty()
;
#line 696
template <class T, class U>
inline bool chmin(T &a, U b) ;
template <class T, class U>
inline bool chmax(T &a, U b) ;
template <class T, enable_if_t<is_same_v<T, i128> || is_same_v<T, u128>, int> = 0>
inline constexpr T abs(T x) ;
namespace template_math_detail
{
template <class T, class U>
constexpr common_type_t<T, U> gcd(T a, U b)
;
#line 720
template <class T, class U>
constexpr common_type_t<T, U> lcm(T a, U b)
;
}
template <class T, enable_if_t<is_integral_ext<T> && !is_same_v<T, bool>, int> = 0>
constexpr common_type_t<i128, T> gcd(i128 a, T b) ;
template <class T, enable_if_t<is_integral_ext<T> && !is_same_v<T, bool> && !is_same_v<T, i128> && !is_same_v<T, u128>, int> = 0>
constexpr common_type_t<T, i128> gcd(T a, i128 b) ;
template <class T, enable_if_t<is_integral_ext<T> && !is_same_v<T, bool>, int> = 0>
constexpr common_type_t<u128, T> gcd(u128 a, T b) ;
template <class T, enable_if_t<is_integral_ext<T> && !is_same_v<T, bool> && !is_same_v<T, i128> && !is_same_v<T, u128>, int> = 0>
constexpr common_type_t<T, u128> gcd(T a, u128 b) ;
template <class T, enable_if_t<is_integral_ext<T> && !is_same_v<T, bool>, int> = 0>
constexpr common_type_t<i128, T> lcm(i128 a, T b) ;
template <class T, enable_if_t<is_integral_ext<T> && !is_same_v<T, bool> && !is_same_v<T, i128> && !is_same_v<T, u128>, int> = 0>
constexpr common_type_t<T, i128> lcm(T a, i128 b) ;
template <class T, enable_if_t<is_integral_ext<T> && !is_same_v<T, bool>, int> = 0>
constexpr common_type_t<u128, T> lcm(u128 a, T b) ;
template <class T, enable_if_t<is_integral_ext<T> && !is_same_v<T, bool> && !is_same_v<T, i128> && !is_same_v<T, u128>, int> = 0>
constexpr common_type_t<T, u128> lcm(T a, u128 b) ;
template <class T = ll, class U, class V, typename = enable_if_t<is_integral_ext<U> && is_integral_ext<V>>>
inline constexpr T divfloor(U a, V b) ;
template <class T = ll, class U, class V, typename = enable_if_t<is_integral_ext<U> && is_integral_ext<V>>>
inline constexpr T divceil(U a, V b) ;
template <class T = ll, class U, class V, typename = enable_if_t<is_integral_ext<U> && is_integral_ext<V>>>
inline constexpr T safemod(U a, V b)
{
const T x = T(a), y = T(b);
if constexpr (is_signed_ext<T>)
if (y == T(-1)) return T(0);
const T r = x % y;
return r != T(0) && is_negative(x) != is_negative(y) ? T(r + y) : r;
}
template <class T = ll, class U, class V>
constexpr T ipow(U a, V b)
;
#line 791
template <class T = ll, class A, class B, class M>
T mul_limited(A a, B b, M m)
;
template <class T = ll, class A, class B>
T mul_limited(A a, B b) ;
template <class T = ll, class A, class B, class M>
T pow_limited(A a, B b, M m)
;
#line 826
template <class T = ll, class A, class B>
T pow_limited(A a, B b) ;
template <class T = ll, class A, class K>
constexpr T iroot(A a, K k)
;
#line 885
template <class T = ll, class A, class K>
constexpr T iroot_ceil(A a, K k)
;
//! https://misawa.github.io/others/avoid_errors/techniques_to_avoid_errors.html
template <class D = decltype(EPS), class A>
int SGN(A a, D eps = EPS) ;
template <class T = ll, class U, class V>
vc<T> base_repr(U val, V base)
;
#line 915
template <class T = ll, class U, class V>
vc<T> base_repr(U val, V base, int n)
;
#line 930
template <const bool use_upper = true, class U>
string base_repr_str(U val, int base)
;
template <const bool use_upper = true, class U>
string base_repr_str(U val, int base, int n)
;
#line 963
#define ALL(a) (a).begin(), (a).end()
template <class T = ll, class V>
inline T SZ(const V &x) { return x.size(); }
#define eb emplace_back
#define LMD(x, fx) ([&](const auto &x) { return fx; })
template <class F>
auto gen_vec(int n, const F &f)
;
#define GEN_VEC(n, i, fi) (gen_vec(n, LMD(i, fi)))
namespace template_vector_detail
{
template <class T, size_t depth>
struct nested_vector
{
using type = vc<typename nested_vector<T, depth - 1>::type>;
};
template <class T>
struct nested_vector<T, 0>
{
using type = T;
};
template <class T, size_t depth, class V>
void populate_dimensions(typename nested_vector<T, depth>::type &result, const V *sizes, const T &init)
;
#line 1008
template <class T, size_t depth, class V>
typename nested_vector<T, depth>::type build_dimensions(const V *sizes, const T &init)
;
}
template <class T, size_t d, size_t i = 0, class V>
auto dvec(const V (&sz)[d], const T &init)
;
template <class T, size_t d, size_t i = 0>
auto dvec(const ll (&sz)[d], const T &init)
;
template <class V, class T>
void fill(V &v, const T &value)
{
if constexpr (is_assignable_v<decltype(*std::begin(v)), const T &>)
std::fill(std::begin(v), std::end(v), value);
else
for (auto &row : v) fill(row, value);
}
template <class T = ll>
optional<T> ctol(const char &c, const string &s)
;
template <class T = ll>
vc<T> stov(const string &s, char first)
;
template <class T = ll>
optional<vc<T>> stov(const string &s, const string &t)
;
template <class T>
string vtos(const vc<T> &v, char first)
;
template <class T>
string vtos(const vc<T> &v, const string &t)
;
template <class T>
vc<T> concat(const vvc<T> &vs)
;
template <class T>
vc<T> concat(const vc<T> &v) ;
template <class T, class... Ts>
vc<T> concat(vc<T> v, const vc<Ts> &...vs)
;
template <class T>
vc<T> merged(const vc<T> &a, const vc<T> &b)
;
template <class T, class I>
T vecget(const vc<T> &v, I i, const T &dflt_negative = -default_infty<T>(), const T &dflt_positive = default_infty<T>())
;
#line 1134
template <class T, auto x, enable_if_t<!is_invocable_v<decltype(x)>, int> = 0>
constexpr T resolved_value()
;
template <class T, auto x, enable_if_t<is_invocable_v<decltype(x)>, long> = 0>
const T &resolved_value()
;
template <class T, auto x = nullptr>
constexpr decltype(auto) resolved_infty()
;
#line 1170
namespace internal
{
template <class M, class I, class = void>
struct HasMonoidPow : false_type {};
template <class M, class I>
struct HasMonoidPow<M, I, void_t<decltype(M::pow(declval<const typename M::S &>(), declval<I>()))>> : true_type {};
}
namespace internal_monoid_power
{
template <class I>
constexpr void check_exponent(I k)
;
template <class M, class I>
constexpr typename M::S binary_inplace(typename M::S &x, I k)
;
template <class M, class I>
constexpr typename M::S binary(typename M::S x, I k)
;
}
template <class M, class I>
constexpr typename M::S pow_monoid(typename M::S x, I k)
;
template <class G, class I>
constexpr typename G::S pow_group(typename G::S x, I k)
;
#line 1238
template <class V>
auto SUM(const V &v)
;
template <class T, class V>
T SUM(const V &v)
;
template <class V>
auto MAX(const V &v) ;
template <class V>
auto MIN(const V &v) ;
template <class I = ll, class V>
optional<I> ARGMAX(const V &v)
;
template <class I = ll, class V>
optional<I> ARGMIN(const V &v)
;
template<class T = ll, class V>
T mex(const V &a)
;
template <class I>
bool is_permutation(const vc<I> &p)
;
#line 1295
template <class T = ll>
vc<T> permid(const int &n, const int &base_index = 0)
;
template <class T>
vc<T> perminv(const vc<T> &p)
;
template <class T, class U>
vc<T> permuted(const vc<T> &a, const vc<U> &p)
;
template <class T, class U, class... Ts>
vc<T> permuted(const vc<T> &p, const vc<U> &q, const vc<Ts> &...rs)
;
template <class V>
V reversed(const V &v) ;
template <class V, class... Args>
V sorted(V v, Args&&... args)
;
template <class V, class Equal = equal_to<>>
void unique(V &v, Equal equal = {}) { v.erase(std::unique(ALL(v), equal), v.end()); }
template <class V, class Equal = equal_to<>>
V uniqued(V v, Equal equal = {}) ;
template <class V, class Compare = less<>, class Equal = equal_to<>>
void sortunique(V &v, Compare comp = {}, Equal equal = {})
;
template <class V, class Compare = less<>, class Equal = equal_to<>>
V sortuniqued(V v, Compare comp = {}, Equal equal = {})
;
template <class V, class U>
void rotate(V &v, U k)
;
template <class V, class U>
V rotated(V v, U k) ;
template <class T>
vvc<T> top(const vvc<T> &a)
{
if (a.empty())
return {};
const int n = a.size(), m = a[0].size();
vvc<T> b(m, vc<T>(n));
repi(i, n)
{
assert(SZ<int>(a[i]) == m);
repi(j, m) b[j][i] = a[i][j];
}
return b;
}
vstr top(const vstr &a)
{
vvc<char> a_(a.size());
repi(i, SZ<int>(a)) a_[i] = {ALL(a[i])};
vvc<char> b_ = top(a_);
vstr b(b_.size());
repi(i, SZ<int>(b)) b[i] = {ALL(b_[i])};
return b;
}
template <class T, class = void>
struct has_e0 : false_type {};
template <class T>
struct has_e0<T, void_t<decltype(T::e0())>> : true_type {};
template <class T>
inline constexpr bool has_e0_v = has_e0<T>::value;
template <class T>
struct MonoidAdd
{
using S = T;
static constexpr S op(S a, S b) ;
static constexpr S e()
;
template <class I, class = decltype(declval<S>() * declval<I>())>
static constexpr S pow(const S &a, I k) { return a * k; }
};
template <class T, auto infty = nullptr>
struct MonoidMin
{
using S = T;
static constexpr S op(S a, S b) ;
static constexpr decltype(auto) e() ;
template <class I>
static constexpr S pow(const S &a, I k) { return k == 0 ? e() : a; }
};
template <class T, auto infty = nullptr>
struct MonoidMax
{
using S = T;
static constexpr S op(S a, S b) ;
static constexpr S e() ;
template <class I>
static constexpr S pow(const S &a, I k) { return k == 0 ? e() : a; }
};
template <class M>
vc<typename M::S> cuml(const vc<typename M::S> &v, int left_index = 0)
;
template <class M>
vc<typename M::S> cumr(const vc<typename M::S> &v, int right_index = 0)
;
template <class T>
vc<T> cumlsum(const vc<T> &v, int left_index = 0)
;
template <class T>
vc<T> cumrsum(const vc<T> &v, int right_index = 0)
;
template <class T>
vc<T> cumlmin(const vc<T> &v, int left_index = 0)
;
template <class T>
vc<T> cumrmin(const vc<T> &v, int right_index = 0)
;
template <class T>
vc<T> cumlmax(const vc<T> &v, int left_index = 0)
;
template <class T>
vc<T> cumrmax(const vc<T> &v, int right_index = 0)
;
template <class T>
vc<T> adjd(const vc<T> &v, int left_index = 0, int right_index = 0)
;
#line 1504
constexpr array<pll, 4> DRULgrid = {{{1, 0}, {0, 1}, {-1, 0}, {0, -1}}};
constexpr array<pll, 4> DRULplane = {{{0, -1}, {1, 0}, {0, 1}, {-1, 0}}};
template <class T>
struct is_random_access_iterator
{
static constexpr bool value = is_same_v<
typename iterator_traits<T>::iterator_category,
random_access_iterator_tag
>;
};
template <class T>
constexpr bool is_random_access_iterator_v = is_random_access_iterator<T>::value;
#line 1569
template <class T = ll, class V, class Value, class Comp = ranges::less, class Proj = identity>
inline T LB(const V &v, const Value &val, Comp comp = {}, Proj proj = {})
;
template <class T = ll, class V, class Value, class Comp = ranges::less, class Proj = identity>
inline T UB(const V &v, const Value &val, Comp comp = {}, Proj proj = {})
;
#define DEFAULT_COMP ranges::less
template <class T = ll, class V, class Value, class Comp = DEFAULT_COMP, class Proj = identity>
inline auto lt_max(const V &v, const Value &val, Comp comp = {}, Proj proj = {})
-> enable_if_t<is_random_access_iterator_v<typename V::iterator>, T>
;
template <class T = ll, class V, class Value, class Comp = DEFAULT_COMP, class Proj = identity>
inline auto leq_max(const V &v, const Value &val, Comp comp = {}, Proj proj = {})
-> enable_if_t<is_random_access_iterator_v<typename V::iterator>, T>
;
template <class T = ll, class V, class Value, class Comp = DEFAULT_COMP, class Proj = identity>
inline auto gt_min(const V &v, const Value &val, Comp comp = {}, Proj proj = {})
-> enable_if_t<is_random_access_iterator_v<typename V::iterator>, T>
;
template <class T = ll, class V, class Value, class Comp = DEFAULT_COMP, class Proj = identity>
inline auto geq_min(const V &v, const Value &val, Comp comp = {}, Proj proj = {})
-> enable_if_t<is_random_access_iterator_v<typename V::iterator>, T>
;
template <class T = ll, class V, class Value, class Comp = DEFAULT_COMP, class Proj = identity>
inline auto lt_cnt(const V &v, const Value &val, Comp comp = {}, Proj proj = {})
-> enable_if_t<is_random_access_iterator_v<typename V::iterator>, T>
;
template <class T = ll, class V, class Value, class Comp = DEFAULT_COMP, class Proj = identity>
inline auto leq_cnt(const V &v, const Value &val, Comp comp = {}, Proj proj = {})
-> enable_if_t<is_random_access_iterator_v<typename V::iterator>, T>
;
template <class T = ll, class V, class Value, class Comp = DEFAULT_COMP, class Proj = identity>
inline auto gt_cnt(const V &v, const Value &val, Comp comp = {}, Proj proj = {})
-> enable_if_t<is_random_access_iterator_v<typename V::iterator>, T>
;
template <class T = ll, class V, class Value, class Comp = DEFAULT_COMP, class Proj = identity>
inline auto geq_cnt(const V &v, const Value &val, Comp comp = {}, Proj proj = {})
-> enable_if_t<is_random_access_iterator_v<typename V::iterator>, T>
;
template <class T = ll, class V, class L, class R, class Comp = DEFAULT_COMP, class Proj = identity>
inline auto in_cnt(const V &v, L l, R r, Comp comp = {}, Proj proj = {})
-> enable_if_t<is_random_access_iterator_v<typename V::iterator>, T>
;
template <class V, class Value>
inline auto lt_max(const V &v, const Value &val)
-> enable_if_t<!is_random_access_iterator_v<typename V::iterator>, optional<typename V::const_iterator>>
;
template <class V, class Value>
inline auto leq_max(const V &v, const Value &val)
-> enable_if_t<!is_random_access_iterator_v<typename V::iterator>, optional<typename V::const_iterator>>
;
template <class V, class Value>
inline auto gt_min(const V &v, const Value &val)
-> enable_if_t<!is_random_access_iterator_v<typename V::iterator>, optional<typename V::const_iterator>>
;
template <class V, class Value>
inline auto geq_min(const V &v, const Value &val)
-> enable_if_t<!is_random_access_iterator_v<typename V::iterator>, optional<typename V::const_iterator>>
;
namespace internal
{
template <class T>
bool binsearch_adjacent(T a, T b)
;
};
template <class T = ll, class Judge, class InitOk, class InitNg>
pair<T, T> binsearch(const Judge &judge, InitOk init_ok, InitNg init_ng, bool check_ok = true, bool check_ng = true)
;
#line 1707
template <class T = ld, class Judge, class InitOk, class InitNg>
T binsearch_real(const Judge &judge, InitOk init_ok, InitNg init_ng, int iteration_count = 100, bool check_ok = true, bool check_ng = true)
;
#line 1723
template <class T = ll, class Judge, class InitVal>
pair<T, T> expsearch(const Judge &judge, InitVal init_val, bool positive = true)
;
#line 1769
template <class T>
inline constexpr ull pow2(T k) ;
template <class T>
inline constexpr ull MASK(T k) { return (1ULL << k) - 1ULL; }
#line 1787
inline constexpr ll bit_width(ll x) { return std::bit_width((ull)x); }
inline constexpr ll bit_floor(ll x) { return std::bit_floor((ull)x); }
inline constexpr ll bit_ceil(ll x) { return std::bit_ceil((ull)x); }
inline constexpr ll countr_zero(ll x) { assert(x != 0); return std::countr_zero((ull)x); }
inline constexpr ll popcount(ll x) { return std::popcount((ull)x); }
inline constexpr bool has_single_bit(ll x) { return std::has_single_bit((ull)x); }
inline constexpr ull lsb_pos(ull x) { assert(x != 0); return countr_zero(x); }
inline constexpr ull msb_pos(ull x) { assert(x != 0); return bit_width(x) - 1; }
inline constexpr ull lsb_mask(ull x) { assert(x != 0); return x & -x; }
inline constexpr ull msb_mask(ull x) { assert(x != 0); return bit_floor(x); }
inline constexpr bool btest(ull x, uint k) { return (x >> k) & 1; }
template <class T>
inline void bset(T &x, uint k, bool b = 1) ;
template <class T>
inline void bflip(T &x, uint k) ;
inline constexpr bool bsubset(ull x, ull y) { return (x & y) == x; }
inline constexpr bool bsupset(ull x, ull y) { return (x & y) == y; }
inline constexpr ull bsetminus(ull x, ull y) { return x & ~y; }
#line 1912
#define CPP_DUMP_DEFINE_DATA(...)
#define dump(...) ((void)0)
#define local(...)
#define oj(...) __VA_ARGS__
#define local_oj(a, b) (b)
template <class T, class Sequence>
vc<T> content(queue<T, Sequence> que)
{
vc<T> res;
while (!que.empty())
{
res.eb(que.front());
que.pop();
}
return res;
}
template <class T, class Sequence, class Compare>
vc<T> content(priority_queue<T, Sequence, Compare> pque)
{
vc<T> res;
while (!pque.empty())
{
res.eb(pque.top());
pque.pop();
}
return res;
}
template <class T>
auto content(const T &obj) { return obj.content(); }
#line 1977
namespace fastio
{
constexpr uint32_t SIZ = 1 << 17;
inline char ibuf[SIZ + 1], obuf[SIZ];
inline uint32_t pil = 0, pir = 0, por = 0;
__attribute__((noinline)) inline void load()
{
const uint32_t remaining = pir - pil;
if (remaining && pil) memmove(ibuf, ibuf + pil, remaining);
pil = 0;
pir = remaining + uint32_t(fread(ibuf + remaining, 1, SIZ - remaining, stdin));
ibuf[pir] = '\x7f';
}
inline bool ascii_space(unsigned char c)
{
return c == ' ' || unsigned(c - '\t') < 5;
}
inline void find_token()
{
for (;;)
{
const char *p = ibuf + pil;
const char *end = ibuf + pir;
while (ascii_space(*p)) ++p;
pil = uint32_t(p - ibuf);
if (p != end) return;
load();
if (!pir) abort();
}
}
__attribute__((noinline)) inline void flush()
{
size_t written = 0;
while (written < por)
{
const size_t n = fwrite(obuf + written, 1, por - written, stdout);
if (!n) abort();
written += n;
}
por = 0;
}
inline void append_bytes(const char *data, size_t size)
{
while (size)
{
const size_t n = min(size, size_t(SIZ - por));
memcpy(obuf + por, data, n);
por += uint32_t(n);
data += n;
size -= n;
if (por == SIZ) flush();
}
}
inline void append_char(char c)
{
if (por == SIZ) flush();
obuf[por++] = c;
}
template <class T>
using integer_magnitude = conditional_t<(sizeof(T) > sizeof(ull)), u128,
conditional_t<(sizeof(T) > sizeof(uint)), ull, uint>>;
inline uint decode_eight(uint64_t columns)
{
columns &= 0x0f0f0f0f0f0f0f0fULL;
const uint64_t pairs = (columns & 0x00ff00ff00ff00ffULL) * 10
+ ((columns >> 8) & 0x00ff00ff00ff00ffULL);
const uint64_t quads = (pairs & 0x0000ffff0000ffffULL) * 100
+ ((pairs >> 16) & 0x0000ffff0000ffffULL);
return uint(quads) * 10000 + uint(quads >> 32);
}
inline constexpr uint powers_of_ten[8] = {1, 10, 100, 1000, 10000, 100000, 1000000, 10000000};
template <class U>
__attribute__((always_inline)) inline void read_wide_integer_tail(
U &value, const char *&p, uint64_t columns, uint64_t stop)
;
#line 2078
template <class T>
__attribute__((always_inline)) inline void rd1_integer(T &x)
{
using U = integer_magnitude<T>;
const char *p = ibuf + pil;
for (;;)
{
while (ascii_space(*p)) ++p;
if (p != ibuf + pir) break;
pil = pir;
load();
if (!pir) abort();
p = ibuf;
}
bool negative = false;
if (*p == '+') ++p;
else if constexpr (numeric_limits<T>::is_signed)
{
if (*p == '-') negative = true, ++p;
}
U value = 0;
for (;;)
{
while (size_t(ibuf + pir - p) >= 8)
{
uint64_t columns;
memcpy(&columns, p, 8);
const uint64_t stop = ((columns + 0x4646464646464646ULL) |
(columns - 0x3030303030303030ULL)) & 0x8080808080808080ULL;
if (stop)
{
if constexpr (sizeof(U) > sizeof(ull))
{
read_wide_integer_tail(value, p, columns, stop);
goto integer_complete;
}
else
{
if (stop & 0x80808080ULL)
{
if (!(stop & 0x808080ULL))
{
value = value * 1000 + uint(columns & 15) * 100
+ uint(columns >> 8 & 15) * 10 + uint(columns >> 16 & 15);
p += 3;
}
else if (!(stop & 0x8080ULL))
{
value = value * 100 + uint(columns & 15) * 10 + uint(columns >> 8 & 15);
p += 2;
}
else if (!(stop & 0x80ULL))
{
value = value * 10 + uint(columns & 15);
++p;
}
goto integer_complete;
}
break;
}
}
columns &= 0x0f0f0f0f0f0f0f0fULL;
const uint64_t pairs = (columns & 0x00ff00ff00ff00ffULL) * 10
+ ((columns >> 8) & 0x00ff00ff00ff00ffULL);
const uint64_t quads = (pairs & 0x0000ffff0000ffffULL) * 100
+ ((pairs >> 16) & 0x0000ffff0000ffffULL);
value = value * 100000000 + uint(quads) * 10000 + uint(quads >> 32);
p += 8;
}
while (unsigned(*p - '0') < 10)
{
value = value * 10 + unsigned(*p - '0');
++p;
}
pil = uint32_t(p - ibuf);
if (pil != pir) break;
load();
p = ibuf;
if (!pir) break;
}
integer_complete:
pil = uint32_t(p - ibuf);
if constexpr (numeric_limits<T>::is_signed)
{
if (negative && value) x = T(-T(value - 1) - 1);
else x = T(value);
}
else x = T(value);
}
template <class T, enable_if_t<numeric_limits<T>::is_integer, int> = 0>
__attribute__((always_inline)) inline void rd1(T &x) { rd1_integer(x); }
inline void rd1(char &x)
;
inline void rd1(string &x)
;
#line 2197
template <class T>
void rd1_real(T &x)
{
string token;
rd1(token);
const char *first = token.data();
const char *last = first + token.size();
if (first != last && *first == '+') ++first;
auto result = from_chars(first, last, x, chars_format::general);
if (result.ec != errc{} || result.ptr != last)
{
char *end;
if constexpr (is_same_v<T, float>) x = strtof(token.c_str(), &end);
else if constexpr (is_same_v<T, double>) x = strtod(token.c_str(), &end);
else x = strtold(token.c_str(), &end);
if (end != token.data() + token.size()) abort();
}
}
inline void rd1(float &x) ;
inline void rd1(double &x) ;
inline void rd1(long double &x) ;
template <class A, class B> void rd1(pair<A, B> &x);
template <class... T> void rd1(tuple<T...> &x);
template <class T, size_t N> void rd1(array<T, N> &x);
template <class T> void rd1(vc<T> &x);
template <class A, class B>
void rd1(pair<A, B> &x) ;
template <class... T>
void rd1(tuple<T...> &x)
;
template <class T, size_t N>
void rd1(array<T, N> &x)
;
template <class T>
void rd1(vc<T> &x)
;
inline constexpr auto decimal_triplets = []
{
array<array<char, 3>, 1000> result{};
for (unsigned i = 0; i < 1000; ++i)
{
result[i][0] = char('0' + i / 100);
result[i][1] = char('0' + i / 10 % 10);
result[i][2] = char('0' + i % 10);
}
return result;
}();
inline void store_nine(char *out, uint x)
{
const uint q = x / 1000, r = q / 1000;
memcpy(out, decimal_triplets[r].data(), 3);
memcpy(out + 3, decimal_triplets[q - r * 1000].data(), 3);
memcpy(out + 6, decimal_triplets[x - q * 1000].data(), 3);
}
inline char *prepend_small(char *end, uint x)
{
while (x >= 1000)
{
const uint q = x / 1000;
end -= 3;
memcpy(end, decimal_triplets[x - q * 1000].data(), 3);
x = q;
}
const unsigned digits = 1 + unsigned(x >= 10) + unsigned(x >= 100);
end -= digits;
memcpy(end, decimal_triplets[x].data() + 3 - digits, digits);
return end;
}
inline char *prepend_wide(char *end, ull x)
{
while (x >= 1000000000ULL)
{
const ull q = x / 1000000000ULL;
end -= 9;
store_nine(end, uint(x - q * 1000000000ULL));
x = q;
}
return prepend_small(end, uint(x));
}
inline char *store_leading_triplet(char *out, uint x)
{
if (x >= 100)
{
memcpy(out, decimal_triplets[x].data(), 3);
return out + 3;
}
if (x >= 10)
{
memcpy(out, decimal_triplets[x].data() + 1, 2);
return out + 2;
}
*out = char('0' + x);
return out + 1;
}
inline char *store_short(char *out, uint x)
{
if (x >= 1000000)
{
const uint q = x / 1000000;
const uint r = x - q * 1000000;
out = store_leading_triplet(out, q);
memcpy(out, decimal_triplets[r / 1000].data(), 3);
memcpy(out + 3, decimal_triplets[r % 1000].data(), 3);
return out + 6;
}
if (x >= 1000)
{
const uint q = x / 1000;
out = store_leading_triplet(out, q);
memcpy(out, decimal_triplets[x - q * 1000].data(), 3);
return out + 3;
}
return store_leading_triplet(out, x);
}
template <class U>
inline char *store_ordinary(char *out, U x)
{
if constexpr (sizeof(U) > sizeof(uint))
{
if (x >= 1000000000000000000ULL)
{
const uint hi = uint(x / 1000000000000000000ULL);
const ull tail = x - ull(hi) * 1000000000000000000ULL;
const uint middle = uint(tail / 1000000000ULL);
out = store_leading_triplet(out, hi);
store_nine(out, middle);
store_nine(out + 9, uint(tail - ull(middle) * 1000000000ULL));
return out + 18;
}
}
if (x >= 1000000000)
{
const uint hi = uint(x / 1000000000);
out = store_short(out, hi);
store_nine(out, uint(x - U(hi) * 1000000000));
return out + 9;
}
return store_short(out, uint(x));
}
template <class U, bool newline>
__attribute__((always_inline)) inline void write_wide_integer(U magnitude, bool negative)
;
#line 2386
template <class T, bool newline = false>
__attribute__((always_inline)) inline void wt1_integer(T x)
{
using U = integer_magnitude<T>;
U magnitude = U(x);
bool negative = false;
if constexpr (numeric_limits<T>::is_signed)
{
negative = x < 0;
const U sign_mask = U(0) - U(negative);
magnitude = (magnitude ^ sign_mask) + U(negative);
}
if constexpr (sizeof(U) <= sizeof(ull))
{
if (por > SIZ - 40) flush();
char *out = obuf + por;
*out = '-';
out += negative;
out = store_ordinary(out, magnitude);
if constexpr (newline) *out++ = '\n';
por = uint32_t(out - obuf);
return;
}
else write_wide_integer<U, newline>(magnitude, negative);
}
template <class T, enable_if_t<numeric_limits<T>::is_integer, int> = 0>
__attribute__((always_inline)) inline void wt1(T x) { wt1_integer(x); }
inline void wt1(char x) ;
inline void wt1(const char *x) ;
inline void wt1(string_view x) ;
inline void wt1(const string &x) ;
template <class T>
void wt1_real(T x)
{
char text[numeric_limits<T>::max_exponent10 + 64];
auto result = to_chars(text, text + sizeof(text), x, chars_format::fixed, 15);
if (result.ec != errc{}) abort();
append_bytes(text, size_t(result.ptr - text));
}
inline void wt1(float x) ;
inline void wt1(double x) ;
inline void wt1(long double x) ;
template <class A, class B> void wt1(const pair<A, B> &x);
template <class... T> void wt1(const tuple<T...> &x);
template <class T, size_t N> void wt1(const array<T, N> &x);
template <class T> void wt1(const vc<T> &x);
template <class A, class B>
void wt1(const pair<A, B> &x) ;
template <class... T>
void wt1(const tuple<T...> &x)
;
template <class T, size_t N>
void wt1(const array<T, N> &x)
;
template <class T>
void wt1(const vc<T> &x)
;
template <class... T>
void read(T &...x) { (rd1(x), ...); }
template <class... T>
void write(T &&...x) { (wt1(forward<T>(x)), ...); }
template <class... T>
void print(T &&...x)
{
if constexpr (sizeof...(T) == 1 &&
((is_integral_v<decay_t<T>> && !is_same_v<decay_t<T>, char>) && ...))
(wt1_integer<decay_t<T>, true>(forward<T>(x)), ...);
else
{
size_t index = 0;
((index++ ? void(append_char(' ')) : void(), wt1(forward<T>(x))), ...);
append_char('\n');
}
}
}
struct Dummy {
Dummy() { atexit(fastio::flush); }
} dummy;
//! https://trap.jp/post/1224/
#line 2532
namespace internal
{
template <class... Ts>
void READnodump(Ts &...a) { fastio::read(a...); }
template <class... T>
void READVECnodump(int n, vc<T> &...v)
;
template <class... T>
void READVEC2nodump(int n, int m, vvc<T> &...v)
;
template <class... T>
void READJAGnodump(int n, vvc<T> &...vs)
;
#line 2579
};
#define READ(...) internal::READnodump(__VA_ARGS__); dump(__VA_ARGS__)
#define LIBRARY_INOUT_TYPE(...) LIBRARY_INOUT_TYPE_IMPL(__VA_ARGS__)(__VA_ARGS__)
#define LIBRARY_INOUT_TYPE_IMPL(T, ...) LIBRARY_INOUT_TYPE_SELECT(LIBRARY_INOUT_TYPE_PROBE T, LIBRARY_INOUT_TYPE_KEEP)
#define LIBRARY_INOUT_TYPE_PROBE(...) _, LIBRARY_INOUT_TYPE_REMOVE
#define LIBRARY_INOUT_TYPE_SELECT(...) LIBRARY_INOUT_TYPE_SECOND(__VA_ARGS__)
#define LIBRARY_INOUT_TYPE_SECOND(_0, _1, ...) _1
#define LIBRARY_INOUT_TYPE_KEEP(...) __VA_ARGS__
#define LIBRARY_INOUT_TYPE_REMOVE(T) LIBRARY_INOUT_TYPE_KEEP T
#define IN(T, ...) LIBRARY_INOUT_TYPE(T) __VA_ARGS__; READ(__VA_ARGS__)
#define CHAR(...) IN(char, __VA_ARGS__)
#define INT(...) IN(int, __VA_ARGS__)
#define LL(...) IN(ll, __VA_ARGS__)
#define STR(...) IN(string, __VA_ARGS__)
#define ARR(T, n, ...) array<LIBRARY_INOUT_TYPE(T), n> __VA_ARGS__; READ(__VA_ARGS__)
#define READVEC(...) internal::READVECnodump(__VA_ARGS__); dump(__VA_ARGS__)
#define READVEC2(...) internal::READVEC2nodump(__VA_ARGS__); dump(__VA_ARGS__)
#define VEC(T, n, ...) vc<LIBRARY_INOUT_TYPE(T)> __VA_ARGS__; READVEC(n, __VA_ARGS__)
#define VEC2(T, n, m, ...) vvc<LIBRARY_INOUT_TYPE(T)> __VA_ARGS__; READVEC2(n, m, __VA_ARGS__)
#define READJAG(...) internal::READJAGnodump(__VA_ARGS__); dump(__VA_ARGS__)
#define JAG(T, n, ...) vvc<LIBRARY_INOUT_TYPE(T)> __VA_ARGS__; READJAG(n, __VA_ARGS__)
#define ENDL '\n'
#line 2688
#define WRITE fastio::write
#define PRINT fastio::print
#define PRINTEXIT(...) do { PRINT(__VA_ARGS__); exit(0); } while (false)
#define PRINTRETURN(...) do { PRINT(__VA_ARGS__); return; } while (false)
template <class T>
void PRINTV(const vc<T> &v) ;
#define PRINTVEXIT(...) do { PRINTV(__VA_ARGS__); exit(0); } while (false)
#define PRINTVRETURN(...) do { PRINTV(__VA_ARGS__); return; } while (false)
template <class T, class Add>
void offset(vc<T> &v, const Add &add) ;
template <class T, class Add>
void offset(vvc<T> &v, const Add &add) ;
template <class T, const size_t m>
array<vc<T>, m> unzip(const vc<array<T, m>> &vt)
;
template <class T, const size_t m>
vc<array<T, m>> zip(const array<vc<T>, m> &tv)
;
#line 2737
template <class T, class U>
pair<vc<T>, vc<U>> unzip(const vc<pair<T, U>> &vt)
;
template <class T, class U>
vc<pair<T, U>> zip(const pair<vc<T>, vc<U>> &tv)
;
namespace internal
{
template <size_t... I, class V, class Tp>
auto vt_to_tv_impl(V &tv, const Tp &t, index_sequence<I...>, size_t index)
;
template <size_t... I, class Tp>
auto tv_to_vt_impl(const Tp &tv, index_sequence<I...>, size_t index)
;
};
template <class... Ts>
auto unzip(const vc<tuple<Ts...>> &vt)
;
template <class... Ts>
auto zip(const tuple<vc<Ts>...> &tv)
;
namespace internal
{
template <class... Ts>
auto zip_vectors(const vc<Ts> &...v)
;
};
#define UNZIP(vt, ...) auto [__VA_ARGS__] = unzip(vt)
#define ZIP(vt, ...) auto vt = internal::zip_vectors(__VA_ARGS__)
mt19937_64 mt;
template <class T = ll, class U1, class U2>
T randint(U1 l, U2 r)
;
template <class T = ll, class U1, class U2>
T randrange(U1 l, U2 r)
;
template <class T = double, class U1, class U2>
T randreal(U1 l, U2 r)
;
bool randbool(double p)
{
assert(0 <= p && p <= 1);
return bernoulli_distribution(p)(mt);
}
namespace internal
{
template <bool does_sort, class V, class T>
void random_sample_range(V &res, T l, T r)
;
#line 2888
};
template <int k, bool does_sort, class T = ll, class U1, class U2>
array<T, k> random_sample_range_array(U1 l, U2 r)
;
template <bool does_sort, class T = ll, class U1, class U2>
vc<T> random_sample_range_vector(U1 l, U2 r, int k)
;
#line 2917
/*!
* Portions corresponding to AtCoder Library (CC0-1.0): modular arithmetic operators.
* Upstream: atcoder/modint.hpp
* https://github.com/atcoder/ac-library/tree/864245a00b00dd008d1abfdc239618fdb7d139da
* License: https://github.com/atcoder/ac-library/blob/864245a00b00dd008d1abfdc239618fdb7d139da/LICENSE
*/
#line 2936
namespace larger_int_detail
{
template <class T>
[[deprecated("larger_int: 128-bit integer is not widened; intermediate arithmetic may overflow")]]
constexpr bool warn_no_wider_integer()
{
return true;
}
}
template <class T>
struct larger_int
{
private:
static constexpr bool check()
;
static_assert(check());
public:
using type = T;
};
#define LARGER_INT(T, U) \
template <> \
struct larger_int<T> \
{ \
using type = U; \
};
LARGER_INT(signed char, short)
LARGER_INT(short, int)
LARGER_INT(int, long long)
LARGER_INT(long, __int128_t)
LARGER_INT(long long, __int128_t)
LARGER_INT(unsigned char, unsigned short)
LARGER_INT(unsigned short, unsigned int)
LARGER_INT(unsigned int, unsigned long long)
LARGER_INT(unsigned long, __uint128_t)
LARGER_INT(unsigned long long, __uint128_t)
#undef LARGER_INT
template <class T>
struct Rational;
template <class T>
struct larger_int<Rational<T>>
{
using type = Rational<typename larger_int<T>::type>;
};
template <class T>
using larger_int_t = typename larger_int<T>::type;
/*!
* Portions corresponding to AtCoder Library (CC0-1.0): primality testing.
* Upstream: atcoder/internal_math.hpp
* https://github.com/atcoder/ac-library/tree/864245a00b00dd008d1abfdc239618fdb7d139da
* License: https://github.com/atcoder/ac-library/blob/864245a00b00dd008d1abfdc239618fdb7d139da/LICENSE
*/
namespace internal
{
template <class T>
constexpr ll powmod_constexpr(ll x, ll n, T m)
{
if (m == 1)
return 0;
using U = make_unsigned_t<T>;
using L = larger_int_t<U>;
U r = 1, y = safemod(x, m);
while (n)
{
if (n & 1)
r = L(r) * y % m;
y = L(y) * y % m;
n >>= 1;
}
return r;
}
template <class T>
constexpr bool isprime_constexpr(T n)
{
if constexpr (sizeof(T) > 4)
{
if (n <= INT_MAX)
return isprime_constexpr<int>(n);
}
if (n <= 1)
return false;
if (n == 2 || n == 7 || n == 61)
return true;
if (n % 2 == 0)
return false;
ll d = n - 1;
while (d % 2 == 0)
d /= 2;
using U = make_unsigned_t<T>;
using L = larger_int_t<U>;
auto miller_rabin = [&](const auto &bases) constexpr
{
for (ll a : bases)
{
ll t = d, y = powmod_constexpr(a, t, n);
while (t != n - 1 && y != 1 && y != n - 1)
{
y = L(y) * y % n;
t <<= 1;
}
if (y != n - 1 && t % 2 == 0)
return false;
}
return true;
};
if constexpr(sizeof(T) <= 4)
{
constexpr ll bases[3] = {2, 7, 61};
return miller_rabin(bases);
}
else
{
constexpr ll bases[7] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
return miller_rabin(bases);
}
}
template <auto n>
constexpr bool isprime = isprime_constexpr(n);
};
namespace internal
{
template <auto M>
struct policy_static
{
using mod_type = decltype(M);
using value_type = make_unsigned_t<mod_type>;
using calc_type = conditional_t<(M <= UINT_MAX), ull, larger_int_t<value_type>>;
static constexpr bool is_prime = isprime_constexpr(M);
static constexpr mod_type mod() { return M; }
static constexpr value_type umod() { return M; }
static constexpr value_type init(value_type v) { return v; }
static constexpr mod_type val(value_type v) { return v; }
static constexpr value_type mul(value_type a, value_type b)
{
return (value_type)((calc_type(a) * b) % M);
}
};
};
/*!
* Portions corresponding to AtCoder Library (CC0-1.0): Barrett reduction.
* Upstream: atcoder/internal_math.hpp
* https://github.com/atcoder/ac-library/tree/864245a00b00dd008d1abfdc239618fdb7d139da
* License: https://github.com/atcoder/ac-library/blob/864245a00b00dd008d1abfdc239618fdb7d139da/LICENSE
*/
namespace internal
{
struct barrett32
{
uint m = 0;
ull im = 0;
barrett32() = default;
explicit barrett32(uint m) : m(m), im((ull)(-1) / m + 1) {}
uint umod() const { return m; }
uint mul(uint a, uint b) const
{
ull z = a;
z *= b;
ull x = ull((u128(z) * im) >> 64);
ull y = x * m;
return uint(z - y + (z < y ? m : 0));
}
};
template <int id>
struct policy_barrett32
{
using value_type = uint;
using calc_type = ull;
using mod_type = int;
static constexpr bool is_prime = false;
static inline barrett32 reducer{};
static void set_mod(mod_type m)
{
assert(m >= 1 && "modint: set_mod(m) requires m >= 1");
reducer = barrett32(m);
}
static mod_type mod() { return reducer.umod(); }
static value_type umod() ;
static value_type init(value_type v) { return v; }
static mod_type val(value_type v) ;
static value_type mul(value_type a, value_type b) { return reducer.mul(a, b); }
};
};
namespace internal
{
inline constexpr ull inv64(ull a)
{
ull x = a;
while (a * x != 1) x *= 2 - a * x;
return x;
}
struct montgomery64odd
{
ull m = 0, im = 0, sq = 0;
montgomery64odd() = default;
explicit montgomery64odd(ull m) : m(m), im(inv64(m)), sq(-u128(m) % m) {}
ull umod() const { return m; }
ull reduce(u128 x) const
{
auto t = (x + u128(m) * (-im * ull(x))) >> 64;
if (t >= m) t -= m;
return (ull)t;
}
ull inv_reduce(i128 v) const { return reduce(u128(v % m + m) * sq); }
};
//! https://www.mathenachia.blog/even-mod-montgomery-impl/
struct montgomery64
{
ull m = 0, mx = 0, imx = 0, d = 0, q = 0;
uint b = 0;
montgomery64() = default;
explicit montgomery64(ull m) : m(m)
{
b = countr_zero(m), mx = m >> b;
imx = inv64(mx);
d = powmod_constexpr((mx + 1) / 2, b, mx);
u128 sq = -u128(mx) % mx;
q = (1 + (((sq - 1) * d) << b)) % m;
}
ull umod() const ;
ull reduce(u128 x) const
;
#line 3222
ull inv_reduce(i128 v) const ;
};
template <int id>
struct policy_montgomery64_odd
{
using value_type = ull;
using calc_type = u128;
using mod_type = ll;
static constexpr bool is_prime = false;
static inline montgomery64odd reducer{};
static void set_mod(mod_type m)
{
assert(m >= 1 && m % 2 == 1 && "modint: set_mod(m) requires a positive odd modulus");
reducer = montgomery64odd(m);
}
static mod_type mod() { return reducer.umod(); }
static value_type umod() ;
static value_type init(value_type v) { return reducer.inv_reduce(v); }
static mod_type val(value_type v) ;
static value_type mul(value_type a, value_type b) { return reducer.reduce((calc_type)a * b); }
};
template <int id>
struct policy_montgomery64
{
using value_type = ull;
using calc_type = u128;
using mod_type = ll;
static constexpr bool is_prime = false;
static inline montgomery64 reducer{};
static void set_mod(mod_type m)
;
static mod_type mod() ;
static value_type umod() ;
static value_type init(value_type v) { return reducer.inv_reduce(v); }
static mod_type val(value_type v) ;
static value_type mul(value_type a, value_type b) ;
};
};
template <class T = ll>
constexpr tuple<T, T, T> extgcd(T a, T b)
{
if (a == 0 && b == 0)
return {0, 0, 0};
T x1 = 1, y1 = 0, z1 = a;
T x2 = 0, y2 = 1, z2 = b;
while (z2 != 0)
{
T q = z1 / z2;
tie(x1, x2) = make_pair(x2, x1 - q * x2);
tie(y1, y2) = make_pair(y2, y1 - q * y2);
tie(z1, z2) = make_pair(z2, z1 - q * z2);
}
if (z1 < 0)
x1 = -x1, y1 = -y1, z1 = -z1;
return {z1, x1, y1};
}
namespace internal
{
template <class Policy, class PowInv = void>
struct modint_impl
{
using V = typename Policy::value_type;
using M = typename Policy::mod_type;
using mint = modint_impl;
private:
V _v;
public:
static constexpr M mod()
{
M m = Policy::mod();
assert(m > 0 && "modint: modulus is not set; call mint::set_mod(m) before use");
return m;
}
template <class T = Policy>
static auto set_mod(M m) -> decltype(T::set_mod(m)) { return T::set_mod(m); }
static void precompute()
;
static mint raw(V v)
{
mint x;
x._v = v;
return x;
}
modint_impl() : _v(0) { (void)mod(); }
template <class T, typename = enable_if_t<is_integral_ext<T>>>
modint_impl(T v)
{
const V m = mod();
V rem;
if constexpr (is_signed_ext<T>)
{
using S = make_signed_t<V>;
S x = v % S(m);
if (x < 0)
x += m;
rem = x;
}
else
rem = V(v % m);
_v = Policy::init(rem);
};
M val() const { return Policy::val(_v); }
mint &operator+=(const mint &rhs)
{
_v += rhs._v;
if (_v >= Policy::umod())
_v -= Policy::umod();
return *this;
}
mint &operator-=(const mint &rhs)
{
_v -= rhs._v;
if (_v >= Policy::umod())
_v += Policy::umod();
return *this;
}
mint &operator*=(const mint &rhs)
{
_v = Policy::mul(_v, rhs._v);
return *this;
}
mint &operator/=(const mint &rhs)
{
return *this *= rhs.inv();
}
mint &operator++()
;
mint &operator--()
;
mint operator++(int)
;
mint operator--(int)
;
mint operator+() const ;
mint operator-() const { return mint() - *this; }
template <class T>
mint pow(T n) const
{
assert(n >= 0);
if constexpr (!is_void_v<PowInv>)
return raw(PowInv::get().pow(val(), n));
mint x = *this, r = 1;
while (n)
{
if (n & 1)
r *= x;
x *= x;
n >>= 1;
}
return r;
}
mint inv() const
{
if constexpr (!is_void_v<PowInv>)
{
return raw(PowInv::get().inv(val()));
}
else if constexpr (Policy::is_prime)
{
return pow(mod() - 2);
}
else
{
auto [g, x, y] = extgcd<M>(val(), mod());
assert(g == 1);
return mint(x);
}
}
friend mint operator+(const mint &lhs, const mint &rhs) { return mint(lhs) += rhs; }
friend mint operator-(const mint &lhs, const mint &rhs) { return mint(lhs) -= rhs; }
friend mint operator*(const mint &lhs, const mint &rhs) { return mint(lhs) *= rhs; }
friend mint operator/(const mint &lhs, const mint &rhs) { return mint(lhs) /= rhs; }
friend bool operator==(const mint &lhs, const mint &rhs) { return lhs._v == rhs._v; }
friend bool operator!=(const mint &lhs, const mint &rhs) { return lhs._v != rhs._v; }
friend M safe_hash_key(const mint &x) { return x.val(); }
#line 3472
friend void rd1(mint &x)
{
long long a;
fastio::rd1(a);
x = a;
}
friend void wt1(const mint &x)
{
fastio::wt1(x.val());
}
};
};
template <int mod>
using static_modint32 = internal::modint_impl<internal::policy_static<mod>>;
template <int id>
using dynamic_modint32 = internal::modint_impl<internal::policy_barrett32<id>>;
template <ll mod>
using static_modint64 = internal::modint_impl<internal::policy_static<mod>>;
template <int id>
using dynamic_modint64_odd = internal::modint_impl<internal::policy_montgomery64_odd<id>>;
template <int id>
using dynamic_modint64 = internal::modint_impl<internal::policy_montgomery64<id>>;
using modint998244353 = static_modint32<998244353>;
using modint1000000007 = static_modint32<1000000007>;
using modint = dynamic_modint32<-1>;
using modint61 = static_modint64<(1LL << 61) - 1>;
using modint64 = dynamic_modint64<-1>;
template <class T>
struct is_modint : std::false_type
{
};
template <class Policy, class PowInv>
struct is_modint<internal::modint_impl<Policy, PowInv>> : std::true_type
{
};
template <class T>
inline constexpr bool is_modint_v = is_modint<T>::value;
template <class T>
struct is_static_modint : false_type {};
template <int m, class PowInv>
struct is_static_modint<internal::modint_impl<internal::policy_static<m>, PowInv>> : true_type {};
template <ll m, class PowInv>
struct is_static_modint<internal::modint_impl<internal::policy_static<m>, PowInv>> : true_type {};
template <class T>
inline constexpr bool is_static_modint_v = is_static_modint<T>::value;
template <class T>
struct is_dynamic_modint : false_type {};
template <int id>
struct is_dynamic_modint<dynamic_modint32<id>> : true_type {};
template <int id>
struct is_dynamic_modint<dynamic_modint64_odd<id>> : true_type {};
template <int id>
struct is_dynamic_modint<dynamic_modint64<id>> : true_type {};
template <class T>
inline constexpr bool is_dynamic_modint_v = is_dynamic_modint<T>::value;
template <typename, typename = void>
struct has_mod : std::false_type
{
};
template <typename T>
struct has_mod<T, std::void_t<decltype(T::mod())>> : std::true_type
{
};
template <class mint>
struct modint_less
{
bool operator()(const mint &a, const mint &b) const
;
};
template <class mint>
struct modint_hash
{
auto operator()(const mint &x) const
;
};
template <class mint>
struct PowerTable
{
private:
decltype(mint::mod()) mod{};
mint base{};
vc<mint> pw;
public:
CPP_DUMP_DEFINE_DATA(mod, base, pw);
PowerTable() {}
PowerTable(const mint &base) : mod(mint::mod()), base(base), pw(1, 1) {}
void reserve(int n)
;
#line 3605
mint pow(int n)
;
};
template <class T>
struct Binomial
{
private:
inline static decltype(T::mod()) mod;
public:
inline static vc<T> fac_, finv_, inv_;
static void reserve(int n)
;
#line 3657
static T fac(int n)
;
static T finv(int n)
;
static T inv(int n)
;
static T P(int n, int k)
;
static T C(int n, int k)
;
static T H(int n, int k)
;
};
template <class mint>
mint stom(string s)
;
#line 3750
template <class T = ll, class U = larger_int_t<T>>
pair<T, T> svp2d(const pair<T, T> &a, const pair<T, T> &b)
;
#line 3788
template <class mint>
pair<decltype(mint(0).val()), decltype(mint(0).val())>
mint_to_rat(const mint &x)
{
auto [p, q] = svp2d({x.val(), 1}, {mint::mod(), 0});
if (q < 0)
p = -p, q = -q;
return {p, q};
}
namespace cpp_dump
{
template <class T>
struct rat_value
{
T p, q;
friend ostream &operator<<(ostream &os, const rat_value &x)
{
os << x.p;
if (x.q != 1)
os << '/' << x.q;
return os;
}
};
struct mint_to_rat_fn
{
template <class T>
requires requires(const T &x)
{
{ T::mod() } -> std::convertible_to<ll>;
{ x.val() } -> std::convertible_to<ll>;
}
constexpr auto operator()(const T &x) const
-> rat_value<typename decltype(mint_to_rat(x))::first_type>
{
auto [p, q] = mint_to_rat(x);
return {p, q};
}
};
template <class T>
struct rat_object
{
T value;
};
template <class T>
inline constexpr bool rat_is_map = false;
template <class... Args>
inline constexpr bool rat_is_map<std::map<Args...>> = true;
template <class... Args>
inline constexpr bool rat_is_map<std::multimap<Args...>> = true;
template <class... Args>
inline constexpr bool rat_is_map<std::unordered_map<Args...>> = true;
template <class... Args>
inline constexpr bool rat_is_map<std::unordered_multimap<Args...>> = true;
struct rat_closure : std::ranges::range_adaptor_closure<rat_closure>
{
template <class T>
static constexpr bool can_convert()
;
#line 3880
template <typename T>
requires(!requires { std::views::all | std::declval<T>(); })
constexpr decltype(auto) operator()(T &&t) const
{
if constexpr (!std::ranges::range<T> && std::invocable<mint_to_rat_fn, T>)
{
return mint_to_rat_fn{}(std::forward<T>(t));
}
else if constexpr (can_convert<T>())
{
if constexpr (!std::ranges::range<T> || rat_is_map<std::remove_cvref_t<T>>)
return rat_object<T>{std::forward<T>(t)};
else
return std::forward<T>(t) | std::views::transform(*this);
}
else
{
return std::forward<T>(t);
}
}
};
#line 3951
template <typename T>
requires(!requires { std::views::all | std::declval<T>(); })
constexpr decltype(auto) operator|(T &&t, rat_closure c)
{
return c(std::forward<T>(t));
}
constexpr rat_closure rat()
{
return rat_closure{};
}
}
using mint = modint998244353;
// using mint = modint1000000007;
// using mint = static_modint<1000000000>;
// using mint = modint;
using bi = Binomial<mint>;
void init()
{
oj(mt.seed(random_device()()));
}
#line 3989
namespace internal_monoid_power
{
template <class T, class I, class = void>
struct HasMemberScalarProduct : false_type {};
template <class T, class I>
struct HasMemberScalarProduct<T, I, void_t<decltype(T(declval<const T &>().operator*({declval<I>()})))>> : true_type {};
template <class T, class I, class = void>
struct HasFreeScalarProduct : false_type {};
template <class T, class I>
struct HasFreeScalarProduct<T, I, void_t<decltype(T(operator*(declval<const T &>(), {declval<I>()})))>> : true_type {};
template <class T, class I>
struct HasScalarProduct : bool_constant<is_arithmetic_v<T> || is_same_v<T, f128> ||
HasMemberScalarProduct<T, I>::value || HasFreeScalarProduct<T, I>::value> {};
template <class T, class I>
constexpr T scalar_product(const T &x, I k)
;
template <class T, class I>
constexpr T scale(const T &x, I k)
;
#line 4038
template <class T, class I>
constexpr T bounded_scale(const T &x, I k, const T &cap)
;
#line 4063
template <class M, class I>
using ConstRefPower = typename M::S (*)(const typename M::S &, I);
template <class M, class I>
using ValuePower = typename M::S (*)(typename M::S, I);
template <class M, class I, class = void>
struct HasConstRefPower : false_type {};
template <class M, class I>
struct HasConstRefPower<M, I, void_t<decltype(static_cast<ConstRefPower<M, I>>(M::pow))>> : true_type {};
template <class M, class I, class = void>
struct HasValuePower : false_type {};
template <class M, class I>
struct HasValuePower<M, I, void_t<decltype(static_cast<ValuePower<M, I>>(M::pow))>> : true_type {};
template <class M, class I>
constexpr typename M::S power(const typename M::S &x, I k)
;
}
template <class S_, auto op_, auto e_>
struct Monoid
{
using S = S_;
static constexpr auto op = op_;
static constexpr auto e = e_;
};
template <class S_, auto op_, auto e_, auto inv_>
struct Group
{
using S = S_;
static constexpr auto op = op_;
static constexpr auto e = e_;
static constexpr auto inv = inv_;
};
template <class S_, auto add_, auto e0_, auto mul_, auto e1_>
struct SemiRing
{
using S = S_;
static constexpr auto add = add_;
static constexpr auto e0 = e0_;
static constexpr auto mul = mul_;
static constexpr auto e1 = e1_;
};
template <class S_, auto add_, auto e0_, auto minus_, auto mul_, auto e1_>
struct Ring
{
using S = S_;
static constexpr auto add = add_;
static constexpr auto e0 = e0_;
static constexpr auto minus = minus_;
static constexpr auto mul = mul_;
static constexpr auto e1 = e1_;
};
template <class S_, auto add_, auto e0_, auto minus_, auto mul_, auto e1_, auto inv_>
struct Field
{
using S = S_;
static constexpr auto add = add_;
static constexpr auto e0 = e0_;
static constexpr auto minus = minus_;
static constexpr auto mul = mul_;
static constexpr auto e1 = e1_;
static constexpr auto inv = inv_;
};
template <class M, class = void>
struct HasInverse : false_type {};
template <class M>
struct HasInverse<M, void_t<decltype(M::inv(declval<const typename M::S &>()))>> : true_type {};
template <class M>
struct OppositeMonoid
{
using S = typename M::S;
static constexpr S op(const S &a, const S &b) ;
static constexpr auto e = M::e;
template <class I>
static constexpr S pow(const S &x, I k) { return internal_monoid_power::power<M>(x, k); }
};
template <class G>
struct OppositeGroup
{
using S = typename G::S;
static constexpr S op(const S &a, const S &b) ;
static constexpr auto e = G::e;
static constexpr auto inv = G::inv;
template <class I>
static constexpr S pow(const S &x, I k) { return internal_monoid_power::power<G>(x, k); }
};
template <class M>
struct NormalAndOppositeMonoid
{
struct S
{
typename M::S normal;
typename M::S opposite;
S() {}
template <class... Args,
std::enable_if_t<std::is_constructible_v<typename M::S, Args...>, std::nullptr_t> = nullptr>
S(Args &&...args)
: normal(std::forward<Args>(args)...), opposite(normal) {}
S rev() ;
S(const typename M::S &normal, const typename M::S &opposite) : normal(normal), opposite(opposite) {}
};
static constexpr S op(const S &a, const S &b) ;
static constexpr S e() ;
template <class I>
static constexpr S pow(const S &x, I k) { return {internal_monoid_power::power<M>(x.normal, k), internal_monoid_power::power<M>(x.opposite, k)}; }
};
template <class G>
struct NormalAndOppositeGroup
{
struct S
{
typename G::S normal;
typename G::S opposite;
S() {}
template <class... Args,
std::enable_if_t<std::is_constructible_v<typename G::S, Args...>, std::nullptr_t> = nullptr>
S(Args &&...args)
: normal(std::forward<Args>(args)...), opposite(normal) {}
S rev() ;
S(const typename G::S &normal, const typename G::S &opposite) : normal(normal), opposite(opposite) {}
};
static constexpr S op(const S &a, const S &b) ;
static constexpr S e() ;
static constexpr S inv(const S &a) { return {G::inv(a.normal), G::inv(a.opposite)}; }
template <class I>
static constexpr S pow(const S &x, I k) { return {internal_monoid_power::power<G>(x.normal, k), internal_monoid_power::power<G>(x.opposite, k)}; }
};
template <class SR>
using MonoidOfSemiRingAdd = Monoid<typename SR::S, SR::add, SR::e0>;
template <class SR>
using MonoidOfSemiRingMul = Monoid<typename SR::S, SR::mul, SR::e1>;
template <class R>
using GroupOfRingAdd = Group<typename R::S, R::add, R::e0, R::minus>;
template <class K>
using GroupOfFieldMul = Group<typename K::S, K::mul, K::e1, K::inv>;
template <class Madd, class Mmul>
struct SemiRingFromMonoidMonoid
{
static_assert(is_same_v<typename Madd::S, typename Mmul::S>, "Madd::S and Mmul::S must be identical");
using S = typename Madd::S;
static constexpr auto add = Madd::op;
static constexpr auto e0 = Madd::e;
static constexpr auto mul = Mmul::op;
static constexpr auto e1 = Mmul::e;
};
template <class Gadd, class Mmul>
struct RingFromGroupMonoid
{
static_assert(is_same_v<typename Gadd::S, typename Mmul::S>, "Gadd::S and Mmul::S must be identical");
using S = typename Gadd::S;
static constexpr auto add = Gadd::op;
static constexpr auto e0 = Gadd::e;
static constexpr auto minus = Gadd::inv;
static constexpr auto mul = Mmul::op;
static constexpr auto e1 = Mmul::e;
};
template <class Gadd, class Gmul>
struct FieldFromGroupGroup
{
static_assert(is_same_v<typename Gadd::S, typename Gmul::S>, "Gadd::S and Gmul::S must be identical");
using S = typename Gadd::S;
static constexpr auto add = Gadd::op;
static constexpr auto e0 = Gadd::e;
static constexpr auto minus = Gadd::inv;
static constexpr auto mul = Gmul::op;
static constexpr auto e1 = Gmul::e;
static constexpr auto inv = Gmul::inv;
};
template <class T, class = void>
struct has_e1 : false_type {};
template <class T>
struct has_e1<T, void_t<decltype(T::e1())>> : true_type {};
template <class T>
inline constexpr bool has_e1_v = has_e1<T>::value;
template <class T>
struct MonoidMul
{
using S = T;
static constexpr S op(S a, S b) ;
static constexpr S e()
;
};
template <class T>
struct GroupAddSub
{
using S = T;
static constexpr S op(S a, S b) ;
static constexpr S e()
;
static constexpr S inv(S a) { return -a; }
template <class I, class = decltype(declval<S>() * declval<I>())>
static constexpr S pow(const S &a, I k) { return a * k; }
};
template <class T>
struct GroupMulDiv
{
using S = T;
static constexpr S op(S a, S b) ;
static constexpr S e()
;
static constexpr S inv(S a) { return e() / a; }
};
template <class T, auto infty = nullptr>
using SemiRingMinPlus = SemiRingFromMonoidMonoid<MonoidMin<T, infty>, MonoidAdd<T>>;
template <class T, auto infty = nullptr>
using SemiRingMaxPlus = SemiRingFromMonoidMonoid<MonoidMax<T, infty>, MonoidAdd<T>>;
template <class T>
using RingAddSubMul = RingFromGroupMonoid<GroupAddSub<T>, MonoidMul<T>>;
template <class T>
using FieldAddSubMulDiv = FieldFromGroupGroup<GroupAddSub<T>, GroupMulDiv<T>>;
namespace internal
{
template <class S>
struct ordinary_mod32 : false_type {};
template <int mod, class PowInv>
struct ordinary_mod32<modint_impl<policy_static<mod>, PowInv>> : true_type {};
template <int id>
struct ordinary_mod32<dynamic_modint32<id>> : true_type {};
}
namespace internal
{
template <class G>
struct ordinary_mod32_add_group : false_type {};
template <class S>
struct ordinary_mod32_add_group<GroupAddSub<S>> : ordinary_mod32<S> {};
template <class S>
struct ordinary_mod32_add_group<Group<S, GroupAddSub<S>::op, GroupAddSub<S>::e, GroupAddSub<S>::inv>>
: ordinary_mod32<S> {};
template <class F>
struct ordinary_mod32_field : false_type {};
template <class S>
struct ordinary_mod32_field<FieldAddSubMulDiv<S>> : ordinary_mod32<S> {};
}
#include <algorithm>
#include <cstddef>
#include <cstdint>
namespace internal
{
template <class NextProduct>
inline std::uint32_t block_dot32(std::size_t n, std::uint32_t mod, std::size_t block,
std::uint32_t seed, NextProduct next)
{
std::uint64_t carry = seed;
while (n)
{
const std::size_t count = std::min(n, block);
std::size_t i = 0;
std::uint64_t s0 = carry, s1 = 0, s2 = 0, s3 = 0;
for (; i + 4 <= count; i += 4)
{
s0 += next();
s1 += next();
s2 += next();
s3 += next();
}
for (; i < count; ++i) s0 += next();
carry = (s0 + s1 + s2 + s3) % mod;
n -= count;
}
return static_cast<std::uint32_t>(carry);
}
}
namespace internal
{
template <class S>
using dot_product_mod32_value = ordinary_mod32<S>;
template <int block, class S, class It1, class It2>
S dot_product_mod32_impl(int n, It1 a, It2 b)
{
return S::raw(block_dot32(n, S::mod(), block, 0, [&]()
{
return ull((a++)->val()) * (b++)->val();
}));
}
template <class S, class It1, class It2>
S dot_product_mod32_values(int n, It1 a, It2 b)
{
static_assert(ordinary_mod32<S>::value);
assert(n >= 0);
if (S::mod() <= (1 << 30)) return dot_product_mod32_impl<16, S>(n, a, b);
return dot_product_mod32_impl<4, S>(n, a, b);
}
}
namespace internal
{
template <class SR>
struct dot_product_mod32 : false_type {};
template <class S>
struct dot_product_mod32<FieldAddSubMulDiv<S>> : dot_product_mod32_value<S> {};
template <class S>
struct dot_product_mod32<RingAddSubMul<S>> : dot_product_mod32_value<S> {};
template <class S>
struct dot_product_mod32<SemiRingFromMonoidMonoid<MonoidAdd<S>, MonoidMul<S>>>
: dot_product_mod32_value<S> {};
}
template <class SR, class It1, class It2>
typename SR::S dot_product(int n, It1 a, It2 b)
;
#line 4474
//! https://codeforces.com/blog/entry/61306
template <class F>
struct BerlekampMassey
{
using S = typename F::S;
vc<S> a, b, c, tmp;
int pos = -1;
S inv_x = F::e0();
CPP_DUMP_DEFINE_DATA(a, c, pos);
int degree() const ;
void append(S value)
;
#line 4522
vc<S> coefficients() const
;
};
template <class F>
vc<typename F::S> berlekamp_massey(const vc<typename F::S> &a)
;
#line 4552
uint32_t bit_reverse32(uint32_t x)
{
x = ((x & 0x55555555) << 1) | ((x & 0xAAAAAAAA) >> 1);
x = ((x & 0x33333333) << 2) | ((x & 0xCCCCCCCC) >> 2);
x = ((x & 0x0F0F0F0F) << 4) | ((x & 0xF0F0F0F0) >> 4);
x = ((x & 0x00FF00FF) << 8) | ((x & 0xFF00FF00) >> 8);
return (x << 16) | (x >> 16);
}
int bitrev(int pw2, int i)
{
assert(pw2 > 0 && has_single_bit((uint)pw2));
assert(0 <= i && i < pw2);
if (pw2 == 1)
return 0;
return bit_reverse32(i) >> (32 - countr_zero((uint)pw2));
}
/*!
* Portions corresponding to AtCoder Library (CC0-1.0): primitive-root search and NTT butterflies.
* Upstream: atcoder/convolution.hpp and atcoder/internal_math.hpp
* https://github.com/atcoder/ac-library/tree/864245a00b00dd008d1abfdc239618fdb7d139da
* License: https://github.com/atcoder/ac-library/blob/864245a00b00dd008d1abfdc239618fdb7d139da/LICENSE
*/
namespace internal
{
template <class T>
constexpr T primitive_root_constexpr(T m)
{
if (m == 2)
return 1;
if (m == 167772161)
return 3;
if (m == 469762049)
return 3;
if (m == 754974721)
return 11;
if (m == 998244353)
return 3;
if (m == 1107296257)
return 10;
if (m == 1711276033)
return 29;
if (m == 1811939329)
return 13;
if (m == 2013265921)
return 31;
if (m == 2113929217)
return 5;
T divs[20] = {};
divs[0] = 2;
int cnt = 1;
T x = (m - 1) / 2;
while (x % 2 == 0)
x /= 2;
for (T i = 3; i <= x / i; i += 2)
{
if (x % i == 0)
{
divs[cnt++] = i;
while (x % i == 0)
{
x /= i;
}
}
}
if (x > 1)
{
divs[cnt++] = x;
}
for (T g = 2;; g++)
{
bool ok = true;
for (int i = 0; i < cnt; i++)
{
if (powmod_constexpr(g, (m - 1) / divs[i], m) == 1)
{
ok = false;
break;
}
}
if (ok)
return g;
}
}
template <auto m>
constexpr auto primitive_root_for_convolution = primitive_root_constexpr(m);
template <class T>
constexpr T ntt_generator_constexpr(T m)
{
if (m <= INT_MAX) return primitive_root_constexpr(m);
if ((m - 1) % 4 == 2) return m - 1;
for (T g = 2; ; ++g)
if (powmod_constexpr(g, (m - 1) / 2, m) == m - 1) return g;
}
template <auto m>
constexpr auto ntt_generator_for_convolution = ntt_generator_constexpr(m);
template <class mint, auto g = internal::ntt_generator_for_convolution<mint::mod()>>
struct fft_info
{
static constexpr int rank2 = countr_zero(mint::mod() - 1);
std::array<mint, std::max(3, rank2 + 1)> root;
std::array<mint, std::max(3, rank2 + 1)> iroot;
std::array<mint, std::max(1, rank2 - 2 + 1)> rate2;
std::array<mint, std::max(1, rank2 - 2 + 1)> irate2;
std::array<mint, std::max(1, rank2 - 3 + 1)> rate3;
std::array<mint, std::max(1, rank2 - 3 + 1)> irate3;
fft_info()
{
root[rank2] = mint(g).pow((mint::mod() - 1) >> rank2);
iroot[rank2] = root[rank2].inv();
for (int i = rank2 - 1; i >= 0; i--)
{
root[i] = root[i + 1] * root[i + 1];
iroot[i] = iroot[i + 1] * iroot[i + 1];
}
auto set_rates = [&](int shift, auto &rate, auto &irate)
{
mint prod = 1, iprod = 1;
for (int i = 0; i <= rank2 - shift; ++i)
{
const int j = i + shift;
rate[i] = root[j] * prod;
irate[i] = iroot[j] * iprod;
prod *= iroot[j];
iprod *= root[j];
}
};
set_rates(2, rate2, irate2);
set_rates(3, rate3, irate3);
}
};
}
template <class mint>
constexpr bool ntt_ok(int n)
{
if (n <= 0)
return false;
if constexpr (is_static_modint_v<mint>)
{
if constexpr (!internal::isprime<mint::mod()>)
return false;
constexpr int rank2 = countr_zero(mint::mod() - 1);
return ull(n) <= (1ULL << rank2);
}
else
return false;
}
template <int id>
void ntt(vc<dynamic_modint32<id>> &) ;
template <int id>
void intt(vc<dynamic_modint32<id>> &) ;
template <auto mod, class PowInv>
void ntt(vc<internal::modint_impl<internal::policy_static<mod>, PowInv>> &a)
{
using mint = internal::modint_impl<internal::policy_static<mod>, PowInv>;
int n = int(a.size());
assert(n > 0);
int h = countr_zero((unsigned int)n);
assert(n == (1 << h));
assert(ntt_ok<mint>(n));
using wide = conditional_t<(mint::mod() <= INT_MAX), ull, u128>;
static const internal::fft_info<mint> info;
int len = 0;
while (len < h)
{
if (h - len == 1)
{
mint rot = 1;
for (int s = 0; s < (1 << len); s++)
{
const int i = 2 * s;
auto l = a[i], r = a[i + 1] * rot;
a[i] = l + r;
a[i + 1] = l - r;
if (s + 1 != (1 << len))
rot *= info.rate2[countr_zero(~(unsigned int)(s))];
}
len++;
}
else
{
int p = 1 << (h - len - 2);
mint rot = 1, imag = info.root[2];
for (int s = 0; s < (1 << len); s++)
{
mint rot2 = rot * rot;
mint rot3 = rot2 * rot;
int offset = s << (h - len);
for (int i = 0; i < p; i++)
{
const int j = i + offset;
auto mod2 = wide(mint::mod()) * mint::mod();
auto a0 = wide(a[j].val());
auto a1 = wide(a[j + p].val()) * rot.val();
auto a2 = wide(a[j + 2 * p].val()) * rot2.val();
auto a3 = wide(a[j + 3 * p].val()) * rot3.val();
auto a1na3imag =
wide(mint(a1 + mod2 - a3).val()) * imag.val();
auto na2 = mod2 - a2;
a[j] = a0 + a2 + a1 + a3;
a[j + p] = a0 + a2 + (2 * mod2 - (a1 + a3));
a[j + 2 * p] = a0 + na2 + a1na3imag;
a[j + 3 * p] = a0 + na2 + (mod2 - a1na3imag);
}
if (s + 1 != (1 << len))
rot *= info.rate3[countr_zero(~(unsigned int)(s))];
}
len += 2;
}
}
}
template <auto mod, class PowInv>
void intt(vc<internal::modint_impl<internal::policy_static<mod>, PowInv>> &a)
{
using mint = internal::modint_impl<internal::policy_static<mod>, PowInv>;
int n = int(a.size());
assert(n > 0);
int h = countr_zero((unsigned int)n);
assert(n == (1 << h));
assert(ntt_ok<mint>(n));
using wide = conditional_t<(mint::mod() <= INT_MAX), ull, u128>;
static const internal::fft_info<mint> info;
int len = h;
while (len)
{
if (len == 1)
{
const int p = 1 << (h - 1);
for (int i = 0; i < p; i++)
{
auto l = a[i], r = a[i + p];
a[i] = l + r;
a[i + p] = wide(mint::mod()) + l.val() - wide(r.val());
}
len--;
}
else
{
int p = 1 << (h - len);
mint irot = 1, iimag = info.iroot[2];
for (int s = 0; s < (1 << (len - 2)); s++)
{
mint irot2 = irot * irot;
mint irot3 = irot2 * irot;
int offset = s << (h - len + 2);
for (int i = 0; i < p; i++)
{
const int j = i + offset;
auto a0 = wide(a[j].val());
auto a1 = wide(a[j + p].val());
auto a2 = wide(a[j + 2 * p].val());
auto a3 = wide(a[j + 3 * p].val());
auto a2na3iimag =
wide(mint((wide(mint::mod()) + a2 - a3) * iimag.val()).val());
a[j] = a0 + a1 + a2 + a3;
a[j + p] =
(a0 + (mint::mod() - a1) + a2na3iimag) * irot.val();
a[j + 2 * p] =
(a0 + a1 + (mint::mod() - a2) + (mint::mod() - a3)) *
irot2.val();
a[j + 3 * p] =
(a0 + (mint::mod() - a1) + (mint::mod() - a2na3iimag)) *
irot3.val();
}
if (s + 1 != (1 << (len - 2)))
irot *= info.irate3[countr_zero(~(unsigned int)(s))];
}
len -= 2;
}
}
}
template <class T>
T convolution_point_get(const vc<T> &a, const vc<T> &b, int p)
{
const int n = a.size(), m = b.size();
if constexpr (internal::dot_product_mod32_value<T>::value)
{
if (p < 0 || ll(p) >= ll(n) + m - 1 || n == 0 || m == 0) return T(0);
const int l = max(0, p - m + 1), r = min(n, p + 1);
return internal::dot_product_mod32_values<T>(r - l, a.begin() + l, b.rbegin() + (m - 1 - p + l));
}
T res = 0;
repi(i, max(0, p - m + 1), min(n, p + 1)) res += a[i] * b[p - i];
return res;
}
namespace internal
{
template <class mint>
vc<mint> convolution_naive(const vc<mint> &a, const vc<mint> &b)
{
const int n = a.size(), m = b.size();
const int cnta = n - count(ALL(a), 0), cntb = m - count(ALL(b), 0);
vc<mint> c(n + m - 1);
if (ll(cnta) * 2 < n && ll(cntb) * 2 < m)
{
vc<pair<int, mint>> an, bn;
repi(i, n) if (a[i] != 0) an.eb(i, a[i]);
repi(j, m) if (b[j] != 0) bn.eb(j, b[j]);
for (auto [i, x] : an) for (auto [j, y] : bn) c[i+j] += x*y;
return c;
}
if constexpr (dot_product_mod32_value<mint>::value)
{
if (min(n, m) >= 16 && ll(cnta) * 2 >= n && ll(cntb) * 2 >= m)
{
repi(p, n + m - 1) c[p] = convolution_point_get(a, b, p);
return c;
}
}
if ((ll)m * cnta > (ll)n * cntb)
{
repi(j, m)
{
if (b[j] == 0)
continue;
repi(i, n) c[i + j] += a[i] * b[j];
}
}
else
{
repi(i, n)
{
if (a[i] == 0)
continue;
repi(j, m) c[i + j] += a[i] * b[j];
}
}
return c;
}
template <class mint>
pair<int, int> convolution_support(const vc<mint> &a)
{
int l = 0, r = a.size();
while (l < r && a[l] == 0) ++l;
while (l < r && a[r - 1] == 0) --r;
return {l, r};
}
inline int convolution_ntt_length(int n, int m)
{
const int size = n + m - 1, z = bit_ceil(size), half = z / 2;
return size - half <= 32 && max(n, m) <= half ? half : z;
}
template <class mint>
bool convolution_sparse_preferred(const vc<mint> &a, const vc<mint> &b, int ca, int cb, int moduli)
{
const bool sparse_a = ll(ca) * 2 < int(a.size()), sparse_b = ll(cb) * 2 < int(b.size());
if (!sparse_a && !sparse_b) return false;
const ll work = sparse_a && sparse_b ? ll(ca) * cb : min(ll(a.size()) * cb, ll(b.size()) * ca);
const auto [al, ar] = convolution_support(a);
const auto [bl, br] = convolution_support(b);
const int z = convolution_ntt_length(ar - al, br - bl);
const ll budget = ll(moduli) * z * countr_zero(unsigned(z)) / 4;
return work <= budget;
}
template <class mint>
vc<mint> convolution_ntt(vc<mint> a, vc<mint> b)
{
if (a.empty() || b.empty()) return {};
const int full_size = int(a.size()) + int(b.size()) - 1;
const auto [al, ar] = convolution_support(a);
const auto [bl, br] = convolution_support(b);
if (al == ar || bl == br) return vc<mint>(full_size);
int z = convolution_ntt_length(a.size(), b.size()), shift = 0;
const int cropped = convolution_ntt_length(ar - al, br - bl);
if (cropped < z)
{
shift = al + bl;
a.resize(ar), b.resize(br);
if (al) a.erase(a.begin(), a.begin() + al);
if (bl) b.erase(b.begin(), b.begin() + bl);
z = cropped;
}
const int size = int(a.size()) + int(b.size()) - 1;
vc<mint> tail(max(0, size - z));
repi(i, tail.size()) tail[i] = convolution_point_get(a, b, z + i);
if (a == b)
{
a.resize(z);
ntt(a);
repi(i, z) a[i] *= a[i];
}
else
{
a.resize(z), b.resize(z);
ntt(a), ntt(b);
repi(i, z) a[i] *= b[i];
}
intt(a);
a.resize(min(z, size));
const mint iz = mint(z).inv();
for (auto &x : a) x *= iz;
a.resize(size);
repi(i, tail.size())
{
a[i] -= tail[i];
a[z + i] = tail[i];
}
if (shift) a.insert(a.begin(), shift, mint(0));
a.resize(full_size);
return a;
}
}
template <class mint>
vc<mint> ntt_convolution(const vc<mint> &a, const vc<mint> &b)
;
#line 5024
/*!
* Portions corresponding to AtCoder Library (CC0-1.0): two-congruence merge.
* Upstream: atcoder/math.hpp
* https://github.com/atcoder/ac-library/tree/864245a00b00dd008d1abfdc239618fdb7d139da
* License: https://github.com/atcoder/ac-library/blob/864245a00b00dd008d1abfdc239618fdb7d139da/LICENSE
*/
template <class T = ll, class R0, class R1, class M0, class M1>
constexpr optional<pair<T, T>> crt2(R0 r0_, R1 r1_, M0 m0_, M1 m1_)
;
#line 5066
template <class T = ll, class V1, class V2>
constexpr optional<pair<T, T>> crt(const V1 &rs, const V2 &ms)
;
#line 5081
namespace crt_internal
{
namespace crt_word_detail
{
template <class U>
constexpr U add_mod(U a, U b, U m)
;
template <class U>
constexpr U sub_mod(U a, U b, U m)
;
template <class U>
constexpr U mul_mod(U a, U b, U m)
;
#line 5117
template <class U>
constexpr U inverse_mod(U a, U m)
;
#line 5137
template <class U, class R>
constexpr U normalize(R x, U m)
;
template <class mint, class U, class V, class A, class B>
constexpr void append(U residue, U modulus, size_t i, const V &ms,
A &remainders, B &coefficients, mint &answer, mint &period)
{
using M = remove_cv_t<typename V::value_type>;
const U digit = mul_mod(sub_mod(residue, remainders[i], modulus),
inverse_mod(coefficients[i], modulus), modulus);
answer = answer + period * mint(M(digit));
period = period * mint(ms[i]);
for (size_t j = i + 1; j < ms.size(); ++j)
{
const U m = U(ms[j]);
if (m == 1) continue;
const U coefficient = coefficients[j];
remainders[j] = add_mod(remainders[j], mul_mod(coefficient, digit % m, m), m);
coefficients[j] = mul_mod(coefficient, modulus % m, m);
}
}
}
template <class mint, class V1, class V2>
pair<mint, mint> crt_mod_word(const V1 &rs, const V2 &ms)
;
#line 5203
template <class mint, class R, class M0, size_t n>
constexpr pair<mint, mint> crt_mod_word_constexpr(const array<R, n> &rs, const array<M0, n> &ms)
;
#line 5238
namespace crt_word_dispatch
{
template <class V>
constexpr bool runtime_requires_word(const V &ms)
;
#line 5255
}
template <class mint, class U1, class U2, size_t n>
constexpr pair<mint, mint> crt_mod_normalized(const array<U1, n> &rs, const array<U2, n> &ms)
{
using T = larger_int_t<U2>;
assert(rs.size() == ms.size());
mint r = 0, m = 1;
array<T, n> rr{}, mm{};
fill(ALL(mm), 1);
repi(i, n)
{
assert(ms[i] >= U2(1));
assert(U1(0) <= rs[i] && U2(rs[i]) < ms[i]);
auto [g, im, _] = extgcd<T>(mm[i], ms[i]);
assert(g == 1);
T t = safemod<T>((rs[i] - rr[i]) * im, ms[i]);
r += t * m, m *= ms[i];
repi(j, i + 1, n)
{
rr[j] += t * mm[j] % ms[j];
if (rr[j] >= ms[j])
rr[j] -= ms[j];
mm[j] *= ms[i], mm[j] %= ms[j];
}
}
return {r, m};
}
}
template <class mint, class V1, class V2>
pair<mint, mint> crt_mod(const V1 &rs, const V2 &ms)
;
#line 5379
template <class mint, class U1, class U2, size_t n>
constexpr pair<mint, mint> crt_mod_constexpr(const array<U1, n> &rs, const array<U2, n> &ms)
;
#line 5405
namespace pre_crt_v2_adl
{
void gcd() = delete;
template <class T>
auto call(const T &a, const T &b) -> decltype(gcd(a, b)) ;
}
namespace pre_crt_v2_detail
{
template <class T, bool = is_integral_ext<T>> struct work_type { using type = T; };
template <class T> struct work_type<T, true> { using type = make_unsigned_ext_t<T>; };
template <class T, class = void> struct has_adl_gcd : false_type {};
template <class T> struct has_adl_gcd<T, void_t<decltype(pre_crt_v2_adl::call(declval<T>(), declval<T>()))>> : true_type {};
template <class T>
constexpr bool unbounded_integer = numeric_limits<T>::is_specialized &&
numeric_limits<T>::is_integer && !numeric_limits<T>::is_bounded;
template <class T>
T gcd(T a, T b)
;
#line 5439
template <class T>
T add_mod(const T &a, const T &b, const T &m)
;
template <class T, class R>
T normalize(R x, const T &m)
;
#line 5488
template <class T>
void distribute(T &a, T &b, T g)
;
}
template <class V1, class V2>
bool pre_crt(const V1 &rs, V2 &ms)
;
#line 5562
namespace internal
{
template <size_t j, int mod, class T, size_t k>
void convolution_crt_helper(const vc<T> &a, const vc<T> &b, vc<array<T, k>> &cs)
{
using mint = static_modint32<mod>;
const int n = a.size(), m = b.size();
auto c = convolution_ntt(vc<mint>(ALL(a)), vc<mint>(ALL(b)));
repi(i, n + m - 1) cs[i][j] = c[i].val();
}
template <int ...ms, class T>
vc<T> convolution_crt(const vc<T> &a, const vc<T> &b)
;
#line 5599
template <class mint, int... ms, class T>
vc<mint> convolution_crt_mod(const vc<T> &a, const vc<T> &b)
{
const int n = a.size(), m = b.size();
constexpr size_t k = sizeof...(ms);
vc<array<T, k>> cs(n + m - 1);
constexpr array<int, k> ms_arr = {ms...};
[&]<size_t... Is>(index_sequence<Is...>)
{
(convolution_crt_helper<Is, ms_arr[Is], T, k>(a, b, cs), ...);
}(make_index_sequence<k>{});
vc<mint> c(n + m - 1);
repi(i, n + m - 1) c[i] = crt_internal::crt_mod_normalized<mint>(cs[i], ms_arr).first;
return c;
}
}
template <class mint, typename = std::enable_if_t<!is_integral_ext<mint>>>
vc<mint> convolution(const vc<mint> &a, const vc<mint> &b)
{
const int n = a.size(), m = b.size();
if (n == 0 || m == 0)
return {};
const int cnta = n - count(ALL(a), 0), cntb = m - count(ALL(b), 0);
if (cnta == 0 || cntb == 0) return vc<mint>(n + m - 1);
if (ntt_ok<mint>(n + m - 1))
{
if (min(cnta, cntb) <= 60 || internal::convolution_sparse_preferred(a, b, cnta, cntb, 1))
return internal::convolution_naive(a, b);
return internal::convolution_ntt(a, b);
}
else
{
if (min(cnta, cntb) <= 300 || internal::convolution_sparse_preferred(a, b, cnta, cntb, 3))
return internal::convolution_naive(a, b);
assert(ntt_ok<static_modint32<469762049>>(n + m - 1) && "|a| + |b| - 1 <= 2^26");
vc<ll> a_(n), b_(m);
repi(i, n) a_[i] = a[i].val();
repi(j, m) b_[j] = b[j].val();
return internal::convolution_crt_mod<mint, 469762049, 1811939329, 2013265921>(a_, b_);
}
}
template <int mod = 998244353, class T, typename = enable_if_t<is_integral_ext<T>>>
vc<T> convolution(const vc<T> &a, const vc<T> &b)
;
/*!
* Portions corresponding to AtCoder Library (CC0-1.0): Miller-Rabin termination logic.
* Upstream: atcoder/internal_math.hpp
* https://github.com/atcoder/ac-library/tree/864245a00b00dd008d1abfdc239618fdb7d139da
* License: https://github.com/atcoder/ac-library/blob/864245a00b00dd008d1abfdc239618fdb7d139da/LICENSE
*/
namespace internal
{
template <class mint, class Array>
bool is_prime_impl(ll n, const Array &bases)
{
if (n <= 1)
return false;
if (n == 2 || n == 7 || n == 61)
return true;
if (n % 2 == 0)
return false;
ll d = (n - 1) >> countr_zero(n - 1);
mint::set_mod(n);
for (ll a : bases)
{
ll t = d;
mint y = mint(a).pow(t);
while (t != n - 1 && y != 1 && y != n - 1)
{
y *= y;
t <<= 1;
}
if (y != n - 1 && t % 2 == 0)
return false;
}
return true;
}
};
bool is_prime(ll n)
{
static constexpr array<ll, 3> bases32 = {2, 7, 61};
static constexpr array<ll, 7> bases64 = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
if (n <= INT_MAX)
{
using mint = dynamic_modint32<INT_MIN>;
return internal::is_prime_impl<mint>(n, bases32);
}
else
{
using mint = dynamic_modint64_odd<INT_MIN>;
return internal::is_prime_impl<mint>(n, bases64);
}
}
template <class mint>
optional<mint> sqrt_mod(mint n)
;
#line 5765
namespace internal
{
enum class FPSSparseOperation
{
inv,
div,
exp,
pow,
sqrt
};
inline ll fps_ntt_weight(int n)
{
return ll(n) * countr_zero(unsigned(n));
}
inline ll fps_inv_weight(int n, bool native)
{
int N = bit_ceil(n);
if (native && n >= 512 && ll(n) * 4 <= ll(N) * 3)
{
int s = min(8, N / 2), m = N / (2 * s), r = (n + m - 1) / m;
return (4 * r + s - 3) * fps_ntt_weight(2 * m) + fps_inv_weight(m, true);
}
ll cost = 0;
for (int m = 1; m < n; m *= 2)
cost += native ? 5 * fps_ntt_weight(2 * m) : 6 * fps_ntt_weight(2 * m) + 9 * fps_ntt_weight(4 * m);
return cost;
}
inline ll fps_exp_bs_weight(int n)
{
int N = bit_ceil(n);
ll cost = 0;
for (int m = 2; m < n; m *= 2)
cost += 5 * fps_ntt_weight(m) + 6 * fps_ntt_weight(2 * m);
if (N >= 4)
cost -= 3 * fps_ntt_weight(N) - 5 * fps_ntt_weight(N / 2);
return cost;
}
inline ll fps_exp_weight(int n, bool native)
{
if (native && n >= 512)
{
int N = bit_ceil(n), s = min(n <= 4096 ? 4 : 16, N / 2), m = N / (2 * s);
int high = (n + m - 1) / m - s;
return (7 * s + 6 * high - 4) * fps_ntt_weight(2 * m) + fps_exp_bs_weight(m) + fps_inv_weight(m, true);
}
if (native)
return fps_exp_bs_weight(n);
ll cost = 0;
for (int m = 2; m < n; m *= 2)
cost += 6 * fps_ntt_weight(m) + 36 * fps_ntt_weight(2 * m);
return cost;
}
inline ll fps_sqrt_weight(int n, bool native)
{
if (native && n >= 512)
{
int N = bit_ceil(n), s = min(16, N), m = N / s, r = (n + m - 1) / m;
return (4 * r - 3) * fps_ntt_weight(2 * m) + fps_sqrt_weight(m, true) + fps_inv_weight(m, true);
}
ll cost = 0;
for (int m = 1; m < n; m *= 2)
cost += (native ? 5 : 15) * fps_ntt_weight(2 * m) + fps_inv_weight(m, native);
return cost;
}
template <FPSSparseOperation op, class F>
bool fps_use_sparse(const F &f, int n, int offset = 0)
{
using mint = typename F::value_type;
if (n <= 64)
return true;
bool native = ntt_ok<mint>(bit_ceil(n));
ll cost;
if constexpr (op == FPSSparseOperation::exp)
cost = fps_exp_weight(n, native);
else if constexpr (op == FPSSparseOperation::sqrt)
cost = fps_sqrt_weight(n, native);
else
cost = fps_inv_weight(n, native);
int scale = 14, repetitions = 1;
if constexpr (op == FPSSparseOperation::div || op == FPSSparseOperation::pow)
cost += (ntt_ok<mint>(2 * n - 1) ? 3 : 9) * fps_ntt_weight(bit_ceil(2 * n - 1));
if constexpr (op == FPSSparseOperation::exp)
scale = 18;
if constexpr (op == FPSSparseOperation::pow)
cost += fps_exp_weight(n, native), repetitions = 2;
if constexpr (op == FPSSparseOperation::sqrt)
scale = 17, repetitions = 2;
if (!native)
scale = 25;
if constexpr (!internal::ordinary_mod32<mint>::value)
scale = max(1, scale / 2);
ll budget = cost * scale / (10 * repetitions), work = 0;
for (int j = 1; j < n && offset + j < f.sz(); ++j)
if (f[offset + j] != 0)
{
work += n - j;
if (work > budget)
return false;
}
return true;
}
template <class F>
optional<F> fps_tail_power(const F &f, int n, typename F::value_type constant, typename F::value_type linear,
typename F::value_type quadratic)
{
int d = 1;
while (d < min(n, f.sz()) && f[d] == 0)
++d;
if (d < min(n, f.sz()) && ll(d) * 3 < n)
return nullopt;
F res = f.resized(n) * linear;
if (n)
res[0] = constant;
if (ll(d) * 2 < n && d < f.sz())
{
int need = n - 2 * d;
F h(need);
repi(i, min(need, f.sz() - d)) h[i] = f[d + i];
F square = (h * h).resized(need) * quadratic;
repi(i, need) res[2 * d + i] += square[i];
}
return res;
}
template <class F> F fps_inv_newton(const F &input, int n)
{
using mint = typename F::value_type;
if (n == 0)
return {};
F f, g2, g{input.get(0).inv()};
for (int m = 1; m < n; m *= 2)
{
if (ntt_ok<mint>(2 * m))
{
f = input.resized(2 * m), g2 = F(g);
ntt(f);
g2.resize(2 * m), ntt(g2);
repi(i, 2 * m) f[i] *= g2[i];
intt(f);
f >>= m;
f.resize(2 * m), ntt(f);
repi(i, 2 * m) f[i] *= g2[i];
intt(f);
mint iz = mint(2 * m).inv();
iz *= -iz;
repi(i, m) f[i] *= iz;
g.insert(g.end(), f.begin(), f.begin() + m);
}
else
g = (g * mint(2) - g * g * input.resized(2 * m)).resized(2 * m);
}
return g.resized(n);
}
template <class F, class Visit>
inline void fps_for_block_pair(const vc<F> &a, const vc<F> &b, int k, Visit visit)
{
for (int i = max(0, k - int(b.size()) + 1); i <= k && i < int(a.size()); ++i)
{
int j = k - i;
if (!a[i].empty() && !b[j].empty()) visit(i, j);
}
}
template <class F> inline void fps_intt_prefix(F &f, int m)
{
intt(f);
f.resize(m);
}
template <class F> inline void fps_intt_block(F &f, int m)
{
using mint = typename F::value_type;
fps_intt_prefix(f, m);
mint iz = mint(2 * m).inv();
for (auto &v : f) v *= iz;
}
template <class F>
inline void fps_ntt_product_prefix(F &f, const F &transformed, int n)
{
ntt(f);
repi(i, 2 * n) f[i] *= transformed[i];
fps_intt_prefix(f, n);
}
template <class F> F fps_block_product(const vc<F> &a, const vc<F> &b, int k, int m, bool delayed = false)
{
using mint = typename F::value_type;
F freq(2 * m);
bool square = &a == &b;
if constexpr (internal::ordinary_mod32<mint>::value)
if (delayed)
{
const ull mod = mint::mod(), max_product = (mod - 1) * (mod - 1);
const ull batch = numeric_limits<ull>::max() / max_product;
square = square && mod <= (1ULL << 31);
vc<ull> pos(2 * m), neg(2 * m);
ull cnt = 0;
auto flush = [&]()
{
repi(j, 2 * m)
{
freq[j] += mint(pos[j] % mod);
freq[j] -= mint(neg[j] % mod);
pos[j] = neg[j] = 0;
}
cnt = 0;
};
for (int carry : {0, 1})
fps_for_block_pair(a, b, k - carry, [&](int i, int j)
{
if (square && i > j)
return;
int times = square && i < j ? 2 : 1;
if (cnt + times > batch)
flush();
repi(t, 2 * m)
{
ull v = ull(a[i][t].val()) * b[j][t].val() * times;
if (carry && t >= m)
neg[t] += v;
else
pos[t] += v;
}
cnt += times;
if (cnt == batch)
flush();
});
if (cnt)
flush();
fps_intt_block(freq, m);
return freq;
}
fps_for_block_pair(a, b, k, [&](int i, int j)
{
repi(t, 2 * m) freq[t] += a[i][t] * b[j][t];
});
fps_for_block_pair(a, b, k - 1, [&](int i, int j)
{
repi(t, 2 * m)
{
mint v = a[i][t] * b[j][t];
if (t < m)
freq[t] += v;
else
freq[t] -= v;
}
});
fps_intt_block(freq, m);
return freq;
}
template <class F> F fps_exp_bs(const F &h, int n, bool split_last = true)
;
#line 6096
template <class F> F fps_inv_block(const F &f, int n, int blocks, bool partial = true)
{
using mint = typename F::value_type;
assert(n >= 0 && f.get(0) != 0);
if (!n)
return {};
if (n == 1)
return {f[0].inv()};
int N = bit_ceil(n), s = min(blocks, max(1, N / 2)), m = N / (2 * s);
int r = partial ? (n + m - 1) / m : 2 * s, h = r - s;
vc<F> ff(r), g(s), gf(s), df(h);
g[0] = f.resized(m).inv(m);
gf[0] = g[0].resized(2 * m);
ntt(gf[0]);
repi(k, r)
{
ff[k].resize(2 * m);
repi(i, m) ff[k][i] = f.get(k * m + i);
ntt(ff[k]);
}
repi(k, 1, s)
{
F v = fps_block_product(ff, gf, k, m, true).resized(2 * m);
fps_ntt_product_prefix(v, gf[0], m);
mint iz = -mint(2 * m).inv();
for (auto &x : v)
x *= iz;
g[k] = std::move(v);
gf[k] = g[k].resized(2 * m);
ntt(gf[k]);
}
repi(k, h)
{
df[k] = -fps_block_product(ff, gf, s + k, m, true);
df[k].resize(2 * m);
ntt(df[k]);
}
F res(r * m);
repi(k, s) copy(g[k].begin(), g[k].end(), res.begin() + k * m);
repi(k, h)
{
F v = fps_block_product(df, gf, k, m, true);
copy(v.begin(), v.end(), res.begin() + (s + k) * m);
}
res.resize(n);
return res;
}
template <class F> F fps_sqrt_newton(const F &f, int n)
;
#line 6160
template <class F> F fps_sqrt_block(const F &f, int n, int blocks, bool partial = true)
;
#line 6193
template <class F>
F fps_exp_block_product(const vc<F> &a, const vc<F> &b, int k, int m, bool delayed,
vc<ull> &pos, vc<ull> &neg)
;
#line 6266
template <class F> F fps_exp_block(const F &f, int n, int blocks, bool delayed = true)
;
#line 6336
}
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <numeric>
#include <type_traits>
#include <utility>
#include <vector>
namespace internal::sparse_fps
{
using std::size_t, std::uint32_t, std::uint64_t;
template <class Mint, class = void>
struct has_canonical32_interface : std::false_type
{
};
template <class Mint>
struct has_canonical32_interface<
Mint, std::void_t<decltype(Mint::mod()),
decltype(std::declval<const Mint &>().val())>>
: std::bool_constant<
std::is_integral_v<decltype(Mint::mod())> &&
std::is_integral_v<decltype(std::declval<const Mint &>().val())> &&
sizeof(decltype(Mint::mod())) <= sizeof(uint32_t) &&
sizeof(decltype(std::declval<const Mint &>().val())) <=
sizeof(uint32_t)>
{
};
template <class Mint>
struct input_shape
{
int stride = 1;
int length = 0;
std::vector<std::pair<int, Mint>> terms;
};
template <class Fps>
auto describe_input(const Fps &f, int n)
;
#line 6405
inline size_t block_size32(uint32_t p)
{
assert(p >= 2);
const uint64_t largest = uint64_t(p) - 1;
const auto safe = (std::numeric_limits<uint64_t>::max() - largest) / (largest * largest);
return static_cast<size_t>(std::min<uint64_t>(16, safe));
}
template <bool UseCanonical32, class Fps>
Fps logarithm(const Fps &f, int n)
;
#line 6500
}
template <class mint>
struct FormalPowerSeries : vc<mint>
{
using F = FormalPowerSeries;
using vc<mint>::vc;
using vc<mint>::operator=;
using vc<mint>::size;
using vc<mint>::empty;
using vc<mint>::back;
using vc<mint>::pop_back;
using vc<mint>::begin;
using vc<mint>::resize;
using vc<mint>::front;
FormalPowerSeries(const vc<mint> &f) : vc<mint>(f) {}
FormalPowerSeries(vc<mint> &&f) : vc<mint>(std::move(f)) {}
int sz() const { return size(); }
void shrink()
{
while (!empty() && back() == 0)
pop_back();
}
mint get(int i) const { return 0 <= i && i < sz() ? (*this)[i] : 0; }
F resized(int n) const
{
assert(n >= 0);
F res(n);
copy_n(begin(), min(sz(), n), res.begin());
return res;
}
F rev(int d = -1) const
{
F res(*this);
if (d >= 0)
res.resize(d);
reverse(ALL(res));
return res;
}
int cnt_nz() const ;
optional<pair<int, mint>> nz_front() const
;
vc<pair<int, mint>> nz() const
;
mint eval(const mint &x) const
;
F operator-() const
{
F res(*this);
fem(a : res) a = -a;
return res;
}
F &operator*=(const mint &k)
{
fem(a : *this) a *= k;
return *this;
}
F operator*(const mint &k) const { F res(*this); res *= k; return res; }
friend F operator*(const mint &k, const F &f) { return f * k; }
F &operator/=(const mint &k)
;
F operator/(const mint &k) const { F res(*this); res /= k; return res; }
F &operator+=(const F &g)
;
F operator+(const F &g) const ;
F &operator-=(const F &g)
{
const int n = size(), m = g.size();
resize(max(n, m));
repi(i, m)(*this)[i] -= g[i];
return *this;
}
F operator-(const F &g) const { F res(*this); res -= g; return res; }
F &operator*=(const F &g) ;
F operator*(const F &g) const { return convolution(*this, g); }
F div_sparse_destructive(const F &g, int n)
{
assert(n >= 0);
assert(g.get(0) != 0);
const mint iv = g.front().inv();
vc<pair<int, mint>> terms;
repi(j, 1, min(n, g.sz())) if (g[j] != 0) terms.eb(j, g[j]);
resize(n);
if constexpr (internal::ordinary_mod32<mint>::value)
{
const auto mod = mint::mod();
const int block = mod <= (1U << 30) ? 16 : 4;
for (auto &[j, b] : terms) b *= -iv;
const int first = terms.empty() ? n : terms.front().first;
repi(i, first) (*this)[i] *= iv;
size_t active = 0;
repi(i, first, n)
{
while (active < terms.size() && terms[active].first <= i) ++active;
size_t k = 0;
const uint seed = ((*this)[i] * iv).val();
const auto sum = internal::block_dot32(active, mod, block, seed, [&]()
{
const auto &[j, b] = terms[k++];
return ull((*this)[i - j].val()) * b.val();
});
(*this)[i] = mint::raw(sum);
}
}
else
repi(i, n)
{
fec([j, b] : terms)
{
if (j > i)
break;
(*this)[i] -= (*this)[i - j] * b;
}
(*this)[i] *= iv;
}
return *this;
}
F div_sparse(const F &g, int n) const { return F(*this).div_sparse_destructive(g, n); }
F inv(int n) const
{
assert(n >= 0);
assert(get(0) != 0);
if (n == 0)
return {};
mint iv = get(0).inv();
if (auto tail = internal::fps_tail_power(*this, n, iv, -iv * iv, iv * iv * iv)) return std::move(*tail);
if (internal::fps_use_sparse<internal::FPSSparseOperation::inv>(*this, n))
return F{1}.div_sparse(*this, n);
const int N = bit_ceil(n);
if (n >= 512 && ll(n) * 4 <= ll(N) * 3 && ntt_ok<mint>(N))
return internal::fps_inv_block(*this, n, 8);
return internal::fps_inv_newton(*this, n);
}
F div(const F &g, int n) const
;
F div_poly(const F &g) const
{
assert(!g.empty() && g.back() != 0);
const int k = sz() - g.sz() + 1;
if (k <= 0)
return {};
F a(k), b(min(k, g.sz()));
copy_n(this->rbegin(), k, a.begin());
copy_n(g.rbegin(), b.sz(), b.begin());
return (a * b.inv(k)).resized(k).rev();
}
pair<F, F> divmod(const F &g) const
{
F q = div_poly(g);
const int l = sz() - q.sz();
F r = resized(l) - (q.resized(min(l, q.sz())) * g.resized(min(l, g.sz()))).resized(l);
r.shrink();
return {q, r};
}
F operator%(const F &g) const ;
F &operator%=(const F &g) ;
F circular_mod(int n) const
;
F operator<<(int k) const
;
F operator>>(int k) const
{
F res(max(0, sz() - k));
repi(i, sz() - k) res[i] = (*this)[i + k];
return res;
}
F &operator<<=(int k) ;
F &operator>>=(int k) { return *this = *this >> k; }
F diff() const
;
F integ() const
;
F log(int n) const
{
assert(n >= 0);
assert(get(0) == 1);
if (n == 0)
return {};
if (internal::fps_use_sparse<internal::FPSSparseOperation::div>(*this, n - 1))
return internal::sparse_fps::logarithm<
internal::ordinary_mod32<mint>::value>(*this, n);
F f = resized(n);
return f.diff().div(f, n - 1).integ();
}
static F diff_eq(const F &a, const F &b, int d)
;
#line 6812
F exp_sparse(int n) const
;
F pow_sparse(ll k, int n) const
;
#line 6851
optional<F> sqrt_sparse(int n) const
;
#line 6874
F exp(int n) const
{
assert(n >= 0);
assert(get(0) == 0);
if (n == 0)
return {};
int first = 1;
while (first < min(n, sz()) && (*this)[first] == 0) ++first;
if (first >= min(n, sz())) return F{1}.resized(n);
if (ll(first) * 2 >= n)
{
F res = resized(n);
res[0] = 1;
return res;
}
if (internal::fps_use_sparse<internal::FPSSparseOperation::exp>(*this, n)) return exp_sparse(n);
if (ntt_ok<mint>(bit_ceil(n)))
{
if (n < 512) return internal::fps_exp_bs(*this, n);
return internal::fps_exp_block(*this, n, n <= 4096 ? 4 : 16);
}
Binomial<mint>::reserve(n - 1);
F f{1, get(1)}, g{1};
for (int m = 2; m < n; m *= 2)
{
F fg2 = f * (g * g);
g.resize(m);
repi(i, m / 2, m) g[i] = -fg2.get(i);
F q(m - 1);
repi(i, 1, min(m, sz())) q[i - 1] = mint(i) * (*this)[i];
F r = (f * q).circular_mod(m), s(m);
repi(i, m)
s[(i + 1) % m] = (i + 1 < m ? mint(i + 1) * f[i + 1] : mint(0)) - r[i];
F t = (g * s).resized(m), u(min(m, n - m));
repi(i, u.sz()) u[i] = get(m + i) - t[i] * Binomial<mint>::inv_[m + i];
F v = (f * u).resized(u.sz());
f.insert(f.end(), v.begin(), v.end());
}
return f.resized(n);
}
F pow(ll k, int n) const
;
#line 6957
optional<F> sqrt(int n) const
{
assert(n >= 0);
const auto first = nz_front();
if (!first) return F(n);
const auto [d0, a0] = *first;
if (d0 % 2 != 0)
return nullopt;
const auto root = sqrt_mod(a0);
if (!root) return nullopt;
const mint r = *root;
if (d0 / 2 >= n)
return F(n);
const int need = n - d0 / 2;
F shifted(need);
copy_n(begin() + d0, min(need, sz() - d0), shifted.begin());
mint iv = a0.inv(), half = Binomial<mint>::inv(2);
if (auto tail = internal::fps_tail_power(shifted, need, r, r * iv * half, -r * iv * iv * half * half * half))
return std::move(*tail) << (d0 / 2);
if (internal::fps_use_sparse<internal::FPSSparseOperation::sqrt>(*this, need, d0))
return sqrt_sparse(n);
F f = std::move(shifted); f *= iv;
F g = need >= 512 && ntt_ok<mint>(bit_ceil(need))
? internal::fps_sqrt_block(f, need, 16)
: internal::fps_sqrt_newton(f, need);
return (g * r) << (d0 / 2);
}
F pow_mod(ll k, const F &g) const
;
F egf() const
;
F ogf() const
;
F mul_bin_destructive(int d, mint c)
;
F mul_bin(int d, mint c) const ;
F div_bin_destructive(int k, mint c, int n)
;
F div_bin(int k, mint c, int n) const ;
};
template <class mint>
mint bostan_mori(const FormalPowerSeries<mint> &p, const FormalPowerSeries<mint> &q, ll k)
{
using F = FormalPowerSeries<mint>;
if (k < 0)
return 0;
static const internal::fft_info<mint> info;
auto [r, u] = p.divmod(q);
mint res = k < r.sz() ? r[(int)k] : 0;
const int d = SZ(q) - 1;
if (d == 0)
return res;
auto finish = [&](const F &a, const F &b) -> mint
{
if (k == 0)
return res + a.get(0) / b[0];
if (k == 1)
{
const mint ib0 = b[0].inv();
return res + (a.get(1) - a.get(0) * b.get(1) * ib0) * ib0;
}
const int n = (int)k + 1;
return res + convolution_point_get<mint>(a, b.resized(n).inv(n), (int)k);
};
if (k <= d)
return finish(u, q);
const int cutoff = d < 64 ? 0 : d;
if (ntt_ok<mint>(2 * d + 1))
{
const int z = bit_ceil(2 * d + 1);
const mint i2 = mint(2).inv();
const mint iz = mint(z / 2).inv();
const mint root = info.root[bit_width(z / 2)];
const mint iroot = info.iroot[bit_width(z / 2)];
vc<mint> ipw(z / 2);
{
mint itmp = 1;
repi(i, z / 2)
{
ipw[bitrev(z / 2, i)] = itmp;
itmp *= iroot;
}
}
F v = q;
u.resize(z / 2), v.resize(z / 2);
F u2 = u, v2 = v;
ntt(u), ntt(v);
while (k > cutoff)
{
{
mint tmp = 1;
repi(i, z / 2)
{
u2[i] = u2[i] * tmp;
v2[i] = v2[i] * tmp;
tmp *= root;
}
}
ntt(u2), ntt(v2);
if (k & 1)
{
repi(i, z / 4)
{
const mint x = v[2 * i], y = v[2 * i + 1];
u[i] = (u[2 * i] * y - u[2 * i + 1] * x) * ipw[i] * i2;
v[i] = x * y;
}
repi(i, z / 4)
{
const mint x = v2[2 * i], y = v2[2 * i + 1];
u[i + z / 4] = (u2[2 * i] * y - u2[2 * i + 1] * x) * ipw[i + z / 4] * i2;
v[i + z / 4] = x * y;
}
}
else
{
repi(i, z / 4)
{
const mint x = v[2 * i], y = v[2 * i + 1];
u[i] = (u[2 * i] * y + u[2 * i + 1] * x) * i2;
v[i] = x * y;
}
repi(i, z / 4)
{
const mint x = v2[2 * i], y = v2[2 * i + 1];
u[i + z / 4] = (u2[2 * i] * y + u2[2 * i + 1] * x) * i2;
v[i + z / 4] = x * y;
}
}
u2 = u, v2 = v;
intt(u2), intt(v2);
repi(i, z / 2) u2[i] *= iz, v2[i] *= iz;
k >>= 1;
}
return finish(u2, v2);
}
else
{
F v = q;
u.resize(d + 1), v.resize(d + 1);
while (k > cutoff)
{
F w = v;
repi(i, 1, d + 1, 2) w[i] = -w[i];
F u2 = u * w, v2 = v * w;
repi(i, d + 1)
{
if (2 * i + (k & 1) < SZ(u2))
u[i] = u2[2 * i + (k & 1)];
if (2 * i < SZ(v2))
v[i] = v2[2 * i];
}
k >>= 1;
}
return finish(u, v);
}
}
template <class mint>
pair<FormalPowerSeries<mint>, FormalPowerSeries<mint>> linear_recurrence_gf(const vc<mint> &a, const vc<mint> &c)
{
using F = FormalPowerSeries<mint>;
const int d = SZ(c) - 1;
assert(d >= 0);
assert(SZ(a) >= d);
F q = -F(c);
q[0] = 1;
F p = (F(a) * q).resized(d);
q.shrink();
return {p, q};
}
template <class mint>
mint bmbm(const vc<mint> &a, ll k, bool show_coefs = true)
;
using fps = FormalPowerSeries<mint>;
void main2()
{
// LL(N, M);
// const ll MAX_N = 5;
// fps f(M + 1);
// rep(i, 1, M + 1) f[i] = 1;
// vc<fps> pw(MAX_N);
// pw[0] = {1};
// rep(i, 1, MAX_N) pw[i] = pw[i - 1] * f;
// // dump(pw | cp::index());
// vc<mint> A(MAX_N);
// rep(n, MAX_N)
// {
// rep(k, n + 1) rep(s, M * n) A[n] += pw[k].get(s) * pw[n - k].get(s);
// }
// // dump(A | cp::index());
// dump(A | cp::index());
// PRINT(bmbm(A, N));
LL(N, M);
fps A = {0, 0, M, M * (M - 1)};
fps C = {-1, M - 1, M, 0};
auto [P, Q] = linear_recurrence_gf(A, C);
PRINT(bostan_mori(P, Q, N));
}
void test()
{
}
template <auto init, auto main2, auto test>
struct Main
{
Main()
{
cauto CERR = [](string val, string color)
{
string s = "\033[" + color + "m" + val + "\033[m";
};
CERR("\n[FAST_IO]\n\n", "32");
cout << fixed << setprecision(20);
init();
#line 7276
CERR("\n[SINGLE_TESTCASE]\n\n", "36");
main2();
#line 7291
}
};
Main<init, main2, test> main_dummy;
int main() {}
// submit-code: hash=130e2dad5b02 source=130e2dad5b02394d7ed184f637af65f292ac4efa22a29e7e2dc63b09ba64117c input=3d9d376e9eeb75b979e7a62bff962a340c0bb02b2c20af82345e0ebd2db9a3bb body=c8477bd8d7cca80c729556841df39b832ac40d70e8afc1268c4f1ef40c798279