#pragma GCC target("avx2") #pragma GCC optimize("O3") #pragma GCC optimize("unroll-loops") // #pragma comment(linker, "/stack:200000000") #include #include #include namespace suisen { // ! utility template using constraints_t = std::enable_if_t, std::nullptr_t>; template constexpr decltype(auto) constexpr_if(Then&& then, OrElse&& or_else) { if constexpr (cond_v) { return std::forward(then); } else { return std::forward(or_else); } } // ! function template using is_same_as_invoke_result = std::is_same, ReturnType>; template using is_uni_op = is_same_as_invoke_result; template using is_bin_op = is_same_as_invoke_result; template using is_comparator = std::is_same, bool>; // ! integral template >> constexpr int bit_num = std::numeric_limits>::digits; template struct is_nbit { static constexpr bool value = bit_num == n; }; template static constexpr bool is_nbit_v = is_nbit::value; // ? template struct safely_multipliable {}; template <> struct safely_multipliable { using type = long long; }; template <> struct safely_multipliable { using type = __int128_t; }; template <> struct safely_multipliable { using type = unsigned long long; }; template <> struct safely_multipliable { using type = __uint128_t; }; template <> struct safely_multipliable { using type = float; }; template <> struct safely_multipliable { using type = double; }; template <> struct safely_multipliable { using type = long double; }; template using safely_multipliable_t = typename safely_multipliable::type; } // namespace suisen // ! type aliases using i128 = __int128_t; using u128 = __uint128_t; using ll = long long; using uint = unsigned int; using ull = unsigned long long; template using vec = std::vector; template using vec2 = vec>; template using vec3 = vec>; template using vec4 = vec>; template using pq_greater = std::priority_queue, std::greater>; template using umap = std::unordered_map; // ! macros (capital: internal macro) #define OVERLOAD2(_1,_2,name,...) name #define OVERLOAD3(_1,_2,_3,name,...) name #define OVERLOAD4(_1,_2,_3,_4,name,...) name #define REP4(i,l,r,s) for(std::remove_reference_t>i=(l);i<(r);i+=(s)) #define REP3(i,l,r) REP4(i,l,r,1) #define REP2(i,n) REP3(i,0,n) #define REPINF3(i,l,s) for(std::remove_reference_t>i=(l);;i+=(s)) #define REPINF2(i,l) REPINF3(i,l,1) #define REPINF1(i) REPINF2(i,0) #define RREP4(i,l,r,s) for(std::remove_reference_t>i=(l)+fld((r)-(l)-1,s)*(s);i>=(l);i-=(s)) #define RREP3(i,l,r) RREP4(i,l,r,1) #define RREP2(i,n) RREP3(i,0,n) #define rep(...) OVERLOAD4(__VA_ARGS__, REP4 , REP3 , REP2 )(__VA_ARGS__) #define rrep(...) OVERLOAD4(__VA_ARGS__, RREP4 , RREP3 , RREP2 )(__VA_ARGS__) #define repinf(...) OVERLOAD3(__VA_ARGS__, REPINF3, REPINF2, REPINF1)(__VA_ARGS__) #define CAT_I(a, b) a##b #define CAT(a, b) CAT_I(a, b) #define UNIQVAR(tag) CAT(tag, __LINE__) #define loop(n) for (std::remove_reference_t> UNIQVAR(loop_variable) = n; UNIQVAR(loop_variable) --> 0;) #define all(iterable) (iterable).begin(), (iterable).end() #define input(type, ...) type __VA_ARGS__; read(__VA_ARGS__) // ! I/O utilities // pair template std::ostream& operator<<(std::ostream& out, const std::pair &a) { return out << a.first << ' ' << a.second; } // tuple template std::ostream& operator<<(std::ostream& out, const std::tuple &a) { if constexpr (N >= std::tuple_size_v>) { return out; } else { out << std::get(a); if constexpr (N + 1 < std::tuple_size_v>) { out << ' '; } return operator<<(out, a); } } // vector template std::ostream& operator<<(std::ostream& out, const std::vector &a) { for (auto it = a.begin(); it != a.end();) { out << *it; if (++it != a.end()) out << ' '; } return out; } // array template std::ostream& operator<<(std::ostream& out, const std::array &a) { for (auto it = a.begin(); it != a.end();) { out << *it; if (++it != a.end()) out << ' '; } return out; } inline void print() { std::cout << '\n'; } template inline void print(const Head &head, const Tail &...tails) { std::cout << head; if (sizeof...(tails)) std::cout << ' '; print(tails...); } template auto print_all(const Iterable& v, std::string sep = " ", std::string end = "\n") -> decltype(std::cout << *v.begin(), void()) { for (auto it = v.begin(); it != v.end();) { std::cout << *it; if (++it != v.end()) std::cout << sep; } std::cout << end; } // pair template std::istream& operator>>(std::istream& in, std::pair &a) { return in >> a.first >> a.second; } // tuple template std::istream& operator>>(std::istream& in, std::tuple &a) { if constexpr (N >= std::tuple_size_v>) { return in; } else { return operator>>(in >> std::get(a), a); } } // vector template std::istream& operator>>(std::istream& in, std::vector &a) { for (auto it = a.begin(); it != a.end(); ++it) in >> *it; return in; } // array template std::istream& operator>>(std::istream& in, std::array &a) { for (auto it = a.begin(); it != a.end(); ++it) in >> *it; return in; } template void read(Args &...args) { ( std::cin >> ... >> args ); } // ! integral utilities // Returns pow(-1, n) template constexpr inline int pow_m1(T n) { return -(n & 1) | 1; } // Returns pow(-1, n) template <> constexpr inline int pow_m1(bool n) { return -int(n) | 1; } // Returns floor(x / y) template constexpr inline T fld(const T x, const T y) { return (x ^ y) >= 0 ? x / y : (x - (y + pow_m1(y >= 0))) / y; } template constexpr inline T cld(const T x, const T y) { return (x ^ y) <= 0 ? x / y : (x + (y + pow_m1(y >= 0))) / y; } template > = nullptr> constexpr inline int popcount(const T x) { return __builtin_popcount(x); } template > = nullptr> constexpr inline int popcount(const T x) { return __builtin_popcount(x); } template > = nullptr> constexpr inline int popcount(const T x) { return __builtin_popcountll(x); } template > = nullptr> constexpr inline int count_lz(const T x) { return x ? __builtin_clz(x) : suisen::bit_num; } template > = nullptr> constexpr inline int count_lz(const T x) { return x ? __builtin_clz(x) : suisen::bit_num; } template > = nullptr> constexpr inline int count_lz(const T x) { return x ? __builtin_clzll(x) : suisen::bit_num; } template > = nullptr> constexpr inline int count_tz(const T x) { return x ? __builtin_ctz(x) : suisen::bit_num; } template > = nullptr> constexpr inline int count_tz(const T x) { return x ? __builtin_ctz(x) : suisen::bit_num; } template > = nullptr> constexpr inline int count_tz(const T x) { return x ? __builtin_ctzll(x) : suisen::bit_num; } template constexpr inline int floor_log2(const T x) { return suisen::bit_num - 1 - count_lz(x); } template constexpr inline int ceil_log2(const T x) { return floor_log2(x) + ((x & -x) != x); } template constexpr inline int kth_bit(const T x, const unsigned int k) { return (x >> k) & 1; } template constexpr inline int parity(const T x) { return popcount(x) & 1; } struct all_subset { struct all_subset_iter { const int s; int t; constexpr all_subset_iter(int s) : s(s), t(s + 1) {} constexpr auto operator*() const { return t; } constexpr auto operator++() {} constexpr auto operator!=(std::nullptr_t) { return t ? (--t &= s, true) : false; } }; int s; constexpr all_subset(int s) : s(s) {} constexpr auto begin() { return all_subset_iter(s); } constexpr auto end() { return nullptr; } }; // ! container template > = nullptr> auto priqueue_comp(const Comparator comparator) { return std::priority_queue, Comparator>(comparator); } template auto isize(const Iterable &iterable) -> decltype(int(iterable.size())) { return iterable.size(); } template > = nullptr> auto generate_vector(int n, Gen generator) { std::vector v(n); for (int i = 0; i < n; ++i) v[i] = generator(i); return v; } template auto generate_range_vector(T l, T r) { return generate_vector(r - l, [l](int i) { return l + i; }); } template auto generate_range_vector(T n) { return generate_range_vector(0, n); } template void sort_unique_erase(std::vector &a) { std::sort(a.begin(), a.end()); a.erase(std::unique(a.begin(), a.end()), a.end()); } template auto foreach_adjacent_values(InputIterator first, InputIterator last, BiConsumer f) -> decltype(f(*first++, *last), void()) { if (first != last) for (auto itr = first, itl = itr++; itr != last; itl = itr++) f(*itl, *itr); } template auto foreach_adjacent_values(Container c, BiConsumer f) -> decltype(c.begin(), c.end(), void()){ foreach_adjacent_values(c.begin(), c.end(), f); } // ! other utilities // x <- min(x, y). returns true iff `x` has chenged. template inline bool chmin(T &x, const T &y) { if (y >= x) return false; x = y; return true; } // x <- max(x, y). returns true iff `x` has chenged. template inline bool chmax(T &x, const T &y) { if (y <= x) return false; x = y; return true; } namespace suisen {} using namespace suisen; using namespace std; struct io_setup { io_setup(int precision = 20) { std::ios::sync_with_stdio(false); std::cin.tie(nullptr); std::cout << std::fixed << std::setprecision(precision); } } io_setup_{}; // ! code from here #include #include #ifdef _MSC_VER #include #endif #include #ifdef _MSC_VER #endif namespace atcoder { namespace internal { // @param m `1 <= m` // @return x mod m constexpr long long safe_mod(long long x, long long m) { x %= m; if (x < 0) x += m; return x; } // Fast modular multiplication by barrett reduction // Reference: https://en.wikipedia.org/wiki/Barrett_reduction // NOTE: reconsider after Ice Lake struct barrett { unsigned int _m; unsigned long long im; // @param m `1 <= m < 2^31` explicit barrett(unsigned int m) : _m(m), im((unsigned long long)(-1) / m + 1) {} // @return m unsigned int umod() const { return _m; } // @param a `0 <= a < m` // @param b `0 <= b < m` // @return `a * b % m` unsigned int mul(unsigned int a, unsigned int b) const { // [1] m = 1 // a = b = im = 0, so okay // [2] m >= 2 // im = ceil(2^64 / m) // -> im * m = 2^64 + r (0 <= r < m) // let z = a*b = c*m + d (0 <= c, d < m) // a*b * im = (c*m + d) * im = c*(im*m) + d*im = c*2^64 + c*r + d*im // c*r + d*im < m * m + m * im < m * m + 2^64 + m <= 2^64 + m * (m + 1) < 2^64 * 2 // ((ab * im) >> 64) == c or c + 1 unsigned long long z = a; z *= b; #ifdef _MSC_VER unsigned long long x; _umul128(z, im, &x); #else unsigned long long x = (unsigned long long)(((unsigned __int128)(z)*im) >> 64); #endif unsigned int v = (unsigned int)(z - x * _m); if (_m <= v) v += _m; return v; } }; // @param n `0 <= n` // @param m `1 <= m` // @return `(x ** n) % m` constexpr long long pow_mod_constexpr(long long x, long long n, int m) { if (m == 1) return 0; unsigned int _m = (unsigned int)(m); unsigned long long r = 1; unsigned long long y = safe_mod(x, m); while (n) { if (n & 1) r = (r * y) % _m; y = (y * y) % _m; n >>= 1; } return r; } // Reference: // M. Forisek and J. Jancina, // Fast Primality Testing for Integers That Fit into a Machine Word // @param n `0 <= n` constexpr bool is_prime_constexpr(int n) { if (n <= 1) return false; if (n == 2 || n == 7 || n == 61) return true; if (n % 2 == 0) return false; long long d = n - 1; while (d % 2 == 0) d /= 2; constexpr long long bases[3] = {2, 7, 61}; for (long long a : bases) { long long t = d; long long y = pow_mod_constexpr(a, t, n); while (t != n - 1 && y != 1 && y != n - 1) { y = y * y % n; t <<= 1; } if (y != n - 1 && t % 2 == 0) { return false; } } return true; } template constexpr bool is_prime = is_prime_constexpr(n); // @param b `1 <= b` // @return pair(g, x) s.t. g = gcd(a, b), xa = g (mod b), 0 <= x < b/g constexpr std::pair inv_gcd(long long a, long long b) { a = safe_mod(a, b); if (a == 0) return {b, 0}; // Contracts: // [1] s - m0 * a = 0 (mod b) // [2] t - m1 * a = 0 (mod b) // [3] s * |m1| + t * |m0| <= b long long s = b, t = a; long long m0 = 0, m1 = 1; while (t) { long long u = s / t; s -= t * u; m0 -= m1 * u; // |m1 * u| <= |m1| * s <= b // [3]: // (s - t * u) * |m1| + t * |m0 - m1 * u| // <= s * |m1| - t * u * |m1| + t * (|m0| + |m1| * u) // = s * |m1| + t * |m0| <= b auto tmp = s; s = t; t = tmp; tmp = m0; m0 = m1; m1 = tmp; } // by [3]: |m0| <= b/g // by g != b: |m0| < b/g if (m0 < 0) m0 += b / s; return {s, m0}; } // Compile time primitive root // @param m must be prime // @return primitive root (and minimum in now) constexpr int primitive_root_constexpr(int 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; int divs[20] = {}; divs[0] = 2; int cnt = 1; int x = (m - 1) / 2; while (x % 2 == 0) x /= 2; for (int i = 3; (long long)(i)*i <= x; i += 2) { if (x % i == 0) { divs[cnt++] = i; while (x % i == 0) { x /= i; } } } if (x > 1) { divs[cnt++] = x; } for (int g = 2;; g++) { bool ok = true; for (int i = 0; i < cnt; i++) { if (pow_mod_constexpr(g, (m - 1) / divs[i], m) == 1) { ok = false; break; } } if (ok) return g; } } template constexpr int primitive_root = primitive_root_constexpr(m); // @param n `n < 2^32` // @param m `1 <= m < 2^32` // @return sum_{i=0}^{n-1} floor((ai + b) / m) (mod 2^64) unsigned long long floor_sum_unsigned(unsigned long long n, unsigned long long m, unsigned long long a, unsigned long long b) { unsigned long long ans = 0; while (true) { if (a >= m) { ans += n * (n - 1) / 2 * (a / m); a %= m; } if (b >= m) { ans += n * (b / m); b %= m; } unsigned long long y_max = a * n + b; if (y_max < m) break; // y_max < m * (n + 1) // floor(y_max / m) <= n n = (unsigned long long)(y_max / m); b = (unsigned long long)(y_max % m); std::swap(m, a); } return ans; } } // namespace internal } // namespace atcoder namespace atcoder { namespace internal { #ifndef _MSC_VER template using is_signed_int128 = typename std::conditional::value || std::is_same::value, std::true_type, std::false_type>::type; template using is_unsigned_int128 = typename std::conditional::value || std::is_same::value, std::true_type, std::false_type>::type; template using make_unsigned_int128 = typename std::conditional::value, __uint128_t, unsigned __int128>; template using is_integral = typename std::conditional::value || is_signed_int128::value || is_unsigned_int128::value, std::true_type, std::false_type>::type; template using is_signed_int = typename std::conditional<(is_integral::value && std::is_signed::value) || is_signed_int128::value, std::true_type, std::false_type>::type; template using is_unsigned_int = typename std::conditional<(is_integral::value && std::is_unsigned::value) || is_unsigned_int128::value, std::true_type, std::false_type>::type; template using to_unsigned = typename std::conditional< is_signed_int128::value, make_unsigned_int128, typename std::conditional::value, std::make_unsigned, std::common_type>::type>::type; #else template using is_integral = typename std::is_integral; template using is_signed_int = typename std::conditional::value && std::is_signed::value, std::true_type, std::false_type>::type; template using is_unsigned_int = typename std::conditional::value && std::is_unsigned::value, std::true_type, std::false_type>::type; template using to_unsigned = typename std::conditional::value, std::make_unsigned, std::common_type>::type; #endif template using is_signed_int_t = std::enable_if_t::value>; template using is_unsigned_int_t = std::enable_if_t::value>; template using to_unsigned_t = typename to_unsigned::type; } // namespace internal } // namespace atcoder namespace atcoder { namespace internal { struct modint_base {}; struct static_modint_base : modint_base {}; template using is_modint = std::is_base_of; template using is_modint_t = std::enable_if_t::value>; } // namespace internal template * = nullptr> struct static_modint : internal::static_modint_base { using mint = static_modint; public: static constexpr int mod() { return m; } static mint raw(int v) { mint x; x._v = v; return x; } static_modint() : _v(0) {} template * = nullptr> static_modint(T v) { long long x = (long long)(v % (long long)(umod())); if (x < 0) x += umod(); _v = (unsigned int)(x); } template * = nullptr> static_modint(T v) { _v = (unsigned int)(v % umod()); } unsigned int val() const { return _v; } mint& operator++() { _v++; if (_v == umod()) _v = 0; return *this; } mint& operator--() { if (_v == 0) _v = umod(); _v--; return *this; } mint operator++(int) { mint result = *this; ++*this; return result; } mint operator--(int) { mint result = *this; --*this; return result; } mint& operator+=(const mint& rhs) { _v += rhs._v; if (_v >= umod()) _v -= umod(); return *this; } mint& operator-=(const mint& rhs) { _v -= rhs._v; if (_v >= umod()) _v += umod(); return *this; } mint& operator*=(const mint& rhs) { unsigned long long z = _v; z *= rhs._v; _v = (unsigned int)(z % umod()); return *this; } mint& operator/=(const mint& rhs) { return *this = *this * rhs.inv(); } mint operator+() const { return *this; } mint operator-() const { return mint() - *this; } mint pow(long long n) const { assert(0 <= n); mint x = *this, r = 1; while (n) { if (n & 1) r *= x; x *= x; n >>= 1; } return r; } mint inv() const { if (prime) { assert(_v); return pow(umod() - 2); } else { auto eg = internal::inv_gcd(_v, m); assert(eg.first == 1); return eg.second; } } 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; } private: unsigned int _v; static constexpr unsigned int umod() { return m; } static constexpr bool prime = internal::is_prime; }; template struct dynamic_modint : internal::modint_base { using mint = dynamic_modint; public: static int mod() { return (int)(bt.umod()); } static void set_mod(int m) { assert(1 <= m); bt = internal::barrett(m); } static mint raw(int v) { mint x; x._v = v; return x; } dynamic_modint() : _v(0) {} template * = nullptr> dynamic_modint(T v) { long long x = (long long)(v % (long long)(mod())); if (x < 0) x += mod(); _v = (unsigned int)(x); } template * = nullptr> dynamic_modint(T v) { _v = (unsigned int)(v % mod()); } unsigned int val() const { return _v; } mint& operator++() { _v++; if (_v == umod()) _v = 0; return *this; } mint& operator--() { if (_v == 0) _v = umod(); _v--; return *this; } mint operator++(int) { mint result = *this; ++*this; return result; } mint operator--(int) { mint result = *this; --*this; return result; } mint& operator+=(const mint& rhs) { _v += rhs._v; if (_v >= umod()) _v -= umod(); return *this; } mint& operator-=(const mint& rhs) { _v += mod() - rhs._v; if (_v >= umod()) _v -= umod(); return *this; } mint& operator*=(const mint& rhs) { _v = bt.mul(_v, rhs._v); return *this; } mint& operator/=(const mint& rhs) { return *this = *this * rhs.inv(); } mint operator+() const { return *this; } mint operator-() const { return mint() - *this; } mint pow(long long n) const { assert(0 <= n); mint x = *this, r = 1; while (n) { if (n & 1) r *= x; x *= x; n >>= 1; } return r; } mint inv() const { auto eg = internal::inv_gcd(_v, mod()); assert(eg.first == 1); return eg.second; } 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; } private: unsigned int _v; static internal::barrett bt; static unsigned int umod() { return bt.umod(); } }; template internal::barrett dynamic_modint::bt(998244353); using modint998244353 = static_modint<998244353>; using modint1000000007 = static_modint<1000000007>; using modint = dynamic_modint<-1>; namespace internal { template using is_static_modint = std::is_base_of; template using is_static_modint_t = std::enable_if_t::value>; template struct is_dynamic_modint : public std::false_type {}; template struct is_dynamic_modint> : public std::true_type {}; template using is_dynamic_modint_t = std::enable_if_t::value>; } // namespace internal } // namespace atcoder using mint = atcoder::modint1000000007; std::istream& operator>>(std::istream& in, mint& a) { long long e; in >> e; a = e; return in; } std::ostream& operator<<(std::ostream& out, const mint& a) { out << a.val(); return out; } #include #include namespace suisen { struct Permutation : public std::vector { using base_type = std::vector; using base_type::base_type; Permutation() : Permutation(0) {} Permutation(std::size_t n) : Permutation(int(n)) {} Permutation(int n) : base_type(n) { std::iota(begin(), end(), 0); } Permutation(const std::vector& a) : std::vector(a) {} Permutation(std::vector&& a) : std::vector(std::move(a)) {} // returns b s.t. b[i] = a[p[i]] template std::vector permute(const std::vector& a) const { const int n = a.size(); std::vector res(n); for (int i = 0; i < n; ++i) res[i] = a[(*this)[i]]; return res; } // returns b s.t. b[p[i]] = a[i] template std::vector inv_permute(const std::vector& a) const { const int n = a.size(); std::vector res(n); for (int i = 0; i < n; ++i) res[(*this)[i]] = a[i]; return res; } // returns p * q s.t. (p * q)[i] = p[q[i]] friend Permutation operator*(const Permutation& p, const Permutation& q) { return q.permute(p); } Permutation& operator*=(const Permutation& q) { return *this = *this * q; } Permutation inv() const { const int n = size(); Permutation res(n); for (int i = 0; i < n; ++i) res[(*this)[i]] = i; return res; } Permutation pow(long long k) const { assert(k >= 0); const int n = size(); std::vector seen(n, false); Permutation res(n); for (int s = 0; s < n; ++s) { if (seen[s]) continue; std::vector cycle { s }; for (int v = (*this)[s]; v != s; v = (*this)[v]) cycle.push_back(v); const int l = cycle.size(); for (int j = 0; j < l; ++j) { int v = cycle[j]; seen[v] = true; res[v] = cycle[(j + k) % l]; } } return res; } template > static Permutation index_sort(const std::vector& a, Comp comp = Comp{}) { Permutation p(a.size()); std::sort(p.begin(), p.end(), [&](int i, int j) { return comp(a[i], a[j]); }); return p; } }; template struct PermutationHash { hash_t operator()(const Permutation &p) const { return hash(p); } /** * minimal perfect hash function for permutations. * complexity: O(n) time, O(n) extra space * reference: https://twitter.com/noshi91/status/1452081886025555972?s=20 */ static hash_t hash(const Permutation &per) { hash_t h = 0; const int n = per.size(); Permutation p = per; Permutation q = per.inv(); for (int i = n - 1; i >= 0; --i) { h = h * (i + 1) + p[i]; p[q[i]] = p[i]; q[p[i]] = q[i]; } return h; } static Permutation unhash(int n, hash_t h) { Permutation p = Permutation(n), q = p; for (int i = 0; i < n; ++i) { p[i] = h % (i + 1), q[i] = q[p[i]]; q[p[i]] = p[q[i]] = i; h /= i + 1; } return p; } }; } // namespace suisen namespace suisen { template std::pair& operator+=(std::pair &p1, const std::pair &p2) { p1.first += p2.first, p1.second += p2.second; return p1; } template std::pair operator+(const std::pair &p1, const std::pair &p2) { return {p1.first + p2.first, p1.second + p2.second}; } template std::pair& operator-=(std::pair &p1, const std::pair &p2) { p1.first -= p2.first, p1.second -= p2.second; return p1; } template std::pair operator-(const std::pair &p1, const std::pair &p2) { return {p1.first - p2.first, p1.second - p2.second}; } template std::pair& operator*=(std::pair &p, const V m) { p.first *= m, p.second *= m; return p; } template std::pair operator*(const std::pair &p, const V m) { return {p.first * m, p.second * m}; } template std::pair operator*(const V m, const std::pair &p) { return {p.first * m, p.second * m}; } } // namespace suisen using S = pair; struct FenwickTree1D { int N; vector BIT; FenwickTree1D(int N, S qdef = S{}) : N(N), BIT(N + 1, qdef) {} template FenwickTree1D(int N, F f, S qdef = S{}) : FenwickTree1D(N, qdef) { for (int i = 1; i <= N; i++) { BIT[i] += f(); int j = i + (i & -i); if (j <= N) BIT[j] += BIT[i]; } } template FenwickTree1D(It st, It en, S qdef = S{}) : FenwickTree1D(en - st, [&] { return *st++; }, qdef) {} vector values() { vector ret(BIT.begin() + 1, BIT.end()); for (int i = N; i >= 1; i--) { int j = i + (i & -i); if (j <= N) ret[j - 1] -= ret[i - 1]; } return ret; } void update(int i, S v) { for (i++; i <= N; i += i & -i) BIT[i] += v; } S query(int r) { S ret = BIT[0]; for (r++; r > 0; r -= r & -r) ret += BIT[r]; return ret; } S query(int l, int r) { return query(r) - query(l - 1); } }; // add (x, y, v) // query [x1, x2] x [y1, y2] struct Query { int t, xlo, xhi, ylo, yhi; S v; Query(int t, int xlo, int xhi, int ylo, int yhi, S v) : t(t), xlo(xlo), xhi(xhi), ylo(ylo), yhi(yhi), v(v) {} }; void add_point(vector& queries, int x, int y, S v) { queries.push_back(Query(0, x, x, y, y, v)); } void add_query(vector& queries, int l, int d, int r, int u) { queries.push_back(Query(1, l, r, d, u, S{})); } struct Event { int t, x, ylo, yhi; S v; Event(int t, int x, int ylo, int yhi, S v) : t(t), x(x), ylo(ylo), yhi(yhi), v(v) {} bool operator<(const Event& e) const { return make_pair(x, t) < make_pair(e.x, e.t); } }; vector solve(vector queries) { vector events; events.reserve(queries.size() * 2); vector YS; YS.reserve(queries.size()); for (auto&& q : queries) if (q.t == 0) YS.push_back(q.ylo); sort(YS.begin(), YS.end()); YS.erase(unique(YS.begin(), YS.end()), YS.end()); int qcnt = 0; for (auto&& q : queries) { if (q.t == 0) q.ylo = lower_bound(YS.begin(), YS.end(), q.ylo) - YS.begin(); else { q.ylo = lower_bound(YS.begin(), YS.end(), q.ylo) - YS.begin(); q.yhi = int(upper_bound(YS.begin(), YS.end(), q.yhi) - YS.begin()) - 1; q.v.first = qcnt++; if (q.ylo > q.yhi) q.t = -1; } } vector ans(qcnt); FenwickTree1D FT(YS.size()); function rec = [&](int l, int r) { if (l == r) return; int m = l + (r - l) / 2; rec(l, m); rec(m + 1, r); events.clear(); for (int i = l; i <= m; i++) if (queries[i].t == 0) events.emplace_back(0, queries[i].xlo, queries[i].ylo, queries[i].ylo, queries[i].v); if (events.empty()) return; int curSize = events.size(); for (int i = m + 1; i <= r; i++) if (queries[i].t == 1) { events.emplace_back(-1, queries[i].xlo, queries[i].ylo, queries[i].yhi, queries[i].v); events.emplace_back(1, queries[i].xhi, queries[i].ylo, queries[i].yhi, queries[i].v); } if (curSize == int(events.size())) return; sort(events.begin(), events.end()); for (auto&& e : events) { if (e.t == 0) FT.update(e.ylo, e.v); else ans[e.v.first] += e.t * FT.query(e.ylo, e.yhi); } for (auto&& e : events) if (e.t == 0) FT.update(e.ylo, { -e.v.first, -e.v.second }); }; if (qcnt > 0) rec(0, int(queries.size()) - 1); return ans; } constexpr int inf = numeric_limits::max() - 10; int main() { input(int, n); vector x(n), y(n); read(x, y); vector answers; loop(2) { auto p = Permutation::index_sort(x); x = p.permute(x); y = p.permute(y); vector queries; mint ans = 0; rep(i, n) { int u = x[i] - y[i]; int v = x[i] + y[i]; add_query(queries, -inf, -inf, u - 1, v - 1); add_point(queries, u, v, { 1, x[i] }); } int i = 0; for (auto [num, sum] : solve(queries)) { ans += mint(num) * x[i++] - sum; } answers.push_back(ans * 2); x.swap(y); } print(answers); return 0; }