Submission #43574330
Source Code Expand
#ifdef INCLUDED_MAIN
auto solve() {
GET(H, W);
GETVSTR(SS, H);
auto G = MaxFlow<ll>(H * W + 2);
ll start = H * W, end = H * W + 1;
auto getpos = [&](ll h, ll w) {
return h * W + w;
};
vll dh = {0, 1, 0, -1}, dw = {1, 0, -1, 0};
rep(h, H) {
rep(w, W) {
if (SS[h][w] == '#') continue;
ll p = getpos(h, w);
if ((h + w) % 2 == 0) {
G.append_edge(start, p, 1);
rep(i, 4) {
ll hh = h + dh[i], ww = w + dw[i];
if (!IN(0, hh, H - 1) || !IN(0, ww, W - 1)) continue;
if (SS[hh][ww] == '#') continue;
ll np = getpos(hh, ww);
G.append_edge(p, np, 1);
}
} else {
G.append_edge(p, end, 1);
}
}
}
auto ans = G.flow(start, end);
auto edges = G.edges();
rep(i, len(edges)) {
if (edges[i].flow == 0) continue; // not used
if (edges[i].from == start) continue;
if (edges[i].to == end) continue;
ll sh = edges[i].from / W, sw = edges[i].from % W;
ll eh = edges[i].to / W, ew = edges[i].to % W;
vll dir = {eh - sh, ew - sw};
if (dir == vll{0, 1}) {
SS[sh][sw] = '>';
SS[eh][ew] = '<';
} else if (dir == vll{1, 0}){
SS[sh][sw] = 'v';
SS[eh][ew] = '^';
} else if (dir == vll{0, -1}) {
SS[sh][sw] = '<';
SS[eh][ew] = '>';
} else {
SS[sh][sw] = '^';
SS[eh][ew] = 'v';
}
}
print(ans);
rep(i, len(SS)) print(SS[i]);
return _0;
}
int main() {
// mint::set_mod(998244353);
// mint::set_mod(1000000007);
mint::set_mod(1); // mint使うときに上のコメントをはずす。
auto ans = solve();
// print(ans);
}
// ラムダ再帰
// auto ff = [&](auto &&f, ll x) {};
// ff(ff, 0);
// 以下は動作確認未実施
#else
#define INCLUDED_MAIN
#ifdef LOCAL
#include "../mytemplate.hpp"
#else
#include <algorithm>
#include <bits/extc++.h>
#include <bitset>
#include <cassert>
#include <cctype>
#include <climits>
#include <cstddef>
#include <deque>
#include <functional>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <map>
#include <queue>
#include <random>
#include <set>
#include <stack>
#include <string_view>
#include <type_traits>
#include <utility>
#include <regex>
#endif
using namespace std;
// clang-format off
/* accelration */
// 高速バイナリ生成
#ifndef LOCAL
#pragma GCC target("avx")
#pragma GCC optimize("O3")
#pragma GCC optimize("unroll-loops")
#endif
// cin cout の結びつけ解除, stdioと同期しない(入出力非同期化)
// cとstdの入出力を混在させるとバグるので注意
struct IOSetting {IOSetting() {std::cin.tie(0); ios::sync_with_stdio(false); cout << fixed << setprecision(15);}} iosetting;
// unordered_mapでpair, vector, tupleをkeyにするためのコード
// (参考文献) https://qiita.com/hamamu/items/4d081751b69aa3bb3557
template<class T> size_t HashCombine(const size_t seed,const T &v){
return seed^(std::hash<T>()(v)+0x9e3779b9+(seed<<6)+(seed>>2));
}
/* pair用 */
template<class T,class S> struct std::hash<std::pair<T,S>>{
size_t operator()(const std::pair<T,S> &keyval) const noexcept {
return HashCombine(std::hash<T>()(keyval.first), keyval.second);
}
};
/* complex用 */
template<class T> struct std::hash<complex<T>>{
size_t operator()(const complex<T> &x) const noexcept {
size_t s=0;
s=HashCombine(s,x.real());
s=HashCombine(s,x.imag());
return s;
}
};
/* vector用 */
template<class T> struct std::hash<std::vector<T>>{
size_t operator()(const std::vector<T> &keyval) const noexcept {
size_t s=0;
for (auto&& v: keyval) s=HashCombine(s,v);
return s;
}
};
/* deque用 */
template<class T> struct std::hash<std::deque<T>>{
size_t operator()(const std::deque<T> &keyval) const noexcept {
size_t s=0;
for (auto&& v: keyval) s=HashCombine(s,v);
return s;
}
};
/* tuple用 */
template<int N> struct HashTupleCore{
template<class Tuple> size_t operator()(const Tuple &keyval) const noexcept{
size_t s=HashTupleCore<N-1>()(keyval);
return HashCombine(s,std::get<N-1>(keyval));
}
};
template <> struct HashTupleCore<0>{
template<class Tuple> size_t operator()(const Tuple &keyval) const noexcept{ return 0; }
};
template<class... Args> struct std::hash<std::tuple<Args...>>{
size_t operator()(const tuple<Args...> &keyval) const noexcept {
return HashTupleCore<tuple_size<tuple<Args...>>::value>()(keyval);
}
};
/* alias */
using ull = __uint128_t;
//using ll = long long; // __int128でTLEするときに切り替える。
using ll = __int128;
using ld = long double;
using vi = vector<int>;
using vl = vector<long>;
using vll = vector<ll>;
using vd = vector<ld>;
using vvi = vector<vi>;
using vvl = vector<vl>;
using vvll = vector<vll>;
using vvvll = vector<vvll>;
using vvd = vector<vd>;
using vvvd = vector<vvd>;
using vc = vector<char>;
using vvc = vector<vc>;
using vs = vector<string>;
using vvs = vector<vs>;
using vvvs = vector<vvs>;
using pii = pair<int, int>;
using pll = pair<ll, ll>;
using umpll = unordered_map<ll, ll>;
using umpsl = unordered_map<string, ll>;
using mpll = map<ll, ll>;
using sll = set<ll>;
using msll = multiset<ll>;
using heapqll = priority_queue<ll, vll, greater<ll>>;
using heapqllrev = priority_queue<ll>;
using dll = deque<ll>;
/* REP macro */
#define _overload4(_1,_2,_3,_4,name,...) name
#define _rep(i,n) reps(i,0,n)
#define reps(i,a,n) for (ll i = (a); i < (ll)(n); ++i)
#define repsp(i,a,n,s) for (ll i = (a); i < (ll)(n); i += s)
#define rep(...) _overload4(__VA_ARGS__,repsp, reps,_rep,)(__VA_ARGS__)
#define _overload4(_1,_2,_3,_4,name,...) name
#define _rrep(i,n) rreps(i,n,0)
#define rreps(i,a,n) for (ll i = (a); i >= (ll)(n); --i)
#define rrepsp(i,a,n,s) for (ll i = (a); i >= (ll)(n); i -= s)
#define rrep(...) _overload4(__VA_ARGS__, rrepsp, rreps, _rrep,)(__VA_ARGS__)
#define repd(i,n) for(ll i=n-1;i>=0;i--)
#define rrepd(i,n) for(ll i=n;i>=1;i--)
#define repdict(key, value, dict) for (const auto& [key, value] : dict)
#define repset(x, st) for(auto x : st)
/* define short */
#define endl "\n"
#define pf emplace_front
#define pb emplace_back
#define popleft pop_front
#define popright pop_back
#define mp make_pair
#define ump unordered_map
#define all(obj) (obj).begin(), (obj).end()
#define rall(obj) (obj).rbegin(), (obj).rend()
#define len(x) (ll)(x.size())
#define MAX(x) *max_element(all(x))
#define MIN(x) *min_element(all(x))
#define ARGMAX(x) distance(x.begin(), max_element(all(x)))
#define ARGMIN(x) distance(x.begin(), min_element(all(x)))
#define CLAMP(L, X, R) min(max(L, X), R)
#define IN(L, X, R) (L <= X && X <= R)
#define SET(x) sll(all(x))
#define VEC(x) vll(all(x))
#define _overload6(_1,_2,_3,_4,_5,_6,name,...) name
#define GET1(x) ll x = in_ll();
#define GET2(x, y) ll x, y; {vll _A_ = in_lls(); x = _A_[0]; y = _A_[1];}
#define GET3(x, y, z) ll x, y, z; {vll _A_ = in_lls(); x = _A_[0]; y = _A_[1]; z = _A_[2];}
#define GET4(x, y, z, a) ll x, y, z, a; {vll _A_ = in_lls(); x = _A_[0]; y = _A_[1]; z = _A_[2]; a = _A_[3];}
#define GET5(x, y, z, a, b) ll x, y, z, a, b; {vll _A_ = in_lls(); x = _A_[0]; y = _A_[1]; z = _A_[2]; a = _A_[3]; b = _A_[4];}
#define GET6(x, y, z, a, b, c) ll x, y, z, a, b, c; {vll _A_ = in_lls(); x = _A_[0]; y = _A_[1]; z = _A_[2]; a = _A_[3]; b = _A_[4]; c = _A_[5];}
#define GET(...) _overload6(__VA_ARGS__,GET6, GET5, GET4, GET3, GET2, GET1)(__VA_ARGS__)
#define GETVLL(x) vll x = in_lls();
#define GETVVLL(x, N) vvll x; rep(i, N) {GETVLL(ab); x.pb(ab);}
#define GET1D(x) ld x = in_d();
#define GET2D(x, y) ld x, y; {vd _A_ = in_ds(); x = _A_[0]; y = _A_[1];}
#define GET3D(x, y, z) ld x, y, z; {vd _A_ = in_ds(); x = _A_[0]; y = _A_[1]; z = _A_[2];}
#define GET4D(x, y, z, a) ld x, y, z, a; {vd _A_ = in_ds(); x = _A_[0]; y = _A_[1]; z = _A_[2]; a = _A_[3];}
#define GET5D(x, y, z, a, b) ld x, y, z, a, b; {vd _A_ = in_ds(); x = _A_[0]; y = _A_[1]; z = _A_[2]; a = _A_[3]; b = _A_[4];}
#define GET6D(x, y, z, a, b, c) ld x, y, z, a, b, c; {vd _A_ = in_ds(); x = _A_[0]; y = _A_[1]; z = _A_[2]; a = _A_[3]; b = _A_[4]; c = _A_[5];}
#define GETD(...) _overload6(__VA_ARGS__,GET6D, GET5D, GET4D, GET3D, GET2D, GET1D)(__VA_ARGS__)
#define GETVD(x) vd x = in_ds();
#define GETVVD(x, N) vvd x; rep(i, N) {GETVD(ab); x.pb(ab);}
#define GETSTR(x) string x = in_str();
#define GETSTRS(x) vs x; x = in_strs();
#define GETVVS(x, N) vvs x; rep(i, N) x.pb(in_strs());
#define GETVSTR(x, N) vs x; rep(i, N) x.pb(in_str());
#define GETPOINT(p) Point p; {GET2(x, y); p = Point{x, y};}
#define GETPOINTS(p, N) vector<Point> p; rep(i, N) {GET2(x, y); p.pb(Point{x, y});}
#define GETCOMPLEX(p) complex<ld> p; {GET2D(x, y); p = complex<ld>{x, y};}
#define GETCOMPLEXS(p, N) vector<complex<ld>> p; rep(i, N) {GET2D(x, y); p.pb(complex<ld>{x, y});}
#define _overload7(_1,_2,_3,_4,_5,_6,_7,name,...) name
#define INI1(x, vec) auto x = vec[0];
#define INI2(x, y, vec) auto x = vec[0], y = vec[1];
#define INI3(x, y, z, vec) auto x = vec[0], y = vec[1], z = vec[2];
#define INI4(x, y, z, a, vec) auto x = vec[0], y = vec[1], z = vec[2], a = vec[3];
#define INI5(x, y, z, a, b, vec) auto x = vec[0], y = vec[1], z = vec[2], a = vec[3], b = vec[4];
#define INI6(x, y, z, a, b, c, vec) auto x = vec[0], y = vec[1], z = vec[2], a = vec[3], b = vec[4], c = vec[5];
#define INI(...) _overload7(__VA_ARGS__,INI6, INI5, INI4, INI3, INI2, INI1)(__VA_ARGS__)
#define SETPERM(x, N) vll x(N); iota(all(x), 0);
#define SETPERMS(x, s, N) vll x(N); iota(all(x), s);
#define UNUSED(x) ((void)x);
#define printF(x) print(x); cout << flush;
// [INT|LLONG|DBL|LDBL]_[MAX|MIN] 最大最小表現
// 型変換
#define CHARSTR(x) (""s + x)
#define STRBIN2LL(x) ((ll)std::stoull(x, nullptr, 2))
#define STRLL(x) ((ll)parse(x))
#define STRD(x) std::stod(x)
#define CHARLL(x) ((ll)std::stoll(CHARSTR(x)))
/* sort */
#define SORT(x) stable_sort(all(x))
#define RSORT(x) stable_sort(rall(x))
#define SORT_IDX(x, idx) stable_sort(all(x), [&](const vll &_a_, const vll &_b_){return _a_[idx] < _b_[idx];})
#define RSORT_IDX(x, idx) stable_sort(all(x), [&](const vll &_a_, const vll &_b_){return _a_[idx] > _b_[idx];})
#define LB_IDX_VEC(c, x) distance((c).begin(), lower_bound(all(c), x)) // O(log N) x未満の最大値についてその右側のidxが求まる
#define UB_IDX_VEC(c, x) distance((c).begin(), upper_bound(all(c), x)) // O(log N) x以上の最小値についてその右側のidxが求まる
#define LB_ITR_VEC(c, x) lower_bound(all(c), x)
#define UB_ITR_VEC(c, x) upper_bound(all(c), x)
// #define LB_IDX_SET(c, x) distance((c).begin(), c.lower_bound(x)) // O(N)
// #define UB_IDX_SET(c, x) distance((c).begin(), c.upper_bound(x)) // O(N)
#define LB_ITR_SET(c, x) c.lower_bound(x)
#define UB_ITR_SET(c, x) c.upper_bound(x)
#define LB_ITR_MAP(c, x) c.lower_bound(x)
#define UB_ITR_MAP(c, x) c.upper_bound(x)
#define KEY_CHANGE(c, k1, k2) { auto i_ = c.extract(k1); i_.key() = k2; c.insert(std::move(i_));}
#define EXIST(key, dict) (dict.find(key) != dict.end())
#define REV(x) reverse(all(x))
// multisetでのerase
#define ERASE(x, s) {auto itr_ = s.find((x)); if (itr_ != s.end()) s.erase(itr_); }
// 第一引数と第二引数を比較し、第一引数(a)をより大きい/小さい値に上書き
template <typename T> inline bool chmin(T& a, const T& b) {bool compare = a > b; if (a > b) a = b; return compare;}
template <typename T> inline bool chmax(T& a, const T& b) {bool compare = a < b; if (a < b) a = b; return compare;}
inline string YESNO(bool cond) {return cond ? "YES" : "NO";}
inline string yesno(bool cond) {return cond ? "yes" : "no";}
inline string YesNo(bool cond) {return cond ? "Yes" : "No";}
namespace // 直値のデフォルトの型をllに。
{
ll _0 = 0;
ll _1 = 1;
ll _2 = 2;
ll _3 = 3;
ll _4 = 4;
ll _5 = 5;
ll _6 = 6;
ll _7 = 7;
ll _8 = 8;
ll _9 = 9;
ll _10 = 10;
ll _11 = 11;
ll _12 = 12;
ll _13 = 13;
ll _14 = 14;
ll _15 = 15;
ll _16 = 16;
ll _17 = 17;
ll _30 = 30;
ll _31 = 31;
ll _32 = 32;
ll _33 = 33;
ll _63 = 63;
ll _64 = 64;
ll _65 = 65;
ll _66 = 66;
ll _126 = 126;
ll _127 = 127;
ll _128 = 128;
ll _129 = 129;
};
void ignore_warning() // ワーニング対策
{
_0 = _0;
_1 = _1;
_2 = _2;
_3 = _3;
_4 = _4;
_5 = _5;
_6 = _6;
_7 = _7;
_8 = _8;
_9 = _9;
_10 = _10;
_11 = _11;
_12 = _12;
_13 = _13;
_14 = _14;
_15 = _15;
_16 = _16;
_17 = _17;
_30 = _30;
_31 = _31;
_32 = _32;
_33 = _33;
_63 = _63;
_64 = _64;
_65 = _65;
_66 = _66;
_126 = _126;
_127 = _127;
_128 = _128;
_129 = _129;
}
/* helper func */
std::ostream &operator<<(std::ostream &dest, __int128 value) {
std::ostream::sentry s(dest);
if (s) {
__uint128_t tmp = value < 0 ? -value : value;
char buffer[128];
char *d = std::end(buffer);
do {
--d;
*d = "0123456789"[tmp % 10];
tmp /= 10;
} while (tmp != 0);
if (value < 0) {
--d;
*d = '-';
}
int len = std::end(buffer) - d;
if (dest.rdbuf()->sputn(d, len) != len) {
dest.setstate(std::ios_base::badbit);
}
}
return dest;
}
ll parse(string &s) {
ll ret = 0;
bool isplus = true;
for (ll i = 0; i < s.length(); i++)
if ('0' <= s[i] && s[i] <= '9')
ret = 10 * ret + s[i] - '0';
else if (s[i] == '-')
isplus ^= isplus;
return isplus ? ret : -ret;
}
string STR(const vector<char> &cs) {
return string(cs.begin(), cs.end());
}
string RSTR(const vector<char> &cs) {
return string(cs.rbegin(), cs.rend());
}
template <typename T>
string STR(T v) {
ostringstream ss;
ss << v;
return ss.str();
}
namespace internal {
template <class T> struct simple_queue {
std::vector<T> payload;
int pos = 0;
void reserve(int n) { payload.reserve(n); }
int size() const { return int(payload.size()) - pos; }
bool empty() const { return pos == int(payload.size()); }
void push(const T& t) { payload.push_back(t); }
T& front() { return payload[pos]; }
void clear() {
payload.clear();
pos = 0;
}
void pop() { pos++; }
};
// @param n `0 <= n`
// @return minimum non-negative `x` s.t. `n <= 2**x`
int ceil_pow2(int n) {
int x = 0;
while ((1U << x) < (unsigned int)(n)) x++;
return x;
}
template <class T>
using is_signed_int128 =
typename std::conditional<std::is_same<T, __int128_t>::value ||
std::is_same<T, __int128>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_unsigned_int128 =
typename std::conditional<std::is_same<T, __uint128_t>::value ||
std::is_same<T, unsigned __int128>::value,
std::true_type,
std::false_type>::type;
template <class T>
using make_unsigned_int128 =
typename std::conditional<std::is_same<T, __int128_t>::value,
__uint128_t,
unsigned __int128>;
template <class T>
using is_integral = typename std::conditional<std::is_integral<T>::value ||
is_signed_int128<T>::value ||
is_unsigned_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_signed_int = typename std::conditional<(is_integral<T>::value &&
std::is_signed<T>::value) ||
is_signed_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using is_unsigned_int =
typename std::conditional<(is_integral<T>::value &&
std::is_unsigned<T>::value) ||
is_unsigned_int128<T>::value,
std::true_type,
std::false_type>::type;
template <class T>
using to_unsigned = typename std::conditional<
is_signed_int128<T>::value,
make_unsigned_int128<T>,
typename std::conditional<std::is_signed<T>::value,
std::make_unsigned<T>,
std::common_type<T>>::type>::type;
template <class T>
using is_signed_int_t = std::enable_if_t<is_signed_int<T>::value>;
template <class T>
using is_unsigned_int_t = std::enable_if_t<is_unsigned_int<T>::value>;
template <class T> using to_unsigned_t = typename to_unsigned<T>::type;
// 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`
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;
}
};
struct modint_base {};
struct static_modint_base : modint_base {};
// @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;
}
// @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 <int n> 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<long long, long long> 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 <int m> constexpr int primitive_root = primitive_root_constexpr(m);
} // namespace internal
template<int m, std::enable_if_t<(1 <= m)> * = 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<class T, internal::is_signed_int_t<T> * = nullptr>
static_modint(T v) {
long long x = (long long)(v % (long long)(umod()));
if (x < 0) x += umod();
_v = (unsigned int)(x);
}
template<class T, internal::is_unsigned_int_t<T> * = nullptr>
static_modint(T v) {
_v = (unsigned int)(v % umod());
}
static_modint(bool v) { _v = ((unsigned int)(v) % umod()); }
ll val() const { return (ll)_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<m>;
};
template<int id>
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<class T, internal::is_signed_int_t<T> * = nullptr>
dynamic_modint(T v) {
long long x = (long long)(v % (long long)(mod()));
if (x < 0) x += mod();
_v = (unsigned int)(x);
}
template<class T, internal::is_unsigned_int_t<T> * = nullptr>
dynamic_modint(T v) {
_v = (unsigned int)(v % mod());
}
dynamic_modint(bool v) { _v = ((unsigned int)(v) % mod()); }
ll val() const { return (ll)_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<int id>
internal::barrett dynamic_modint<id>::bt = 998244353;
using modint998244353 = static_modint<998244353>;
using modint1000000007 = static_modint<1000000007>;
using modint = dynamic_modint<-1>;
namespace internal {
template<class T>
using is_static_modint = std::is_base_of<internal::static_modint_base, T>;
template<class T>
using is_static_modint_t = std::enable_if_t<is_static_modint<T>::value>;
template<class>
struct is_dynamic_modint : public std::false_type {};
template<int id>
struct is_dynamic_modint<dynamic_modint<id>> : public std::true_type {};
template<class T>
using is_dynamic_modint_t = std::enable_if_t<is_dynamic_modint<T>::value>;
} // namespace internal
using mint = modint;
using vm = vector<modint>;
using vvm = vector<vm>;
template <typename T>
T SUM(const vector<T> &v) {
T total = 0;
rep(i, len(v)) {
total += v[i];
}
return total;
}
ll POPLEFT(dll &q) {
const ll v = q.front();
q.pop_front();
return v;
}
ll POP(dll &q) {
const ll v = q.back();
q.pop_back();
return v;
}
void HEAPPUSH(heapqll &q, ll v) {
q.push(v);
}
ll HEAPPOP(heapqll &q) {
const ll v = q.top();
q.pop();
return v;
}
void HEAPPUSH(heapqllrev &q, ll v) {
q.push(v);
}
ll HEAPPOP(heapqllrev &q) {
const ll v = q.top();
q.pop();
return v;
}
// 配列の回転O(N)
// SETPERM(X, 10);
// ROTATE_LEFT(X, 3);
// -> X: [ 3 4 5 6 7 8 9 0 1 2 ]
template<typename T>
void ROTATE_LEFT(vector<T> &vec, ll N) {
N = N % len(vec);
rotate(vec.begin(), vec.begin() + N, vec.end());
}
// 配列の回転O(N)
// SETPERM(X, 10);
// ROTATE_LEFT(X, 3);
// -> X: [ 7 8 9 0 1 2 3 4 5 6 ]
template<typename T>
void ROTATE_RIGHT(vector<T> &vec, ll N) {
N = N % len(vec);
rotate(vec.rbegin(), vec.rbegin() + N, vec.rend());
}
// 文字列区間swap[L, R)
string rangeswap(const string &S, ll L, ll R) {
string T = S;
ll cnt = (R - L) >> 1;
rep (i, cnt) swap(T[L + i], T[R - i - 1]);
return T;
}
template<class... T>
constexpr auto min(T... a){
return min(initializer_list<common_type_t<T...>>{a...});
}
template<class... T>
constexpr auto max(T... a){
return max(initializer_list<common_type_t<T...>>{a...});
}
// 幾何関連データ構造
constexpr ld eps = 1e-9;
// ラジアン->度
ld rad2Deg(ld rad) { return rad * 180.0 / M_PI; }
// 度->ラジアン
ld deg2Rad(ld deg) { return deg * M_PI / 180.0; }
// 小数座標の点p(x, y)
struct PointD {
ld x = 0.0;
ld y = 0.0;
explicit PointD() : x(0.0), y(0.0) {}
explicit PointD(ld ldx, ld ldy) : x(ldx), y(ldy) {}
explicit PointD(const vd &xs) : x(xs[0]), y(xs[1]) {}
explicit PointD(ll lx, ll ly) : x(static_cast<ld>(lx)), y(static_cast<ld>(ly)) {}
explicit PointD(const vll &xs) : x(xs[0]), y(xs[1]) {}
ld norm() const {
return hypot(x, y);
}
ld arg() const {
return fmodl(atan2l(y, x) + 2.0 * M_PI, 2.0 * M_PI);
}
PointD u() const {
ld L = norm();
return PointD(x / L, y / L);
}
// 1 : left, 0 : same, -1 : right
constexpr int toleft(const PointD &r) const {
auto t = x * r.y - y * r.x;
return t > eps ? 1 : t < -eps ? -1 : 0;
}
};
bool operator==(const PointD &a, const PointD &b) { return abs(a.x - b.x) < eps && abs(a.y - b.y) < eps;}
bool operator!=(const PointD &a, const PointD &b) { return !(a == b);}
PointD operator+(const PointD &a, const PointD &b) { return PointD(a.x + b.x, a.y + b.y); }
PointD operator-(const PointD &a, const PointD &b) { return PointD(a.x - b.x, a.y - b.y); }
PointD operator-(const PointD &a) { return PointD(-a.x, -a.y); }
PointD operator*(ld t, const PointD &b) { return PointD(t * b.x, t * b.y); }
PointD operator/(const PointD &a, ld t) { return PointD(a.x / t, a.y / t); }
ld operator*(const PointD &a, const PointD &b) { return a.x * b.x + a.y * b.y; } // dot product
ld operator%(const PointD &a, const PointD &b) { return a.x * b.y - a.y * b.x; } // cross product
bool operator<(const PointD &a, const PointD &b) { // lexicographical compare
if (abs(a.x - b.x) > eps) return a.x < b.x;
return a.y + eps < b.y;
}
// 整数座標の点p(x, y)
struct Point {
ll x = 0;
ll y = 0;
ld norm() const {
return hypot(x, y);
}
ld arg() const {
return fmodl(atan2l(y, x) + 2.0 * M_PI, 2.0 * M_PI);
}
PointD u() const {
ld L = norm();
return PointD(x / L, y / L);
}
// 1 : left, 0 : same, -1 : right
constexpr ll toleft(const Point &r) const {
auto t = x * r.y - y * r.x;
return t > 0 ? 1 : t < 0 ? -1 : 0;
}
};
bool operator==(const Point &p, const Point &q) { return p.x == q.x && p.y == q.y; }
bool operator!=(const Point &p, const Point &q) { return !(p == q); }
Point operator+(const Point &a, const Point &b) { return {a.x + b.x, a.y + b.y}; }
Point operator-(const Point &a, const Point &b) { return {a.x - b.x, a.y - b.y}; }
Point operator-(const Point &a) { return {-a.x, -a.y}; }
Point operator*(ll t, const Point &b) { return {t * b.x, t * b.y}; }
Point operator/(const Point &a, ll t) { return {a.x / t, a.y / t}; }
PointD operator/(const Point &a, ld t) { return PointD{(ld)a.x / t, (ld)a.y / t}; }
ll operator*(const Point &a, const Point &b) { return a.x * b.x + a.y * b.y; } // dot product
ll operator%(const Point &a, const Point &b) { return a.x * b.y - a.y * b.x; } // cross product
bool operator<(const Point &a, const Point &b) { // lexicographical compare
if (abs(a.x - b.x) > 0) return a.x < b.x;
return a.y < b.y;
}
// 小数係数の直線ax + by + c = 0
struct LineD {
ld a;
ld b;
ld c;
explicit LineD() : a(0.0), b(0.0), c(0.0) {}
explicit LineD(ld lda, ld ldb, ld ldc) : a(lda), b(ldb), c(ldc) {}
explicit LineD(const vd &xs) : a(xs[0]), b(xs[1]), c(xs[2]) {}
explicit LineD(ll la, ll lb, ll lc) : a(static_cast<ld>(la)), b(static_cast<ld>(lb)), c(static_cast<ld>(lc)){}
explicit LineD(const vll &xs) : a(xs[0]), b(xs[1]), c(xs[2]) {}
// 直線の方向ベクトルの大きさ
ld norm() const {
return hypot(a, b);
}
PointD u() {
ld L = norm();
return PointD(b / L, -a / L);
}
};
bool operator==(const LineD &l, const LineD &m) { return abs(l.a - m.a) < eps && abs(l.b - m.b) < eps && abs(l.c - m.c) < eps;}
bool operator!=(const LineD &l, const LineD &m) { return !(l == m);}
// 整数係数の直線ax + by + c = 0
struct Line {
ll a;
ll b;
ll c;
// 直線の方向ベクトルの大きさ
ld norm() const {
return hypot(a, b);
}
PointD u() {
ld L = norm();
return PointD((ld)b / L, -(ld)a / L);
}
};
bool operator==(const Line &l, const Line &m) { return l.a == m.a && l.b == m.b && l.c == m.c;}
bool operator!=(const Line &l, const Line &m) { return !(l == m);}
// 小数の2点p, qを結ぶ線分
struct SegmentLineD {
PointD p;
PointD q;
PointD dir() const { return q - p; }
ld norm() const {
auto d = dir();
return hypot(d.x, d.y);
}
PointD u() const {
assert(p != q);
return dir() / norm();
}
// pから長さtの位置
PointD vec(ld t) const {
return t * u() + p;
}
PointD midpos() const {
return PointD{(p.x + q.x) / 2.0, (p.y + q.y) / 2.0};
}
};
bool operator==(const SegmentLineD &l, const SegmentLineD &m) {
return abs(l.p.x - m.p.x) < eps &&
abs(l.p.y - m.p.y) < eps &&
abs(l.q.x - m.q.x) < eps &&
abs(l.q.y - m.q.y) < eps;
}
bool operator!=(const SegmentLineD &l, const SegmentLineD &m) { return !(l == m); }
LineD getLineD(const PointD &P, const PointD &Q) {
assert(P.x != Q.x || P.y != Q.y);
if (P.x == Q.x) {
return LineD(1.0, 0.0, -P.x);
} else if (P.y == Q.y) {
return LineD(0.0, 1.0, -P.y);
} else {
return LineD(P.y - Q.y, Q.x - P.x, P % Q);
}
}
// 小数の2点p, qを結ぶ線分
struct SegmentLine {
Point p;
Point q;
Point dir() const { return q - p; }
ld norm() const {
auto d = dir();
return hypot(d.x, d.y);
}
PointD u() const {
assert(p != q);
return dir() / norm();
}
// pから長さtの位置
PointD vec(ld t) const {
PointD s = PointD{p.x, p.y};
return t * u() + s;
}
PointD midpos() const {
return PointD{(ld)(p.x + q.x) / 2.0, (ld)(p.y + q.y) / 2.0};
}
};
bool operator==(const SegmentLine &l, const SegmentLine &m) {
return l.p.x == m.p.x && l.p.y == m.p.y && l.q.x == m.q.x && l.q.y == m.q.y;
}
bool operator!=(const SegmentLine &l, const SegmentLine &m) { return !(l == m); }
Line getLine(const Point &P, const Point &Q) {
assert(P.x != Q.x || P.y != Q.y);
if (P.x == Q.x) {
return {1, 0, -P.x};
} else if (P.y == Q.y) {
return {0, 1, -P.y};
} else {
ll a = P.y - Q.y;
ll b = Q.x - P.x;
ll c = P % Q;
ll g = gcd(gcd(a, b), c);
a /= g, b /=g, c /= g;
if (a > 0) {
return {a, b, c};
} else {
return {-a, -b, -c};
}
}
}
// 垂直二等分線
LineD getPerpendicularBisector(const SegmentLineD &L) {
ld a = 2.0 * (L.q.x - L.p.x);
ld b = 2.0 * (L.q.y - L.p.y);
ld c = L.p.x * L.p.x - L.q.x * L.q.x + L.p.y * L.p.y - L.q.y * L.q.y;
return LineD(a, b, c);
}
LineD getPerpendicularBisector(const PointD &p, const PointD &q) {
ld a = 2.0 * (q.x - p.x);
ld b = 2.0 * (q.y - p.y);
ld c = p.x * p.x - q.x * q.x + p.y * p.y - q.y * q.y;
return LineD(a, b, c);
}
// 2直線の交点算出
PointD intersectingPoint(const LineD &L1, const LineD &L2) {
ld D = L1.a * L2.b - L1.b * L2.a;
assert(abs(D) > eps);
return PointD{(-L1.c * L2.b + L1.b * L2.c) / D, (-L1.a * L2.c + L1.c * L2.a) / D};
}
// 外接円の中心点と半径を取得
pair<PointD, ld> circumcircle(const PointD &p1, const PointD &p2, const PointD &p3) {
LineD L1 = getPerpendicularBisector(p1, p2);
LineD L2 = getPerpendicularBisector(p2, p3);
PointD center = intersectingPoint(L1, L2);
ld r = (center - p1).norm();
return make_pair(center, r);
}
class Affine {
public:
vvd T_;
Affine()
: T_(3, vd(3)) {}
Affine(const vvd &a)
: T_(a) {}
Affine(ld v)
: T_(3, vd(3)) {
rep (i, 3) T_[i][i] = v;
}
ld det() const {
ld a = T_[1][1] * T_[2][2] - T_[1][2] * T_[2][1];
ld b = T_[1][0] * T_[2][2] - T_[1][2] * T_[2][0];
ld c = T_[1][0] * T_[2][1] - T_[1][1] * T_[2][0];
return a * T_[0][0] - b * T_[0][1] + c * T_[0][2];
}
bool isRegular() const { return det() != 0.0; }
Affine transpose() const {
Affine c;
rep (i, 3)
rep (j, 3) c.T_[i][j] = T_[j][i];
return c;
}
Affine inv() const {
ld D = det();
assert(abs(D) > eps);
Affine c;
c.T_[0][0] = (T_[1][1] * T_[2][2] - T_[2][1] * T_[1][2]) / D;
c.T_[1][0] = -(T_[1][0] * T_[2][2] - T_[2][0] * T_[1][2]) / D;
c.T_[2][0] = (T_[1][0] * T_[2][1] - T_[2][0] * T_[1][1]) / D;
c.T_[0][1] = -(T_[0][1] * T_[2][2] - T_[2][1] * T_[0][2]) / D;
c.T_[1][1] = (T_[0][0] * T_[2][2] - T_[2][0] * T_[0][2]) / D;
c.T_[2][1] = -(T_[0][0] * T_[2][1] - T_[2][0] * T_[0][1]) / D;
c.T_[0][2] = (T_[0][1] * T_[1][2] - T_[1][1] * T_[0][2]) / D;
c.T_[1][2] = -(T_[0][0] * T_[1][2] - T_[1][0] * T_[0][2]) / D;
c.T_[2][2] = (T_[0][0] * T_[1][1] - T_[1][0] * T_[0][1]) / D;
return c;
}
Affine rot(ld deg) const {
ld d = fmodl(fmodl(deg, 360.0) + 360.0, 360.0);
Affine T;
if (d == 0.0) {
return *this;
} else if (d == 90.0) {
T = Affine({
{0.0, -1.0, 0.0},
{1.0, 0.0, 0.0},
{0.0, 0.0, 1.0}
});
} else if (d == 180.0) {
T = Affine({
{-1.0, 0.0, 0.0},
{ 0.0, -1.0, 0.0},
{ 0.0, 0.0, 1.0}
});
} else if (d == 270.0) {
T = Affine({
{ 0.0, 1.0, 0.0},
{-1.0, 0.0, 0.0},
{ 0.0, 0.0, 1.0}
});
} else {
ld rad = deg2Rad(deg);
T = Affine({
{cosl(rad), -sinl(rad), 0.0},
{ sin(rad), cosl(rad), 0.0},
{ 0.0, 0.0, 1.0}
});
}
Affine c;
rep (i, 3)
rep (j, 3)
rep (k, 3) c.T_[i][j] += T.T_[i][k] * T_[k][j];
return c;
}
Affine move(ld tx, ld ty) const {
Affine c = *this;
c.T_[0][2] += tx;
c.T_[1][2] += ty;
return c;
}
Affine move(const Point &t) const {
Affine c = *this;
c.T_[0][2] += t.x;
c.T_[1][2] += t.y;
return c;
}
Affine move(const PointD &t) const {
Affine c = *this;
c.T_[0][2] += t.x;
c.T_[1][2] += t.y;
return c;
}
};
bool operator==(const Affine &a, const Affine &b) {
Affine c;
rep (i, 3)
rep (j, 3)
if (a.T_[i][j] != b.T_[i][j]) return false;
return true;
}
Affine operator+(const Affine &a, const Affine &b) {
Affine c;
rep (i, 3)
rep (j, 3) c.T_[i][j] = a.T_[i][j] + b.T_[i][j];
return c;
}
Affine operator-(const Affine &a, const Affine &b) {
Affine c;
rep (i, 3)
rep (j, 3) c.T_[i][j] = a.T_[i][j] - b.T_[i][j];
return c;
}
Affine operator-(const Affine &a) {
Affine c;
rep (i, 3)
rep (j, 3) c.T_[i][j] = -a.T_[i][j];
return c;
}
Affine operator*(ld t, const Affine &b) {
Affine c;
rep (i, 3)
rep (j, 3) c.T_[i][j] = t * b.T_[i][j];
return c;
}
PointD operator*(Affine T, const PointD &p) {
PointD q;
q.x = T.T_[0][0] * p.x + T.T_[0][1] * p.y + T.T_[0][2];
q.y = T.T_[1][0] * p.x + T.T_[1][1] * p.y + T.T_[1][2];
return q;
}
Affine operator/(const Affine &a, ld t) {
assert(t != 0.0);
Affine c;
rep (i, 3)
rep (j, 3) c.T_[i][j] = a.T_[i][j] / t;
return c;
}
Affine operator*(const Affine &a, const Affine &b) {
Affine c;
rep (i, 3)
rep (j, 3)
rep (k, 3) c.T_[i][j] += a.T_[i][k] * b.T_[k][j];
return c;
}
Affine operator/(const Affine &a, const Affine &b) {
return a * b.inv();
}
ld norm(const PointD &d) { return d.norm(); }
ld norm(const Point &d) { return d.norm(); }
// 直線L1, L2の交差判定
bool JudgeLineIntersection(const LineD &L1, const LineD &L2) {
return abs(L1.a * L2.b - L1.b * L2.a) < eps;
}
bool JudgeLineIntersection(const Line &L1, const Line &L2) {
return abs(L1.a * L2.b - L1.b * L2.a) < eps;
}
// 直線L1, 線分L2の交差判定
bool JudgeLineIntersection(const LineD &L1, const SegmentLineD &L2) {
PointD d = PointD(L1.a, L1.b);
ld v1 = d * L2.p, v2 = d * L2.q;
ld tmp = v2 - v1;
if (tmp < 0.0) {
return IN(tmp, v1 + L1.c, 0.0);
} else {
return IN(0.0, v1 + L1.c, tmp);
}
}
// 直線L1, 線分L2の交差判定
bool JudgeLineIntersection(const Line &L1, const SegmentLine &L2) {
Point d = {L1.a, L1.b};
ll v1 = d * L2.p, v2 = d * L2.q;
ll tmp = v2 - v1;
if (tmp <= 0) {
return IN(tmp, v1 + L1.c, 0);
} else {
return IN(0, v1 + L1.c, tmp);
}
}
// 点pと線分Lとの最短距離
ld SegmentlinePointDist(const PointD &p, const SegmentLineD &L) {
PointD d = L.dir();
PointD pp = p - L.p;
PointD pq = p - L.q;
PointD rd = -d;
if (d * pp < 0) {
return sqrt(pp * pp);
} else if (rd * pq < 0) {
return sqrt(pq * pq);
} else {
return abs(d % pp) / d.norm();
}
}
// 点pと線分Lとの最短距離
ld SegmentlinePointDist(const Point &p, const SegmentLine &L) {
Point d = L.dir();
Point pp = p - L.p;
Point pq = p - L.q;
Point rd = -d;
if (d * pp < 0) {
return sqrt(pp * pp);
} else if (rd * pq < 0) {
return sqrt(pq * pq);
} else {
return abs(d % pp) / d.norm();
}
}
// 線分L1, L2の交差判定
bool JudgeSegmentLineIntersection(const SegmentLineD &L1, const SegmentLineD &L2) {
PointD vec1 = L2.p - L1.p;
PointD vec2 = L2.q - L1.p;
ld c1 = L1.dir() % vec1;
ld c2 = L1.dir() % vec2;
if (c1 * c2 > 0.0) return false;
PointD vec3 = L1.p - L2.p;
PointD vec4 = L1.q - L2.p;
ld c3 = L2.dir() % vec3;
ld c4 = L2.dir() % vec4;
if (c3 * c4 > 0.0) return false;
if (c1 == 0 || c2 == 0 || c3 == 0 || c4 == 0) {
return SegmentlinePointDist(L2.p, L1) < eps || SegmentlinePointDist(L2.q, L1) < eps
|| SegmentlinePointDist(L1.p, L2) < eps || SegmentlinePointDist(L1.q, L2) < eps;
} else {
return true;
}
}
// 線分L1, L2の交差判定
bool JudgeSegmentLineIntersection(const SegmentLine &L1, const SegmentLine &L2) {
Point vec1 = L2.p - L1.p;
Point vec2 = L2.q - L1.p;
ll c1 = L1.dir() % vec1;
ll c2 = L1.dir() % vec2;
if (c1 * c2 > 0) return false;
Point vec3 = L1.p - L2.p;
Point vec4 = L1.q - L2.p;
ll c3 = L2.dir() % vec3;
ll c4 = L2.dir() % vec4;
if (c3 * c4 > 0) return false;
if (c1 == 0 || c2 == 0 || c3 == 0 || c4 == 0) {
return SegmentlinePointDist(L2.p, L1) < eps || SegmentlinePointDist(L2.q, L1) < eps
|| SegmentlinePointDist(L1.p, L2) < eps || SegmentlinePointDist(L1.q, L2) < eps;
} else {
return true;
}
}
// 三角形をなしているか
bool isTriangle(const PointD &P1, const PointD &P2, const PointD &P3) {
SegmentLineD L12 = {P1, P2};
SegmentLineD L13 = {P1, P3};
return abs(L12.dir() % L13.dir()) > 0;
}
// 3点が一直線上に並んでいるかどうか
bool isSameLine(const PointD &P1, const PointD &P2, const PointD &P3) {
return !isTriangle(P1, P2, P3);
}
// 三角形をなしているか
bool isTriangle(const Point &P1, const Point &P2, const Point &P3) {
SegmentLine L12 = {P1, P2};
SegmentLine L13 = {P1, P3};
return abs(L12.dir() % L13.dir()) > 0;
}
// 3点が一直線上に並んでいるかどうか
bool isSameLine(const Point &P1, const Point &P2, const Point &P3) {
return !isTriangle(P1, P2, P3);
}
bool existPointOnLine(const Point &P, const Line &L) {
return L.a * P.x + L.b * P.y + L.c == 0;
}
// 凸多角形か。頂点集合は反時計周りで設定
bool isConvexPolygon(const vector<Point> &ps) {
ll N = len(ps);
reps (i, 1, N + 1) {
Point p = ps[i - 1] - ps[i % N];
Point q = ps[(i + 1) % N] - ps[i % N];
if (q % p < 0) return false;
}
return true;
}
// 凸多角形か。頂点集合は反時計周りで設定
bool isConvexPolygon(const vector<PointD> &ps) {
ll N = len(ps);
reps (i, 1, len(ps) + 1) {
PointD p = ps[i - 1] - ps[i % N];
PointD q = ps[(i + 1) % N] - ps[i % N];
if (q % p < 0) return false;
}
return true;
}
enum class IsIn : ll {
ON = (ll)-1,
OUT = (ll)0,
IN = (ll)1
};
// 凸多角形の中に点が存在しているか
// O(logN)
// (参考https://atcoder.jp/contests/abc296/submissions/40217933)
IsIn IsPointInConvexPolygon(const vector<Point> &poly, const Point &pt) {
if (poly.size() == 1) return pt == poly[0] ? IsIn::ON : IsIn::OUT;
if (poly.size() == 2) {
SegmentLine l;
l.p = poly[0], l.q = poly[1];
return SegmentlinePointDist(pt, l) < eps ? IsIn::ON : IsIn::OUT;
}
if (pt == poly[0]) return IsIn::ON;
if ((poly[1] - poly[0]).toleft(pt - poly[0]) == -_1 || (poly.back() - poly[0]).toleft(pt - poly[0]) == _1) return IsIn::OUT;
const auto cmp = [&](const Point &u, const Point &v) { return (u - poly[0]).toleft(v - poly[0]) == _1; };
const size_t i = lower_bound(poly.begin() + 1, poly.end(), pt, cmp) - poly.begin();
SegmentLine l;
l.p = poly[0], l.q = poly[i];
if (i == 1) return SegmentlinePointDist(pt, l) < eps ? IsIn::ON : IsIn::OUT;
if (i == poly.size() - 1 && SegmentlinePointDist(pt, l) < eps) return IsIn::ON;
SegmentLine l1;
l1.p = poly[i - 1], l1.q = poly[i];
if (SegmentlinePointDist(pt, l1) < eps) return IsIn::ON;
return (poly[i] - poly[i - 1]).toleft(pt - poly[i - 1]) > 0 ? IsIn::IN : IsIn::OUT;
return IsIn::ON;
}
// 凸でない多角形に対して、点の内外判定を行なう。
// O(len(poly))でpolyは線を描く順に点を並べる。終点は閉じてなくてよくて
// 反時計周りでも時計周りでもOK。
IsIn IsPointInPolygon(const vector<Point> &poly, const Point &pt) {
ll N = len(poly);
if (N == 1) {
return poly[0] == pt ? IsIn::ON : IsIn::OUT;
}
if (N == 2) {
auto l = SegmentLine{poly[0], poly[1]};
return SegmentlinePointDist(pt, l) < eps ? IsIn::ON : IsIn::OUT;
}
complex<ld> AB = complex<ld>{(ld)pt.x, (ld)pt.y};
vector<complex<ld>> XY;
rep(i, len(poly)) XY.pb(complex<ld>{(ld)poly[i].x, (ld)poly[i].y});
complex<ld> r = complex<ld>{0.1, 0.01};
vector<PointD> pts;
rep(i, N) {
XY[i] -= AB;
XY[i] *= r;
pts.pb(PointD{XY[i].real(), XY[i].imag()});
}
ll xmcnt = 0;
rep(i, N) {
if (pts[i].x * pts[(i + 1) % N].x > 0.0) continue;
auto l = getLineD(pts[i], pts[(i + 1) % N]);
ld b = -l.c / l.b;
if (abs(b) < eps * eps) return IsIn::ON; // epsで判定すると誤判定するので、もっと0に近いときに線上に載ってる判定する。
if (b < 0.0) xmcnt++;
}
return xmcnt % 2 == 0 ? IsIn::OUT : IsIn::IN;
}
// 2,3,4が2円交差してる。(中心x, y, 半径rを入力する)
// [1] 一方の円が他方の円を完全に含み、2 つの円は接していない
// [2] 一方の円が他方の円を完全に含み、2 つの円は接している
// [3] 2 つの円が互いに交差する
// [4] 2 つの円の内部に共通部分は存在しないが、2 つの円は接している
// [5] 2 つの円の内部に共通部分は存在せず、2 つの円は接していない
ll twoCirclesState(const vll &XYR1, const vll &XYR2) {
INI3(x1, y1, r1, XYR1);
INI3(x2, y2, r2, XYR2);
if ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) < (r1 - r2) * (r1 - r2)) return _1;
if ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) == (r1 - r2) * (r1 - r2)) return _2;
if ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) < (r1 + r2) * (r1 + r2)) return _3;
if ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) == (r1 + r2) * (r1 + r2)) return _4;
if ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) > (r1 + r2) * (r1 + r2)) return _5;
return _0;
}
// 偏角ソートのcomp関数。(0, 0)は未定義
// p1の偏角 < p2の偏角
bool LessInArg(const PointD &p1, const PointD &p2) {
if (p1.y == 0.0 && p2.y == 0.0) return p1.x > 0.0 && p2.x < 0.0;
else if (p1.y >= 0.0 && p2.y < 0.0) return true;
else if (p1.y < 0.0 && p2.y >= 0.0) return false;
else return p1.x * p2.y > p2.x * p1.y;
}
// 偏角ソートのcomp関数。(0, 0)は未定義
// p1の偏角 < p2の偏角
bool LessInArg(const Point &p1, const Point &p2) {
if (p1.y == 0 && p2.y == 0) return p1.x > 0 && p2.x < 0;
else if (p1.y >= 0 && p2.y < 0) return true;
else if (p1.y < 0 && p2.y >= 0) return false;
else if (p1.x * p2.y > p2.x * p1.y) return true;
else if (p1.x * p2.y < p2.x * p1.y) return false;
else return p1.x * p1.x + p1.y * p1.y >= p2.x * p2.x + p2.y * p2.y;
}
// 偏角ソートのcomp関数。(0, 0)は未定義
// p1の偏角 < p2の偏角
template <typename T>
bool LessInArg(const complex<T> &p1, const complex<T> &p2) {
if (p1.imag() == 0 && p2.imag() == 0) return p1.real() > 0 && p2.real() < 0;
else if (p1.imag() >= 0 && p2.imag() < 0) return true;
else if (p1.imag() < 0 && p2.imag() >= 0) return false;
else if (p1.real() * p2.imag() > p2.real() * p1.imag()) return true;
else if (p1.real() * p2.imag() < p2.real() * p1.imag()) return false;
else return p1.real() * p1.real() + p1.imag() * p1.imag() >= p2.real() * p2.real() + p2.imag() * p2.imag();
}
// 偏角ソート(引数を直接ソートする)
void SORTARG(vector<Point> &ps, bool reverse = false) {
auto cmp = [&](Point &x, Point &y) {
return LessInArg(x, y) ^ reverse;
};
sort(all(ps), cmp);
}
// 偏角ソート(引数を直接ソートする)
void SORTARG(vector<PointD> &ps, bool reverse = false) {
auto cmp = [&](PointD &x, PointD &y) {
return LessInArg(x, y) ^ reverse;
};
sort(all(ps), cmp);
}
// 偏角ソート(引数を直接ソートする)
template <typename T>
void SORTARG(vector<complex<T>> &ps, bool reverse = false) {
auto cmp = [&](complex<T> &x, complex<T> &y) {
return LessInArg(x, y) ^ reverse;
};
sort(all(ps), cmp);
}
template <typename T>
bool isSamePointsAfterRotationMove(const vector<complex<T>> &ps, const vector<complex<T>> &qs) {
assert(len(ps) == len(qs));
if (len(ps) == 1) return true;
vector<complex<T>> PS, PT;
const ll N = len(ps);
rep(i, N) PS.pb(complex<T>{ps[i].real() * N, ps[i].imag() * N}), PT.pb(complex<T>{qs[i].real() * N, qs[i].imag() * N});
complex<ld> cs = SUM(PS), ct = SUM(PT);
cs /= N, ct /= N;
rep(i, N) PS[i] -= cs, PT[i] -= ct;
SORTARG(PT);
ll baseidx = 0;
ld val = 0.0;
rep(i, N) if (val < abs(PT[i])) val = abs(PT[i]), baseidx = i;
rep(i, N) {
if (abs(abs(PS[i]) - abs(PT[baseidx])) > eps) continue;
auto z = PT[baseidx] / PS[i];
vector<complex<ld>> tmp;
rep(j, N) tmp.pb(PS[j] * z);
SORTARG(tmp);
bool isok = true;
rep(j, N) if (abs(tmp[i] - PT[i]) > eps) isok = false;
if (isok) return true;
}
return false;
}
bool isSamePointsAfterRotationMove(const vector<Point> &ps, const vector<Point> &qs) {
assert(len(ps) == len(qs));
if (len(ps) == 1) return true;
vector<complex<ld>> cps, cqs;
rep(i, len(ps)) {
cps.pb(complex<ld>{(ld)ps[i].x, (ld)ps[i].y});
cqs.pb(complex<ld>{(ld)qs[i].x, (ld)qs[i].y});
}
return isSamePointsAfterRotationMove(cps, cqs);
}
// マトリックスのある文字cを点群の座標に変換
vector<Point> convMatrixToPoints(const vs &mat, char c) {
vector<Point> ps;
rep(i, len(mat)) {
rep(j, len(mat[i])) {
if (mat[i][j] == c) {
ps.pb(Point{i, j});
}
}
}
return ps;
}
// マトリックスのある値vを点群の座標に変換
vector<Point> convMatrixToPoints(const vvll &mat, ll v) {
vector<Point> ps;
rep(i, len(mat)) {
rep(j, len(mat[i])) {
if (mat[i][j] == v) {
ps.pb(Point{i, j});
}
}
}
return ps;
}
// ps, qsを平行移動させて点群をマージする。
vector<Point> mergePoints(const vector<Point> &ps, const Point &offset_ps, const vector<Point> &qs, const Point &offset_qs) {
set<Point> tbl;
rep(i, len(ps)) {
tbl.insert(ps[i] + offset_ps);
}
rep(i, len(qs)) {
tbl.insert(qs[i] + offset_qs);
}
return vector<Point>(all(tbl));
}
// 平行移動のみで点群が一致するか
bool isSimilarPoints(const vector<Point> &ps, const vector<Point> &qs) {
Point psminxy = {LLONG_MAX, LLONG_MAX}, qsminxy = {LLONG_MAX, LLONG_MAX};
rep(i, len(ps)) {
chmin(psminxy.x, ps[i].x);
chmin(psminxy.y, ps[i].y);
}
rep(i, len(qs)) {
chmin(qsminxy.x, qs[i].x);
chmin(qsminxy.y, qs[i].y);
}
set<Point> moveps, moveqs;
rep(i, len(ps)) moveps.insert(Point{ps[i].x - psminxy.x, ps[i].y - psminxy.y});
rep(i, len(qs)) moveqs.insert(Point{qs[i].x - qsminxy.x, qs[i].y - qsminxy.y});
return moveps == moveqs;
}
// unordered_mapでPoint, PointD, Line, LineD, SegmentLine, SegmentLineDをkeyにするためのコード
/* PointD用 */
template <>
struct std::hash<PointD>{
size_t operator()(const PointD &keyval) const noexcept {
size_t s = HashCombine(0, keyval.x);
return HashCombine(s, keyval.y);
}
};
/* Point用 */
template <>
struct std::hash<Point>{
size_t operator()(const Point &keyval) const noexcept {
size_t s = HashCombine(0, keyval.x);
return HashCombine(s, keyval.y);
}
};
/* LineD用 */
template <>
struct std::hash<LineD>{
size_t operator()(const LineD &keyval) const noexcept {
size_t s = HashCombine(0, keyval.a);
s = HashCombine(s, keyval.b);
return HashCombine(s, keyval.c);
}
};
/* Line用 */
template <>
struct std::hash<Line>{
size_t operator()(const Line &keyval) const noexcept {
size_t s = HashCombine(0, keyval.a);
s = HashCombine(s, keyval.b);
return HashCombine(s, keyval.c);
}
};
/* SegmentLineD用 */
template <>
struct std::hash<SegmentLineD>{
size_t operator()(const SegmentLineD &keyval) const noexcept {
size_t s = HashCombine(0, keyval.p.x);
s = HashCombine(s, keyval.p.y);
s = HashCombine(s, keyval.q.x);
return HashCombine(s, keyval.q.y);
}
};
/* SegmentLine用 */
template <>
struct std::hash<SegmentLine>{
size_t operator()(const SegmentLine &keyval) const noexcept {
size_t s = HashCombine(0, keyval.p.x);
s = HashCombine(s, keyval.p.y);
s = HashCombine(s, keyval.q.x);
return HashCombine(s, keyval.q.y);
}
};
/* func */
inline ll in_ll() {string s; getline(cin, s); return STRLL(s);}
inline ld in_d() {string s; getline(cin, s); return STRD(s);}
inline string in_str() {string s; getline(cin, s); return s;}
/* debug */
namespace debug_print_func {
std::ostream& os = std::cout;
template <class Tp> auto has_cbegin(int) -> decltype(std::cbegin(std::declval<Tp>()), std::true_type {});
template <class Tp> auto has_cbegin(...) -> std::false_type;
template <class Tp> auto has_value_type(int) -> decltype(std::declval<typename Tp::value_type>(), std::true_type {});
template <class Tp> auto has_value_type(...) -> std::false_type;
template <class Tp>[[maybe_unused]] constexpr bool is_iteratable_container_v = decltype(has_cbegin<Tp>(int {}))::value;
template <class Tp>[[maybe_unused]] constexpr bool is_container_v = decltype(has_value_type<Tp>(int {}))::value
|| is_iteratable_container_v<Tp>;
template <> [[maybe_unused]] constexpr bool is_iteratable_container_v<std::string_view> = false;
template <> [[maybe_unused]] constexpr bool is_container_v<std::string_view> = false;
#if (defined _GLIBCXX_STRING) || (defined _LIBCPP_STRING)
template <> [[maybe_unused]] constexpr bool is_iteratable_container_v<std::string> = false;
template <> [[maybe_unused]] constexpr bool is_container_v<std::string> = false;
#endif
template <class Tp, class... Ts> struct first_element { using type = Tp; };
template <class... Ts> using first_t = typename first_element<Ts...>::type;
template <class Tp, std::enable_if_t<!decltype(has_value_type<Tp>(int {}))::value, std::nullptr_t> = nullptr>
auto check_elem(int) -> decltype(*std::cbegin(std::declval<Tp>()));
template <class Tp, std::enable_if_t<decltype(has_value_type<Tp>(int {}))::value, std::nullptr_t> = nullptr>
auto check_elem(int) -> typename Tp::value_type;
template <class Tp>
auto check_elem(...) -> void;
template <class Tp> using elem_t = decltype(check_elem<Tp>(int {}));
template <class Tp> [[maybe_unused]] constexpr bool is_multidim_container_v = is_container_v<Tp>
&& is_container_v<elem_t<Tp>>;
template <class Tp> std::enable_if_t<!is_container_v<Tp>> out(const Tp&);
void out(const char&);
void out(const char*);
void out(const std::string_view&);
#if (defined _GLIBCXX_STRING) || (defined _LIBCPP_STRING)
void out(const std::string&);
#endif
#ifdef __SIZEOF_INT128__
void out(const __int128&);
void out(const unsigned __int128&);
#endif
template <class Tp1, class Tp2> void out(const std::pair<Tp1, Tp2>&);
#if (defined _GLIBCXX_TUPLE) || (defined _LIBCPP_TUPLE)
template <class... Ts> void out(const std::tuple<Ts...>&);
#endif
#if (defined _GLIBCXX_STACK) || (defined _LIBCPP_STACK)
template <class... Ts> void out(std::stack<Ts...>);
#endif
#if (defined _GLIBCXX_QUEUE) || (defined _LIBCPP_QUEUE)
template <class... Ts> void out(std::queue<Ts...>);
template <class... Ts> void out(std::priority_queue<Ts...>);
#endif
template <class C>
std::enable_if_t<is_iteratable_container_v<C>> out(const C&);
template <class Tp> std::enable_if_t<!is_container_v<Tp>> out(const Tp& arg) {
os << arg;
}
void out(const char& arg) {
os << arg;
}
void out(const char* arg) {
os << arg;
}
void out(const ld arg) {
if (arg == LDBL_MAX) {
#ifdef LOCAL
os << "∞";
#else
os << arg;
#endif
} else if (arg == -LDBL_MAX) {
#ifdef LOCAL
os << "-∞";
#else
os << arg;
#endif
} else {
os << arg;
}
}
template <typename T>
void out(const std::complex<T>& arg) {
os << arg.real() << " + " << arg.imag() << "i";
}
void out(const std::string_view& arg) {
os << arg;
}
#if (defined _GLIBCXX_STRING) || (defined _LIBCPP_STRING)
void out(const std::string& arg) {
os << arg;
}
#endif
#if (defined _GLIBCXX_STRING) || (defined _LIBCPP_STRING)
void out(const vs& arg) {
rep(i, len(arg)) {
out(arg[i]);
if (i != len(arg) - 1) cout << endl;
}
}
#endif
#ifdef __SIZEOF_INT128__
void out(const __int128& arg) {
if (arg == ULLONG_MAX) {
#ifdef LOCAL
os << "∞";
#else
int sign = (arg < 0) ? (-1) : 1;
if (sign == -1) os << '-';
__int128 base = sign;
while (sign * arg >= sign * base * 10) base *= 10;
while (base) {
os << static_cast<char>('0' + (arg / base % 10));
base /= 10;
}
#endif
} else {
int sign = (arg < 0) ? (-1) : 1;
if (sign == -1) os << '-';
__int128 base = sign;
while (sign * arg >= sign * base * 10) base *= 10;
while (base) {
os << static_cast<char>('0' + (arg / base % 10));
base /= 10;
}
}
}
void out(const unsigned __int128& arg) {
if (arg == ULLONG_MAX) {
#ifdef LOCAL
os << "∞";
#else
unsigned __int128 base = 1;
while (arg >= base * 10) base *= 10;
while (base) {
os << static_cast<char>('0' + (arg / base % 10));
base /= 10;
}
#endif
} else {
unsigned __int128 base = 1;
while (arg >= base * 10) base *= 10;
while (base) {
os << static_cast<char>('0' + (arg / base % 10));
base /= 10;
}
}
}
#endif
void out(const mint &arg) {
out(arg.val());
}
void out(const Point& arg) {
out(arg.x);
os << " ";
out(arg.y);
}
void out(const SegmentLine& arg) {
out(arg.p);
os << " ";
out(arg.q);
}
template <class Tp1, class Tp2> void out(const std::pair<Tp1, Tp2>& arg) {
out(arg.first);
os << " ";
out(arg.second);
}
#if (defined _GLIBCXX_TUPLE) || (defined _LIBCPP_TUPLE)
template <class T, std::size_t... Is> void print_tuple(const T& arg, std::index_sequence<Is...>) {
static_cast<void>(((os << (Is == 0 ? "" : ", "), out(std::get<Is>(arg))), ...));
}
template <class... Ts> void out(const std::tuple<Ts...>& arg) {
print_tuple(arg, std::make_index_sequence<sizeof...(Ts)>());
}
#endif
#if (defined _GLIBCXX_STACK) || (defined _LIBCPP_STACK)
template <class... Ts> void out(std::stack<Ts...> arg) {
if (arg.empty()) {
os << "<empty stack>";
return;
}
while (!arg.empty()) {
out(arg.top());
os << ' ';
arg.pop();
}
}
#endif
#if (defined _GLIBCXX_QUEUE) || (defined _LIBCPP_QUEUE)
template <class... Ts> void out(std::queue<Ts...> arg) {
if (arg.empty()) {
os << "<empty queue>";
return;
}
while (!arg.empty()) {
out(arg.front());
os << ' ';
arg.pop();
}
}
template <class... Ts> void out(std::priority_queue<Ts...> arg) {
if (arg.empty()) {
os << "<empty priority_queue>";
return;
}
while (!arg.empty()) {
out(arg.top());
os << ' ';
arg.pop();
}
}
#endif
template <class Container>
std::enable_if_t<is_iteratable_container_v<Container>> out(const Container& arg) {
if (std::distance(std::cbegin(arg), std::cend(arg)) == 0) {
os << "<empty container>";
return;
}
std::for_each(std::cbegin(arg), std::cend(arg), [](const elem_t<Container>& elem) {
out(elem);
os << ' ';
});
}
template <class Tp> std::enable_if_t<!is_multidim_container_v<Tp>>
print(std::string_view name, const Tp& arg) {
os << name << ": ";
out(arg);
if constexpr (is_container_v<Tp>)
os << '\n';
}
template <class Tp> std::enable_if_t<is_multidim_container_v<Tp>>
print(std::string_view name, const Tp& arg) {
os << name << ": ";
if (std::distance(std::cbegin(arg), std::cend(arg)) == 0) {
os << "<empty multidimensional container>\n";
return;
}
std::for_each(std::cbegin(arg), std::cend(arg),
[&name, is_first_elem = true](const elem_t<Tp>& elem) mutable {
if (is_first_elem)
is_first_elem = false;
else
for (std::size_t i = 0; i < name.length() + 2; i++)
os << ' ';
out(elem);
os << '\n';
});
}
template <class Tp, class... Ts> void multi_print(std::string_view names, const Tp& arg, const Ts&... args) {
if constexpr (sizeof...(Ts) == 0) {
names.remove_suffix(
std::distance(
names.crbegin(),
std::find_if_not(names.crbegin(), names.crend(),
[](const char c) { return std::isspace(c); })
)
);
print(names, arg);
if constexpr (!is_container_v<Tp>)
os << '\n';
} else {
std::size_t comma_pos = 0;
for (std::size_t i = 0, paren_depth = 0, inside_quote = false; i < names.length(); i++) {
if (!inside_quote && paren_depth == 0 && i > 0 && names[i - 1] != '\'' && names[i] == ',') {
comma_pos = i;
break;
}
if (names[i] == '\"') {
if (i > 0 && names[i - 1] == '\\') continue;
inside_quote ^= true;
}
if (!inside_quote && names[i] == '(' && (i == 0 || names[i - 1] != '\''))
paren_depth++;
if (!inside_quote && names[i] == ')' && (i == 0 || names[i - 1] != '\''))
paren_depth--;
}
const std::size_t first_varname_length = comma_pos - std::distance(
names.crend() - comma_pos,
std::find_if_not(
names.crend() - comma_pos, names.crend(),
[](const char c) { return std::isspace(c); }
)
);
print(names.substr(0, first_varname_length), arg);
if constexpr (!is_container_v<Tp>) {
if constexpr (is_container_v<first_t<Ts...>>)
os << '\n';
else
os << " | ";
}
const std::size_t next_varname_begins_at = std::distance(
names.cbegin(),
std::find_if_not(
names.cbegin() + comma_pos + 1, names.cend(),
[](const char c) { return std::isspace(c); }
)
);
names.remove_prefix(next_varname_begins_at);
multi_print(names, args...);
}
}
template <class Tp> std::enable_if_t<!is_multidim_container_v<Tp>>
print_out(std::string_view name, const Tp& arg) {
UNUSED(name);
out(arg);
if constexpr (is_container_v<Tp>)
os << '\n';
}
template <class Tp> std::enable_if_t<is_multidim_container_v<Tp>>
print_out(std::string_view name, const Tp& arg) {
if (std::distance(std::cbegin(arg), std::cend(arg)) == 0) {
os << "<empty multidimensional container>\n";
return;
}
std::for_each(std::cbegin(arg), std::cend(arg),
[&name, is_first_elem = true](const elem_t<Tp>& elem) mutable {
if (is_first_elem)
is_first_elem = false;
else
for (std::size_t i = 0; i < name.length() + 2; i++)
os << ' ';
out(elem);
os << '\n';
});
}
template <class Tp, class... Ts> void _print(std::string_view names, const Tp& arg, const Ts&... args) {
if constexpr (sizeof...(Ts) == 0) {
names.remove_suffix(
std::distance(
names.crbegin(),
std::find_if_not(names.crbegin(), names.crend(),
[](const char c) { return std::isspace(c); })
)
);
print_out(names, arg);
if constexpr (!is_container_v<Tp>)
os << '\n';
} else {
std::size_t comma_pos = 0;
for (std::size_t i = 0, paren_depth = 0, inside_quote = false; i < names.length(); i++) {
if (!inside_quote && paren_depth == 0 && i > 0 && names[i - 1] != '\'' && names[i] == ',') {
comma_pos = i;
break;
}
if (names[i] == '\"') {
if (i > 0 && names[i - 1] == '\\') continue;
inside_quote ^= true;
}
if (!inside_quote && names[i] == '(' && (i == 0 || names[i - 1] != '\''))
paren_depth++;
if (!inside_quote && names[i] == ')' && (i == 0 || names[i - 1] != '\''))
paren_depth--;
}
const std::size_t first_varname_length = comma_pos - std::distance(
names.crend() - comma_pos,
std::find_if_not(
names.crend() - comma_pos, names.crend(),
[](const char c) { return std::isspace(c); }
)
);
print_out(names.substr(0, first_varname_length), arg);
if constexpr (!is_container_v<Tp>) {
if constexpr (is_container_v<first_t<Ts...>>)
os << '\n';
else
os << " ";
}
const std::size_t next_varname_begins_at = std::distance(
names.cbegin(),
std::find_if_not(
names.cbegin() + comma_pos + 1, names.cend(),
[](const char c) { return std::isspace(c); }
)
);
names.remove_prefix(next_varname_begins_at);
_print(names, args...);
}
}
} // namespace debug_print
#define print(...) debug_print_func::_print(#__VA_ARGS__, __VA_ARGS__);
#ifdef LOCAL
# define debug(...) cerr << "\033[33m(line:" << __LINE__ << ") " << endl; debug_print_func::multi_print(#__VA_ARGS__, __VA_ARGS__); cerr << "\033[m";
#else
# define debug(...) ;
#endif
/* 標準入力 */
vs in_strs(const string &delimiter = " ")
{
string s;
getline(cin, s);
vs output;
bitset<255> delims;
for (unsigned char c: delimiter)
{
delims[c] = true;
}
string::const_iterator beg;
bool in_token = false;
for( string::const_iterator it = s.cbegin(), end = s.cend(); it != end; ++it )
{
if( delims[*it] )
{
if( in_token )
{
output.pb(beg, it);
in_token = false;
}
}
else if( !in_token )
{
beg = it;
in_token = true;
}
}
if( in_token )
output.pb(beg, s.cend());
return output;
}
inline vll in_lls()
{
vll vals;
vs tokens = in_strs();
for (string i: tokens) vals.pb(STRLL(i));
return vals;
}
inline vd in_ds()
{
vd vals;
vs tokens = in_strs();
for (string i: tokens) vals.pb(STRD(i));
return vals;
}
inline vvll in_llss(ll line) // 複数行文字列解析
{
vvll valss;
rep(i, line) valss.pb(in_lls());
return valss;
}
inline vs in_vs(ll line) // 複数行文字列解析
{
vs vecs;
rep(i, line) {
vecs.pb(in_str());
}
return vecs;
}
inline ll popcnt(ll x) { return __builtin_popcountll(x); }
template <typename T, typename U>
T ceil(T x, U y) {
return (x > 0 ? (x + y - 1) / y : x / y);
}
template <typename T, typename U>
T floor(T x, U y) {
return (x > 0 ? x / y : (x - y + 1) / y);
}
template <typename T>
vector<T> accsum(const vector<T> &vec, bool need0 = true, ll ansmod = 0) {
if (len(vec) == 0) return vector<T>();
vector<T> acc = {0};
ll idx = 0;
if (!need0) {
acc[0] = ansmod == 0 ? vec[0] : vec[0] % ansmod;
idx = 1;
}
if (ansmod != 0) {
rep (i, idx, len(vec)) acc.pb((acc[len(acc) - 1] + vec[i]) % ansmod);
} else {
rep (i, idx, len(vec)) acc.pb(acc[len(acc) - 1] + vec[i]);
}
return acc;
}
template <typename T>
vector<T> accmax(const vector<T> &vec) {
if (len(vec) == 0) return vector<T>();
vector<T> acc = {vec[0]};
reps (i, 1, len(vec)) acc.pb(max(acc[i - 1], vec[i]));
return acc;
}
template <typename T>
vector<T> accmin(const vector<T> &vec) {
if (len(vec) == 0) return vector<T>();
vector<T> acc = {vec[0]};
reps (i, 1, len(vec)) acc.pb(min(acc[i - 1], vec[i]));
return acc;
}
template <typename T, typename U>
pair<T, T> divmod(T x, U y) {
T q = floor(x, y);
return {q, x - q * y};
}
inline ll sumk(ll n)
{
return n > 0 ? n * (n + 1) / 2 : 0;
}
inline ll sumk2(ll n)
{
return n > 0 ? n * (n + 1) * (2 * n + 1) / 6 : 0;
}
inline mint sumk(mint n)
{
return n * (n + 1) / 2;
}
inline mint sumk2(mint n)
{
return n * (n + 1) * (2 * n + 1) / 6;
}
ll accxor(ll n)
{
ll v = (n + 1) % 4;
if (v == 0)
return 0;
else if (v == 1)
return n;
else if (v == 2)
return 1;
else
return 1 ^ n;
}
inline string alpha()
{
return "abcdefghijklmnopqrstuvwxyz";
}
inline ll alpha_num(char c)
{
return ll(c) - ll('a');
}
inline string alpha_big()
{
return "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
}
inline ll alpha_big_num(char c)
{
return ll(c) - ll('A');
}
inline char alpha_big2small(char C) {
static string s = alpha();
ll idx = alpha_big_num(C);
return IN(0, idx, 25) ? s[idx] : C;
}
inline char alpha_small2big(char c) {
static string s = alpha_big();
ll idx = alpha_num(c);
return IN(0, idx, 25) ? s[idx] : c;
}
// 10進数の値Nをb進数で表したときの桁和。
ll digitsum(ll N, ll b) {
ll ret = 0;
while (N) {
ret += N % b;
N /= b;
}
return ret;
}
// 10進数文字列の各桁和
ll digitsum(const string &s) {
ll val = 0;
rep (i, len(s)) {
val += CHARLL(s[i]);
}
return val;
}
string zerofill(ll v, ll outputlen)
{
string s = STR(v);
string zerostr(outputlen - len(s), '0');
return zerostr + s;
}
// ランレングス圧縮
// auto rle = RunLengthEncoding(S);
// rep(i, len(rle)) {
// auto &[c, cnt] = rle[i];
// }
vector<pair<char, ll>> RunLengthEncoding(const string &s) {
vector<pair<char, ll>> tbl;
char c = s[0];
ll cnt = 1;
ll N = len(s);
reps (i, 1, N) {
if (c == s[i]) {
cnt++;
}
else {
tbl.pb(mp(c, cnt));
c = s[i];
cnt = 1;
}
}
tbl.pb(mp(c, cnt));
return tbl;
}
// ランレングス圧縮
// auto rle = RunLengthEncoding(S);
// rep(i, len(rle)) {
// auto &[c, cnt] = rle[i];
// }
template <typename T>
vector<pair<T, ll>> RunLengthEncoding(const vector<T> &s) {
vector<pair<T, ll>> tbl;
T c = s[0];
ll cnt = 1;
ll N = len(s);
reps (i, 1, N) {
if (c == s[i]) {
cnt++;
}
else {
tbl.pb(mp(c, cnt));
c = s[i];
cnt = 1;
}
}
tbl.pb(mp(c, cnt));
return tbl;
}
// 文字列連結
string repeatstr(const string &s, ll num) {
if (num == 0) return "";
string ret = "";
bitset<128> tmp = num;
bool isok = false;
repd (i, 128) {
if (!isok && tmp[i]) isok = true;
if (!isok) continue;
ret += ret;
if (tmp[i]) {
ret += s;
}
}
return ret;
}
// [lidx, ridx)の区間の文字列を取得 substr("0123456789", 2, 6) -> "2345"
// 第3引数は文字数ではない
string substr(const string &s, ll lidx, ll ridx) {
if (ridx <= lidx) return "";
return s.substr(lidx, ridx - lidx);
}
// 入力した文字列から入力文字列の順序を並び替えずに
// 辞書順最小のK文字の部分文字列を取得
string GetDictOrderMinimumSubStr(const string &S, ll K) {
const ll N = len(S);
const string alp = alpha();
vvll tbl(26, vll(N + 1, LLONG_MAX));
rep(i, 26) {
repd(j, N) {
if (S[j] == alp[i]) tbl[i][j] = j + 1;
else tbl[i][j] = tbl[i][j + 1];
}
}
ll neednum = K, idx = 0;
vector<char> ans;
while(neednum) {
rep(i, 26) {
if (neednum - 1 <= N - tbl[i][idx]) {
ans.pb(alp[i]);
idx = tbl[i][idx];
neednum--;
break;
}
}
}
return string(all(ans));
}
// // N = 200000, Q = 200000でchecknum = 2のとき200msほど
// // N = 200000, Q = 200000でchecknum = 16のとき1200msほど
// // checknumは1から16までで大きければ大きいほど一致しているの誤判定の確率が低くなる。
// class RollingHash {
// const vector<__int128> mods_ = {100000000000645093, 100000000000913821, 100000000001253899, 100000000001318983,
// 100000000001387039, 100000000001437331, 100000000001536913, 100000000001603269,
// 100000000001653931, 100000000001705759, 100000000001756759, 100000000001809399,
// 100000000001878349, 100000000001930267, 100000000001997103, 100000000002030859};
// vector<vector<__int128>> hashcodes_;
// const __int128 b_ = ((__int128)1 << (__int128)61) - 1;
// vector<vector<__int128>> bases_;
// const __int128 checknum_ = 2; // WAして衝突してそうなら、ここの数字を増やす。
// void build_(const vector<__int128> &v) {
// rep(i, checknum_) {
// vector<__int128> tmp(len(v) + 1);
// hashcodes_.pb(tmp);
// tmp[0] = 1;
// bases_.pb(tmp);
// }
// rep(i, len(v)) {
// rep(j, checknum_) {
// hashcodes_[j][i + 1] = (hashcodes_[j][i] * b_ + v[i]) % mods_[j];
// bases_[j][i + 1] = (bases_[j][i] * b_) % mods_[j];
// }
// }
// }
// __int128 getHashCode_(__int128 L, __int128 R, __int128 checkidx) {
// assert(IN(0, checkidx, checknum_ - 1));
// __int128 v = hashcodes_[checkidx][R] - (hashcodes_[checkidx][L] * bases_[checkidx][R - L] % mods_[checkidx]);
// if (v < 0) v = (v + mods_[checkidx]) % mods_[checkidx];
// return v;
// }
// __int128 connectionHashCode_(__int128 lHash, __int128 rHash, __int128 rstrlength, __int128 checkidx) {
// return (rHash + lHash * bases_[checkidx][rstrlength]) % mods_[checkidx];
// }
// public:
// RollingHash(const string &s) {
// vector<__int128> v;
// rep(i, len(s)) v.pb((__int128)s[i]);
// build_(v);
// }
// RollingHash(const vector<__int128> &v) {
// build_(v);
// }
// vector<__int128> getHashCode(__int128 L, __int128 R) {
// vector<__int128> hashcodes;
// rep(i, checknum_) {
// hashcodes.pb(getHashCode_(L, R, i));
// }
// return hashcodes;
// }
// vector<__int128> connectionHashCode(const vector<__int128> &lHashs, const vector<__int128> &rHashs, __int128 rstrlength) {
// vector<__int128> hashcodes;
// rep(i, len(lHashs)) {
// hashcodes.pb(connectionHashCode_(lHashs[i], rHashs[i], rstrlength, i));
// }
// return hashcodes;
// }
// // 区間[L1, R1), [L2, R2)が一致しているか
// bool equal(__int128 L1, __int128 R1, __int128 L2, __int128 R2) {
// return getHashCodes(L1, R1) == getHashCodes(L2, R2);
// }
// };
// 衝突するなら上を使う
class RollingHash {
using hash = pair<modint998244353, modint1000000007>;
vector<hash> hashcodes_;
const hash b_ = mp((1UL << 31) - 1, (1UL << 31) - 1);
vector<hash> bases_;
void build_(const vector<hash> &v) {
hashcodes_.reserve(len(v) + 1);
bases_.reserve(len(v) + 1);
bases_[0] = mp(1, 1);
rep(i, len(v)) {
hashcodes_[i + 1].first = hashcodes_[i].first * b_.first + v[i].first;
hashcodes_[i + 1].second = hashcodes_[i].second * b_.second + v[i].second;
bases_[i + 1].first = bases_[i].first * b_.first;
bases_[i + 1].second = bases_[i].second * b_.second;
}
}
public:
RollingHash(const string &s) {
vector<hash> v;
rep(i, len(s)) v.pb(mp(s[i], s[i]));
build_(v);
}
RollingHash(const vector<ll> &arr) {
vector<hash> v;
rep(i, len(arr)) v.pb(mp(arr[i], arr[i]));
build_(v);
}
hash getHashCode(ll L, ll R) {
return mp(hashcodes_[R].first - hashcodes_[L].first * bases_[R - L].first, hashcodes_[R].second - hashcodes_[L].second * bases_[R - L].second);
}
hash connectionHashCode(hash lHash, hash rHash, ll rlength) {
auto h1 = rHash.first + lHash.first * bases_[rlength].first;
auto h2 = rHash.second + lHash.second * bases_[rlength].second;
return mp(h1, h2);
}
// 区間[L1, R1), [L2, R2)が一致しているか
bool equal(ll L1, ll R1, ll L2, ll R2) {
return getHashCode(L1, R1) == getHashCode(L2, R2);
}
};
// 入力した文字列からK個異なる文字列で辞書順最小のものを取得する。O(len(S))
// ABC009-C参照
string GetDictOrderMinimumKthDiff(const string &S, ll K) {
ll N = len(S);
vll cnts(26);
rep(i, N) {
cnts[alpha_num(S[i])]++;
}
vll nowcnts = cnts;
string ans = "";
string al = alpha();
rep(i, N) {
rep(j, 26) {
if (nowcnts[j] > 0) {
cnts[alpha_num(S[i])]--, nowcnts[j]--;
ll samecnt = j == alpha_num(S[i]) ? 1 : 0;
rep(k, 26) {
samecnt += min(cnts[k], nowcnts[k]);
}
ll diffcnt = N - i - samecnt;
if (diffcnt <= K) {
ans += STR(al[j]);
if (j != alpha_num(S[i])) K--;
break;
}
cnts[alpha_num(S[i])]++, nowcnts[j]++;
}
}
}
return ans;
}
// 連続でない部分列xを用いて、両端から順に該当する文字を交換したときの
// 辞書順最小の文字列を返す。
// "dcab"のとき"acdb"が解。x={0,2}のときs[0]とs[2]が入れ替わる。
string GetDictOrderMinimumExchangeStr(const string &str) {
ll N = len(str);
string S = str;
vll hist(26);
rep(i, N) hist[alpha_num(S[i])]++;
ll histidx = 0;
ll l = 0, r = N - 1;
while (l < r) {
while(histidx < 26 && !hist[histidx]) histidx++;
if (histidx == 26) break;
if (alpha_num(S[l]) <= histidx) {
hist[alpha_num(S[l])]--;
l++;
continue;
}
while (l < r && alpha_num(S[r]) != histidx) {
hist[alpha_num(S[r])]--;
r--;
}
if (r <= l) break;
hist[alpha_num(S[l])]--;
hist[alpha_num(S[r])]--;
swap(S[l], S[r]);
l++, r--;
}
return S;
}
// join_strの文字列でvectorの要素を結合していく
template<class T>
string join(const vector<T> &vec, const string &join_str) {
string s = "";
if (len(vec) == 0) return s;
s += STR(vec[0]);
reps (i, 1, len(vec)) {
s += join_str + STR(vec[i]);
}
return s;
}
// 文字列検索(正規表現OK) O(|S|)っぽい
bool search_string(const string &S, const string ®ex_strkey) {
return std::regex_search(S, std::regex(regex_strkey));
}
// 文字列置換(正規表現OK)
string replace(const string &S, const string ®ex_strkey, const string &replacestr) {
return regex_replace(S, regex(regex_strkey), replacestr);
}
vll sa_naive(const vll& s) {
ll n = s.size();
vll sa(n);
std::iota(sa.begin(), sa.end(), 0);
std::sort(sa.begin(), sa.end(), [&](ll l, ll r) {
if (l == r) return false;
while (l < n && r < n) {
if (s[l] != s[r]) return s[l] < s[r];
l++;
r++;
}
return l == n;
});
return sa;
}
vll sa_doubling(const vll& s) {
ll n = s.size();
vll sa(n), rnk = s, tmp(n);
std::iota(sa.begin(), sa.end(), 0);
for (ll k = 1; k < n; k *= 2) {
auto cmp = [&](ll x, ll y) {
if (rnk[x] != rnk[y]) return rnk[x] < rnk[y];
ll rx = x + k < n ? rnk[x + k] : -1;
ll ry = y + k < n ? rnk[y + k] : -1;
return rx < ry;
};
std::sort(sa.begin(), sa.end(), cmp);
tmp[sa[0]] = 0;
for (ll i = 1; i < n; i++) {
tmp[sa[i]] = tmp[sa[i - 1]] + (cmp(sa[i - 1], sa[i]) ? 1 : 0);
}
std::swap(tmp, rnk);
}
return sa;
}
// SA-IS, linear-time suffix array construction
// Reference:
// G. Nong, S. Zhang, and W. H. Chan,
// Two Efficient Algorithms for Linear Time Suffix Array Construction
template <int THRESHOLD_NAIVE = 10, int THRESHOLD_DOUBLING = 40>
vll sa_is(const vll& s, ll upper) {
ll n = s.size();
if (n == 0) return {};
if (n == 1) return {0};
if (n == 2) {
if (s[0] < s[1]) {
return {0, 1};
} else {
return {1, 0};
}
}
if (n < THRESHOLD_NAIVE) {
return sa_naive(s);
}
if (n < THRESHOLD_DOUBLING) {
return sa_doubling(s);
}
vll sa(n);
std::vector<bool> ls(n);
for (int i = n - 2; i >= 0; i--) {
ls[i] = (s[i] == s[i + 1]) ? ls[i + 1] : (s[i] < s[i + 1]);
}
vll sum_l(upper + 1), sum_s(upper + 1);
rep (i, n) {
if (!ls[i]) {
sum_s[s[i]]++;
} else {
sum_l[s[i] + 1]++;
}
}
rep (i, upper + 1) {
sum_s[i] += sum_l[i];
if (i < upper) sum_l[i + 1] += sum_s[i];
}
auto induce = [&](const vll& lms) {
std::fill(sa.begin(), sa.end(), -1);
vll buf(upper + 1);
std::copy(sum_s.begin(), sum_s.end(), buf.begin());
for (auto d : lms) {
if (d == n) continue;
sa[buf[s[d]]++] = d;
}
std::copy(sum_l.begin(), sum_l.end(), buf.begin());
sa[buf[s[n - 1]]++] = n - 1;
for (int i = 0; i < n; i++) {
int v = sa[i];
if (v >= 1 && !ls[v - 1]) {
sa[buf[s[v - 1]]++] = v - 1;
}
}
std::copy(sum_l.begin(), sum_l.end(), buf.begin());
for (int i = n - 1; i >= 0; i--) {
int v = sa[i];
if (v >= 1 && ls[v - 1]) {
sa[--buf[s[v - 1] + 1]] = v - 1;
}
}
};
vll lms_map(n + 1, -1);
ll m = 0;
for (ll i = 1; i < n; i++) {
if (!ls[i - 1] && ls[i]) {
lms_map[i] = m++;
}
}
vll lms;
lms.reserve(m);
for (ll i = 1; i < n; i++) {
if (!ls[i - 1] && ls[i]) {
lms.push_back(i);
}
}
induce(lms);
if (m) {
vll sorted_lms;
sorted_lms.reserve(m);
for (ll v : sa) {
if (lms_map[v] != -1) sorted_lms.push_back(v);
}
vll rec_s(m);
ll rec_upper = 0;
rec_s[lms_map[sorted_lms[0]]] = 0;
for (ll i = 1; i < m; i++) {
ll l = sorted_lms[i - 1], r = sorted_lms[i];
ll end_l = (lms_map[l] + 1 < m) ? lms[lms_map[l] + 1] : n;
ll end_r = (lms_map[r] + 1 < m) ? lms[lms_map[r] + 1] : n;
bool same = true;
if (end_l - l != end_r - r) {
same = false;
} else {
while (l < end_l) {
if (s[l] != s[r]) {
break;
}
l++;
r++;
}
if (l == n || s[l] != s[r]) same = false;
}
if (!same) rec_upper++;
rec_s[lms_map[sorted_lms[i]]] = rec_upper;
}
auto rec_sa =
sa_is<THRESHOLD_NAIVE, THRESHOLD_DOUBLING>(rec_s, rec_upper);
rep(i, m) {
sorted_lms[i] = lms[rec_sa[i]];
}
induce(sorted_lms);
}
return sa;
}
// O(n + upper)
vll suffix_array(const vll& s, ll upper) {
assert(0 <= upper);
for (ll d : s) {
assert(0 <= d && d <= upper);
}
auto sa = sa_is(s, upper);
return sa;
}
// 時間O(nlogn), 空間O(n)
template <class T>
vll suffix_array(const std::vector<T> &s) {
ll n = s.size();
vll idx(n);
iota(idx.begin(), idx.end(), 0);
sort(idx.begin(), idx.end(), [&](int l, int r) { return s[l] < s[r]; });
vll s2(n);
ll now = 0;
rep(i, n) {
if (i && s[idx[i - 1]] != s[idx[i]]) now++;
s2[idx[i]] = now;
}
return sa_is(s2, now);
}
// O(n)
// "abracadabra"のときsuffix_array("abracadabra");は
// [10, 7, 0, 3, 5, 8, 1, 4, 6, 9, 2]
// これは接尾辞をソートして、その接尾辞の開始位置はどこかを表す。
// abracadabra, bracadabra, racadabra, acadabra, cadabra,
// adabra, dabra, abra, bra, ra, a
// 11個の接尾辞を持ち、それをソートすると
// a, abra, abracadabra, acadabra, adabra
// bra, bracadabra, cadabra, dabra, ra, racadabra
// となる。この並びで接尾辞の開始位置を求めると上記の
// [10, 7, 0, 3, 5, 8, 1, 4, 6, 9, 2]となる。
vll suffix_array(const std::string& s) {
ll n = s.size();
vll s2(n);
rep(i, n) {
s2[i] = s[i];
}
return sa_is(s2, 255);
}
// Reference:
// T. Kasai, G. Lee, H. Arimura, S. Arikawa, and K. Park,
// Linear-Time Longest-Common-Prefix Computation in Suffix Arrays and Its
// Applications
template <class T>
vll lcp_array(const std::vector<T>& s, const vll& sa) {
ll n = s.size();
assert(n >= 1);
vll rnk(n);
rep(i, n) {
rnk[sa[i]] = i;
}
vll lcp(n - 1);
ll h = 0;
rep(i, n) {
if (h > 0) h--;
if (rnk[i] == 0) continue;
ll j = sa[rnk[i] - 1];
for (; j + h < n && i + h < n; h++) {
if (s[j + h] != s[i + h]) break;
}
lcp[rnk[i] - 1] = h;
}
return lcp;
}
// string s = "abracadabra";
// vll sa = suffix_array(s);
// vll la = lcp_array(s, sa);
// print(sa);
// print(la);
// とすると
// saは[10, 7, 0, 3, 5, 8, 1, 4, 6, 9, 2]
// laは[1, 4, 1, 1, 0, 3, 0, 0, 0, 2]
// となる。
// laは、接尾辞をソートしたものに対して、
// 1個前の接尾辞とどれだけ共通部分があるのかを計算した結果となる。
// a, abra, abracadabra, acadabra, adabra, bra, bracadabra, cadabra, dabra, ra, racadabra
// aとabraは1, abraとabracadabraは4, abracadabraとacadabraは1というような感じで。
// O(n)
vll lcp_array(const string& s, const vll& sa) {
ll n = s.size();
vll s2(n);
rep(i, n) {
s2[i] = s[i];
}
return lcp_array(s2, sa);
}
template <class T> vll z_algorithm(const std::vector<T>& s) {
ll n = (ll)(s.size());
if (n == 0) return {};
vll z(n);
z[0] = 0;
for (ll i = 1, j = 0; i < n; i++) {
ll& k = z[i];
k = (j + z[j] <= i) ? 0 : std::min(j + z[j] - i, z[i - j]);
while (i + k < n && s[k] == s[i + k]) k++;
if (j + z[j] < i + z[i]) j = i;
}
z[0] = n;
return z;
}
// 入力の長さを n として、長さ n の配列を返す。
// i 番目の要素は s[0..n)とs[i..n)のLCP(Longest Common Prefix)の長さ。
// 文字列全体と、i番目から開始した文字列との共通文字列の長さを返す。
// z_algorithm("abcdeabcde");とすると
// [10, 0, 0, 0, 0, 5, 0, 0, 0, 0]
// O(n)
vll z_algorithm(const std::string& s) {
ll n = (ll)(s.size());
vll s2(n);
for (ll i = 0; i < n; i++) {
s2[i] = s[i];
}
return z_algorithm(s2);
}
// valの値をbase_stringの文字列にN進数のように変換する。(ABC171-C)
// 例えば、base_string=alpha();のとき
// valが1のときはa, 26のときはz, 27のときはaa, 702のときはzz、703のときはaaaなど。
string convValueToString(ll val, const string &base_string) {
ll length = len(base_string);
string s;
while (val) {
val--;
s += STR(base_string[val % length]);
val /= length;
}
REV(s);
return s;
}
// sの値をbase_stringの文字列にN進数のように変換する。(ABC285-C)
// 例えば、base_string={{'A', 1}, {'B', 2}, {'C', 3}};のとき
// valが1のときはA, 2のときはB, 4のときはAA, など。
// 引数にconstついてないけど、両方入力。
ll convStringToValue(string &s, ump<char, ll> &base_string) {
ll length = len(base_string);
ll val = 0;
for (auto c : s) {
val = val * length + base_string[c];
}
return val;
}
// O(NM)
// 2つの文字列S, Tの編集距離を求める。
// 操作1：S 中の文字を 1 つ選び、削除する。
// 操作2：S 中の文字を 1 つ選び、別の文字に変更する。
// 操作3：S 中の適当な位置に、文字を 1 つ挿入する。
// tokyoとkyotoなら距離は4。
class EditDistance {
string s_;
string t_;
public:
EditDistance(const string &s, const string &t)
: s_(s)
, t_(t) {}
// O(NM)
ll get_distance() {
ll slen = len(s_), tlen = len(t_);
vvll dp(slen + 1, vll(tlen + 1, LLONG_MAX));
dp[0][0] = 0;
rep (i, slen) {
rep (j, tlen) {
if (s_[i] == t_[j]) {
chmin(dp[i + 1][j], dp[i][j] + 1);
chmin(dp[i][j + 1], dp[i][j] + 1);
chmin(dp[i + 1][j + 1], dp[i][j]);
} else {
chmin(dp[i + 1][j], dp[i][j] + 1);
chmin(dp[i][j + 1], dp[i][j] + 1);
chmin(dp[i + 1][j + 1], dp[i][j] + 1);
}
}
}
rep (i, slen) {
chmin(dp[i + 1][tlen], dp[i][tlen] + 1);
}
rep (i, tlen) {
chmin(dp[slen][i + 1], dp[slen][i] + 1);
}
return dp[len(s_)][len(t_)];
}
};
class EditDistanceVector {
vll s_;
vll t_;
public:
EditDistanceVector(const vll &s, const vll &t)
: s_(s)
, t_(t) {}
// O(NM)
ll get_distance() {
ll slen = len(s_), tlen = len(t_);
vvll dp(slen + 1, vll(tlen + 1, LLONG_MAX));
dp[0][0] = 0;
rep (i, slen) {
rep (j, tlen) {
if (s_[i] == t_[j]) {
chmin(dp[i + 1][j], dp[i][j] + 1);
chmin(dp[i][j + 1], dp[i][j] + 1);
chmin(dp[i + 1][j + 1], dp[i][j]);
} else {
chmin(dp[i + 1][j], dp[i][j] + 1);
chmin(dp[i][j + 1], dp[i][j] + 1);
chmin(dp[i + 1][j + 1], dp[i][j] + 1);
}
}
}
rep (i, slen) {
chmin(dp[i + 1][tlen], dp[i][tlen] + 1);
}
rep (i, tlen) {
chmin(dp[slen][i + 1], dp[slen][i] + 1);
}
return dp[len(s_)][len(t_)];
}
};
// O(NM)
// 最長共通部分列の大きさを求める。
// abbcedとbcdeならbcdの長さ3。
class LCS {
string s_;
string t_;
vvll dp_;
void build() {
ll slen = len(s_), tlen = len(t_);
rep (i, slen) {
rep (j, tlen) {
if (s_[i] == t_[j]) {
chmax(dp_[i + 1][j + 1], dp_[i + 1][j]);
chmax(dp_[i + 1][j + 1], dp_[i][j + 1]);
chmax(dp_[i + 1][j + 1], dp_[i][j] + 1);
} else {
chmax(dp_[i + 1][j + 1], dp_[i + 1][j]);
chmax(dp_[i + 1][j + 1], dp_[i][j + 1]);
}
}
}
}
public:
LCS(const string &s, const string &t)
: s_(s)
, t_(t)
, dp_(len(s) + 1, vll(len(t) + 1)) {
build();
}
ll get_length() {
return dp_[len(s_)][len(t_)];
}
// O(max(N, M)) 最長共通部分文字列の取得
string get_str() {
vector<char> revans;
ll h = len(s_);
ll w = len(t_);
while (dp_[h][w] != 0)
{
if (h == 0) break;
if (w == 0) break;
if (dp_[h][w - 1] == dp_[h][w]) {
w--;
} else if (dp_[h - 1][w] == dp_[h][w]) {
h--;
} else if (dp_[h - 1][w - 1] + 1 == dp_[h][w])
{
revans.pb(t_[w - 1]);
h--, w--;
}
}
return std::string(revans.rbegin(), revans.rend());
}
};
string ll2str(ll x, ll base) {
if(x == 0) return "0";
stringstream ss;
string ret;
auto ll2base = [&]() {
stringstream tmp;
string cs = "0123456789" + alpha() + alpha_big();
while (x > 0) {
tmp << cs[(x % base)];
x /= base;
}
ret = tmp.str();
REV(ret);
};
switch(base) {
case 8:
ss << std::oct << x; ret = ss.str(); break;
case 10:
ss << std::dec << x; ret = ss.str(); break;
case 16:
ss << std::hex << x; ret = ss.str(); break;
default:
ll2base();
break;
}
return ret;
}
inline vvll init_tbl(int h, int w, ll val = 0)
{
vvll tbl(h, vll(w, val));
return tbl;
}
constexpr ll safe_mod(ll x, ll m) {
x %= m;
if (x < 0) x += m;
return x;
}
constexpr std::pair<ll, ll> inv_gcd(ll a, ll b) {
a = safe_mod(a, b);
if (a == 0) return {b, 0};
ll s = b, t = a;
ll m0 = 0, m1 = 1;
while (t) {
ll u = s / t;
s -= t * u;
m0 -= m1 * u;
auto tmp = s;
s = t;
t = tmp;
tmp = m0;
m0 = m1;
m1 = tmp;
}
if (m0 < 0) m0 += b / s;
return {s, m0};
}
ll inv_mod(ll x, ll m) {
assert(1 <= m);
auto z = inv_gcd(x, m);
assert(z.first == 1);
return z.second;
}
/*
vll r = {2, 3, 2};
vll m = {3, 5, 7};
INI2(a, b, crt(r, m));
debug(a, b);
としたとき
a: 23 | b: 105
これは、x≡2 (mod 3), x≡3 (mod 5), x≡2 (mod 7)の解が
x = 105k + 23ということを示す
*/
vll crt(const std::vector<ll>& r, const std::vector<ll>& m) {
assert(r.size() == m.size());
ll n = ll(r.size());
// Contracts: 0 <= r0 < m0
ll r0 = 0, m0 = 1;
for (ll i = 0; i < n; i++) {
assert(1 <= m[i]);
ll r1 = safe_mod(r[i], m[i]), m1 = m[i];
if (m0 < m1) {
std::swap(r0, r1);
std::swap(m0, m1);
}
if (m0 % m1 == 0) {
if (r0 % m1 != r1) return {0, 0};
continue;
}
ll g, im;
std::tie(g, im) = inv_gcd(m0, m1);
ll u1 = (m1 / g);
if ((r1 - r0) % g) return {0, 0};
ll x = (r1 - r0) / g % u1 * im % u1;
r0 += x * m0;
m0 *= u1;
if (r0 < 0) r0 += m0;
}
return vll{r0, m0};
}
// rep(i, n) floor((a * i + B), M)を高速に求める。
ll floor_sum(ll n, ll m, ll a, ll b) {
ll ans = 0;
if (a >= m) {
ans += (n - 1) * n * (a / m) / 2;
a %= m;
}
if (b >= m) {
ans += n * (b / m);
b %= m;
}
ll y_max = (a * n + b) / m, x_max = (y_max * m - b);
if (y_max == 0) return ans;
ans += (n - (x_max + a - 1) / a) * y_max;
ans += floor_sum(y_max, a, m, (a - x_max % a) % a);
return ans;
}
// dはaとbの最小公倍数、xとyはax + by = gcd(a, b)を満たすx, y
ll extgcd(ll a, ll b, ll &x, ll &y)
{
ll d = a;
if(b != 0) {
d = extgcd(b, a % b, y, x);
y -= (a / b) * x;
} else {
x = 1;
y = 0;
}
return d;
}
ll POW(ll n, ll r)
{
if (r == 0) return 1;
else if (r % 2 == 0) return POW(n * n, (ll)(r / 2));
else return n * POW(n, r - 1);
}
#define mod_m2p(a, m) (((m) + (a)) % (m))
#define mod_add(a, b, m) (((a) + (b)) % (m))
#define mod_sub(a, b, m) (((m) + (a) - (b)) % (m))
#define mod_mul(a, b, m) (mod_m2p(((a) % (m)) * ((b) % (m)), (m)))
ll mod_bipow_(ll x, ll y, ll m) { // x^y by bisection method
if (y == 0) return 1 % m;
else if (y == 1) return x % m;
else if (y % 2 == 0) {
ll val = mod_bipow_(x, (ll)(y / 2), m);
return mod_mul(val, val, m);
} else {
ll val = mod_bipow_(x, (ll)(y / 2), m);
return mod_mul(mod_mul(val, val, m), x, m);
}
}
ll mod_inv(ll x, ll pm) { return mod_bipow_(mod_m2p(x, pm), pm - 2, pm); } // x^{-1} = x^{MOD-2} (MOD: prime number)
ll mod_div(ll a, ll b, ll m) { return mod_mul(mod_m2p(a, m), mod_inv(mod_m2p(b, m), m), m); } // a/b = a*b^{-1}
ll mod_bipow(ll x, ll y, ll m) {
if (y < 0) {
ll xx = mod_div((ll)1, x, m);
return mod_bipow_(xx, -y, m);
}
return mod_bipow_(x, y, m);
}
// a x ≡ b^k + c (mod m)を求める。(aとmが合成数のとき)
// このときは、b^kをm*aの法として求めてaで割ればいい。
// bがaで割り切れる前提。
// p / q (mod N) = p (mod N * q) / q
ll mod_div_composite_number(ll a, ll b, ll k, ll c, ll m) {
return safe_mod((mod_bipow(b, k, m * a) + c) / a, m);
}
// ax ≡ b (mod m)の合同式を解く。解は正の値
// 解が無ければ-1,
ll congruent_expression(ll a, ll b, ll m) {
ll g = gcd(gcd(m, a), b);
a /= g, b /= g, m /= g;
if (gcd(a, m) != 1) return -_1;
ll x = 0, y = 0;
extgcd(a, m, x, y);
return safe_mod(mod_mul(x, b, m), m);
}
// 三分探索(下に凸関数を定義する。)
// 関数f(x)が最小値となるxを返す。
// // 目的関数(最小化したい)
// auto f = [a, b](ld x) {
// return a / sqrt(x + 1.0) + x * b;
// };
// ld minimalx = TernarySearch(f);
ld TernarySearch(const function<ld(ld)> &func, ld low = 0.0, ld high = 1e18) {
ld l = low, h = high;
int cnt = 1000;
while (cnt--) {
ld c1 = (l * 2 + h) / 3;
ld c2 = (l + h * 2) / 3;
if (func(c1) > func(c2)) l = c1;
else h = c2;
}
return l;
}
// 小数を表す文字列を小数部分が整数で表せるように数値をオフセットして
// 整数値にして返す。
// 例えば、dblstr2ll("123.456", 3)は123456
// 例えば、dblstr2ll("123.0456", 5)は12304560
// LLONG_MAXを超えないように注意
ll dblstr2ll(const string &dblstr, ll fractional_part_cnt) {
ll idx = 0;
string X = "", Y = "";
while(idx != len(dblstr) && dblstr[idx] != '.') {
X += dblstr[idx];
idx++;
}
idx++;
while(idx < len(dblstr)) {
Y += dblstr[idx];
idx++;
}
return STRLL(X) * POW(10, fractional_part_cnt) + STRLL(Y) * POW(10, fractional_part_cnt - len(Y));
}
// 尺取りのベース関数
auto syakutori(const vll &vec) {
ll N = len(vec);
ll val = 0, cnt = 0;
ll l = 0, r = 0;
while (l < N)
{
// if (val == N) ++cnt; カウントする条件を記載
bool cond = N < val; // 判定範囲を超えたかどうか
if ((r == N) || cond)
{
// val -= vec[l]; // 範囲を縮小するときに演算する処理
++l;
}
else // 範囲を拡大するときに演算する処理
{
// val += vec[r];
++r;
}
}
return cnt;
}
// 昇順ソートされたvecの[l, r]の区間の差がV以下の組み合わせ数を求める。
// 例えば[-1, 0, 2]でV=1のときは1, V=2のときは2, V=3のときは3
auto RangeCount(const vll &vec, ll V) {
ll N = len(vec);
ll cnt = 0;
ll l = 0, r = 1;
while (l < N)
{
bool cond = r == N || V < vec[r] - vec[l]; // 判定範囲を超えたかどうか
if (cond)
{
if (vec[r - 1] - vec[l] <= V) cnt += r - 1 - l;
++l;
}
else ++r;
}
return cnt;
}
// 連続した要素の部分和がsvalと等しくなる個数をカウント。O(N)
auto subsum(const vll &vec, const ll sval) {
ll N = len(vec);
ll val = 0, cnt = 0;
ll l = 0, r = 0;
while (l < N)
{
if (val == sval)
++cnt;
if ((r == N) || sval < val)
{
val -= vec[l];
++l;
}
else
{
val += vec[r];
++r;
}
}
return cnt;
}
// 全区間[L, R]の区間和の総和
// [2, 3, 5]のとき、ret = 2 + (2 + 3) + (2 + 3 + 5) + 3 + (3 + 5) + 5が戻り値
ll rangesum(const vll &v, ll ansmod) {
ll N = len(v);
ll ret = 0;
rep (i, N) {
ret += v[i] * (i + 1) * (N - i);
ret %= ansmod;
}
return ret % ansmod;
}
// 全区間[L, R]の区間和の総和
// [2, 3, 5]のとき、ret = 2 + (2 + 3) + (2 + 3 + 5) + 3 + (3 + 5) + 5が戻り値
ll rangesum(const vll &v) {
ll N = len(v);
ll ret = 0;
rep (i, N) {
ret += v[i] * (i + 1) * (N - i);
}
return ret;
}
template <typename T>
class Counter
{
public:
unordered_map<T, ll> tbl_;
ll max_cnt_ = 0;
T max_key_;
ll min_cnt_ = -1;
T min_key_;
Counter(const vector<T> &vec) {
rep(i, len(vec)) {
ll v = ++tbl_[vec[i]];
if (max_cnt_ < v) {
max_cnt_ = v;
max_key_ = vec[i];
}
}
}
ll count(T key) {
return tbl_[key];
}
T maxkey() {
return max_key_;
}
ll maxcnt() {
return max_cnt_;
}
T minkey() {
if (min_cnt_ == -1) {
mincnt();
}
return min_key_;
}
ll mincnt() {
if (min_cnt_ == -1) {
min_key_ = tbl_.begin()->first;
min_cnt_ = tbl_.begin()->second;
for(auto itr = tbl_.begin(); itr != tbl_.end(); itr++){
if(min_cnt_ > itr->second) {
min_key_ = itr->first;
min_cnt_ = itr->second;
}
}
}
return min_cnt_;
}
};
// 簡易二項係数 // 100C50=10 ** 29くらいまでなら計算できる。
ll comb(ll n, ll k)
{
static bool iscalc = false;
static vvll C(101, vll(101));
if (!iscalc) {
C[0][0] = 1;
reps (i, 1, 101) {
rep (j, i + 1) {
if (j == 0) C[i][0] = 1;
else if (j == i) C[i][i] = 1;
else C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
}
}
}
return C[n][k];
}
// nが大きくてkが小さいときO(k) nが大きい時はllを使う。
ll combNbig(ll n, ll k) {
if(k < 0 || n < k) return 0;
ll ret = 1;
for(ll i = 1; i <= k; ++i) {
ret *= n--;
ret /= i;
}
return ret;
}
ll combNbigMod(ll n, ll k) {
if(k < 0 || n < k) return 0;
mint numerator = 1;
mint denominator = 1;
for(ll i = 1; i <= k; ++i) {
numerator *= n--;
denominator *= i;
}
return (numerator / denominator).val();
}
ll mysqrt(ll n) {
ll ok = 0, ng = n + 1;
while (ng - ok > 1) {
ll mid = (ng + ok) >> 1;
if (mid * mid <= n) {
ok = mid;
} else {
ng = mid;
}
}
return ok;
}
ll det2(const vvll &mat) {
return mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0];
}
// 等差数列の和。初項a0, 公差d, 0番目が初項のときのn番目までの総和。n番目を含む
ll sumArithmeticSequence(ll a0, ll d, ll n) {
if (n < 0) return _0;
return (n + 1) * (2 * a0 + n * d) / 2;
}
// 等差数列の和。初項a0, 公差d, 0番目が初項のときのn番目までの総和。n番目を含む
mint sumArithmeticSequence(mint a0, mint d, mint n) {
return (n + 1) * (2 * a0 + n * d) / 2;
}
ll gcds(const vll &v) {
if (len(v) == 0) return _1;
ll g = v[0];
rep(i, len(v)) g = gcd(g, v[i]);
return g;
}
// 約数個数列挙(MAXNまで)
vll divisors_count(ll MAXN = 1000000)
{
vll div = vll(MAXN + 1, 0);
reps(n, 1, MAXN + 1) repsp(m, n, MAXN + 1, n) {
div[m] += 1;
}
return div;
}
// Nの階乗がkで何回割れるかを返す。(ルジャンドルの公式)
ll LegendresFormula(ll N, ll k) {
ll v = k;
ll cnt = 0;
while(v <= N) {
cnt += floor(N, v);
v *= k;
}
return cnt;
}
// 約数列挙
// 約数の個数は大体n^(1/3)個
vll make_divisors(ll n) {
vll divisors;
for(ll i = 1; i * i <= n; ++i) {
if(n % i == 0) {
divisors.pb(i);
if(i != n / i) {
divisors.pb(n / i);
}
}
}
return divisors;
}
// N以下の素数列挙(Nはせいぜい10^7くらいまで)
inline vll eratosthenes(ll N) {
vll ps;
vector<bool> primes(N + 1);
rep(i, len(primes)) primes[i] = true;
primes[0] = primes[1] = false;
rep(i, len(primes)) {
if (primes[i]) {
ps.pb(i);
for (ull j = i + i; j < len(primes); j += i) {
primes[j] = false;
}
}
}
return ps;
}
// 高速素数判定
bool suspect(ll a, ll s, ll d, ll md){
ll x = 1, xx = a, one = x, minusone = md - 1;
while(d > 0){
if(d & 1) x = mod_mul(x, xx, md);
xx = mod_mul(xx, xx, md);
d >>= 1;
}
if ((x % md) == (one % md)) return true;
for (ll r = 0; r < s; ++r) {
if((x % md) == (minusone % md)) return true;
x = mod_mul(x, x, md);
}
return false;
}
// 64ビット整数までの高速素数判定
bool is_prime_miller_rabin(ll n){
if (n <= 1 || (n > 2 && n % 2 == 0)) return false;
ll d = n - 1, s = 0;
while (!(d&1)) {++s; d >>= 1;}
vll v = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
if(n < 4759123141LL) v = {2, 7, 61};
for (auto &&p : v) {
if(p >= n) break;
if(!suspect(p, s, d, n)) return false;
}
return true;
}
// 高速素数判定
bool is_prime(ll p, int k = 50)
{
ll q = abs(p);
if (q == 2) return true;
if ((q < 2) || ((q & 1) == 0)) return false;
ll d = (q - 1) >> 1;
while ((d & 1) == 0) d >>= 1;
static std::random_device rd;
static std::mt19937_64 engine(rd());
std::uniform_int_distribution<ll> distr(1, q - 1);
rep(i, k){
ll a = distr(engine);
ll t = d;
ll y = mod_bipow(a, t, q);
while ((t != q - 1) && (y != 1) && (y != q - 1))
{
y = mod_bipow(y, 2, q);
t <<= 1;
}
if ((y != q - 1) && ((t & 1) == 0)) return false;
}
return true;
}
// Pollard のロー素因数分解法で使う関数。
ll pollard(ll N) {
if (N % 2 == 0) return _2;
if (is_prime_miller_rabin(N)) return N;
auto f = [&](ll x) -> ll {
return (__int128_t(x) * x + 1) % N;
};
ll step = 0;
while (true) {
++step;
ll x = step, y = f(x);
while (true) {
ll p = gcd(abs(y - x + N), N);
if (p == 0 || p == N) break;
if (p != 1) return p;
x = f(x);
y = f(f(y));
}
}
}
// Pollard のロー素因数分解法
// 1秒以内に10^18以下の整数の素因数分解を100回くらいできる
umpll prime(ll N) {
umpll primes;
if (N == 1) return primes;
ll p = pollard(N);
if (p == N) {
primes[p]++;
return primes;
}
umpll left = prime(p);
umpll right = prime(N / p);
repdict (k, v, left) primes[k] += v;
repdict (k, v, right) primes[k] += v;
return primes;
}
// 高速素因数分解(MAXNの範囲まで)
class Prime
{
public:
static vll sieve; // 何回もPrimeのインスタンスを生成するときはstaticをはずして、下の実体もコメントアウトする。
Prime(ll MAXN = 10000000) {
rep(i, MAXN + 1) sieve.pb(i);
ll p = 2;
while (p * p <= MAXN) {
if (sieve[p] == p) {
repsp(q, 2 * p, MAXN + 1, p) {
if (sieve[q] == q) sieve[q] = p;
}
}
p += 1;
}
}
Counter<ll> prime(ll n) {
if (n == 1) return Counter<ll>(vll());
vll tmp;
while (n > 1) {
tmp.pb(sieve[n]);
n = (ll)(n / sieve[n]);
}
return Counter<ll>(tmp);
}
};
vll Prime::sieve = vll();
// 素数と使用出来る個数が与えられたものから、N以下を満たす値を求める。
// 例えば2が4個, 3が2個から35以下の値を列挙する。
// そのときpsにはps = vvll{{2, 4}, {3, 2}};とする。
vll makeDivisorsFromPrimeNOrLess(vvll ps, ll N = LLONG_MAX) {
vll divs;
auto f = [&](auto &&f_, ll idx, ll L, ll V) {
if (idx == L) {
divs.pb(V);
return;
}
ll tmp = 1;
rep (i, ps[idx][1] + 1) {
if (N < V * tmp) break;
f_(f_, idx + 1, L, V * tmp);
tmp *= ps[idx][0];
}
};
f(f, 0, len(ps), 1);
return divs;
}
// Nの2乗の約数を求める。
// Nを2乗してから約数求めるより速い。
vll makePrime2Divisors(ll N) {
static Prime P = Prime();
auto p = P.prime(N);
vvll ps;
repdict (k, v, p.tbl_) {
ps.pb(vll{k, v * 2});
}
return makeDivisorsFromPrimeNOrLess(ps);
}
ll LIS(const vector<ll> &a, bool strict = true) // trueのとき厳密に増加する列
{
vll lis;
for(auto &p : a) {
vll::iterator it;
if(strict) it = lower_bound(all(lis), p);
else it = upper_bound(all(lis), p);
if(end(lis) == it) lis.emplace_back(p);
else *it = p;
}
return len(lis);
}
// 各長さでの最長部分増加列が求まる
vll LISList(const vll &vec) {
ll N = len(vec);
vll ans(N);
vll lis = {vec[0]};
rep (i, N) {
if (vec[i] > lis[len(lis) - 1]) {
lis.pb(vec[i]);
} else {
ll idx = LB_IDX_VEC(lis, vec[i]);
lis[idx] = vec[i];
}
ans[i] = len(lis);
}
return ans;
}
// [[xi, yi], [xj, yj]...]のLISを求める
// xi < xj, yi < yjを満たす最大長
ll LIS2D(const vvll &vec, bool strict = true) {
vvll tmp = vec;
std::sort(std::begin(tmp), std::end(tmp),
[](const vll& a, const vll& b) {
if (a[0] == b[0]) { // 第一要素同じなら降順ソート
return a[1] > b[1];
}
return a[0] < b[0]; // 第一要素は昇順
}
);
vll v;
rep(i, len(tmp)) v.pb(tmp[i][1]);
return LIS(v, strict);
}
ll LDS(const vector<ll> &a, bool strict = true) // trueのとき厳密に減少する列
{
vll b;
rep(i, len(a)) {
b.pb(-a[i]);
}
return LIS(b, strict);
}
ll randint(ll l, ll r)
{
static random_device rnd;
static mt19937_64 mt(rnd());
uniform_int_distribution<> rand(l, r - 1);
return rand(mt);
}
// 「#」がもの、「.」が何もないとすると、一回り周辺を「outside」で埋める
// マップを作成する。「.」が0で「#」が1を表す
inline vvll create_map(int h, int w, ll outside = 0)
{
vs mp = in_vs(h);
vvll tbl = init_tbl(h + 2, w + 2, outside);
reps(i, 1, h + 1) reps(j, 1, w + 1) {
tbl[i][j] = mp[i - 1][j - 1] == '.' ? 0 : 1;
}
return tbl;
}
inline bool substrcheck(const string &in, const string &ptn)
{
boyer_moore_searcher searcher(all(ptn));
auto it = search(all(ptn), searcher);
return it != in.end();
}
// nCkをmで割った余りを求める
class Combination {
const ll n_;
vector<mint> facts_;
vector<mint> inv_facts_;
public:
Combination(ll N) : n_(2 * N), facts_(n_ + 1), inv_facts_(n_ + 1) {
rep(i, n_ + 1) facts_[i] = i == 0 ? 1 : facts_[i - 1] * i;
for (ll i = n_; i >= 0; i--) inv_facts_[i] = i == n_ ? facts_[n_].inv() : inv_facts_[i + 1] * (i + 1); // (i!)^{-1}=((i+1)!)^{-1}*(i+1)
}
mint nPr(ll n, ll r) {
if (n < r) return 0;
return facts_[n] * inv_facts_[n - r];
}
mint nCr(ll n, ll r) {
if (n < r) return 0;
return facts_[n] * inv_facts_[r] * inv_facts_[n - r];
}
mint nHr(ll n, ll r) {
return nCr(n + r - 1, r);
}
mint catalan(ll n) { // https://ja.wikipedia.org/wiki/%E3%82%AB%E3%82%BF%E3%83%A9%E3%83%B3%E6%95%B0
return nCr(2 * n, n) / (n + 1);
}
// カタラン数・・・(2n C n)/(n + 1) = (2n C n) - (2n C n-1)
// c0 = 1, c_n = rep(i, n) c[i] * c[n - i - 1]
// c0から順に1,1,2,5,14,42,132,429,1430,...
// 1 <= a1 <= a2 <= a3 <= a4 <= ... <= ak <= Nの組み合わせの数。
// CombinationのコンストラクタにN + Kを入れておくこと。
mint loopcnt(ll n, ll k) {
assert(n + k <= n_);
return nCr(n - 1 + k, n - 1);
}
};
// Lに開始し、Rに終了し、終了と同時刻に区間[L, R]を割り当てれる。
// 重ねることはできない。
ll max_cnt_interval_scheduling(vvll l) // コピーが実行される
{
SORT_IDX(l, 1);
ll r1 = l[0][1];
ll cnt = 1;
reps(i, 1, len(l)) {
ll l1 = l[i][0], r2 = l[i][1];
if (l1 >= r1) {
cnt += 1;
r1 = r2;
}
}
return cnt;
}
// 自作 imos 2D
template <class T>
class Imos2D
{
private:
vector<vector<T>> m_Data;
int m_W, m_H;
bool m_Built;
public:
Imos2D(int W, int H)
: m_Data(H + 1, vector<T>(W + 1, 0)), m_W(W), m_H(H), m_Built(false)
{
}
// [(sx, sy), (tx, ty)) に add を加算する(半開区間)
void RangeAdd(int sx, int sy, int tx, int ty, T add)
{
assert(!m_Built);
m_Data[sy][sx] += add; // 左上
m_Data[sy][tx] -= add; // 右上
m_Data[ty][sx] -= add; // 左下
m_Data[ty][tx] += add; // 右下
}
void Build()
{
assert(!m_Built);
// 横方向の累積和
for (int y = 0; y < m_H; y++)
{
for (int x = 1; x < m_W; x++)
{
m_Data[y][x] += m_Data[y][x - 1];
}
}
// 縦方向の累積和
for (int y = 1; y < m_H; y++)
{
for (int x = 0; x < m_W; x++)
{
m_Data[y][x] += m_Data[y - 1][x];
}
}
m_Built = true;
}
int Get(int x, int y)
{
assert(m_Built);
return m_Data[y][x];
}
};
class IntegralImage {
public:
vvll data_;
IntegralImage(ull H, ull W) : data_(H + 1, vll(W + 1, 0)) {}
void add(ull h, ull w, ll z) {
++h, ++w;
if(h >= len(data_) || w >= len(data_[0])) return;
data_[h][w] += z;
}
void build() {
for(ull i = 1; i < len(data_); i++) {
for(ull j = 1; j < len(data_[i]); j++) {
data_[i][j] += data_[i][j - 1] + data_[i - 1][j] - data_[i - 1][j - 1];
}
}
}
// matrixの升目の添え字で考えるのではなく
// matrixの格子点で左上、右下で囲われた領域の総和を求める
ll get(int sh, int sw, int gh, int gw) {
return (data_[gh][gw] - data_[sh][gw] - data_[gh][sw] + data_[sh][sw]);
}
};
// class UnionFind {
// ll n_;
// vll size_;
// vll par_;
// vll link_;
// public:
// UnionFind(ll n) : n_(n), size_(n, 1), par_(n), link_(n) {
// iota(all(par_), 0);
// iota(all(link_), 0);
// }
// ll find(ll x) {
// while (par_[x] != x) {
// par_[x] = par_[par_[x]];
// x = par_[x];
// }
// return x;
// }
// ll operator[](ll x) { return find(x); }
// bool unite(ll x, ll y, bool fixed = false) {
// x = find(x);
// y = find(y);
// if (x == y) { return false; }
// if (!fixed && y < x) swap(x, y);
// size_[x] += size_[y];
// size_[y] = 0;
// par_[y] = x;
// swap(link_[x], link_[y]);
// return true;
// }
// vll find_all() {
// vll A(n_);
// rep(i, n_) A[i] = find(i);
// return A;
// }
// vll members(ll x) {
// vll mems = vll{x};
// for (ll y = link_[x]; y != x; y = link_[y]) mems.pb(y);
// return mems;
// }
// ll size(ll x) {
// return size_[find(x)];
// }
// bool same(ll x, ll y) {
// return find(x) == find(y);
// }
// vll roots() {
// vll rs;
// rep(i, n_) if (size_[i] > 0) rs.pb(i);
// return rs;
// }
// ll group_count() {
// return len(roots());
// }
// unordered_map<ll, vll> all_group_members() {
// unordered_map<ll, vll> group_members;
// rep(member, n_) group_members[find(member)].pb(member);
// return group_members;
// }
// };
// 重みつきUnionFind
// 1-originで扱う。
class UnionFind {
ll n_;
vll size_;
vll par_;
vll link_;
vll rank_;
vll par_diff_;
public:
// コストが∞となるサイクルがあった場合、超頂点0番と連結する。
UnionFind(ll n) : n_(n + 1), size_(n_, 1), par_(n_), link_(n_), rank_(n_), par_diff_(n_)
{
iota(all(par_), 0);
iota(all(link_), 0);
}
// 要素xが属する木の根を再帰的に見つける
ll find(ll x) {
if (par_[x] == x) return x;
else { // 経路圧縮 + 累積和
ll ret = find(par_[x]);
if (par_diff_[par_[x]] == LLONG_MAX) par_diff_[x] = LLONG_MAX;
else par_diff_[x] += par_diff_[par_[x]];
return par_[x] = ret;
}
}
ll operator[](ll x) { return find(x); }
bool unite(ll x, ll y, ll w = 0) {
if (x != 0 && same(x, y) && diff(x, y) != w) unite(0, y);
ll rx = find(x);
ll ry = find(y);
ll wt = w;
wt += weight(x);
wt -= weight(y);
if (rx == ry) {
return false;
}
if (ry < rx) {
swap(rx, ry);
wt = -wt;
}
size_[rx] += size_[ry];
if (rank_[rx] == rank_[ry]) rank_[rx]++;
size_[ry] = 0;
par_[ry] = rx;
par_diff_[ry] = wt;
swap(link_[rx], link_[ry]);
return true;
}
vll find_all() {
vll A(n_);
rep(i, n_) A[i] = find(i);
return A;
}
vll members(ll x) {
vll mems = vll{x};
for (ll y = link_[x]; y != x; y = link_[y]) mems.pb(y);
return mems;
}
ll size(ll x) {
return size_[find(x)];
}
bool same(ll x, ll y) {
return find(x) == find(y);
}
vll roots() {
vll rs;
reps(i, 1, n_) if (size_[i] > 0) rs.pb(i);
return rs;
}
ll group_count() {
return len(roots());
}
// 超頂点0番の情報は含めない
unordered_map<ll, vll> all_group_members() {
unordered_map<ll, vll> group_members;
reps(member, 1, n_) group_members[find(member)].pb(member);
return group_members;
}
// 経路圧縮 + costを返す
ll weight(ll x) {
find(x);
return par_diff_[x];
}
// yのcost - xのcost
ll diff(ll x, ll y) {
if (same(0, x) || same(0, y)) return LLONG_MAX;
return weight(y) - weight(x);
}
};
class Graph
{
private:
const ll V_;
const bool dir_;
const ll ansmod_;
// dist(-1で初期化), cntは呼び出し元でN個分の配列を与えること。connect_vtxsは空のvectorでよい。
void bfs_(ll sv, vll &dist, vll &connect_vtxs, vll &cnt, vll &root, ll search_depth = 1000000000)
{
if (dist[sv] != -1) return;
dll q = dll();
q.pb(sv);
dist[sv] = 0;
connect_vtxs.pb(sv);
cnt[sv] = 1;
if (search_depth == 0) return;
while (len(q) > 0) {
ll p = q[0];
q.popleft();
if (!EXIST(p, edges_)) continue;
vector<pll> &evw = edges_[p];
for (const auto& [ev, w] : evw) {
bool isignore = false;
rep(i, len(ignore_edges_[p])) {
const auto& [igev, igw] = ignore_edges_[p][i];
if (ev == igev && w == igw) {
isignore = true;
}
}
if (isignore) continue;
if (dist[ev] != -1) {
if (dist[ev] == dist[p] + 1) {
cnt[ev] += cnt[p];
cnt[ev] %= ansmod_;
}
continue;
}
dist[ev] = dist[p] + 1;
root[ev] = p;
cnt[ev] = cnt[p];
connect_vtxs.pb(ev);
if (dist[ev] < search_depth)
{
if (w == 0) q.pf(ev);
else q.pb(ev);
}
}
}
}
// 木で無向辺のみ対応
void dfs_(ll sv, ll pre = -1) {
connect_vtxs_.pb(sv);
root_[sv] = pre;
traverse_preorder_path_.pb(sv);
vector<pll> &evw = edges_[sv];
ll child_cnt = 0;
for (const auto& [ev, w] : evw) {
if (ev == pre) continue;
bool isignore = false; // 無効化されている辺は無視する。
rep(i, len(ignore_edges_[sv])) {
const auto& [igev, igw] = ignore_edges_[sv][i];
if (ev == igev && w == igw) {
isignore = true;
}
}
if (isignore) continue;
dfs_(ev, sv);
child_cnt++;
}
if (child_cnt == 0) { // 辺が1つのときで根でないとき
traverse_reached_[sv] = true;
traverse_inorder_path_.pb(sv);
}
if (pre != -1 && !traverse_reached_[pre]) {
traverse_reached_[pre] = true;
traverse_inorder_path_.pb(pre);
}
traverse_postorder_path_.pb(sv);
}
public:
unordered_map<ll, vector<pll>> edges_;
unordered_map<ll, vector<pll>> ignore_edges_; // 多重辺で同じ重みの辺があるとダメ
vll outdeg_;
vll indeg_;
vll path_; // dfsでたどった経路
vll traverse_path_; // dfsでたどり着けるノード順
vll traverse_preorder_path_; // dfsでたどり着けるノード順 (行きがけ順)
vll traverse_inorder_path_; // dfsでたどり着けるノード順 (通りがけ順)
vll traverse_postorder_path_; // dfsでたどり着けるノード順 (帰りがけ順) (木DPで使えそう)
vector<bool> traverse_reached_; // dfsでたどり着いたチェック
vll dist_;
vll connect_vtxs_;
vll cnt_;
vll root_;
Graph(ll V, bool dir, ll ansmod = 1000000007) : V_(V), dir_(dir), ansmod_(ansmod), traverse_reached_(V_, false), dist_(V_, -1), cnt_(V_), root_(V_, -1){
outdeg_ = vll(V);
indeg_ = vll(V);
}
void append_edge(ll sv, ll ev, ll weight = 1)
{
vector<pll> &u = edges_[sv];
pll v = make_pair(ev, weight);
u.pb(v);
outdeg_[sv]++;
indeg_[ev]++;
if (!dir_) {
swap(sv, ev);
outdeg_[sv]++;
indeg_[ev]++;
vector<pll> &ru = edges_[sv];
pll rv = make_pair(ev, weight);
ru.pb(rv);
}
}
void ignore_edge(ll sv, ll ev, ll weight = 1) {
vector<pll> &u = ignore_edges_[sv];
pll v = make_pair(ev, weight);
u.pb(v);
if (!dir_) {
swap(sv, ev);
vector<pll> &ru = ignore_edges_[sv];
pll rv = make_pair(ev, weight);
ru.pb(rv);
}
}
void ignore_edge_clear() {
ignore_edges_.clear();
}
void bfs(ll sv, bool reset = true, ll search_depth = 1000000000) {
if (reset) {
rep(i, len(connect_vtxs_)) {
dist_[connect_vtxs_[i]] = -1;
cnt_[connect_vtxs_[i]] = 0;
root_[connect_vtxs_[i]] = -1;
}
connect_vtxs_.clear();
}
bfs_(sv, dist_, connect_vtxs_, cnt_, root_, search_depth);
}
void dfs(ll sv, bool reset = false) {
if (traverse_reached_[sv]) return; // すでに到達済みの点はdfsしない
if (reset) {
traverse_preorder_path_.clear();
traverse_inorder_path_.clear();
traverse_postorder_path_.clear();
rep(i, len(connect_vtxs_)) {
traverse_reached_[connect_vtxs_[i]] = false; // 「他の始点から到達してたときは無視したい」場合はこれをコメントアウト
root_[connect_vtxs_[i]] = -1;
}
connect_vtxs_.clear();
}
dfs_(sv);
}
void dfs_path(ll sv, ll N)
{
path_.clear();
traverse_path_.clear();
vector<bool> reached = vector<bool>(N, false);
dll q = dll();
q.pb(sv);
while (len(q) != 0 ) {
ll p = q.back();
q.popright();
if (p >= 0) {
reached[p] = true;
path_.pb(p);
traverse_path_.pb(p);
if (EXIST(p, edges_)) {
vector<pll> &evw = edges_[p];
for (const auto& [ev, w] : evw) {
if (reached[ev]) continue;
bool isignore = false;
rep(i, len(ignore_edges_[p])) {
const auto& [igev, igw] = ignore_edges_[p][i];
if (ev == igev && w == igw) {
isignore = true;
}
}
if (isignore) continue;
q.pb(~p);
q.pb(ev);
}
}
}
else path_.pb(~p);
}
}
vll topo_sort()
{
heapqll q;
vll to = vll(V_);
vll topo_vertex_list;
repdict(u, vtxs, edges_) {
for (const auto& [ev, w] : vtxs) {
++to[ev];
}
}
rep(i, V_) {
if (to[i] != 0) continue;
HEAPPUSH(q, i);
}
while (len(q) != 0) {
ll v = HEAPPOP(q);
topo_vertex_list.pb(v);
for (const auto [ev, w] : edges_[v]) {
--to[ev];
if (to[ev]) continue;
HEAPPUSH(q, ev);
}
}
return topo_vertex_list;
}
ll longest_path() { // 有向グラフで非連結グラフの最大パスの長さを求める。
dll q;
vll dist(V_);
rep (i, V_) {
if (indeg_[i] == 0) q.pb(i);
}
while (len(q) != 0) {
ll u = POPLEFT(q);
vector<pll> &v = edges_[u];
rep (i, len(v)) {
chmax(dist[v[i].first], dist[u] + 1);
indeg_[v[i].first]--;
if (indeg_[v[i].first] == 0) q.pb(v[i].first);
}
}
return MAX(dist);
}
// 無向グラフのみ対応。
// start->endのパスが無ければ空を返す
vll get_path(ll start, ll end, ll vertexidx_offset = 0) {
bfs(start);
if (dist_[end] == -1) return vll{};
ll pos = end;
vll path = {end + vertexidx_offset};
while(pos != start) {
pos = root_[pos];
path.pb(pos + vertexidx_offset);
}
REV(path);
return path;
}
// 先にbfsを実行しておく。
vll parts_tree_size() {
vvll tmp;
rep (i, V_) {
tmp.pb(vll{dist_[i], i});
}
SORT(tmp);
REV(tmp);
vll ans(V_);
rep (i, V_) {
INI2(d, idx, tmp[i]);
UNUSED(d);
if (root_[idx] != -1) ans[root_[idx]] += ans[idx] + 1;
}
return ans;
}
// 二部グラフ判定 O(NlogN)
bool is_bipartite() {
auto U = UnionFind(2 * V_);
repdict(sv, evw, edges_) {
for (const auto& [ev, w] : evw) {
bool isignore = false; // 無効化されている辺は無視する。
rep(i, len(ignore_edges_[sv])) {
const auto& [igev, igw] = ignore_edges_[sv][i];
if (ev == igev && w == igw) {
isignore = true;
}
}
if (isignore) continue;
U.unite(sv, ev + V_);
U.unite(sv + V_, ev);
}
}
rep(i, V_) {
if (U.same(i, i + V_)) return false;
}
return true;
}
// 奇数閉路のグラフは二部グラフチェックで確認できる。
bool exist_odd_cycle() {
return !is_bipartite();
}
};
// ベルマンフォード法
// 負の辺が混ざっててもOK、負のループ検出ができる。
class GraphBellmanFord {
using Weight = ll;
using Flow = ll;
struct Edge
{
ll src, dst;
ll rev;
Weight weight;
Flow cap;
Edge() : src(0), dst(0), weight(0) {}
Edge(ll s, ll d, Weight w) : src(s), dst(d), weight(w) {}
};
using Edges = std::vector<Edge>;
const ll V_;
const bool dir_;
Edges edges_;
public:
GraphBellmanFord(ll V, bool dir) : V_(V), dir_(dir), edges_(Edges(V)) {}
void append_edge(ll sv, ll ev, ll weight = 1)
{
edges_.pb(Edge(sv, ev, weight));
if (!dir_) {
edges_.pb(Edge(ev, sv, weight));
}
}
std::pair<std::vector<Weight>, bool> bellmanFord(ll s)
{
const Weight inf = LLONG_MAX;
std::vector<Weight> dist(V_, inf);
dist[s] = 0;
bool negCycle = false;
for (ll i = 0;; i++)
{
bool update = false;
//すべての辺について、その辺をとおった場合に最短経路が更新できる場合は更新する
for (auto &e : edges_)
{
if (dist[e.src] != inf && dist[e.dst] > dist[e.src] + e.weight)
{
dist[e.dst] = dist[e.src] + e.weight;
update = true;
}
}
//更新がなくなったらおはり
if (!update)
break;
// n回以上更新されてたら負閉路がある
if (i > V_)
{
negCycle = true;
break;
}
}
return std::make_pair(dist, !negCycle);
}
//ゴールを指定して、それまでのパスに負閉路がなかったらOK
std::pair<std::vector<Weight>, bool> bellmanFord(ll s, ll d)
{
const Weight inf = std::numeric_limits<Weight>::max() / 8;
//初期化、スタート地点以外の距離は無限大
std::vector<Weight> dist(V_, inf);
dist[s] = 0;
bool negCycle = false;
for (int i = 0; i < V_ * 2; i++)
{
bool update = false;
//すべての辺について、その辺をとおった場合に最短経路が更新できる場合は更新する
for (auto &e : edges_)
{
if (dist[e.src] != inf && dist[e.dst] > dist[e.src] + e.weight)
{
// n回目の更新で d が更新されてたら問答無用で負閉路ありとしてNG
if (i >= V_ - 1 && e.dst == d)
{
negCycle = true;
}
// 終点以外に負閉路がある場合はそこの距離を十分小さい値に置き換える
else if (i >= V_ - 1)
{
dist[e.dst] = -inf;
update = true;
}
else
{
dist[e.dst] = dist[e.src] + e.weight;
update = true;
}
}
}
//更新がなくなったらおはり
if (!update)
break;
}
return std::make_pair(dist, !negCycle);
}
};
// GraphSCCで使用。
template<class E>
struct csr {
vll start;
std::vector<E> elist;
csr(ll n, const std::vector<std::pair<ll, E>> &edges)
: start(n + 1)
, elist(edges.size()) {
for (auto e : edges) {
start[e.first + 1]++;
}
for (ll i = 1; i <= n; i++) {
start[i] += start[i - 1];
}
auto counter = start;
for (auto e : edges) {
elist[counter[e.first]++] = e.second;
}
}
};
// 有向グラフのみ対応
struct GraphSCC {
public:
GraphSCC(ll n)
: _n(n) {}
ll num_vertices() { return _n; }
void append_edge(ll from, ll to) { edges.push_back({from, {to}}); }
// グループの数と各頂点がどの今日連結成分に属しているかのidが設定されたものが返る。
std::pair<ll, std::vector<ll>> scc_ids() {
auto g = csr<edge>(_n, edges);
ll now_ord = 0, group_num = 0;
std::vector<ll> visited, low(_n), ord(_n, -1), ids(_n);
visited.reserve(_n);
auto dfs = [&](auto self, ll v) -> void {
low[v] = ord[v] = now_ord++;
visited.push_back(v);
for (ll i = g.start[v]; i < g.start[v + 1]; i++) {
auto to = g.elist[i].to;
if (ord[to] == -1) {
self(self, to);
low[v] = std::min(low[v], low[to]);
} else {
low[v] = std::min(low[v], ord[to]);
}
}
if (low[v] == ord[v]) {
while (true) {
ll u = visited.back();
visited.pop_back();
ord[u] = _n;
ids[u] = group_num;
if (u == v) break;
}
group_num++;
}
};
for (ll i = 0; i < _n; i++) {
if (ord[i] == -1) dfs(dfs, i);
}
for (auto &x : ids) {
x = group_num - 1 - x;
}
return {group_num, ids};
}
// 強連結成分ごとにvectorが作られて、その中に頂点番号が含まれている。
vvll scc() {
auto ids = scc_ids();
ll group_num = ids.first;
vll counts(group_num);
for (auto x : ids.second) counts[x]++;
vvll groups(ids.first);
for (ll i = 0; i < group_num; i++) {
groups[i].reserve(counts[i]);
}
for (ll i = 0; i < _n; i++) {
groups[ids.second[i]].push_back(i);
}
return groups;
}
// cycleに入る頂点の数を求める。
ll getReachCycleCnt() {
auto sccgroups = scc();
umpll convedges;
vll cyclecheck(len(sccgroups));
rep(i, len(sccgroups)) {
rep(j, len(sccgroups[i])) {
convedges[sccgroups[i][j]] = i;
}
if (len(sccgroups[i]) > 1) cyclecheck[i] = 1;
}
auto Gr = Graph(len(sccgroups), true);
rep(i, len(edges)) {
ll U = edges[i].first, V = edges[i].second.to;
Gr.append_edge(convedges[V], convedges[U]);
}
rep(i, len(sccgroups)) {
if (cyclecheck[i]) Gr.bfs(i, false);
}
ll ans = 0;
rep(i, len(Gr.connect_vtxs_)) {
ans += len(sccgroups[Gr.connect_vtxs_[i]]);
}
return ans;
}
private:
ll _n;
struct edge {
ll to;
};
std::vector<std::pair<ll, edge>> edges;
};
class GraphLCA
{
private:
const ll V_;
const bool dir_;
bool is_build_exe_;
vector<vector<pll>> edges_;
public:
vll path; // 蟻本の vs、オイラーツアーを保持
vll depth; // 蟻本の depth、path[i] であるノードの深さを保持
vll in_order; // 蟻本の id、ノードiがオイラーツアーで最初に出てくるインデックスを保持
vector<pll> dat;
ll segn;
const pll INF = std::make_pair(LLONG_MAX, LLONG_MAX);
GraphLCA(ll V, bool dir) : V_(V), dir_(dir), is_build_exe_(false), edges_(V_), path(V * 2 - 1), depth(V * 2 - 1), in_order(V) {}
void append_edge(ll sv, ll ev, ll weight = 1)
{
vector<pll> &u = edges_[sv];
pll v = make_pair(ev, weight);
u.pb(v);
if (!dir_) {
swap(sv, ev);
vector<pll> &ru = edges_[sv];
pll rv = make_pair(ev, weight);
ru.pb(rv);
}
}
// 木を作るための根を入力
void build(ll root){
ll k = 0;
dfs(root, -1, 0, k);
// セグ木を構築、持つのはpair(depth, index) => depth が最小となる index がわかる
for (segn = 1; segn < V_ * 2 - 1; segn <<= 1)
;
dat.assign(segn * 2, INF);
rep(i, depth.size())
dat[segn + i] = std::make_pair(depth[i], i);
rrepd(i, segn - 1)
dat[i] = min(dat[i * 2], dat[i * 2 + 1]);
is_build_exe_ = true;
}
ll get(ll u, ll v) const
{
if (!is_build_exe_) {
debug("GraphLCA Build run required.");
assert(is_build_exe_);
}
ll l = std::min(in_order[u], in_order[v]);
ll r = std::max(in_order[u], in_order[v]) + 1;
return path[range_min(1, segn, l, r).second];
}
void dfs(ll sv, ll pre, ll d, ll &k)
{
// k: オイラーツアーの何番目かを保持する変数
in_order[sv] = k;
path[k] = sv;
depth[k++] = d;
for (const auto& [ev, w] : edges_[sv])
{
if (ev == pre) continue;
dfs(ev, sv, d + 1, k);
// ここに来た時はノードvの子であるe.dstを展開し終わってvに戻ってきたときなので、再度 path と depth に記録する
path[k] = sv;
depth[k++] = d;
}
}
// v : いまみてるノード、w: 今見てるノードに対応する区間長 l: ? r: ?
pll range_min(ll v, ll w, ll l, ll r) const
{
if (r <= l || w == 0)
return INF;
if (r - l == w)
return dat[v];
ll m = w / 2;
auto rmin = range_min(v * 2, m, l, min(r, m));
auto lmin = range_min(v * 2 + 1, m, max(_0, l - m), r - m);
return min(rmin, lmin);
}
};
class GridMap
{
private:
unordered_map<ll, vector<pll>> edges_;
const ll H_;
const ll W_;
const ll ansmod_;
vvll gridmap_;
// const vvll dirHWWeight_ = {{1, 1, 1}, {-1, -1, 1}, {1, -1, 1}, {-1, 1, 1}}; // 十字
const vvll dirHWWeight_ = {{1, 0, 1}, {-1, 0, 1}, {0, 1, 1}, {0, -1, 1}}; // 十字
// const vvll dirHWWeight_ = {{1, 0, 1}, {-1, 0, 1}, {0, 1, 1}, {0, -1, 1},
// {1, 1, 1}, {-1, 1, 1}, {1, -1, 1}, {-1, -1, 1}}; // 全方向
vvll build(ll H, ll W) {
vvll tbl = vvll(H, vll(W, -1));
rep(i, H) {
string s = in_str();
rep(j, W) {
tbl[i][j] = s[j] != '#' ? 0 : -1;
s[j] != '#' ? notwall_++ : wall_++;
if (s[j] == 'S' || s[j] == 's') startpos_ = {i, j};
if (s[j] == 'G' || s[j] == 'g') endpos_ = {i, j};
if (s[j] == 'o') po.pb(vll{i, j});
}
}
return tbl;
}
public:
vvll dist_;
vvvll dirdist_;
vvll cnt_;
vvll connect_vtxs_;
vll startpos_;
vll endpos_;
ll wall_, notwall_;
vvll po;
GridMap(ll H, ll W, ll ansmod = 1000000007) : H_(H), W_(W), ansmod_(ansmod), wall_(0), notwall_(0){
gridmap_ = build(H, W);
reset();
}
void reset() {
dist_ = vvll(H_, vll(W_, LLONG_MAX));
dirdist_ = vvvll(H_, vvll(W_, vll(len(dirHWWeight_), LLONG_MAX)));
cnt_ = vvll(H_, vll(W_));
connect_vtxs_.clear();
}
void bfs01(ll sh, ll sw, ll search_depth = LLONG_MAX) {
if (gridmap_[sh][sw] == -1) return;
if (dist_[sh][sw] != LLONG_MAX) return;
dist_[sh][sw] = 0;
deque<vll> q = deque<vll>();
rep(i, len(dirHWWeight_)) {
q.pb(vll{sh, sw, i, 1});
dirdist_[sh][sw][i] = 1;
}
connect_vtxs_.pb(vll{sh, sw});
cnt_[sh][sw] = 1;
while (len(q) > 0) {
auto qf = q[0];
vll ps = {qf[0], qf[1]};
ll pidx = qf[2];
q.popleft();
rep(i, len(dirHWWeight_)) {
const auto &dir = *(dirHWWeight_.cbegin() + i);
ll dh = dir[0], dw = dir[1], dweight = i == pidx ? 0 : dir[2];
ll posh = ps[0] + dh, posw = ps[1] + dw;
if (!IN(0, posh, H_ - 1) || !IN(0, posw, W_ - 1)) continue;
if (gridmap_[posh][posw] == -1) continue;
ll td = qf[2] + dweight;
connect_vtxs_.pb(vll{posh, posw});
if (dirdist_[posh][posw][i] == dirdist_[ps[0]][ps[1]][i] + dweight) {
cnt_[posh][posw] += cnt_[ps[0]][ps[1]];
cnt_[posh][posw] %= ansmod_;
} else {
cnt_[posh][posw] = cnt_[ps[0]][ps[1]];
}
if (dirdist_[posh][posw][i] == LLONG_MAX || td < dirdist_[posh][posw][i]) {
dirdist_[posh][posw][i] = td;
if (td < search_depth) {
if (i == pidx) {
q.pf(vll{posh, posw, i, td});
} else {
q.pb(vll{posh, posw, i, td});
}
}
}
}
}
rep(i, H_) {
rep(j, W_) {
chmin(dist_[i][j], MIN(dirdist_[i][j]));
}
}
}
void bfs(ll sh, ll sw, ll search_depth = 1000000000) {
if (gridmap_[sh][sw] == -1) return;
if (dist_[sh][sw] != LLONG_MAX) return;
deque<vll> q;
q.pb(vll{sh, sw});
dist_[sh][sw] = 0;
connect_vtxs_.pb(vll{sh, sw});
cnt_[sh][sw] = 1;
while (len(q) > 0) {
vll ps = q[0];
q.popleft();
vvll edges;
for(const auto &dir: dirHWWeight_) {
ll dh = dir[0], dw = dir[1], dweight = dir[2];
ll posh = ps[0] + dh, posw = ps[1] + dw;
if (IN(0, posh, H_ -1) && IN(0, posw, W_ - 1)) {
if (gridmap_[posh][posw] == -1) continue;
if (dist_[posh][posw] >= dist_[ps[0]][ps[1]] + dweight) edges.pb(vll{posh, posw, dweight});
}
}
for (const auto &edge: edges) {
ll ph = edge[0], pw = edge[1], w = edge[2];
if (dist_[ph][pw] == LLONG_MAX || dist_[ph][pw] == dist_[ps[0]][ps[1]] + w) {
cnt_[ph][pw] += cnt_[ps[0]][ps[1]];
cnt_[ph][pw] %= ansmod_;
}
if (dist_[ph][pw] <= dist_[ps[0]][ps[1]] + w) continue;
dist_[ph][pw] = dist_[ps[0]][ps[1]] + w;
cnt_[ph][pw] = cnt_[ps[0]][ps[1]];
connect_vtxs_.pb(vll{ph, pw});
if (dist_[ph][pw] < search_depth)
{
if (w == 0) q.pf(vll{ph, pw});
else q.pb(vll{ph, pw});
}
}
}
}
};
// 巡回セールスマン問題(TSP)
// 始点0, 中継点1～N, 終点N + 1という番号付けで対応。
// 2点間の最短距離を保存した行列の要素があるdistsを引数に与える。
// 各頂点のコスト(どの頂点を尋ねたかの種類ごと)と、始点から終点まで最大何個経由したかの組みを返す。
// iscycleは始点に戻ってきたときの情報を返す。
template <typename T>
pair<vector<vector<T>>, ll> tsp(ll startidx, const vector<vector<T>> &dists, bool iscycle = false, ll costlimit = 1000000007){
ll ans = 0;
ll N = len(dists) - 1; // 終点を除く個数
ll allp = _1 << N;
vector<vector<T>> dp(N + 1, vll(allp, LLONG_MAX));
deque<vll> q;
if (!iscycle) {
q.pb(vll{startidx, _1 << startidx}); // [pos, snum]
dp[startidx][_1 << startidx] = 0;
} else {
q.pb(vll{startidx, 0});
dp[startidx][0] = 0;
}
while(len(q)) {
deque<vll> newq;
while (len(q)) {
vll data = q[0];
INI2(nowpos, nowsnum, data);
q.pop_front();
rep (i, N) {
if (nowpos == i) continue;
ll nextsnum = nowsnum | (_1 << i);
T dist = dists[nowpos][i];
if (dp[i][nextsnum] <= dp[nowpos][nowsnum] + dist) continue;
chmin(dp[i][nextsnum], dp[nowpos][nowsnum] + dist);
if (costlimit < dp[nowpos][nowsnum] + dist + dists[i][N]) continue;
chmin(dp[N][nextsnum], dp[nowpos][nowsnum] + dist + dists[i][N]);
newq.pb(vll{i, nextsnum});
chmax(ans, popcnt(nextsnum) - 1);
}
}
q = newq;
}
return mp(dp, ans);
}
class Dijkstra{
class Edge{
public:
ll to_;
ll cost_;
Edge(ll to, ll cost) : to_(to), cost_(cost) {}
};
vector<vector<Edge>> G_;
ll V_;
bool dir_;
const ll INF = LLONG_MAX;
public:
vll dist_;
vll roots_;
Dijkstra(ll V, bool dir = false) : G_(vector<vector<Edge>>(V)), V_(V), dir_(dir), dist_(vll(V, INF)), roots_(vll(V, -1)) {}
void append_edge(ll from, ll to, ll cost = 1) {
G_[from].pb(Edge(to, cost));
if (!dir_) {
G_[to].pb(Edge(from, cost));
}
}
void shortest_path(ll s) {
priority_queue<pll, vector<pll>, greater<pll>> que;
dist_[s] = 0;
que.push(make_pair(0, s));
while (true) {
if (len(que) == 0) break;
const auto[cost, v] = que.top(); que.pop();
if (dist_[v] < cost) continue;
rep (i, len(G_[v])) {
auto e = G_[v][i];
if (dist_[e.to_] > dist_[v] + e.cost_) {
dist_[e.to_] = dist_[v] + e.cost_;
roots_[e.to_] = v;
que.push(make_pair(dist_[e.to_], e.to_));
}
}
}
}
void reset() {
dist_ = vll(V_, INF);
roots_ = vll(V_, -1);
}
// 無向グラフのみ対応。
// start->endのパスが無ければ空を返す
vll get_path(ll start, ll end, ll vertexidx_offset = 0) {
shortest_path(start);
if (dist_[end] == INF) return vll{};
ll pos = end;
vll path = {end + vertexidx_offset};
while(pos != start) {
pos = roots_[pos];
path.pb(pos + vertexidx_offset);
}
REV(path);
return path;
}
};
// ワーシャルフロイド
class WarshallFloyd {
private:
unordered_map<ll, vvll> edges_;
unordered_map<ll, unordered_map<ll, sll>> ignore_edges_;
const ll V_;
const bool dir_;
const ll ansmod_;
public:
WarshallFloyd(ll V, bool dir, ll ansmod = LLONG_MAX) : V_(V), dir_(dir), ansmod_(ansmod){}
void append_edge(ll sv, ll ev, ll weight = 1)
{
vvll &u = edges_[sv];
vll v = {ev, weight};
u.pb(v);
if (!dir_) {
swap(sv, ev);
vvll &ru = edges_[sv];
vll rv = {ev, weight};
ru.pb(rv);
}
}
void ignore_edge(ll sv, ll ev, ll weight = 1) {
sll &uv = ignore_edges_[sv][ev];
uv.insert(weight);
if (!dir_) {
swap(sv, ev);
sll &ruv = ignore_edges_[sv][ev];
ruv.insert(weight);
}
}
void remove_ignore_edge(ll sv, ll ev, ll weight = 1) {
if (!EXIST(sv, ignore_edges_))
return;
if (!EXIST(ev, ignore_edges_[sv]))
return;
sll &uv = ignore_edges_[sv][ev];
uv.erase(weight);
if (!dir_) {
swap(sv, ev);
sll &ruv = ignore_edges_[sv][ev];
ruv.erase(weight);
}
}
vvll calcdist()
{
auto const INF = LLONG_MAX / 2;
vvll d(V_, vll(V_, INF));
rep(i, V_) d[i][i] = 0;
repdict(k, v, edges_) rep(i, len(v)) {
if (EXIST(k, ignore_edges_)) if (EXIST(v[i][0], ignore_edges_[k])) if (EXIST(v[i][1], ignore_edges_[k][v[i][0]])) continue;
d[k][v[i][0]] = std::min(d[k][v[i][0]], v[i][1]);
}
rep(k, V_) rep(i, V_) rep(j, V_)
{
if (d[i][k] != INF && d[k][j] != INF)
d[i][j] = std::min(d[i][j], d[i][k] + d[k][j]);
}
return d;
}
};
// https://atcoder.github.io/ac-library/production/document_ja/maxflow.html
template <class Cap> struct MaxFlow {
public:
MaxFlow() : _n(0) {}
MaxFlow(ll n) : _n(n), g(n) {}
ll append_edge(ll from, ll to, Cap cap) {
assert(0 <= from && from < _n);
assert(0 <= to && to < _n);
assert(0 <= cap);
ll m = ll(pos.size());
pos.push_back({from, ll(g[from].size())});
ll from_id = ll(g[from].size());
ll to_id = ll(g[to].size());
if (from == to) to_id++;
g[from].push_back(_edge{to, to_id, cap});
g[to].push_back(_edge{from, from_id, 0});
return m;
}
struct edge {
ll from, to;
Cap cap, flow;
};
edge get_edge(ll i) {
ll m = ll(pos.size());
assert(0 <= i && i < m);
auto _e = g[pos[i].first][pos[i].second];
auto _re = g[_e.to][_e.rev];
return edge{pos[i].first, _e.to, _e.cap + _re.cap, _re.cap};
}
std::vector<edge> edges() {
ll m = ll(pos.size());
std::vector<edge> result;
for (ll i = 0; i < m; i++) {
result.push_back(get_edge(i));
}
return result;
}
void change_edge(ll i, Cap new_cap, Cap new_flow) {
ll m = ll(pos.size());
assert(0 <= i && i < m);
assert(0 <= new_flow && new_flow <= new_cap);
auto& _e = g[pos[i].first][pos[i].second];
auto& _re = g[_e.to][_e.rev];
_e.cap = new_cap - new_flow;
_re.cap = new_flow;
}
Cap flow(ll s, ll t) {
return flow(s, t, std::numeric_limits<Cap>::max());
}
Cap flow(ll s, ll t, Cap flow_limit) {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);
std::vector<ll> level(_n), iter(_n);
internal::simple_queue<ll> que;
auto bfs = [&]() {
std::fill(level.begin(), level.end(), -1);
level[s] = 0;
que.clear();
que.push(s);
while (!que.empty()) {
ll v = que.front();
que.pop();
for (auto e : g[v]) {
if (e.cap == 0 || level[e.to] >= 0) continue;
level[e.to] = level[v] + 1;
if (e.to == t) return;
que.push(e.to);
}
}
};
auto dfs = [&](auto self, ll v, Cap up) {
if (v == s) return up;
Cap res = 0;
ll level_v = level[v];
for (ll& i = iter[v]; i < ll(g[v].size()); i++) {
_edge& e = g[v][i];
if (level_v <= level[e.to] || g[e.to][e.rev].cap == 0) continue;
Cap d =
self(self, e.to, std::min(up - res, g[e.to][e.rev].cap));
if (d <= 0) continue;
g[v][i].cap += d;
g[e.to][e.rev].cap -= d;
res += d;
if (res == up) break;
}
return res;
};
Cap flow = 0;
while (flow < flow_limit) {
bfs();
if (level[t] == -1) break;
std::fill(iter.begin(), iter.end(), 0);
while (flow < flow_limit) {
Cap f = dfs(dfs, t, flow_limit - flow);
if (!f) break;
flow += f;
}
}
return flow;
}
std::vector<bool> min_cut(ll s) {
std::vector<bool> visited(_n);
internal::simple_queue<ll> que;
que.push(s);
while (!que.empty()) {
ll p = que.front();
que.pop();
visited[p] = true;
for (auto e : g[p]) {
if (e.cap && !visited[e.to]) {
visited[e.to] = true;
que.push(e.to);
}
}
}
return visited;
}
private:
ll _n;
struct _edge {
ll to, rev;
Cap cap;
};
std::vector<std::pair<ll, ll>> pos;
std::vector<std::vector<_edge>> g;
};
template<class Cap, class Cost>
struct MinCostFlow {
public:
MinCostFlow() {}
MinCostFlow(ll n)
: _n(n)
, g(n) {}
ll append_edge(ll from, ll to, Cap cap, Cost cost) {
assert(0 <= from && from < _n);
assert(0 <= to && to < _n);
ll m = ll(pos.size());
pos.push_back({from, ll(g[from].size())});
ll from_id = ll(g[from].size());
ll to_id = ll(g[to].size());
if (from == to) to_id++;
g[from].push_back(_edge{to, to_id, cap, cost});
g[to].push_back(_edge{from, from_id, 0, -cost});
return m;
}
struct edge {
ll from, to;
Cap cap, flow;
Cost cost;
};
edge get_edge(ll i) {
ll m = ll(pos.size());
assert(0 <= i && i < m);
auto _e = g[pos[i].first][pos[i].second];
auto _re = g[_e.to][_e.rev];
return edge{
pos[i].first,
_e.to,
_e.cap + _re.cap,
_re.cap,
_e.cost,
};
}
vector<edge> edges() {
ll m = ll(pos.size());
vector<edge> result(m);
rep (i, m) result[i] = get_edge(i);
return result;
}
pair<Cap, Cost> flow(ll s, ll t) {
return flow(s, t, std::numeric_limits<Cap>::max());
}
pair<Cap, Cost> flow(ll s, ll t, Cap flow_limit) {
return slope(s, t, flow_limit).back();
}
vector<std::pair<Cap, Cost>> slope(ll s, ll t) {
return slope(s, t, std::numeric_limits<Cap>::max());
}
vector<std::pair<Cap, Cost>> slope(ll s, ll t, Cap flow_limit) {
assert(0 <= s && s < _n);
assert(0 <= t && t < _n);
assert(s != t);
vector<Cost> dual(_n, 0), dist(_n);
vll pv(_n), pe(_n);
vector<bool> vis(_n);
auto dual_ref = [&]() {
std::fill(dist.begin(), dist.end(),
std::numeric_limits<Cost>::max());
std::fill(pv.begin(), pv.end(), -1);
std::fill(pe.begin(), pe.end(), -1);
std::fill(vis.begin(), vis.end(), false);
struct Q {
Cost key;
ll to;
bool operator<(Q r) const { return key > r.key; }
};
std::priority_queue<Q> que;
dist[s] = 0;
que.push(Q{0, s});
while (!que.empty()) {
ll v = que.top().to;
que.pop();
if (vis[v]) continue;
vis[v] = true;
if (v == t) break;
rep (i, len(g[v])) {
auto e = g[v][i];
if (vis[e.to] || !e.cap) continue;
Cost cost = e.cost - dual[e.to] + dual[v];
if (dist[e.to] - dist[v] > cost) {
dist[e.to] = dist[v] + cost;
pv[e.to] = v;
pe[e.to] = i;
que.push(Q{dist[e.to], e.to});
}
}
}
if (!vis[t]) {
return false;
}
rep (v, _n) {
if (!vis[v]) continue;
dual[v] -= dist[t] - dist[v];
}
return true;
};
Cap flow = 0;
Cost cost = 0, prev_cost_per_flow = -1;
vector<pair<Cap, Cost>> result;
result.push_back({flow, cost});
while (flow < flow_limit) {
if (!dual_ref()) break;
Cap c = flow_limit - flow;
for (ll v = t; v != s; v = pv[v]) {
c = min(c, g[pv[v]][pe[v]].cap);
}
for (ll v = t; v != s; v = pv[v]) {
auto &e = g[pv[v]][pe[v]];
e.cap -= c;
g[v][e.rev].cap += c;
}
Cost d = -dual[s];
flow += c;
cost += c * d;
if (prev_cost_per_flow == d) {
result.pop_back();
}
result.push_back({flow, cost});
prev_cost_per_flow = d;
}
return result;
}
private:
ll _n;
struct _edge {
ll to, rev;
Cap cap;
Cost cost;
};
std::vector<std::pair<ll, ll>> pos;
std::vector<std::vector<_edge>> g;
};
template <class T> struct fenwick_tree {
public:
fenwick_tree() : _n(0) {}
fenwick_tree(int n) : _n(n), data(n) {}
void add(int p, T x) {
assert(0 <= p && p < _n);
p++;
while (p <= _n) {
data[p - 1] += x;
p += p & -p;
}
}
T sum(int l, int r) { // [l, r)
assert(0 <= l && l <= r && r <= _n);
return l == 0 ? sum(r) : sum(r) - sum(l);
}
private:
int _n;
std::vector<T> data;
T sum(int r) {
T s = 0;
while (r > 0) {
s += data[r - 1];
r -= r & -r;
}
return s;
}
};
template<class T>
class Bit
{
public:
fenwick_tree<T> ft;
ll n_;
Bit(ll n = 0) : ft(n), n_(n) {}
T sum(ll i) { // [0, i)
return ft.sum(0, i);
}
T sum(ll l, ll r) { // [l, r)
return ft.sum(l, r);
}
void add(ll i, T x) {
ft.add(i, x);
}
void set(ll p, T x) {
ft.add(p, -ft.sum(p, p + 1) + x);
}
// Bitの使い方として個数をaddしていったときにidx番目の値を取得できる。
T getval(ll idx) { // 0-origin, idx番目の値を求める
ll l = 0, r = n_;
while(r - l > 1) {
ll mid = (r + l) >> 1;
ll cnt = sum(mid);
if (cnt <= idx) {
l = mid;
} else {
r = mid;
}
}
return l;
}
T inversion(const vector<T> &vec) {
ft = fenwick_tree<T>(MAX(vec) + 1);
T val = 0;
rep(i, len(vec)) {
ft.add(vec[i], 1);
val += i + 1 - sum(vec[i] + 1);
}
return val;
}
T big_inversion(const vector<T> &vec) {
auto [d, revd] = compcoord(vec);
vector<T> v;
rep(i, len(vec)) {
v.pb(d[vec[i]]);
}
return inversion(v);
}
};
// マンハッタン
vll xy_Manhattan(const vll &p) {
return vll{p[0] - p[1], p[0] + p[1]};
}
template <typename T>
pair<unordered_map<T, ll>, vector<T>> compcoord(const vector<T> &vec)
{
set<T> s = set<T>(all(vec));
unordered_map<T, ll> d;
ll idx = 0;
repset (v, s) d[v] = idx++;
vector<T> revd = vector<T>(len(s));
repdict(k, v, d) revd[v] = k;
return make_pair(d, revd);
}
// https://howardhinnant.github.io/combinations.html
namespace howardhinnantdetail
{
// Rotates two discontinuous ranges to put *first2 where *first1 is.
// If last1 == first2 this would be equivalent to rotate(first1, first2, last2),
// but instead the rotate "jumps" over the discontinuity [last1, first2) -
// which need not be a valid range.
// In order to make it faster, the length of [first1, last1) is passed in as d1,
// and d2 must be the length of [first2, last2).
// In a perfect world the d1 > d2 case would have used swap_ranges and
// reverse_iterator, but reverse_iterator is too inefficient.
template <class BidirIter>
void
rotate_discontinuous(BidirIter first1, BidirIter last1,
typename std::iterator_traits<BidirIter>::difference_type d1,
BidirIter first2, BidirIter last2,
typename std::iterator_traits<BidirIter>::difference_type d2)
{
using std::swap;
if (d1 <= d2)
std::rotate(first2, std::swap_ranges(first1, last1, first2), last2);
else
{
BidirIter i1 = last1;
while (first2 != last2)
swap(*--i1, *--last2);
std::rotate(first1, i1, last1);
}
}
// Rotates the three discontinuous ranges to put *first2 where *first1 is.
// Just like rotate_discontinuous, except the second range is now represented by
// two discontinuous ranges: [first2, last2) + [first3, last3).
template <class BidirIter>
void
rotate_discontinuous3(BidirIter first1, BidirIter last1,
typename std::iterator_traits<BidirIter>::difference_type d1,
BidirIter first2, BidirIter last2,
typename std::iterator_traits<BidirIter>::difference_type d2,
BidirIter first3, BidirIter last3,
typename std::iterator_traits<BidirIter>::difference_type d3)
{
rotate_discontinuous(first1, last1, d1, first2, last2, d2);
if (d1 <= d2)
rotate_discontinuous(std::next(first2, d2 - d1), last2, d1, first3, last3, d3);
else
{
rotate_discontinuous(std::next(first1, d2), last1, d1 - d2, first3, last3, d3);
rotate_discontinuous(first2, last2, d2, first3, last3, d3);
}
}
// Call f() for each combination of the elements [first1, last1) + [first2, last2)
// swapped/rotated into the range [first1, last1). As long as f() returns
// false, continue for every combination and then return [first1, last1) and
// [first2, last2) to their original state. If f() returns true, return
// immediately.
// Does the absolute mininum amount of swapping to accomplish its task.
// If f() always returns false it will be called (d1+d2)!/(d1!*d2!) times.
template <class BidirIter, class Function>
bool
combine_discontinuous(BidirIter first1, BidirIter last1,
typename std::iterator_traits<BidirIter>::difference_type d1,
BidirIter first2, BidirIter last2,
typename std::iterator_traits<BidirIter>::difference_type d2,
Function& f,
typename std::iterator_traits<BidirIter>::difference_type d = 0)
{
typedef typename std::iterator_traits<BidirIter>::difference_type D;
using std::swap;
if (d1 == 0 || d2 == 0)
return f();
if (d1 == 1)
{
for (BidirIter i2 = first2; i2 != last2; ++i2)
{
if (f())
return true;
swap(*first1, *i2);
}
}
else
{
BidirIter f1p = std::next(first1);
BidirIter i2 = first2;
for (D d22 = d2; i2 != last2; ++i2, --d22)
{
if (combine_discontinuous(f1p, last1, d1-1, i2, last2, d22, f, d+1))
return true;
swap(*first1, *i2);
}
}
if (f())
return true;
if (d != 0)
rotate_discontinuous(first1, last1, d1, std::next(first2), last2, d2-1);
else
rotate_discontinuous(first1, last1, d1, first2, last2, d2);
return false;
}
// A binder for binding arguments to call combine_discontinuous
template <class Function, class BidirIter>
class call_combine_discontinuous
{
typedef typename std::iterator_traits<BidirIter>::difference_type D;
Function f_;
BidirIter first1_;
BidirIter last1_;
D d1_;
BidirIter first2_;
BidirIter last2_;
D d2_;
public:
call_combine_discontinuous(
BidirIter first1, BidirIter last1,
D d1,
BidirIter first2, BidirIter last2,
D d2,
Function& f)
: f_(f), first1_(first1), last1_(last1), d1_(d1),
first2_(first2), last2_(last2), d2_(d2) {}
bool operator()()
{
return combine_discontinuous(first1_, last1_, d1_, first2_, last2_, d2_, f_);
}
};
// See combine_discontinuous3
template <class BidirIter, class Function>
bool
combine_discontinuous3_(BidirIter first1, BidirIter last1,
typename std::iterator_traits<BidirIter>::difference_type d1,
BidirIter first2, BidirIter last2,
typename std::iterator_traits<BidirIter>::difference_type d2,
BidirIter first3, BidirIter last3,
typename std::iterator_traits<BidirIter>::difference_type d3,
Function& f,
typename std::iterator_traits<BidirIter>::difference_type d = 0)
{
typedef typename std::iterator_traits<BidirIter>::difference_type D;
using std::swap;
if (d1 == 1)
{
for (BidirIter i2 = first2; i2 != last2; ++i2)
{
if (f())
return true;
swap(*first1, *i2);
}
if (f())
return true;
swap(*first1, *std::prev(last2));
swap(*first1, *first3);
for (BidirIter i2 = std::next(first3); i2 != last3; ++i2)
{
if (f())
return true;
swap(*first1, *i2);
}
}
else
{
BidirIter f1p = std::next(first1);
BidirIter i2 = first2;
for (D d22 = d2; i2 != last2; ++i2, --d22)
{
if (combine_discontinuous3_(f1p, last1, d1-1, i2, last2, d22, first3,
last3, d3, f, d+1))
return true;
swap(*first1, *i2);
}
i2 = first3;
for (D d22 = d3; i2 != last3; ++i2, --d22)
{
if (combine_discontinuous(f1p, last1, d1-1, i2, last3, d22, f, d+1))
return true;
swap(*first1, *i2);
}
}
if (f())
return true;
if (d1 == 1)
swap(*std::prev(last2), *first3);
if (d != 0)
{
if (d2 > 1)
rotate_discontinuous3(first1, last1, d1, std::next(first2), last2, d2-1, first3, last3, d3);
else
rotate_discontinuous(first1, last1, d1, first3, last3, d3);
}
else
rotate_discontinuous3(first1, last1, d1, first2, last2, d2, first3, last3, d3);
return false;
}
// Like combine_discontinuous, but swaps/rotates each combination out of
// [first1, last1) + [first2, last2) + [first3, last3) into [first1, last1).
// If f() always returns false, it is called (d1+d2+d3)!/(d1!*(d2+d3)!) times.
template <class BidirIter, class Function>
bool
combine_discontinuous3(BidirIter first1, BidirIter last1,
typename std::iterator_traits<BidirIter>::difference_type d1,
BidirIter first2, BidirIter last2,
typename std::iterator_traits<BidirIter>::difference_type d2,
BidirIter first3, BidirIter last3,
typename std::iterator_traits<BidirIter>::difference_type d3,
Function& f)
{
typedef call_combine_discontinuous<Function&, BidirIter> F;
F fbc(first2, last2, d2, first3, last3, d3, f); // BC
return combine_discontinuous3_(first1, last1, d1, first2, last2, d2, first3, last3, d3, fbc);
}
// See permute
template <class BidirIter, class Function>
bool
permute_(BidirIter first1, BidirIter last1,
typename std::iterator_traits<BidirIter>::difference_type d1,
Function& f)
{
using std::swap;
switch (d1)
{
case 0:
case 1:
return f();
case 2:
if (f())
return true;
swap(*first1, *std::next(first1));
return f();
case 3:
{
if (f())
return true;
BidirIter f2 = std::next(first1);
BidirIter f3 = std::next(f2);
swap(*f2, *f3);
if (f())
return true;
swap(*first1, *f3);
swap(*f2, *f3);
if (f())
return true;
swap(*f2, *f3);
if (f())
return true;
swap(*first1, *f2);
swap(*f2, *f3);
if (f())
return true;
swap(*f2, *f3);
return f();
}
}
BidirIter fp1 = std::next(first1);
for (BidirIter p = fp1; p != last1; ++p)
{
if (permute_(fp1, last1, d1-1, f))
return true;
std::reverse(fp1, last1);
swap(*first1, *p);
}
return permute_(fp1, last1, d1-1, f);
}
// Calls f() for each permutation of [first1, last1)
// Divided into permute and permute_ in a (perhaps futile) attempt to
// squeeze a little more performance out of it.
template <class BidirIter, class Function>
bool
permute(BidirIter first1, BidirIter last1,
typename std::iterator_traits<BidirIter>::difference_type d1,
Function& f)
{
using std::swap;
switch (d1)
{
case 0:
case 1:
return f();
case 2:
{
if (f())
return true;
BidirIter i = std::next(first1);
swap(*first1, *i);
if (f())
return true;
swap(*first1, *i);
}
break;
case 3:
{
if (f())
return true;
BidirIter f2 = std::next(first1);
BidirIter f3 = std::next(f2);
swap(*f2, *f3);
if (f())
return true;
swap(*first1, *f3);
swap(*f2, *f3);
if (f())
return true;
swap(*f2, *f3);
if (f())
return true;
swap(*first1, *f2);
swap(*f2, *f3);
if (f())
return true;
swap(*f2, *f3);
if (f())
return true;
swap(*first1, *f3);
}
break;
default:
BidirIter fp1 = std::next(first1);
for (BidirIter p = fp1; p != last1; ++p)
{
if (permute_(fp1, last1, d1-1, f))
return true;
std::reverse(fp1, last1);
swap(*first1, *p);
}
if (permute_(fp1, last1, d1-1, f))
return true;
std::reverse(first1, last1);
break;
}
return false;
}
// Creates a functor with no arguments which calls f_(first_, last_).
// Also has a variant that takes two It and ignores them.
template <class Function, class It>
class bound_range
{
Function f_;
It first_;
It last_;
public:
bound_range(Function f, It first, It last)
: f_(f), first_(first), last_(last) {}
bool
operator()()
{
return f_(first_, last_);
}
bool
operator()(It, It)
{
return f_(first_, last_);
}
};
// A binder for binding arguments to call permute
template <class Function, class It>
class call_permute
{
typedef typename std::iterator_traits<It>::difference_type D;
Function f_;
It first_;
It last_;
D d_;
public:
call_permute(Function f, It first, It last, D d)
: f_(f), first_(first), last_(last), d_(d) {}
bool
operator()()
{
return permute(first_, last_, d_, f_);
}
};
} // detail
template <class BidirIter, class Function>
Function
for_each_combination(BidirIter first, BidirIter last,
BidirIter mid, Function f)
{
howardhinnantdetail::bound_range<Function&, BidirIter> wfunc(f, first, mid);
howardhinnantdetail::combine_discontinuous(first, mid, std::distance(first, mid),
mid, last, std::distance(mid, last),
wfunc);
return f;
}
class BitwiseFullSearch
{
public:
ll count_;
std::function<void(const vll &ptn, ll &count)> checkcount_; // カウントするロジックのラムダ式を突っ込む。
BitwiseFullSearch(std::function<void(const vll &ptn, ll &count)> f) : count_(0), checkcount_(f) {}
// ここは触らなくてよい(パターンを列挙しているだけ)
bool operator()(vll::iterator it1, vll::iterator it2)
{
vll ptn;
while (it1 != it2)
{
ptn.pb(*it1);
++it1;
}
checkcount_(ptn, count_);
return false;
}
};
ll _comb_ptn_count(ll R, const std::function<void(const vll &ptn, ll &count)> &f, vll &_A_) {
auto B = BitwiseFullSearch(f);
vll::iterator _R_ = _A_.begin() + R;
B = for_each_combination(all(_A_), _R_, B);
return B.count_;
}
ll comb_ptn_count(ll N, ll R, const std::function<void(const vll &ptn, ll &count)> &f) {
SETPERM(_A_, N);
return _comb_ptn_count(R, f, _A_);
}
vvll get_comb_ptn(ll N, ll R) {
vvll cb;
auto f = [&](const vll &ptn, ll &count)
{
UNUSED(count);
cb.pb(ptn);
};
comb_ptn_count(N, R, f);
return cb;
}
ll comb_allptn_count(ll N, const std::function<void(const vll &ptn, ll &count)> &f) {
ll cnt = 0;
SETPERM(_A_, N);
rep(r, N + 1) {
cnt += _comb_ptn_count(r, f, _A_);
}
return cnt;
}
// N が 20とかだとSegmentation faultが発生する。スタック領域が足りなくなるか。
vvll get_comb_allptn(ll N) {
vvll cb;
auto f = [&](const vll &ptn, ll &count)
{
UNUSED(count);
cb.pb(ptn);
};
comb_allptn_count(N, f);
return cb;
}
vvll get_perm_ptn(ll N) {
vvll ptn;
SETPERM(_A_, N);
do {
ptn.pb(_A_);
} while(next_permutation(all(_A_)));
return ptn;
}
// n個からr個有効にしたbitパターンの数値を列挙
// 下1桁目から0, 1, 2番目が有効無効を表す。(Nは60くらいまで)
vll get_bitptn(ll N, ll R) {
vll bitptns;
auto ptns = get_comb_ptn(N, R);
rep (i, len(ptns)) {
ll bitptn = 0;
rep (j, R) {
bitptn += _1 << ptns[i][j];
}
bitptns.pb(bitptn);
}
return bitptns;
}
void _ptn_distribute_kval(ll N, ll K, vvll &addptn, vll &ptn, ll idx = 0) {
if (N == idx) {
addptn.pb(ptn);
return;
}
rep(i, K + 1) {
ptn[idx] = i;
_ptn_distribute_kval(N, K - i, addptn, ptn, idx + 1);
}
return;
}
// N個の配列に整数K以下を各要素に分配するパターンを全通り列挙
// N = 2, K = 3のとき
// ptn_dist = {{0, 0}, {0, 1}, {0, 2}, {0, 3}, }
// {1, 0}, {1, 1}, {1, 2},
// {2, 0}, {2, 1},
// {3, 0}}
vvll ptn_distribute_kval(ll N, ll K) {
vvll ptn_dist;
vll ptn(N);
_ptn_distribute_kval(N, K, ptn_dist, ptn, 0);
return ptn_dist;
}
// vecに含まれる数列をk個取り出したときの総和の全列挙
// vec = {1, 3, 4, 6}, K = 3のとき
// ret = {8, 10, 11, 13}
vll ptn_sumval(const vll &vec, ll k) {
auto ptns = get_comb_ptn(len(vec), k);
vll ret;
rep (i, len(ptns)) {
ll tmp = 0;
rep (j, len(ptns[i])) {
tmp += vec[ptns[i][j]];
}
ret.pb(tmp);
}
SORT(ret);
return ret;
}
template <class S, S (*op)(S, S), S (*e)()> struct segtree {
public:
segtree() : segtree(0) {}
segtree(int n) : segtree(std::vector<S>(n, e())) {}
segtree(const std::vector<S>& v) : _n(int(v.size())) {
log = internal::ceil_pow2(_n);
size = 1 << log;
d = std::vector<S>(2 * size, e());
for (int i = 0; i < _n; i++) d[size + i] = v[i];
for (int i = size - 1; i >= 1; i--) {
update(i);
}
}
void set(int p, S x) { // 0-origin
assert(0 <= p && p < _n);
p += size;
d[p] = x;
for (int i = 1; i <= log; i++) update(p >> i);
}
S get(int p) { // 0-origin
assert(0 <= p && p < _n);
return d[p + size];
}
S prod(int l, int r) { // 0-origin [l, r)
assert(0 <= l && l <= r && r <= _n);
S sml = e(), smr = e();
l += size;
r += size;
while (l < r) {
if (l & 1) sml = op(sml, d[l++]);
if (r & 1) smr = op(d[--r], smr);
l >>= 1;
r >>= 1;
}
return op(sml, smr);
}
S all_prod() { return d[1]; }
template <bool (*f)(S)> int max_right(int l) { // 引数も戻り値も0-origin
return max_right(l, [](S x) { return f(x); });
}
// lより右側でfを満たす最大の右位置を取得
template <class F> int max_right(int l, F f) { // 引数も戻り値も0-origin
assert(0 <= l && l <= _n);
assert(f(e()));
if (l == _n) return _n;
l += size;
S sm = e();
do {
while (l % 2 == 0) l >>= 1;
if (!f(op(sm, d[l]))) {
while (l < size) {
l = (2 * l);
if (f(op(sm, d[l]))) {
sm = op(sm, d[l]);
l++;
}
}
return l - size;
}
sm = op(sm, d[l]);
l++;
} while ((l & -l) != l);
return _n;
}
template <bool (*f)(S)> int min_left(int r) { // 引数も戻り値も0-origin
return min_left(r, [](S x) { return f(x); });
}
// rより左側でfを満たす最小の左位置を取得
template <class F> int min_left(int r, F f) { // 引数も戻り値も0-origin
assert(0 <= r && r <= _n);
assert(f(e()));
if (r == 0) return 0;
r += size;
S sm = e();
do {
r--;
while (r > 1 && (r % 2)) r >>= 1;
if (!f(op(d[r], sm))) {
while (r < size) {
r = (2 * r + 1);
if (f(op(d[r], sm))) {
sm = op(d[r], sm);
r--;
}
}
return r + 1 - size;
}
sm = op(d[r], sm);
} while ((r & -r) != r);
return 0;
}
private:
int _n, size, log;
std::vector<S> d;
void update(int k) { d[k] = op(d[2 * k], d[2 * k + 1]); }
};
template <class S,
S (*op)(S, S),
S (*e)(),
class F,
S (*mapping)(F, S),
F (*composition)(F, F),
F (*id)()>
struct lazy_segtree {
public:
lazy_segtree() : lazy_segtree(0) {}
lazy_segtree(int n) : lazy_segtree(std::vector<S>(n, e())) {}
lazy_segtree(const std::vector<S>& v) : _n(int(v.size())) {
log = internal::ceil_pow2(_n);
size = 1 << log;
d = std::vector<S>(2 * size, e());
lz = std::vector<F>(size, id());
for (int i = 0; i < _n; i++) d[size + i] = v[i];
for (int i = size - 1; i >= 1; i--) {
update(i);
}
}
void set(int p, S x) {
assert(0 <= p && p < _n);
p += size;
for (int i = log; i >= 1; i--) push(p >> i);
d[p] = x;
for (int i = 1; i <= log; i++) update(p >> i);
}
S get(int p) {
assert(0 <= p && p < _n);
p += size;
for (int i = log; i >= 1; i--) push(p >> i);
return d[p];
}
S prod(int l, int r) {
assert(0 <= l && l <= r && r <= _n);
if (l == r) return e();
l += size;
r += size;
for (int i = log; i >= 1; i--) {
if (((l >> i) << i) != l) push(l >> i);
if (((r >> i) << i) != r) push(r >> i);
}
S sml = e(), smr = e();
while (l < r) {
if (l & 1) sml = op(sml, d[l++]);
if (r & 1) smr = op(d[--r], smr);
l >>= 1;
r >>= 1;
}
return op(sml, smr);
}
S all_prod() { return d[1]; }
void apply(int p, F f) {
assert(0 <= p && p < _n);
p += size;
for (int i = log; i >= 1; i--) push(p >> i);
d[p] = mapping(f, d[p]);
for (int i = 1; i <= log; i++) update(p >> i);
}
void apply(int l, int r, F f) {
assert(0 <= l && l <= r && r <= _n);
if (l == r) return;
l += size;
r += size;
for (int i = log; i >= 1; i--) {
if (((l >> i) << i) != l) push(l >> i);
if (((r >> i) << i) != r) push((r - 1) >> i);
}
{
int l2 = l, r2 = r;
while (l < r) {
if (l & 1) all_apply(l++, f);
if (r & 1) all_apply(--r, f);
l >>= 1;
r >>= 1;
}
l = l2;
r = r2;
}
for (int i = 1; i <= log; i++) {
if (((l >> i) << i) != l) update(l >> i);
if (((r >> i) << i) != r) update((r - 1) >> i);
}
}
template <bool (*g)(S)> int max_right(int l) {
return max_right(l, [](S x) { return g(x); });
}
template <class G> int max_right(int l, G g) {
assert(0 <= l && l <= _n);
assert(g(e()));
if (l == _n) return _n;
l += size;
for (int i = log; i >= 1; i--) push(l >> i);
S sm = e();
do {
while (l % 2 == 0) l >>= 1;
if (!g(op(sm, d[l]))) {
while (l < size) {
push(l);
l = (2 * l);
if (g(op(sm, d[l]))) {
sm = op(sm, d[l]);
l++;
}
}
return l - size;
}
sm = op(sm, d[l]);
l++;
} while ((l & -l) != l);
return _n;
}
template <bool (*g)(S)> int min_left(int r) {
return min_left(r, [](S x) { return g(x); });
}
template <class G> int min_left(int r, G g) {
assert(0 <= r && r <= _n);
assert(g(e()));
if (r == 0) return 0;
r += size;
for (int i = log; i >= 1; i--) push((r - 1) >> i);
S sm = e();
do {
r--;
while (r > 1 && (r % 2)) r >>= 1;
if (!g(op(d[r], sm))) {
while (r < size) {
push(r);
r = (2 * r + 1);
if (g(op(d[r], sm))) {
sm = op(d[r], sm);
r--;
}
}
return r + 1 - size;
}
sm = op(d[r], sm);
} while ((r & -r) != r);
return 0;
}
private:
int _n, size, log;
std::vector<S> d;
std::vector<F> lz;
void update(int k) { d[k] = op(d[2 * k], d[2 * k + 1]); }
void all_apply(int k, F f) {
d[k] = mapping(f, d[k]);
if (k < size) lz[k] = composition(f, lz[k]);
}
void push(int k) {
all_apply(2 * k, lz[k]);
all_apply(2 * k + 1, lz[k]);
lz[k] = id();
}
};
// セグ木パターン
ll E0() { return (ll)0; }
ll E1() { return (ll)0; }
ll EM1() { return (ll)-1; }
ll EMAX() { return LLONG_MAX; }
ll EMIN() { return _0; }
ll EOR() { return (ll)0; }
ll EXOR() { return (ll)0; }
ll EAND() { return (ll)(((ll)1 << (ll)62) - (ll)1); }
ll ESUM() { return (ll)0; }
ll segmax(ll a, ll b) { return max(a, b); }
ll segmin(ll a, ll b) { return min(a, b); }
ll segand(ll a, ll b) { return a & b; }
ll segor(ll a, ll b) { return a | b; }
ll segxor(ll a, ll b) { return a ^ b; }
ll segsum(ll a, ll b) { return a + b; }
using SegTreeGCD = segtree<ll, gcd, E0>;
using SegTreeMaxE0 = segtree<ll, segmax, E0>; // 設定する要素の最小値が1のとき
using SegTreeMaxE1 = segtree<ll, segmax, E1>; // 設定する要素の最小値が2のとき
using SegTreeMaxEM1 = segtree<ll, segmax, EM1>; // 設定する要素の最小値が0のとき
using SegTreeMaxEMIN = segtree<ll, segmax, EMIN>;
using SegTreeMinEMAX = segtree<ll, segmin, EMAX>;
using SegTreeAnd = segtree<ll, segand, EAND>;
using SegTreeOr = segtree<ll, segor, EOR>;
using SegTreeXor = segtree<ll, segxor, EXOR>;
using SegTreeSum = segtree<ll, segsum, ESUM>;
// // 遅延セグ木パターン
// ll MPSUM(ll x, ll y) { return x + y; } // 区間加算用
// ll COMPSUM(ll x, ll y) { return x + y; } // 区間加算用
// ll IDSUM() { return 0LL; }
// // 最小値取得、区間加算(区間オフセット)
// using LazySegTreeGetMINChangeAdd = lazy_segtree<ll, segmin, EMAX, ll, MPSUM, COMPSUM, IDSUM>;
// // 最大値取得、区間加算(区間オフセット)
// using LazySegTreeGetMAXChangeAdd = lazy_segtree<ll, segmax, EMIN, ll, MPSUM, COMPSUM, IDSUM>;
// struct RangeSum { ll val; ll size; };
// vector<RangeSum> setTinit(vll &v) {
// vector<RangeSum> r;
// rep(i, len(v)) r.pb(RangeSum{v[i], 1});
// return r;
// }
// RangeSum MPRANGESUM(RangeSum a, RangeSum b) { return {a.val + b.val, a.size + b.size};}
// RangeSum RANGESUME0() { return {0, 0}; }
// RangeSum MPADD(ll f, RangeSum x) { return {x.val + f * x.size, x.size}; }
// ll RANGEID() { return 0; }
// // 区間和取得、区間加算(区間オフセット)
// using LazySegTreeGetSUMChangeAdd = lazy_segtree<RangeSum, MPRANGESUM, RANGESUME0, ll, MPADD, COMPSUM, RANGEID>;
// // 区間和取得、区間affine変換
// const ll ANSMOD = 998244353;
// // const ll ANSMOD = 1000000007;
// struct RangeSum { ll val; ll size; };
// RangeSum MPRANGESUM(RangeSum a, RangeSum b) { return {(a.val + b.val) % ANSMOD, a.size + b.size};}
// RangeSum RANGESUME0() { return {0, 0}; }
// struct AffineParam { ll b; ll c; };
// RangeSum MPADD(AffineParam f, RangeSum x) { return {((f.b * x.val) % ANSMOD + (f.c * x.size) % ANSMOD) % ANSMOD, x.size}; }
// AffineParam COMPSUM(AffineParam x, AffineParam y) { return AffineParam{(x.b * y.b) % ANSMOD, (y.c * x.b + x.c) % ANSMOD}; } // 区間加算用
// AffineParam RANGEID() { return {1, 0}; }
// using LazySegTreeGetSUMChangeAffine = lazy_segtree<RangeSum, MPRANGESUM, RANGESUME0, AffineParam, MPADD, COMPSUM, RANGEID>;
// const ll SEGID = LLONG_MAX - 14235;
// ll MPRANGEUPDATE(ll f, ll x) { return (f == SEGID) ? x : f; }
// RangeSum MPUPDATEADD(ll f, RangeSum x) {
// if (f != SEGID) x.val = f * x.size;
// return x;
// }
// ll COMPUPDATESUM(ll f, ll g) { return (f == SEGID ? g : f); }
// ll INITID() { return SEGID; }
// // 最小値取得、区間更新
// using LazySegTreeGetMINRangeUpdate = lazy_segtree<ll, segmin, EMAX, ll, MPRANGEUPDATE, COMPUPDATESUM, INITID>;
// // 最大値取得、区間更新
// using LazySegTreeGetMAXRangeUpdate = lazy_segtree<ll, segmax, EMIN, ll, MPRANGEUPDATE, COMPUPDATESUM, INITID>;
// // 区間和取得、区間更新
// using LazySegTreeGetSUMRangeUpdate = lazy_segtree<RangeSum, MPRANGESUM, RANGESUME0, ll, MPUPDATEADD, COMPUPDATESUM, INITID>;
// 複数要素を保持はできない
// 複数保持するときは、第二要素にインデックスを持たせる。
template <typename T>
class BinTree {
__gnu_pbds::tree<T, __gnu_pbds::null_type, less<T>, __gnu_pbds::rb_tree_tag, __gnu_pbds::tree_order_statistics_node_update> s;
public:
BinTree() : s() {}
BinTree(const vector<T> &v) { rep(i, len(v)) s.insert(v[i]);}
void add(const T &v) {
s.insert(v);
}
void del(const T &v) {
s.erase(v);
}
T getval(ll idx) {
return *s.find_by_order(idx);
}
T beginval() {
return *s.begin();
}
T lastval() {
return *s.rbegin();
}
ll getidx(const T &val) {
return s.order_of_key(val); // val未満の数が何個あるかわかる。(lower_bound)
}
void modify(const T &prevval, const T &newval) {
del(prevval);
add(newval);
}
};
// https://snuke.hatenablog.com/entry/2016/07/01/000000
// Mo's Algorithm・・・クエリ先読み可能で、区間の値を求めるときに使う。
// 区間クエリ[l, r]両側閉区間、1-originのクエリを設定
// 区間を伸ばしたときに答えを計算するAddFunc(ll idx)
// 区間を縮めたときに答えを計算するDelFunc(ll idx)
// 答えを出力するAnsFunc()
// のラムダ関数を設定する。
template<typename AddFunc, typename DelFunc, typename AnsFunc>
class MosAlgorithm {
int N_;
int Q_;
vector<int> l_, r_, qidx_;
AddFunc addfunc_;
DelFunc delfunc_;
AnsFunc ansfunc_;
ll maxn_;
inline ll hilbertorder(ll x, ll y) {
ll rx, ry, d = 0;
for (ll s = maxn_ >> 1; s; s >>= 1) {
rx = (x & s) > 0, ry = (y & s) > 0;
d += s * s * ((rx * 3) ^ ry);
if (ry) continue;
if (rx) {
x = maxn_ - 1 - x;
y = maxn_ - 1 - y;
}
swap(x, y);
}
return d;
}
public:
vll anslist_;
// 入力される区間クエリ[l, r]は1-originとする。
MosAlgorithm(int N, const vvll &offline_query, AddFunc addf, DelFunc delf, AnsFunc ansf)
: N_(N), Q_(len(offline_query)), addfunc_(addf), delfunc_(delf), ansfunc_(ansf) {
ll n = N_, cnt = 0;
while(n) { n >>= 1, cnt++;}
maxn_ = _1 << cnt;
anslist_.resize(Q_);
rep(i, Q_) {
l_.pb(offline_query[i][0] - 1); // 0-originにする。
r_.pb(offline_query[i][1]); // 0-originにして、境界をr]からr)に片側開区間にする。
qidx_.pb(i);
}
vll eval(Q_);
rep (i, Q_) eval[i] = hilbertorder(l_[i], r_[i]);
sort(all(qidx_),[&](int i, int j) { return eval[i] < eval[j]; });
debug(qidx_);
}
void solve() {
int nl = 0, nr = 0; // [nl, nr)
for (int i : qidx_) {
while (nl > l_[i]) --nl, addfunc_(nl);
while (nr < r_[i]) addfunc_(nr), ++nr;
while (nl < l_[i]) delfunc_(nl), ++nl;
while (nr > r_[i]) --nr, delfunc_(nr);
anslist_[i] = ansfunc_();
}
}
};
template <typename T>
T _Pivot(vector<T> &array, ll start, ll end);
template <typename T>
ll _Partition(vector<T> &array, ll start, ll end, const T &pivot)
{
vector<T> lt, eq, mt;
reps(i, start, end)
{
if (array[i] < pivot)
{
lt.pb(array[i]);
}
else if (array[i] == pivot)
{
eq.pb(array[i]);
}
else
{
mt.pb(array[i]);
}
}
reps(i, start, start + len(lt))
{
array[i] = lt[i - start];
}
start += len(lt);
ll ret = start;
reps(i, start, start + len(eq))
{
array[i] = eq[i - start];
}
start += len(eq);
reps(i, start, start + len(mt))
{
array[i] = mt[i - start];
}
return ret;
}
// [start, end)範囲のk番目(0-index)の値を取得。arrayは破壊的
template <typename T>
T get_kth_value(vector<T> &array, ll k, ll start, ll end)
{
if (end - start <= 1)
return array[start];
ll pivotIndex = -1;
do
{
T pivot = _Pivot(array, start, end);
pivotIndex = _Partition(array, start, end, pivot);
if (pivotIndex < k)
{
start = pivotIndex + 1;
}
else if (pivotIndex > k)
{
end = pivotIndex;
}
} while (pivotIndex != k);
return array[k + start - 1];
}
template <typename T>
T _Median5(vector<T> array, ll start, ll end)
{
reps(i, start, end)
{
reps(j, i + 1, end)
{
if (array[i] > array[j])
swap(array[i], array[j]);
}
}
return array[(end + start) / 2];
}
template <typename T>
T _Pivot(vector<T> &array, ll start, ll end)
{
vector<T> medians;
for (ll i = start; i < end; i += 5)
{
ll subStart = i;
ll subEnd = min(i + 5, end);
medians.pb(_Median5(array, subStart, subEnd));
}
ll n = len(medians);
ll st = 0;
ll ed = n;
ll newk = n / 2;
return get_kth_value(medians, newk, st, ed);
}
// vecの要素が整数(負の数から0含む正の数)の配列にたいして、
// k個選択したときの最大値を求める。
// O(N log N) (vecのsortが最悪計算量)
ll productMaxvalSelectK(const vll &vec, ll K, ll ANSMOD = 1000000007) {
vll A = vec;
const ll N = len(A);
SORT(A);
ll ans = 1;
if (N == K) {
rep(i, N) {
ans *= A[i];
ans %= ANSMOD;
}
return ans;
}
vll minus, plus, zeros;
rep(i, N) {
if (A[i] < 0) minus.pb(A[i]);
else if (A[i] == 0) zeros.pb(A[i]);
else plus.pb(A[i]);
}
REV(plus);
if (len(plus) + len(minus) < K) return _0;
if (len(plus) == 0) {
if (K % 2 == 0) {
rep(i, K) {
ans *= minus[i];
ans %= ANSMOD;
}
return ans;
} else {
if (len(zeros) > 0) return _0;
else {
REV(minus);
rep(i, K) {
ans *= minus[i];
ans %= ANSMOD;
}
return ans;
}
}
}
if (len(plus) + len(minus) == K) {
if (len(minus) % 2 == 1) return _0;
}
if (len(plus) < K) {
vll tmp = plus;
if (((len(tmp) % 2) + (K % 2)) % 2 != 0) tmp.pop_back();
vll minustmp;
rep(i, K - len(tmp)) {
minustmp.pb(minus[i]);
}
ll idx = K - len(tmp);
while(len(tmp) > 1 && len(minus) > idx + 1) {
if (minus[idx] * minus[idx + 1] > tmp[len(tmp) - 1] * tmp[len(tmp) - 2]) {
tmp.pop_back();
tmp.pop_back();
minustmp.pb(minus[idx]);
minustmp.pb(minus[idx + 1]);
idx += 2;
} else break;
}
rep(i, len(tmp)) {
ans *= tmp[i];
ans %= ANSMOD;
}
rep(i, len(minustmp)) {
ans *= minustmp[i];
ans %= ANSMOD;
}
return ans;
}
if (len(plus) >= K) {
vll tmp;
rep(i, K) {
tmp.pb(plus[i]);
}
ll idx = 0;
vll minustmp;
while(len(tmp) > 1 && len(minus) > idx + 1) {
if (minus[idx] * minus[idx + 1] > tmp[len(tmp) - 1] * tmp[len(tmp) - 2]) {
tmp.pop_back();
tmp.pop_back();
minustmp.pb(minus[idx]);
minustmp.pb(minus[idx + 1]);
idx += 2;
} else break;
}
rep(i, len(tmp)) {
ans *= tmp[i];
ans %= ANSMOD;
}
rep(i, len(minustmp)) {
ans *= minustmp[i];
ans %= ANSMOD;
}
return ans;
}
assert(false);
return _0;
}
enum rotOP {
ROTATE_90, // 反時計周り回転
ROTATE_180,
ROTATE_270,
MIRROR_LR,
MIRROR_UD,
TRANSPOSE
};
template<typename T>
vvll rotateMat(const vector<vector<T>> &mat, rotOP op) {
if (len(mat) == 0) return vvll();
ll H = len(mat), W = len(mat[0]);
if (op == ROTATE_90) {
vector<vector<T>> ret(W, vector<T>(H));
rep (i, H) {
rep (j, W) {
ret[W - 1 - j][i] = mat[i][j];
}
}
return ret;
} else if (op == ROTATE_180) {
vector<vector<T>> ret(H, vector<T>(W));
rep (i, H) {
rep (j, W) {
ret[H - 1 - i][W - 1 - j] = mat[i][j];
}
}
return ret;
} else if (op == ROTATE_270) {
vector<vector<T>> ret(W, vector<T>(H));
rep (i, H) {
rep (j, W) {
ret[j][H - 1 - i] = mat[i][j];
}
}
return ret;
} else if (op == MIRROR_LR) {
vector<vector<T>> ret(H, vector<T>(W));
rep (i, H) {
rep (j, W) {
ret[i][W - 1 - j] = mat[i][j];
}
}
return ret;
} else if (op == MIRROR_UD) {
vector<vector<T>> ret(H, vector<T>(W));
rep (i, H) {
rep (j, W) {
ret[H - 1 - i][j] = mat[i][j];
}
}
return ret;
} else if (op == TRANSPOSE) {
vector<vector<T>> ret(W, vector<T>(H));
rep (i, H) {
rep (j, W) {
ret[j][i] = mat[i][j];
}
}
return ret;
}
return mat;
}
// 木構造を用いて履歴管理。ABC273-E
template <class T>
class Tree {
public: // 競プロなので利便性優先
unordered_map<ll, T> elems_;
umpll history_;
umpll parentidx_;
ll curidx_; // 0番をrootとしてrootの要素は-1とする。
ll addidx_;
const ll ROOT_ = -1;
Tree(const T &rootelem) : elems_(), history_(), parentidx_(), curidx_(-1), addidx_(0) {
elems_[curidx_] = rootelem;
parentidx_[curidx_] = ROOT_;
}
void add(const T &elem) { // elemを作成し、現在位置からそこに移動。
parentidx_[addidx_] = curidx_;
curidx_= addidx_;
elems_[curidx_] = elem;
addidx_++;
}
void move(ll idx) { // 指定した頂点に移動
curidx_ = idx;
}
void prev() { // 親に移動
curidx_ = parentidx_[curidx_];
}
T nowgetelem() { // 現在位置の要素を取得
return elems_[curidx_];
}
// historyidxは頂点番号とは別。履歴番号historyidxに現在位置を登録
void save(ll historyidx) {
history_[historyidx] = curidx_;
}
// historyidxは頂点番号とは別。
// 履歴番号historyidxに登録されている位置に現在位置を移動。
// 未登録のidxのときはrootに移動。
void load(ll historyidx) {
if (EXIST(historyidx, history_)) curidx_ = history_[historyidx];
else curidx_ = ROOT_;
}
};
// こんな感じで現在の位置から次の遷移の値を返すラムダ関数を突っ込む。
// auto D = DoublingQuery(N, [&](ll x) { return A[x]; });
// 後はD.getval(0, K);でOK。0の位置からK回遷移した後の位置を取得できる。
template<typename Func>
class DoublingQuery {
vvll tbl_;
const ll n_; // データ数
const ll CNT_ = 100; // 2^100回遷移後まで求めれる。
Func func_;
// 1回処理する遷移関数群(この関数は問題によって書き換え。)
// どの値がどの値に遷移するのかを実装する。
// ll transition_func(ll x) { // 整数xの各桁和 + 元の数の下5桁を返す。
// return (x + digitsum(STR(x))) % 100000;
// }
// 1回遷移後, 2回遷移後, 4回遷移後, 8回遷移後...のデータをtbl_に構築
void build() {
tbl_.clear();
vll tmp(n_);
rep (i, n_) {
tmp[i] = func_(i);
}
tbl_.pb(tmp);
rep (i, CNT_) {
vll newtbl(n_);
rep (j, n_) {
newtbl[j] = tbl_[i][tbl_[i][j]];
}
tbl_.pb(newtbl);
}
}
public:
DoublingQuery(ll N, Func func) // Nには遷移後に可能性のある値の最大値を設定
: tbl_()
, n_(N)
, func_(func) {
build();
}
ll getval(ll pos, ll cnt) {
ll p = pos;
ll num = 0;
while (cnt) {
if (cnt & _1) {
p = tbl_[num][p];
}
cnt >>= 1;
num++;
}
return p;
}
};
// こんな感じで現在の位置から次の遷移の値を返すラムダ関数を突っ込む。
// auto D = CyclePos(N, [&](ll x) { return A[x]; });
// 後はD.getval(0, K);でOK。0の位置からK回遷移した後の位置を取得できる。
template<typename Func>
class CyclePos {
const ll n_; // データ数
Func func_;
public:
CyclePos(ll N, Func func) // Nには遷移後に可能性のある値の最大値を設定
: n_(N)
, func_(func) {
}
// posは0-origin, cntは超巨大な10進数の文字列。10^100000オーダー。
// 計算量はデータ数Nで遷移回数10^Mのとき、大体O(N + M)
ll getval(ll pos, const string &cnt) {
vll checkin(n_, -1);
ll p = pos;
checkin[p] = 0;
ll cyclecnt = 0, notcyclecnt = 0;
ll cyclestartpos = -1;
rep(i, n_) {
ll nextp = func_(p);
if (checkin[nextp] != -1) {
notcyclecnt = checkin[nextp];
cyclecnt = i + 1 - notcyclecnt;
cyclestartpos = nextp;
break;
} else {
checkin[nextp] = i + 1;
}
p = nextp;
}
ll val = 0;
rep(i, len(cnt)) {
val = 10 * val + CHARLL(cnt[i]);
if (cyclecnt != 0) {
if (val >= notcyclecnt + cyclecnt) {
val = (val - notcyclecnt) % cyclecnt + notcyclecnt;
}
}
}
ll ret = pos;
if (val >= notcyclecnt) {
ret = cyclestartpos;
val -= notcyclecnt;
}
rep(i, val) {
ret = func_(ret);
}
return ret;
}
};
// insertした値の中から、大きいほう、または小さい方のK個の総和を求める。
template<typename T, typename Compare = less<T>, typename RCompare = greater<T>>
class PrioritySumStructure {
private:
ll k_;
T sumval_;
priority_queue<T, vector<T>, Compare> in, d_in;
priority_queue<T, vector<T>, RCompare> out, d_out;
private:
void modify() {
while(len(in) - len(d_in) < k_ && !out.empty()) {
auto p = out.top();
out.pop();
if(!d_out.empty() && p == d_out.top()) {
d_out.pop();
} else {
sumval_ += p;
in.emplace(p);
}
}
while(len(in) - len(d_in) > k_) {
auto p = in.top();
in.pop();
if(!d_in.empty() && p == d_in.top()) {
d_in.pop();
} else {
sumval_ -= p;
out.emplace(p);
}
}
while(!d_in.empty() && in.top() == d_in.top()) {
in.pop();
d_in.pop();
}
}
public:
PrioritySumStructure(ll k) : k_(k), sumval_(0) {}
T getval() const {
return sumval_;
}
void insert(T x) {
in.emplace(x);
sumval_ += x;
modify();
}
void erase(T x) {
assert(size());
if(!in.empty() && in.top() == x) {
sumval_ -= x;
in.pop();
} else if(!in.empty() && RCompare()(in.top(), x)) {
sumval_ -= x;
d_in.emplace(x);
} else {
d_out.emplace(x);
}
modify();
}
void set_k(ll kk) {
k_ = kk;
modify();
}
ll get_k() const {
return k_;
}
ll size() const {
return len(in) + out.size() - len(d_in) - d_out.size();
}
};
template< typename T >
using MaximumSum = PrioritySumStructure<T, greater<T>, less<T>>;
// auto maxsum = MaximumSum<ll>(N); // 大きい方N個の合計求める
template< typename T >
using MinimumSum = PrioritySumStructure<T, less<T>, greater<T>>;
// auto minsum = MinimumSum<ll>(N); // 小さい方N個の合計求める
vvll getdir4() {
return {{-1, 0}, {1, 0}, {0, -1}, {0, 1}};
}
vvll getdir8() {
return {{-1, 0}, {1, 0}, {0, -1}, {0, 1}, {-1, 1}, {1, 1}, {1, -1}, {-1, -1}};
}
// 2次元配列の十字方向探索の基本コード
// rep (i, H) {
// rep(j, W) {
// vvll d = getdir4();
// rep(dd, len(d)) {
// INI2(dh, dw, d[dd]);
// ll h = i + dh, w = j + dw;
// if (IN(0, h, H - 1) && IN(0, w, W - 1)) {
// }
// }
// }
// }
template<class T>
vvll manhattan_mst_(const vector<pair<T, T>> &rps, T inf = numeric_limits<T>::max()) {
vector<pair<T, T>> ps = rps;
vector<pair<ll, ll>> edges;
ll n = ll(ps.size());
SETPERM(ids, n);
rep (ph, 4) {
stable_sort(all(ids), [&](auto i, auto j) {
T ixy = (ps[i].first + ps[i].second), jxy = (ps[j].first + ps[j].second);
return tie(ixy, ps[i].second) > tie(jxy, ps[j].second);
});
vector<T> xv;
rep (i, n) xv.push_back(ps[i].first);
stable_sort(all(xv));
xv.erase(unique(all(xv)), xv.end());
using P = pair<T, ll>;
vector<P> fen(n, P(-inf, -1));
for (auto id : ids) {
auto xi = int(lower_bound(all(xv), ps[id].first) - xv.begin());
P ma = P(-inf, -1);
{
ll i = xi + 1;
while (i > 0) {
if (ma.first <= fen[i - 1].first) ma = fen[i - 1];
i -= i & -i;
}
}
if (ma.second != -1) edges.push_back({id, ma.second});
{
T x = ps[id].first - ps[id].second;
ll i = xi + 1;
while (i <= n) {
if (fen[i - 1].first <= x) fen[i - 1] = P(x, id);
i += i & -i;
}
}
}
for (auto &p : ps) {
swap(p.first, p.second);
}
if (ph == 1) {
for (auto &p : ps) {
p.second *= -1;
}
}
}
auto dist = [&](ll i, ll j) {
return abs(ps[i].first - ps[j].first) + abs(ps[i].second - ps[j].second);
};
stable_sort(all(edges), [&](auto x, auto y) {
return dist(x.first, x.second) < dist(y.first, y.second);
});
auto U = UnionFind(n);
vvll res;
for (auto p : edges) {
if (U.same(p.first + 1, p.second + 1)) continue;
res.pb(vll{p.first, p.second});
U.unite(p.first + 1, p.second + 1);
}
return res;
}
// マンハッタンMST O(N log N)
// 木の頂点番号(0-origin)の組の配列を返す
template<class T>
vvll manhattan_mst(const vector<vector<T>> &xy) {
vector<pair<T, T>> ps;
rep (i, len(xy)) {
T x = xy[i][0], y = xy[i][1];
auto p = mp(x, y);
ps.pb(p);
}
return manhattan_mst_(ps);
}
template <typename T>
class Matrix {
public:
// mod = 0のとき、余り計算しない。
static vector<vector<T>> multi(const vector<vector<T>> &A, const vector<vector<T>> &B) {
ll Ah = len(A), Aw = len(A[0]), Bw = len(B[0]);
vector<vector<T>> C(Ah, vector<T>(Bw));
rep (i, Ah)
rep (j, Bw)
rep (k, Aw) {
C[i][j] += A[i][k] * B[k][j];
}
return C;
}
static vector<vector<T>> init_matrix(ll N) {
vector<vector<T>> C(N, vector<T>(N));
rep (i, N) C[i][i] = (T)1;
return C;
}
// 行列の累乗。行列の縦横が少ないときに高速に求めれる。
// O(log n)
// mod = 0のとき、余り計算しない。
static vector<vector<T>> power(const vector<vector<T>> &A, ll r) {
if (r == 0) {
return init_matrix(len(A));
} else if (r % 2 == 0) {
return power(multi(A, A), r >> 1);
} else {
return multi(A, power(A, r - 1));
}
}
};
// 2択で片方しか選択出来ない場合の最大種類数
ll getMaximumCntType(const vvll &cards) {
ll N = len(cards);
vll card;
rep (i, N) {
INI2(A, B, cards[i]);
card.pb(A);
card.pb(B);
}
auto [d, revd] = compcoord(card);
auto U = UnionFind(len(d));
vector<ll> treechecks(len(d) + 1, 1);
rep (i, N) {
INI2(A, B, cards[i]);
ll u = d[A], v = d[B];
u++, v++;
if (U.same(u, v)) treechecks[U.find(u)] = 0;
else {
bool istree = treechecks[U.find(u)] && treechecks[U.find(v)];
U.unite(u, v);
treechecks[U.find(u)] &= istree;
}
}
vll roots = U.roots();
ll ans = 0;
rep (i, len(roots)) {
ans += treechecks[roots[i]] ? U.size(roots[i]) - 1 : U.size(roots[i]);
}
return ans;
}
// https://github.com/xuzijian629/library2/blob/master/treap/implicit_treap.cpp 参照
// T0: 元の配列のモノイド
// T1: T0に対する作用素モノイド
template<class T0, class T1>
class BaseImplicitTreap {
// T0上の演算、単位元
virtual T0 f0(T0, T0) = 0;
const T0 u0;
// T1上の演算、単位元
virtual T1 f1(T1, T1) = 0;
const T1 u1;
// T0に対するT1の作用
virtual T0 g(T0, T1) = 0;
// 多数のt1(T1)に対するf1の合成
virtual T1 p(T1, ll) = 0;
class xorshift {
uint64_t x;
public:
xorshift() {
mt19937 rnd(chrono::steady_clock::now().time_since_epoch().count());
x = rnd();
for (ll i = 0; i < 100; i++) {
random();
}
}
uint64_t random() {
x = x ^ (x << 7);
return x = x ^ (x >> 9);
}
} rnd;
struct Node {
T0 value, acc;
T1 lazy;
ll priority, cnt;
bool rev;
Node *l, *r;
Node(T0 value_, ll priority_, T0 u0_, T1 u1_)
: value(value_)
, acc(u0_)
, lazy(u1_)
, priority(priority_)
, cnt(1)
, rev(false)
, l(nullptr)
, r(nullptr) {}
} *root = nullptr;
using Tree = Node *;
ll cnt(Tree t) { return t ? t->cnt : 0; }
T0 acc(Tree t) { return t ? t->acc : u0; }
void update_cnt(Tree t) {
if (t) {
t->cnt = 1 + cnt(t->l) + cnt(t->r);
}
}
void update_acc(Tree t) {
if (t) {
t->acc = f0(acc(t->l), f0(t->value, acc(t->r)));
}
}
void pushup(Tree t) { update_cnt(t), update_acc(t); }
void pushdown(Tree t) {
if (t && t->rev) {
t->rev = false;
swap(t->l, t->r);
if (t->l) t->l->rev ^= 1;
if (t->r) t->r->rev ^= 1;
}
if (t && t->lazy != u1) {
if (t->l) {
t->l->lazy = f1(t->l->lazy, t->lazy);
t->l->acc = g(t->l->acc, p(t->lazy, cnt(t->l)));
}
if (t->r) {
t->r->lazy = f1(t->r->lazy, t->lazy);
t->r->acc = g(t->r->acc, p(t->lazy, cnt(t->r)));
}
t->value = g(t->value, p(t->lazy, 1));
t->lazy = u1;
}
pushup(t);
}
void split(Tree t, ll key, Tree &l, Tree &r) {
if (!t) {
l = r = nullptr;
return;
}
pushdown(t);
ll implicit_key = cnt(t->l) + 1;
if (key < implicit_key) {
split(t->l, key, l, t->l), r = t;
} else {
split(t->r, key - implicit_key, t->r, r), l = t;
}
pushup(t);
}
void insert(Tree &t, ll key, Tree item) {
Tree t1, t2;
split(t, key, t1, t2);
merge(t1, t1, item);
merge(t, t1, t2);
}
void merge(Tree &t, Tree l, Tree r) {
pushdown(l);
pushdown(r);
if (!l || !r) {
t = l ? l : r;
} else if (l->priority > r->priority) {
merge(l->r, l->r, r), t = l;
} else {
merge(r->l, l, r->l), t = r;
}
pushup(t);
}
void erase(Tree &t, ll key) {
Tree t1, t2, t3;
split(t, key + 1, t1, t2);
split(t1, key, t1, t3);
merge(t, t1, t2);
}
void update(Tree t, ll l, ll r, T1 x) {
if (l >= r) return;
Tree t1, t2, t3;
split(t, l, t1, t2);
split(t2, r - l, t2, t3);
t2->lazy = f1(t2->lazy, x);
t2->acc = g(t2->acc, p(x, cnt(t2)));
merge(t2, t2, t3);
merge(t, t1, t2);
}
T0 query(Tree t, ll l, ll r) {
if (l == r) return u0;
Tree t1, t2, t3;
split(t, l, t1, t2);
split(t2, r - l, t2, t3);
T0 ret = t2->acc;
merge(t2, t2, t3);
merge(t, t1, t2);
return ret;
}
// [l, r)の中で左から何番目か
ll find(Tree t, T0 x, ll offset, bool left = true) {
if (f0(t->acc, x) == x) {
return -1;
} else {
if (left) {
if (t->l && f0(t->l->acc, x) != x) {
return find(t->l, x, offset, left);
} else {
return (f0(t->value, x) != x) ? offset + cnt(t->l) : find(t->r, x, offset + cnt(t->l) + 1, left);
}
} else {
if (t->r && f0(t->r->acc, x) != x) {
return find(t->r, x, offset + cnt(t->l) + 1, left);
} else {
return (f0(t->value, x) != x) ? offset + cnt(t->l) : find(t->l, x, offset, left);
}
}
}
}
void reverse(Tree t, ll l, ll r) {
if (l > r) return;
Tree t1, t2, t3;
split(t, l, t1, t2);
split(t2, r - l, t2, t3);
t2->rev ^= 1;
merge(t2, t2, t3);
merge(t, t1, t2);
}
// [l, r)の先頭がmになるようにシフトさせる。std::rotateと同じ仕様
void rotate(Tree t, ll l, ll m, ll r) {
reverse(t, l, r);
reverse(t, l, l + r - m);
reverse(t, l + r - m, r);
}
void dump(Tree t) {
if (!t) return;
pushdown(t);
dump(t->l);
cout << t->value << " ";
dump(t->r);
}
public:
BaseImplicitTreap(T0 u0_, T1 u1_)
: u0(u0_)
, u1(u1_) {}
void set_by_vector(const vector<T0> &a) {
for (ll i = 0; i < a.size(); i++) {
insert(i, a[i]);
}
}
ll size() { return cnt(root); }
void insert(ll pos, T0 x) { insert(root, pos, new Node(x, rnd.random(), u0, u1)); }
void update(ll l, ll r, T1 x) { update(root, l, r, x); }
T0 query(ll l, ll r) { return query(root, l, r); }
// 二分探索。[l, r)内のkでf0(tr[k], x) != xとなる最左/最右のもの。存在しない場合は-1
// たとえばMinMonoidの場合、x未満の最左/最右の要素の位置を返す
ll binary_search(ll l, ll r, T0 x, bool left = true) {
if (l >= r) return -1;
Tree t1, t2, t3;
split(root, l, t1, t2);
split(t2, r - l, t2, t3);
ll ret = find(t2, x, l, left);
merge(t2, t2, t3);
merge(root, t1, t2);
return ret;
}
void erase(ll pos) { erase(root, pos); }
void reverse(ll l, ll r) { reverse(root, l, r); }
void rotate(ll l, ll m, ll r) { rotate(root, l, m, r); }
void dump() {
dump(root);
cout << endl;
}
T0 operator[](ll pos) { return query(pos, pos + 1); }
};
template<class T0, class T1>
struct MinUpdateQuery : public BaseImplicitTreap<T0, T1> {
using BaseImplicitTreap<T0, T1>::BaseImplicitTreap;
MinUpdateQuery()
: MinUpdateQuery(numeric_limits<T0>::max(), numeric_limits<T1>::min()) {}
T0 f0(T0 x, T0 y) override { return min(x, y); }
T1 f1(T1 x, T1 y) override { return y == numeric_limits<T1>::min() ? x : y; }
T0 g(T0 x, T1 y) override { return y == numeric_limits<T1>::min() ? x : y; }
T1 p(T1 x, ll L) override { return x; }
};
template<class T0, class T1>
struct SumAddQuery : public BaseImplicitTreap<T0, T1> {
using BaseImplicitTreap<T0, T1>::BaseImplicitTreap;
SumAddQuery()
: SumAddQuery(0, 0) {}
T0 f0(T0 x, T0 y) override { return x + y; }
T1 f1(T1 x, T1 y) override { return x + y; }
T0 g(T0 x, T1 y) override { return x + y; }
T1 p(T1 x, ll L) override { return x * L; }
};
template<class T0, class T1>
struct MinAddQuery : public BaseImplicitTreap<T0, T1> {
using BaseImplicitTreap<T0, T1>::BaseImplicitTreap;
MinAddQuery()
: MinAddQuery(numeric_limits<T0>::max(), 0) {}
T0 f0(T0 x, T0 y) override { return min(x, y); }
T1 f1(T1 x, T1 y) override { return x + y; }
T0 g(T0 x, T1 y) override { return x + y; }
T1 p(T1 x, ll L) override { return x; }
};
template<class T0, class T1>
struct SumUpdateQuery : public BaseImplicitTreap<T0, T1> {
using BaseImplicitTreap<T0, T1>::BaseImplicitTreap;
SumUpdateQuery()
: SumUpdateQuery(0, numeric_limits<T1>::min()) {}
T0 f0(T0 x, T0 y) override { return x + y; }
T1 f1(T1 x, T1 y) override { return y == numeric_limits<T1>::min() ? x : y; }
T0 g(T0 x, T1 y) override { return y == numeric_limits<T1>::min() ? x : y; }
T1 p(T1 x, ll L) override { return x == numeric_limits<T1>::min() ? numeric_limits<T1>::min() : x * L; }
};
template<class T0>
struct SumAffineQuery : public BaseImplicitTreap<T0, pair<T0, T0>> {
using T1 = pair<T0, T0>; // first * x + second
using BaseImplicitTreap<T0, T1>::BaseImplicitTreap;
SumAffineQuery()
: SumAffineQuery(0, {1, 0}) {}
T0 f0(T0 x, T0 y) override { return x + y; }
T1 f1(T1 x, T1 y) override { return {x.first * y.first, x.second * y.first + y.second}; }
T0 g(T0 x, T1 y) override { return y.first * x + y.second; }
T1 p(T1 x, ll L) override { return {x.first, x.second * L}; }
// update(i, j, {a, b}); // [i, j)にax + bを作用
// update(i, j, {0, a}); // update
// update(i, j, {1, a}); // 加算
// update(i, j, {a, 0}); // 倍
};
// 余りを使用するパターンを追加。
template<class T0>
struct SumModAffineQuery : public BaseImplicitTreap<T0, pair<T0, T0>> {
ll ANSMOD;
using T1 = pair<T0, T0>; // first * x + second
using BaseImplicitTreap<T0, T1>::BaseImplicitTreap;
SumModAffineQuery()
: SumModAffineQuery(0, {1, 0}) {}
void set_mod(ll ansmod) {ANSMOD = ansmod;}
T0 f0(T0 x, T0 y) override { return mod_add(x, y, ANSMOD); }
T1 f1(T1 x, T1 y) override { return {mod_mul(x.first, y.first, ANSMOD), mod_add(mod_mul(x.second, y.first, ANSMOD), y.second, ANSMOD)}; }
T0 g(T0 x, T1 y) override { return mod_add(mod_mul(y.first, x, ANSMOD), y.second, ANSMOD); }
T1 p(T1 x, ll L) override { return {x.first, mod_mul(x.second, L, ANSMOD)}; }
// update(i, j, {a, b}); // [i, j)にax + bを作用
// update(i, j, {0, a}); // update
// update(i, j, {1, a}); // 加算
// update(i, j, {a, 0}); // 倍
};
template<class T>
struct MinmaxAffineQuery : public BaseImplicitTreap<pair<T, T>, pair<T, T>> {
using T0 = pair<T, T>; // {min, max}
using T1 = pair<T, T>; // first * x + second
using BaseImplicitTreap<T0, T1>::BaseImplicitTreap;
MinmaxAffineQuery()
: MinmaxAffineQuery({numeric_limits<T>::max(), -numeric_limits<T>::max()}, {1, 0}) {
} // TODO: _u1を使うとコンパイル通らない原因不明
T0 f0(T0 x, T0 y) override { return {min(x.first, y.first), max(x.second, y.second)}; }
T1 f1(T1 x, T1 y) override { return {x.first * y.first, x.second * y.first + y.second}; }
T0 g(T0 x, T1 y) override {
T0 ret = {x.first * y.first + y.second, x.second * y.first + y.second};
if (y.first < 0) swap(ret.first, ret.second);
return ret;
}
T1 p(T1 x, ll L) override { return x; }
// update(i, j, {a, b}); // [i, j)にax + bを作用
// update(i, j, {0, a}); // update
// update(i, j, {1, a}); // 加算
// update(i, j, {a, 0}); // 倍
};
// 年と月から存在する最大の日付を返す。
ll yearmonth2day(ll y, ll m)
{
if (m == 2) return y % 4 ? 28 : y % 100 ? 29 : y % 400 ? 28 : 29;
else if (m == 4 || m == 6 || m == 9 || m == 11) return 30;
else return 31;
}
// O(N log N)
// SORT(v)
// rep(i, N) rep(j, i + 1, N) total += v[j] - v[i];
// 降順に並べたときに異なる要素の差の合計を出す。
template<class T>
T sumOfDifference(const vector<T> &v) {
vector<T> A = v;
RSORT(A);
ll ans = 0;
ll N = len(v);
ll base = N - 1;
rep(i, N) {
ans += A[i] * base;
base -= 2;
}
return ans;
}
// 異なる2組(i < j)の座標について、マンハッタン距離の総和を求める。
// O(N log N)
template<class T>
T manhattanDistanceSum(const vector<vector<T>> &XY) {
vll Xs, Ys;
ll N = len(XY);
rep(i, N) {
INI(X, Y, XY[i]);
Xs.pb(X), Ys.pb(Y);
}
return sumOfDifference(Xs) + sumOfDifference(Ys);
}
#include __FILE__
#endif
Submission Info
Submission Time
2023-07-15 15:53:28+0900
Task
D - Maxflow
User
mind_cpp
Language
C++ (GCC 9.2.1)
Score
100
Code Size
234499 Byte
Status
AC
Exec Time
22 ms
Memory
10896 KiB
Compile Error
In file included from ./Main.cpp:7754:
././Main.cpp: In function ‘int main()’:
././Main.cpp:67:10: warning: unused variable ‘ans’ [-Wunused-variable]
67 | auto ans = solve();
| ^~~
Judge Result
Set Name
Sample
All
Score / Max Score
0 / 0
100 / 100
Status
Set Name
Test Cases
Sample
00-sample-01.txt
All
00-sample-01.txt, 01-01.txt, 01-02.txt, 01-03.txt, 01-04.txt, 01-05.txt, 01-06.txt, 01-07.txt, 01-08.txt, 01-09.txt, 01-10.txt, 01-11.txt, 01-12.txt, 01-13.txt, 01-14.txt, 01-15.txt, 01-16.txt, 01-17.txt, 01-18.txt
Case Name
Status
Exec Time
Memory
00-sample-01.txt
AC 6 ms
3936 KiB
01-01.txt
AC 2 ms
3756 KiB
01-02.txt
AC 2 ms
3768 KiB
01-03.txt
AC 2 ms
3948 KiB
01-04.txt
AC 6 ms
3836 KiB
01-05.txt
AC 4 ms
3872 KiB
01-06.txt
AC 2 ms
3752 KiB
01-07.txt
AC 15 ms
9744 KiB
01-08.txt
AC 3 ms
3964 KiB
01-09.txt
AC 3 ms
4044 KiB
01-10.txt
AC 4 ms
4040 KiB
01-11.txt
AC 3 ms
3852 KiB
01-12.txt
AC 22 ms
10896 KiB
01-13.txt
AC 18 ms
8632 KiB
01-14.txt
AC 13 ms
6192 KiB
01-15.txt
AC 6 ms
4684 KiB
01-16.txt
AC 4 ms
3980 KiB
01-17.txt
AC 16 ms
6148 KiB
01-18.txt
AC 12 ms
6272 KiB