// #define MULTI_TESTCASE // #define AOJ_TESTCASE // #define FAST_CIO // #define INTERACTIVE #define INF 4'000'000'000'000'000'037LL #define EPS 1e-11 // https://github.com/miscalculation53/library/tree/wip/template/template_all.hpp // https://github.com/miscalculation53/library/tree/wip/template/template_all_but_modint.hpp // https://github.com/miscalculation53/library/tree/wip/template/template_types.hpp #include using namespace std; // https://github.com/miscalculation53/library/tree/wip/template/template_float128.hpp using namespace std; #if __cplusplus > 202002L && defined(__STDCPP_FLOAT128_T__) #include 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 #if __cplusplus > 202002L && defined(__STDCPP_FLOAT128_T__) && \ (defined(_GLIBCXX_HAVE_FLOAT128_MATH) || defined(_GLIBCXX_LDOUBLE_IS_IEEE_BINARY128)) #define LIBRARY_USE_QUADMATH 0 #else #include #define LIBRARY_USE_QUADMATH 1 //! https://gcc.gnu.org/onlinedocs/libquadmath/Math-Library-Routines.html #define _p_F128_UNARY(name) \ inline f128 name(f128 x) { return f128(name##q(__float128(x))); } _p_F128_UNARY(acos) _p_F128_UNARY(acosh) _p_F128_UNARY(asin) _p_F128_UNARY(asinh) _p_F128_UNARY(atan) _p_F128_UNARY(atanh) _p_F128_UNARY(cbrt) _p_F128_UNARY(ceil) _p_F128_UNARY(cos) _p_F128_UNARY(cosh) _p_F128_UNARY(erf) _p_F128_UNARY(erfc) _p_F128_UNARY(exp) _p_F128_UNARY(exp2) _p_F128_UNARY(expm1) _p_F128_UNARY(fabs) _p_F128_UNARY(floor) _p_F128_UNARY(lgamma) _p_F128_UNARY(log) _p_F128_UNARY(log10) _p_F128_UNARY(log1p) _p_F128_UNARY(log2) _p_F128_UNARY(logb) _p_F128_UNARY(nearbyint) _p_F128_UNARY(rint) _p_F128_UNARY(round) _p_F128_UNARY(sin) _p_F128_UNARY(sinh) _p_F128_UNARY(sqrt) _p_F128_UNARY(tan) _p_F128_UNARY(tanh) _p_F128_UNARY(tgamma) _p_F128_UNARY(trunc) #undef _p_F128_UNARY #define _p_F128_BINARY(name) \ inline f128 name(f128 x, f128 y) { return f128(name##q(__float128(x), __float128(y))); } \ template , int> = 0> \ inline f128 name(f128 x, T y) { return name(x, f128(y)); } \ template , int> = 0> \ inline f128 name(T x, f128 y) { return name(f128(x), y); } _p_F128_BINARY(atan2) _p_F128_BINARY(copysign) _p_F128_BINARY(fdim) _p_F128_BINARY(fmax) _p_F128_BINARY(fmin) _p_F128_BINARY(fmod) _p_F128_BINARY(hypot) _p_F128_BINARY(nextafter) _p_F128_BINARY(pow) _p_F128_BINARY(remainder) #undef _p_F128_BINARY inline f128 fma(f128 x, f128 y, f128 z) { return f128(fmaq(__float128(x), __float128(y), __float128(z))); } inline f128 hypot(f128 x, f128 y, f128 z) { return hypot(hypot(x, y), z); } namespace template_float128_detail { template constexpr bool arithmetic = is_arithmetic_v || is_same_v; } #define _p_F128_TERNARY_MIXED(name) \ template && template_float128_detail::arithmetic, int> = 0> \ inline f128 name(f128 x, T y, U z) { return name(x, f128(y), f128(z)); } \ template && !is_same_v && template_float128_detail::arithmetic, int> = 0> \ inline f128 name(T x, f128 y, U z) { return name(f128(x), y, f128(z)); } \ template && !is_same_v && template_float128_detail::arithmetic && !is_same_v, int> = 0> \ inline f128 name(T x, U y, f128 z) { return name(f128(x), f128(y), z); } inline f128 ldexp(f128 x, int exponent) { return f128(ldexpq(__float128(x), exponent)); } inline f128 scalbn(f128 x, int exponent) { return f128(scalbnq(__float128(x), exponent)); } inline f128 scalbln(f128 x, long exponent) { return f128(scalblnq(__float128(x), exponent)); } inline f128 modf(f128 x, f128 *integer) { __float128 value; const auto fraction = modfq(__float128(x), &value); *integer = f128(value); return f128(fraction); } inline f128 remquo(f128 x, f128 y, int *quotient) { return f128(remquoq(__float128(x), __float128(y), quotient)); } inline long lrint(f128 x) { return lrintq(__float128(x)); } inline long lround(f128 x) { return lroundq(__float128(x)); } inline long long llrint(f128 x) { return llrintq(__float128(x)); } inline long long llround(f128 x) { return llroundq(__float128(x)); } inline int ilogb(f128 x) { return ilogbq(__float128(x)); } #endif inline f128 nexttoward(f128 x, long double y) { return nextafter(x, f128(y)); } namespace template_float128_detail { inline bool parse(string_view s, f128 &x) { if (!s.empty() && s.front() == '+') { s.remove_prefix(1); if (s.empty() || s.front() == '-') return false; } #if LIBRARY_USE_QUADMATH string token(s); char *end; errno = 0; const __float128 value = strtoflt128(token.c_str(), &end); if (end == token.data() || end != token.data() + token.size()) return false; if (errno == ERANGE && (value == 0 || isinfq(value))) return false; x = f128(value); return true; #else auto [p, ec] = from_chars(s.data(), s.data() + s.size(), x); return ec == errc{} && p == s.data() + s.size(); #endif } inline string format(f128 x, chars_format style, int precision) { string s(128, '\0'); while (true) { #if LIBRARY_USE_QUADMATH const char *spec = style == chars_format::fixed ? "%.*Qf" : style == chars_format::scientific ? "%.*Qe" : "%.*Qg"; const int n = quadmath_snprintf(s.data(), s.size(), spec, precision, __float128(x)); if (n < 0) abort(); if (size_t(n) < s.size()) { s.resize(size_t(n)); return s; } s.resize(size_t(n) + 1); #else auto [p, ec] = to_chars(s.data(), s.data() + s.size(), x, style, precision); if (ec == errc{}) { s.resize(size_t(p - s.data())); return s; } s.resize(s.size() * 2); #endif } } } using ld = decltype(EPS); using ll = long long; using uint = unsigned int; using ull = unsigned long long; using pll = pair; using tlll = tuple; using tllll = tuple; #define vc vector template using vvc = vc>; template using vvvc = vc>>; using vb = vc; using vl = vc; using vpll = vc; using vtlll = vc; using vtllll = vc; using vstr = vc; using vvb = vvc; using vvl = vvc; template using pql = priority_queue, greater>; template using pqg = priority_queue; using i128 = __int128_t; using u128 = __uint128_t; 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); } string i128tos(i128 x) { if (x == 0) return "0"; string sign = "", res = ""; u128 ux; if (x < 0) ux = u128(-(x + 1)) + 1, sign = "-"; else ux = u128(x); while (ux > 0) { res += char('0' + ux % 10); ux /= 10; } reverse(res.begin(), res.end()); return sign + res; } 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; } #define cauto const auto // https://github.com/miscalculation53/library/tree/wip/template/template_rep.hpp #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(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__) // https://github.com/miscalculation53/library/tree/wip/template/template_math.hpp // https://github.com/miscalculation53/library/tree/wip/template/template_arithmetic.hpp template constexpr pair &operator+=(pair &a, const P &b); template constexpr pair operator+(pair a, const P &b); template constexpr pair &operator-=(pair &a, const P &b); template constexpr pair operator-(pair a, const P &b); template constexpr pair operator-(pair a); template constexpr array &operator+=(array &a, const A &b); template constexpr array operator+(array a, const A &b); template constexpr array &operator-=(array &a, const A &b); template constexpr array operator-(array a, const A &b); template constexpr array operator-(array a); template constexpr tuple &operator+=(tuple &a, const Tp &b); template constexpr tuple operator+(tuple a, const Tp &b); template constexpr tuple &operator-=(tuple &a, const Tp &b); template constexpr tuple operator-(tuple a, const Tp &b); template constexpr tuple operator-(tuple a); namespace internal { }; // https://github.com/miscalculation53/library/tree/wip/utils/is_integral_ext.hpp template constexpr bool is_integral_ext = is_integral_v || is_same_v || is_same_v; template constexpr bool is_signed_ext = is_signed_v || is_same_v; template constexpr bool is_unsigned_ext = is_unsigned_v || is_same_v; // https://github.com/miscalculation53/library/tree/wip/utils/integer_arithmetic.hpp // https://github.com/miscalculation53/library/tree/wip/utils/make_unsigned_ext.hpp template struct make_unsigned_ext { using type = make_unsigned_t; }; template <> struct make_unsigned_ext { using type = u128; }; template <> struct make_unsigned_ext { using type = u128; }; template using make_unsigned_ext_t = typename make_unsigned_ext::type; template struct make_signed_ext { using type = make_signed_t; }; template <> struct make_signed_ext { using type = i128; }; template using make_signed_ext_t = typename make_signed_ext::type; template constexpr bool is_negative(T x) { static_assert(is_integral_ext); if constexpr (is_signed_ext) return x < 0; return false; } template > struct unsigned_arithmetic { using type = T; }; template struct unsigned_arithmetic { static_assert(!is_same_v); using type = common_type_t>; }; template using unsigned_arithmetic_t = typename unsigned_arithmetic::type; // https://github.com/miscalculation53/library/tree/wip/utils/default_infty.hpp template constexpr decltype(auto) default_infty(); namespace default_infty_detail { template struct is_product : false_type {}; template struct is_product> : true_type {}; template struct is_product> : true_type {}; template struct is_product> : true_type {}; template struct has_infty : false_type {}; template struct has_infty> : true_type {}; template inline constexpr bool unsupported = false; } namespace template_math_detail { } template && is_integral_ext>> inline constexpr T divceil(U a, V b) { return T(a) / T(b) + (T(a) % T(b) && (T(a) ^ T(b)) >= 0); } template && is_integral_ext>> inline constexpr T safemod(U a, V b) { const T x = T(a), y = T(b); if constexpr (is_signed_ext) 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 constexpr T ipow(U a, V b) { assert(b >= 0); if (b == 0) return 1; if (a == 0 || a == 1) return a; if (a < 0 && a == -1) return b & 1 ? -1 : 1; T res = 1, tmp = a; while (true) { if (b & 1) res *= tmp; b >>= 1; if (b == 0) break; tmp *= tmp; } return res; } //! https://misawa.github.io/others/avoid_errors/techniques_to_avoid_errors.html // https://github.com/miscalculation53/library/tree/wip/template/template_vector.hpp #define ALL(a) (a).begin(), (a).end() template inline T SZ(const V &x) { return x.size(); } #define eb emplace_back #define LMD(x,fx) ([&](const auto &x) { return fx; }) #define GEN_VEC(n,i,fi) (gen_vec(n, LMD(i, fi))) namespace template_vector_detail { template struct nested_vector { using type = vc::type>; }; template struct nested_vector { using type = T; }; } // https://github.com/miscalculation53/library/tree/wip/template/template_algo.hpp // https://github.com/miscalculation53/library/tree/wip/utils/resolved_infty.hpp // https://github.com/miscalculation53/library/tree/wip/utils/resolved_value.hpp template vvc top(const vvc &a) { if (a.empty()) return {}; const int n = a.size(), m = a[0].size(); vvc b(m, vc(n)); repi(i, n) { assert(SZ(a[i]) == m); repi(j, m) b[j][i] = a[i][j]; } return b; } vstr top(const vstr &a) { vvc a_(a.size()); repi(i, SZ(a)) a_[i] = {ALL(a[i])}; vvc b_ = top(a_); vstr b(b_.size()); repi(i, SZ(b)) b[i] = {ALL(b_[i])}; return b; } template struct has_e0 : false_type {}; template struct has_e0> : true_type {}; template inline constexpr bool has_e0_v = has_e0::value; template struct MonoidAdd { using S = T; static constexpr S e(); }; template struct MonoidMin { using S = T; static constexpr decltype(auto) e(); }; template struct MonoidMax { using S = T; static constexpr S e(); }; namespace internal { template struct HasMonoidPow : false_type { }; template struct HasMonoidPow(), declval()))>> : true_type { }; } constexpr array DRULgrid = {{{1, 0}, {0, 1}, {-1, 0}, {0, -1}}}; constexpr array DRULplane = {{{0, -1}, {1, 0}, {0, 1}, {-1, 0}}}; // https://github.com/miscalculation53/library/tree/wip/template/template_binsearch.hpp template struct is_random_access_iterator { static constexpr bool value = is_same_v< typename iterator_traits::iterator_category, random_access_iterator_tag >; }; template constexpr bool is_random_access_iterator_v = is_random_access_iterator::value; #define DEFAULT_COMP ranges::less namespace internal { }; // https://github.com/miscalculation53/library/tree/wip/template/template_bit.hpp template inline constexpr ull MASK(T k) { return (1ULL << k) - 1ULL; } 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; } 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; } // https://github.com/miscalculation53/library/tree/wip/template/template_inout.hpp // https://github.com/miscalculation53/library/tree/wip/template/template_dump.hpp #if defined(LOCAL) && LIBRARY_HAS_FLOAT128 namespace cpp_dump::_detail { struct export_command; inline string export_var(const f128 &, const string &, size_t, size_t, bool, const export_command &); } #endif // https://github.com/miscalculation53/library/tree/wip/template/template_dump_map.hpp #define CPP_DUMP_DEFINE_DATA(...) #define dump(...) ((void)0) #define local(...) #define oj(...) __VA_ARGS__ #define local_oj(a,b) (b) #if LIBRARY_HAS_FLOAT128 inline istream &operator>>(istream &is, f128 &x) { string s; if (is >> s) { f128 value; if (template_float128_detail::parse(s, value)) x = value; else is.setstate(ios::failbit); } return is; } inline ostream &operator<<(ostream &os, f128 x) { const auto field = os.flags() & ios::floatfield; const auto format = field == ios::fixed ? chars_format::fixed : field == ios::scientific ? chars_format::scientific : chars_format::general; string s = template_float128_detail::format(x, format, int(os.precision())); if ((os.flags() & ios::showpos) && !signbit(x)) s.insert(s.begin(), '+'); if (os.flags() & ios::uppercase) for (char &c : s) c = char(toupper(static_cast(c))); return os << s; } #endif 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 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 __attribute__((always_inline)) inline void rd1_integer(T &x) { using U = integer_magnitude; 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::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)) { const unsigned digits = unsigned(__builtin_ctzll(stop)) / 8; if (digits) { uint tail; if (digits <= 3) { tail = unsigned(p[0] - '0'); if (digits >= 2) tail = tail * 10 + unsigned(p[1] - '0'); if (digits == 3) tail = tail * 10 + unsigned(p[2] - '0'); } else tail = decode_eight(columns << (8 * (8 - digits))); value = value * powers_of_ten[digits] + tail; p += digits; } 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::is_signed) { if (negative && value) x = T(-T(value - 1) - 1); else x = T(value); } else x = T(value); } template ::is_integer, int> = 0> __attribute__((always_inline)) inline void rd1(T &x) { rd1_integer(x); } inline void rd1(char &x) { find_token(); x = ibuf[pil++]; } inline void rd1(string &x) { find_token(); x.clear(); for (;;) { const char *p = ibuf + pil; const char *end = ibuf + pir; const char *first = p; while (p != end && !ascii_space(*p)) ++p; x.append(first, size_t(p - first)); pil = uint32_t(p - ibuf); if (p != end) return; load(); if (!pir) return; } } template 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) x = strtof(token.c_str(), &end); else if constexpr (is_same_v) 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) { rd1_real(x); } inline void rd1(double &x) { rd1_real(x); } inline void rd1(long double &x) { rd1_real(x); } #if LIBRARY_HAS_FLOAT128 inline void rd1(f128 &x) { string token; rd1(token); if (!template_float128_detail::parse(token, x)) abort(); } #endif template void rd1(pair &x); template void rd1(tuple &x); template void rd1(array &x); template void rd1(vc &x); inline constexpr auto decimal_triplets = [] { array, 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 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 __attribute__((always_inline)) inline void wt1_integer(T x) { using U = integer_magnitude; U magnitude = U(x); bool negative = false; if constexpr (numeric_limits::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; } char text[48]; char *const end = text + sizeof(text); char *first = end; if constexpr (sizeof(U) > sizeof(ull)) { constexpr ull base = 1000000000000000000ULL; while (magnitude >= base) { const U q = magnitude / base; const ull tail = ull(magnitude - q * base); const uint hi = uint(tail / 1000000000ULL); first -= 18; store_nine(first, hi); store_nine(first + 9, uint(tail - ull(hi) * 1000000000ULL)); magnitude = q; } first = prepend_wide(first, ull(magnitude)); } else if constexpr (sizeof(U) > sizeof(uint)) first = prepend_wide(first, magnitude); else first = prepend_small(first, magnitude); if (negative) *--first = '-'; append_bytes(first, size_t(end - first)); if constexpr (newline) append_char('\n'); } inline void wt1(char x) { append_char(x); } inline void wt1(const char *x) { append_bytes(x, strlen(x)); } inline void wt1(string_view x) { append_bytes(x.data(), x.size()); } inline void wt1(const string &x) { append_bytes(x.data(), x.size()); } template void wt1_real(T x) { char text[numeric_limits::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) { wt1_real(x); } inline void wt1(double x) { wt1_real(x); } inline void wt1(long double x) { wt1_real(x); } #if LIBRARY_HAS_FLOAT128 inline void wt1(f128 x) { wt1(template_float128_detail::format(x, chars_format::fixed, 35)); } #endif template void wt1(const pair &x); template void wt1(const tuple &x); template void wt1(const array &x); template void wt1(const vc &x); template void read(T &...x) { (rd1(x), ...); } template void write(T &&...x) { (wt1(forward(x)), ...); } template void print(T &&...x) { if constexpr (sizeof...(T) == 1 && ((is_integral_v> && !is_same_v, char>) && ...)) (wt1_integer, true>(forward(x)), ...); else { size_t index = 0; ((index++ ? void(append_char(' ')) : void(), wt1(forward(x))), ...); append_char('\n'); } } } struct Dummy { Dummy() { atexit(fastio::flush); } } dummy; //! https://trap.jp/post/1224/ namespace internal { template void READnodump(Ts &...a) { fastio::read(a...); } }; #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 __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 __VA_ARGS__; READVEC(n, __VA_ARGS__) #define VEC2(T,n,m,...) vvc __VA_ARGS__; READVEC2(n, m, __VA_ARGS__) #define READJAG(...) internal::READJAGnodump(__VA_ARGS__); dump(__VA_ARGS__) #define JAG(T,n,...) vvc __VA_ARGS__; READJAG(n, __VA_ARGS__) #define ENDL '\n' #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) #define PRINTVEXIT(...) do { PRINTV(__VA_ARGS__); exit(0); } while (false) #define PRINTVRETURN(...) do { PRINTV(__VA_ARGS__); return; } while (false) namespace internal { }; namespace internal { }; #define UNZIP(vt,...) auto [__VA_ARGS__] = unzip(vt) #define ZIP(vt,...) auto vt = internal::zip_vectors(__VA_ARGS__) // https://github.com/miscalculation53/library/tree/wip/template/template_random.hpp mt19937_64 mt; template T randrange(U1 l, U2 r) { assert(T(l) < T(r)); return uniform_int_distribution(T(l), T(r) - 1)(mt); } bool randbool(double p) { assert(0 <= p && p <= 1); return bernoulli_distribution(p)(mt); } namespace internal { }; // https://github.com/miscalculation53/library/tree/wip/math/modint/template_modint.hpp // https://github.com/miscalculation53/library/tree/wip/math/modint/modint.hpp /*! * 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 */ // https://github.com/miscalculation53/library/tree/wip/math/modint/modint_internal_static.hpp // https://github.com/miscalculation53/library/tree/wip/utils/larger_int.hpp namespace larger_int_detail { } template struct larger_int { private: static constexpr bool check(); static_assert(check()); public: using type = T; }; #define LARGER_INT(T,U) template <> struct larger_int { 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 struct Rational; template struct larger_int> { using type = Rational::type>; }; template using larger_int_t = typename larger_int::type; // https://github.com/miscalculation53/library/tree/wip/math/modint/modint_internal_isprime.hpp /*! * 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 constexpr ll powmod_constexpr(ll x, ll n, T m) { if (m == 1) return 0; using U = make_unsigned_t; using L = larger_int_t; 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 constexpr bool isprime_constexpr(T n) { if constexpr (sizeof(T) > 4) { if (n <= INT_MAX) return isprime_constexpr(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; using L = larger_int_t; 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 constexpr bool isprime = isprime_constexpr(n); }; namespace internal { template struct policy_static { using mod_type = decltype(M); using value_type = make_unsigned_t; using calc_type = conditional_t<(M <= UINT_MAX), ull, larger_int_t>; 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); } }; }; // https://github.com/miscalculation53/library/tree/wip/math/modint/modint_internal_barrett32.hpp /*! * 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 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() { return reducer.umod(); } static value_type init(value_type v) { return v; } static mod_type val(value_type v) { return v; } static value_type mul(value_type a, value_type b) { return reducer.mul(a, b); } }; }; // https://github.com/miscalculation53/library/tree/wip/math/modint/modint_internal_montgomery64.hpp 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 { return m; } ull reduce(u128 x) const { if (b == 0) { auto t = (x + u128(mx) * (-imx * ull(x))) >> 64; if (t >= m) t -= m; return (ull)t; } ull p = x & MASK(b); x = (x >> b) + p * d; ull y = p << (64 - b); auto t = (x + u128(mx) * (imx * (y - ull(x)))) >> (64 - b); if (t >= m) { t -= m; if (t >= m) t -= m; } return (ull)t; } ull inv_reduce(i128 v) const { return reduce(u128(v % m + m) * q); } }; template 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() { return reducer.umod(); } static value_type init(value_type v) { return reducer.inv_reduce(v); } static mod_type val(value_type v) { return reducer.reduce(v); } static value_type mul(value_type a, value_type b) { return reducer.reduce((calc_type)a * b); } }; template 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 mod_type mod() { return reducer.umod(); } static value_type init(value_type v) { return reducer.inv_reduce(v); } }; }; // https://github.com/miscalculation53/library/tree/wip/math/extgcd.hpp namespace internal { template 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 static auto set_mod(M m) -> decltype(T::set_mod(m)) { return T::set_mod(m); } modint_impl() : _v(0) { (void)mod(); } template >> modint_impl(T v) { const V m = mod(); V rem; if constexpr (is_signed_ext) { using S = make_signed_t; 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; } template mint pow(T n) const { assert(n >= 0); if constexpr (!is_void_v) 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) { return raw(PowInv::get().inv(val())); } else if constexpr (Policy::is_prime) { return pow(mod() - 2); } else { auto [g, x, y] = extgcd(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(); } friend void rd1(mint &x) { long long a; fastio::rd1(a); x = a; } friend void wt1(const mint &x) { fastio::wt1(x.val()); } }; }; template using static_modint32 = internal::modint_impl>; template using dynamic_modint32 = internal::modint_impl>; template using static_modint64 = internal::modint_impl>; template using dynamic_modint64_odd = internal::modint_impl>; template using dynamic_modint64 = internal::modint_impl>; 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 struct is_modint : std::false_type { }; template struct is_modint> : std::true_type { }; template inline constexpr bool is_modint_v = is_modint::value; template struct is_static_modint : false_type {}; template struct is_static_modint, PowInv>> : true_type {}; template struct is_static_modint, PowInv>> : true_type {}; template inline constexpr bool is_static_modint_v = is_static_modint::value; template struct is_dynamic_modint : false_type {}; template struct is_dynamic_modint> : true_type {}; template struct is_dynamic_modint> : true_type {}; template struct is_dynamic_modint> : true_type {}; template inline constexpr bool is_dynamic_modint_v = is_dynamic_modint::value; template struct has_mod : std::false_type { }; template struct has_mod> : std::true_type { }; template struct modint_less { }; template struct modint_hash { }; // https://github.com/miscalculation53/library/tree/wip/math/modint/power_table.hpp template struct PowerTable { private: decltype(mint::mod()) mod; mint base; vc pw; public: void reserve(int n); }; // https://github.com/miscalculation53/library/tree/wip/math/modint/binomial.hpp template struct Binomial { private: inline static decltype(T::mod()) mod; public: inline static vc fac_, finv_, inv_; static void reserve(int n); static T inv(int n) { n %= T::mod(); if (n < 0) n += T::mod(); assert(n != 0); reserve(n); return inv_[n]; } }; // https://github.com/miscalculation53/library/tree/wip/math/modint/stom.hpp // https://github.com/miscalculation53/library/tree/wip/math/modint/to_rational.hpp // https://github.com/miscalculation53/library/tree/wip/math/svp2d.hpp namespace cpp_dump { template 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 struct rat_object { T value; }; template inline constexpr bool rat_is_map = false; template inline constexpr bool rat_is_map> = true; template inline constexpr bool rat_is_map> = true; template inline constexpr bool rat_is_map> = true; template inline constexpr bool rat_is_map> = true; struct rat_closure : std::ranges::range_adaptor_closure { }; constexpr rat_closure rat() { return rat_closure{}; } } using mint998 = modint998244353; // using mint = modint1000000007; using mint103 = static_modint32<100003>; // using mint = modint; using bi = Binomial; void init() { oj(mt.seed(random_device()())); } unsigned xor128_x = 123456789, xor128_y = 362436069, xor128_z = 521288629, xor128_w = 88675123; unsigned xor128() { unsigned t = xor128_x ^ (xor128_x << 11); xor128_x = xor128_y; xor128_y = xor128_z; xor128_z = xor128_w; return xor128_w = xor128_w ^ (xor128_w >> 19) ^ (t ^ (t >> 8)); } // https://github.com/miscalculation53/library/tree/wip/math/prime/large/order_primitive_root.hpp // https://github.com/miscalculation53/library/tree/wip/math/prime/prime_power.hpp template struct PrimePower { P p; int e; P pe; PrimePower(P p, int e = 1) : p(p), e(e), pe(ipow(p, e)) {} PrimePower(P p, int e, P pe) : p(p), e(e), pe(pe) {} }; tuple ord_pow_div(ll n, ll m) { assert(m >= 2); if (m == 2) { int e = countr_zero(n); return {e, 1LL << e, n >> e}; } if (n % m != 0) return {0, 1, n}; n /= m; if (n % m != 0) return {1, m, n}; n /= m; ll m2 = m * m; auto [f, m2f, nn] = ord_pow_div(n, m2); int e = 2 + 2 * f; ll me = m2f * m2; if (nn % m == 0) e++, me *= m, nn /= m; return {e, me, nn}; } template vc> factorized_mul (const vc> &fac1, const vc> &fac2) { const int n = fac1.size(), m = fac2.size(); vc> fac; fac.reserve(n + m); int i = 0, j = 0; while (i < n && j < m) { if (fac1[i].p < fac2[j].p) fac.emplace_back(fac1[i++]); else if (fac1[i].p > fac2[j].p) fac.emplace_back(fac2[j++]); else { using U = larger_int_t

; fac.emplace_back(fac1[i].p, fac1[i].e + fac2[j].e, U(fac1[i].pe) * U(fac2[j].pe)); i++, j++; } } fac.insert(fac.end(), fac1.begin() + i, fac1.end()); fac.insert(fac.end(), fac2.begin() + j, fac2.end()); return fac; } // https://github.com/miscalculation53/library/tree/wip/math/prime/large/primality_test.hpp /*! * 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 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 bases32 = {2, 7, 61}; static constexpr array bases64 = {2, 325, 9375, 28178, 450775, 9780504, 1795265022}; if (n <= INT_MAX) { using mint = dynamic_modint32; return internal::is_prime_impl(n, bases32); } else { using mint = dynamic_modint64_odd; return internal::is_prime_impl(n, bases64); } } // https://github.com/miscalculation53/library/tree/wip/math/prime/large/factorize.hpp namespace internal { template ll get_prime_factor_impl(ll n) { mint::set_mod(n); int m = pow(n, .125); while (true) { int c = randrange(1, 100); mint x = 2, y = 2, prod = 1; ll g = 1; while (g == 1) { repi(i, m) { x = x * x + c; y = y * y + c, y = y * y + c; prod *= x - y; } g = gcd(prod.val(), n); } if (g == n) continue; if (is_prime(g)) return g; else if (is_prime(n / g)) return n / g; else return get_prime_factor_impl(g); } } ll get_prime_factor(ll n) { if (n <= INT_MAX) { using mint = dynamic_modint32; return get_prime_factor_impl(n); } else { using mint = dynamic_modint64_odd; return get_prime_factor_impl(n); } } }; vc> factorize(ll n) { vc> res; repi(p, 2, 100) { if (n % p == 0) { auto [e, pe, nn] = ord_pow_div(n, p); res.emplace_back(PrimePower(p, e, pe)); n = nn; } } while (n > 1) { if (is_prime(n)) { res.emplace_back(n); break; } ll p = internal::get_prime_factor(n); auto [e, pe, nn] = ord_pow_div(n, p); res.emplace_back(PrimePower(p, e, pe)); n = nn; } sort(ALL(res), [&](const PrimePower &pp1, const PrimePower &pp2) { return pp1.p < pp2.p; }); return res; } // https://github.com/miscalculation53/library/tree/wip/math/prime/large/euler_phi_carmichael.hpp //! https://37zigen.com/primitive-root/#i-4 //! https://x.com/noshi91/status/1317025404700405760 namespace internal { template vc lagrange_basis(const vc &a, mint x) { const int n = a.size(); vc res(n); if (n == 0) return res; auto dfs = [&](auto dfs, mint v, int l, int r) { if (r - l == 1) { res[l] = v; return; } int m = (l + r) / 2; mint v_l = v, v_r = v; repi(i, l, m) v_r = v_r.pow(a[i]); repi(i, m, r) v_l = v_l.pow(a[i]); dfs(dfs, v_l, l, m); dfs(dfs, v_r, m, r); }; dfs(dfs, x, 0, n); return res; } template bool internal_is_primitive_root(ll x, ll p, const vc> &fac) { const int k = fac.size(); mint::set_mod(p); vc a(k); repi(i, k) a[i] = fac[i].pe; auto b = lagrange_basis(a, mint(x)); repi(i, k) { if (b[i].pow(fac[i].pe / fac[i].p) == 1) return false; } return true; } template ll internal_primitive_root(ll p, const vc> &fac) { if (p == 2) return 1; while (true) { ll g = randrange(1, p); if (internal_is_primitive_root(g, p, fac)) return g; } } }; template ll primitive_root(ll p, const vc> &fac) { if (p <= INT_MAX) { using mint = dynamic_modint32; return internal::internal_primitive_root(p, fac); } else { using mint = dynamic_modint64_odd; return internal::internal_primitive_root(p, fac); } } // https://github.com/miscalculation53/library/tree/wip/math/convolution/convolution.hpp /*! * 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 */ // https://github.com/miscalculation53/library/tree/wip/math/crt.hpp /*! * 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 */ namespace crt_internal { namespace crt_word_detail { } namespace crt_word_dispatch { constexpr bool constant_evaluation() { return __builtin_is_constant_evaluated(); } } } namespace pre_crt_v2_adl { void gcd() = delete; template auto call(const T &a, const T &b) -> decltype(gcd(a, b)); } namespace pre_crt_v2_detail { template > struct work_type { using type = T; }; template struct work_type { using type = make_unsigned_ext_t; }; template struct has_adl_gcd : false_type {}; template struct has_adl_gcd(), declval()))>> : true_type {}; template constexpr bool unbounded_integer = numeric_limits::is_specialized && numeric_limits::is_integer && !numeric_limits::is_bounded; } // https://github.com/miscalculation53/library/tree/wip/math/dot_product.hpp // https://github.com/miscalculation53/library/tree/wip/algebra/algebra_basic_ops.hpp // https://github.com/miscalculation53/library/tree/wip/algebra/algebra_base.hpp // https://github.com/miscalculation53/library/tree/wip/algebra/internal_power.hpp //! SPDX-License-Identifier: CC0-1.0 namespace internal_monoid_power { template struct HasMemberScalarProduct : false_type {}; template struct HasMemberScalarProduct().operator*({declval()})))>> : true_type {}; template struct HasFreeScalarProduct : false_type {}; template struct HasFreeScalarProduct(), {declval()})))>> : true_type {}; template struct HasScalarProduct : bool_constant || is_same_v || HasMemberScalarProduct::value || HasFreeScalarProduct::value> {}; template using ConstRefPower = typename M::S (*)(const typename M::S &, I); template using ValuePower = typename M::S (*)(typename M::S, I); template struct HasConstRefPower : false_type {}; template struct HasConstRefPower>(M::pow))>> : true_type {}; template struct HasValuePower : false_type {}; template struct HasValuePower>(M::pow))>> : true_type {}; } template struct Monoid { using S = S_; static constexpr auto op = op_; static constexpr auto e = e_; }; template struct Group { using S = S_; static constexpr auto op = op_; static constexpr auto e = e_; static constexpr auto inv = inv_; }; template 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 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 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 struct HasInverse : false_type {}; template struct HasInverse()))>> : true_type {}; template struct OppositeMonoid { using S = typename M::S; static constexpr auto e = M::e; }; template struct OppositeGroup { using S = typename G::S; static constexpr auto e = G::e; static constexpr auto inv = G::inv; }; template struct NormalAndOppositeMonoid { struct S { typename M::S normal; typename M::S opposite; S(); template , std::nullptr_t> = nullptr> S(Args &&...args); S(const typename M::S &normal, const typename M::S &opposite); }; static constexpr S e(); }; template struct NormalAndOppositeGroup { struct S { typename G::S normal; typename G::S opposite; S(); template , std::nullptr_t> = nullptr> S(Args &&...args); S(const typename G::S &normal, const typename G::S &opposite); }; static constexpr S e(); static constexpr S inv(const S &a) { return {G::inv(a.normal), G::inv(a.opposite)}; } }; template using MonoidOfSemiRingAdd = Monoid; template using MonoidOfSemiRingMul = Monoid; template using GroupOfRingAdd = Group; template using GroupOfFieldMul = Group; template struct SemiRingFromMonoidMonoid { static_assert(is_same_v, "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 struct RingFromGroupMonoid { static_assert(is_same_v, "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 struct FieldFromGroupGroup { static_assert(is_same_v, "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 struct has_e1 : false_type {}; template struct has_e1> : true_type {}; template inline constexpr bool has_e1_v = has_e1::value; template struct MonoidMul { using S = T; static constexpr S e(); }; template struct GroupAddSub { using S = T; static constexpr S e(); static constexpr S inv(S a) { return -a; } }; template struct GroupMulDiv { using S = T; static constexpr S e(); static constexpr S inv(S a) { return e() / a; } }; template using SemiRingMinPlus = SemiRingFromMonoidMonoid, MonoidAdd>; template using SemiRingMaxPlus = SemiRingFromMonoidMonoid, MonoidAdd>; template using RingAddSubMul = RingFromGroupMonoid, MonoidMul>; template using FieldAddSubMulDiv = FieldFromGroupGroup, GroupMulDiv>; // https://github.com/miscalculation53/library/tree/wip/math/modint/internal_mod32.hpp namespace internal { template struct ordinary_mod32 : false_type {}; template struct ordinary_mod32, PowInv>> : true_type {}; template struct ordinary_mod32> : true_type {}; template struct ordinary_mod32_add_group : false_type {}; template struct ordinary_mod32_add_group> : ordinary_mod32 {}; template struct ordinary_mod32_add_group::op, GroupAddSub::e, GroupAddSub::inv>> : ordinary_mod32 {}; template struct ordinary_mod32_field : false_type {}; template struct ordinary_mod32_field> : ordinary_mod32 {}; } namespace internal { template using dot_product_mod32_value = ordinary_mod32; template struct dot_product_mod32 : false_type {}; template struct dot_product_mod32> : dot_product_mod32_value {}; template struct dot_product_mod32> : dot_product_mod32_value {}; template struct dot_product_mod32, MonoidMul>> : dot_product_mod32_value {}; } namespace internal { template 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 constexpr auto primitive_root_for_convolution = primitive_root_constexpr(m); template 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 constexpr auto ntt_generator_for_convolution = ntt_generator_constexpr(m); template > struct fft_info { static constexpr int rank2 = countr_zero(mint::mod() - 1); std::array root; std::array iroot; std::array rate2; std::array irate2; std::array rate3; std::array 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]; } { mint prod = 1, iprod = 1; for (int i = 0; i <= rank2 - 2; i++) { rate2[i] = root[i + 2] * prod; irate2[i] = iroot[i + 2] * iprod; prod *= iroot[i + 2]; iprod *= root[i + 2]; } } { mint prod = 1, iprod = 1; for (int i = 0; i <= rank2 - 3; i++) { rate3[i] = root[i + 3] * prod; irate3[i] = iroot[i + 3] * iprod; prod *= iroot[i + 3]; iprod *= root[i + 3]; } } } }; } template bool ntt_ok(int n) { if (n <= 0) return false; if constexpr (is_static_modint_v) { if constexpr (!internal::isprime) return false; static constexpr int rank2 = countr_zero(mint::mod() - 1); return ull(n) <= (1ULL << rank2); } else return false; } template void ntt(vc> &) { assert(false); } template void ntt(vc, PowInv>> &a) { using mint = internal::modint_impl, PowInv>; int n = int(a.size()); assert(n > 0); int h = countr_zero((unsigned int)n); assert(n == (1 << h)); assert(ntt_ok(n)); using wide = conditional_t<(mint::mod() <= INT_MAX), ull, u128>; static const internal::fft_info info; int len = 0; while (len < h) { if (h - len == 1) { int p = 1 << (h - len - 1); mint rot = 1; for (int s = 0; s < (1 << len); s++) { int offset = s << (h - len); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p] * rot; a[i + offset] = l + r; a[i + offset + p] = 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++) { auto mod2 = wide(mint::mod()) * mint::mod(); auto a0 = wide(a[i + offset].val()); auto a1 = wide(a[i + offset + p].val()) * rot.val(); auto a2 = wide(a[i + offset + 2 * p].val()) * rot2.val(); auto a3 = wide(a[i + offset + 3 * p].val()) * rot3.val(); auto a1na3imag = wide(mint(a1 + mod2 - a3).val()) * imag.val(); auto na2 = mod2 - a2; a[i + offset] = a0 + a2 + a1 + a3; a[i + offset + 1 * p] = a0 + a2 + (2 * mod2 - (a1 + a3)); a[i + offset + 2 * p] = a0 + na2 + a1na3imag; a[i + offset + 3 * p] = a0 + na2 + (mod2 - a1na3imag); } if (s + 1 != (1 << len)) rot *= info.rate3[countr_zero(~(unsigned int)(s))]; } len += 2; } } } template void intt(vc, PowInv>> &a) { using mint = internal::modint_impl, PowInv>; int n = int(a.size()); assert(n > 0); int h = countr_zero((unsigned int)n); assert(n == (1 << h)); assert(ntt_ok(n)); using wide = conditional_t<(mint::mod() <= INT_MAX), ull, u128>; static const internal::fft_info info; int len = h; while (len) { if (len == 1) { int p = 1 << (h - len); mint irot = 1; for (int s = 0; s < (1 << (len - 1)); s++) { int offset = s << (h - len + 1); for (int i = 0; i < p; i++) { auto l = a[i + offset]; auto r = a[i + offset + p]; a[i + offset] = l + r; a[i + offset + p] = (wide(mint::mod()) + l.val() - wide(r.val())) * irot.val(); ; } if (s + 1 != (1 << (len - 1))) irot *= info.irate2[countr_zero(~(unsigned int)(s))]; } 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++) { auto a0 = wide(a[i + offset + 0 * p].val()); auto a1 = wide(a[i + offset + 1 * p].val()); auto a2 = wide(a[i + offset + 2 * p].val()); auto a3 = wide(a[i + offset + 3 * p].val()); auto a2na3iimag = wide(mint((wide(mint::mod()) + a2 - a3) * iimag.val()).val()); a[i + offset] = a0 + a1 + a2 + a3; a[i + offset + 1 * p] = (a0 + (mint::mod() - a1) + a2na3iimag) * irot.val(); a[i + offset + 2 * p] = (a0 + a1 + (mint::mod() - a2) + (mint::mod() - a3)) * irot2.val(); a[i + offset + 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; } } } namespace internal { } void main2() { const ll P = mint103::mod(); const ll g = primitive_root(P, factorize(P - 1)); LL(N, Q); vc A(N); rep(i, N) A[i] = xor128(); dump(A); vc pw(P - 1); pw[0] = 1; rep(i, 1, P - 1) pw[i] = pw[i - 1] * g; vl lg(P); rep(i, P - 1) lg[pw[i].val()] = i; vl B; rep(i, N) if (A[i] != 0) B.eb(lg[A[i].val()]); vl ipw(P); rep(i, P - 1) ipw[pw[i].val()] = i; const ll K = 2048; vc f(P - 1); fe(b : B) f[safemod(-b, P - 1)] = 1; auto nf = f; nf.resize(2 * bit_ceil(SZ(nf))); ntt(nf); vvc hs; rep(block, divceil(P, K)) { vc g(P - 1); rep(ans, block * K, min(P, (block + 1) * K)) { ll z = ipw[ans]; g[z] = 1; } g.resize(2 * bit_ceil(SZ(g))); ntt(g); rep(y, SZ(g)) g[y] *= nf[y]; intt(g); const mint998 iv = mint998(g.size()).inv(); rep(y, SZ(g)) g[y] *= iv; rep(y, P - 1, SZ(g)) g[y % (P - 1)] += g[y]; g.resize(P - 1); hs.eb(g); } rep(_, Q) { LL(x); if (x == 0) PRINT(0); else { ll y = lg[x]; // ll ans = 0; // fe(b : B) chmax(ans, pw[(y + b) % (P - 1)].val()); // PRINT(ans); rep(block, divceil(P, K) - 1, -1, -1) { if (hs[block][y] != 0) { rep(ans, min(P, (block + 1) * K) - 1, block * K - 1, -1) { ll z = ipw[ans]; ll b = safemod(z - y, P - 1); if (f[safemod(-b, P - 1)] != 0) { PRINT(ans); break; } } break; } } } } } void test() { } // https://github.com/miscalculation53/library/tree/wip/template/template_main.hpp template 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(); CERR("\n[SINGLE_TESTCASE]\n\n", "36"); main2(); } }; Main main_dummy; int main() {} // submit-code: hash=09a7acc3bf18 source=09a7acc3bf1830af566dbeac1ed8b9da1799f36d5f0f02d20760e0b8a372cb4e input=df2908542d42c6758df5264f214259cfce58289af866b226b93ea0bf698c522b body=53164d901de5611b4a26a6ac9339c8f5e96891333c59101e7e70fc6e17a0f0ed