#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define allof(obj) (obj).begin(), (obj).end() #define range(i, l, r) for(int i=l;i>1)|y_bit)) #define bit_kth(i, k) ((i >> k)&1) #define bit_highest(i) (i?63-__builtin_clzll(i):-1) #define bit_lowest(i) (i?__builtin_ctzll(i):-1) #define sleepms(t) std::this_thread::sleep_for(std::chrono::milliseconds(t)) using ll = long long; using ld = long double; using ul = uint64_t; using pi = std::pair; using pl = std::pair; using namespace std; template std::ostream &operator<<(std::ostream &dest, const std::pair &p){ dest << p.first << ' ' << p.second; return dest; } template std::ostream &operator<<(std::ostream &dest, const std::vector> &v){ int sz = v.size(); if(sz==0) return dest; for(int i=0;i std::ostream &operator<<(std::ostream &dest, const std::vector &v){ int sz = v.size(); if(sz==0) return dest; for(int i=0;i std::ostream &operator<<(std::ostream &dest, const std::array &v){ if(sz==0) return dest; for(int i=0;i std::ostream &operator<<(std::ostream &dest, const std::set &v){ for(auto itr=v.begin();itr!=v.end();){ dest << *itr; itr++; if(itr!=v.end()) dest << ' '; } return dest; } template std::ostream &operator<<(std::ostream &dest, const std::map &v){ for(auto itr=v.begin();itr!=v.end();){ dest << '(' << itr->first << ", " << itr->second << ')'; itr++; if(itr!=v.end()) dest << '\n'; } return dest; } std::ostream &operator<<(std::ostream &dest, __int128_t 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; } template vector make_vec(size_t sz, T val){return std::vector(sz, val);} template auto make_vec(size_t sz, Tail ...tail){ return std::vector(tail...))>(sz, make_vec(tail...)); } template vector read_vec(size_t sz){ std::vector v(sz); for(int i=0;i<(int)sz;i++) std::cin >> v[i]; return v; } template auto read_vec(size_t sz, Tail ...tail){ auto v = std::vector(tail...))>(sz); for(int i=0;i<(int)sz;i++) v[i] = read_vec(tail...); return v; } long long max(long long a, int b){return std::max(a, (long long)b);} long long max(int a, long long b){return std::max((long long)a, b);} long long min(long long a, int b){return std::min(a, (long long)b);} long long min(int a, long long b){return std::min((long long)a, b);} long long modulo(long long a, long long m){a %= m; return a < 0 ? a + m : a;} template struct safe_vector : std::vector{ using std::vector::vector; T& operator [](size_t i){return this->at(i);} }; template struct safe_array : std::array{ using std::array::array; T& operator [](size_t i){return this->at(i);} }; ll ceil_div(ll x, ll y){ assert(y > 0); return (x + (x > 0 ? y - 1 : 0)) / y; } ll floor_div(ll x, ll y){ assert(y > 0); return (x + (x > 0 ? 0 : -y + 1)) / y; } void io_init(){ std::cin.tie(nullptr); std::ios::sync_with_stdio(false); } #include struct barrett_reduction{ unsigned int mod; unsigned long long m; barrett_reduction(unsigned int _mod) : mod(_mod){ m = ((__uint128_t)1 << 64) / mod; } unsigned int reduce(unsigned int a){ unsigned long long q = ((__uint128_t)a * m) >> 64; a -= q * mod; // 0 <= a < 2 * mod // return a; return a >= mod ? a - mod : a; } unsigned int mul(unsigned int a, unsigned int b){ return reduce((unsigned long long)a * b); } // {gcd(a, mod), x}, such that a * x ≡ gcd(a, mod) std::pair inv(unsigned int a){ if(a >= mod) a = reduce(a); if(a == 0) return {mod, 0}; unsigned int s = mod, t = a; long long m0 = 0, m1 = 1; while(t){ int u = s / t; s -= t * u; m0 -= m1 * u; std::swap(m0, m1); std::swap(s, t); } if(m0 < 0) m0 += mod / s; return {s, m0}; } }; // 64bit mod対応 // R = 2^64 // 偶数modだと壊れる struct montgomery_reduction_64bit{ private: // [0, 2 * MOD) inline uint64_t reduce_unsafe(__uint128_t x) const{ x = (x + ((uint64_t)x * (uint64_t)NEG_INV) * MOD) >> 64; return x; } void _set_mod(uint64_t mod){ assert((mod < (1ULL << 63)) && (mod & 1)); MOD = mod; NEG_INV = 0; __uint128_t s = 1, t = 0; for(int i = 0; i < 64; i++){ if (~t & 1) { t += MOD; NEG_INV += s; } t >>= 1; s <<= 1; } R2 = ((__uint128_t)1 << 64) % MOD; R2 = R2 * R2 % MOD; ONE = generate(1); } __uint128_t MOD, NEG_INV, R2; uint64_t ONE; public: montgomery_reduction_64bit(){} montgomery_reduction_64bit(uint64_t mod){_set_mod(mod);} void set_mod(uint64_t mod){ _set_mod(mod); } uint64_t mod()const{ return MOD; } inline uint64_t one()const{ return ONE; } inline uint64_t generate(uint64_t x)const{ return reduce((__uint128_t)x * R2); } inline uint64_t reduce(__uint128_t x)const{ x = (x + ((uint64_t)x * (uint64_t)NEG_INV) * MOD) >> 64; return x < MOD ? x : x - MOD; } inline uint64_t fix(uint64_t x)const{ return x < MOD ? x : x - MOD; } // [0, 2MOD) inline uint64_t mul(uint64_t mx, uint64_t my)const{ return reduce_unsafe((__uint128_t)mx * my); } inline uint64_t mul_safe(uint64_t mx, uint64_t my)const{ return reduce((__uint128_t)mx * my); } // [0, 2MOD) inline uint64_t add(uint64_t mx, uint64_t my)const{ return (mx >= MOD ? mx - MOD : mx) + (my >= MOD ? my - MOD : my); } inline uint64_t add_safe(uint64_t mx, uint64_t my)const{ uint64_t res = (mx >= MOD ? mx - MOD : mx) + (my >= MOD ? my - MOD : my); return res >= MOD ? res - MOD : res; } // [0, 2MOD) inline uint64_t sub(uint64_t mx, uint64_t my)const{ if(my >= MOD) my -= MOD; return mx >= my ? mx - my : mx + MOD - my; } inline uint64_t sub_safe(uint64_t mx, uint64_t my)const{ if(my >= MOD) my -= MOD; uint64_t res = mx >= my ? mx - my : mx + MOD - my; return res >= MOD ? res - MOD : res; } inline uint64_t pow(uint64_t ma, uint64_t b)const{ uint64_t ret = one(); while(b){ if(b & 1) ret = mul(ret, ma); ma = mul(ma, ma); b >>= 1; } return ret; } inline uint64_t pow_safe(uint64_t ma, uint64_t b)const{ return fix(pow(ma, b)); } }; unsigned long long mod_pow_mr(unsigned long long a, unsigned long long b, unsigned long long m){ montgomery_reduction_64bit mr(m); return mr.reduce(mr.pow(mr.generate(a), b)); } namespace prime_sieve{ std::vector primes, min_factor;// 素数, 各数を割り切る最小の素数 // O(MAX_N loglog MAX_N) // [1, MAX_N]を扱えるように初期化 void init(int MAX_N){ min_factor.resize(MAX_N + 1, -1); for(int i = 2; i <= MAX_N; i++){ if(min_factor[i] == -1){ primes.push_back(i); min_factor[i] = i; } for(int p : primes){ if((long long)p * i > MAX_N || p > min_factor[i]) break; min_factor[p * i] = p; } } } bool is_prime(int n){ assert(n < min_factor.size()); return n == min_factor[n]; } // {{素因数, 数}}, O(log n) std::vector> factorize(int n){ assert(n < min_factor.size()); std::vector> ret; while(n > 1){ int cnt = 0, f = min_factor[n]; while(n % f == 0){ n /= f; cnt++; } ret.push_back({f, cnt}); } return ret; } // 約数列挙, O(√n) std::vector divisor(int n){ auto p = factorize(n); std::vector> x; for(int i = 0; i < p.size(); i++){ x.push_back(std::vector{1}); for(int j = 0; j < p[i].second; j++) x[i].push_back(x[i][j] * p[i].first); } int l = 0, r = 1; std::vector ret{1}; for(int i = 0; i < x.size(); i++){ for(auto e : x[i]){ for(int j = l; j < r; j++){ ret.push_back(ret[j] * e); } } l = r; r = ret.size(); } return std::vector(ret.begin() + l, ret.end()); } // O(logN) unsigned long long totient_function(unsigned long long n){ unsigned long long res = n; int prev = -1; while(n > 1){ if(min_factor[n] > prev){ res -= res / min_factor[n]; prev = min_factor[n]; } n /= min_factor[n]; } return res; } int mobius_function(int x){ int pcnt = 0; while(x > 1){ int y = x / min_factor[x]; if(min_factor[x] == min_factor[y]) return 0; x = y; pcnt++; } return pcnt % 2 == 0 ? 1 : -1; } }; bool _miller_rabin_mr(unsigned long long n, const montgomery_reduction_64bit &mr){ static std::vector small_p{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47}; static std::vector A{2, 325, 9375, 28178, 450775, 9780504, 1795265022}; static std::vector B{2, 7, 61}; if(n <= 1) return false; if(n <= 50){ for(int l = n < 20 ? 0 : 8, r = n < 20 ? 8 : 15; l < r; l++) if(small_p[l] == n) return true; return false; } if(!(n & 1)) return false; unsigned long long d = n - 1; unsigned long long one = mr.one(), mone = mr.generate(n - 1); d >>= __builtin_ctzll(d); for(unsigned long long a : (n >> 32) ? A : B){ if(a % n == 0) continue; unsigned long long d2s = d; // d * 2^s, 0 <= s <= (n - 1)が2で割れる回数 unsigned long long y = mr.pow_safe(mr.generate(a), d); while(d2s != n - 1 && y != one && y != mone){ y = mr.mul_safe(y, y); d2s <<= 1; } if(y != mone && !(d2s & 1)) return false; } return true; } bool miller_rabin_mr(unsigned long long n){ if(n % 2 == 0) return n == 2 ? true : false; montgomery_reduction_64bit mr(n); return _miller_rabin_mr(n, mr); } // https://en.wikipedia.org/wiki/Binary_GCD_algorithm unsigned long long binary_gcd(unsigned long long a, unsigned long long b){ if(!a || !b) return !a ? b : a; int shift = __builtin_ctzll(a | b); // [1] gcd(2a', 2b') = 2 * gcd(a', b') a >>= __builtin_ctzll(a); do{ // if b is odd // gcd(2a', b) = gcd(a', b), if a = 2a'(a is even) // gcd(a, b) = gcd(|a - b|, min(a, b)), if a is odd b >>= __builtin_ctzll(b); // make b odd if(a > b) std::swap(a, b); b -= a; }while(b); // gcd(a, 0) = a return a << shift; // [1] } unsigned long long generate_random_prime(unsigned long long min_n = 2, unsigned long long max_n = ~0ULL){ std::random_device seed_gen; std::mt19937_64 engine(seed_gen()); __uint128_t len = max_n - min_n + 1; // https://en.wikipedia.org/wiki/Prime_number_theorem while(true){ unsigned long long a = engine() % len + min_n; if(miller_rabin_mr(a)){ return a; } } } namespace rho_factorization{ unsigned long long rho(unsigned long long n){ static std::vector small_p{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47}; for(int sp : small_p) if(n % sp == 0) return sp; // n < 50 montgomery_reduction_64bit mr(n); if(_miller_rabin_mr(n, mr)) return n; auto try_factorize = [n, mr](unsigned long long c){ c = mr.generate(c); auto f = [mr, c](unsigned long long mx){ return mr.add(mr.mul(mx, mx), c); }; unsigned long long m = 1LL << ((64 - __builtin_clzll(n)) / 8); unsigned long long y = n, r = 1, q = 1; unsigned long long x, g, k, ys; do{ x = y; y = mr.generate(y); for(int i = 0; i < r; i++) y = f(y); y = mr.reduce(y); k = 0; while(k < r && g == 1){ q = mr.generate(q); y = mr.generate(y); ys = y; for(int i = 0; i < std::min(m, r - k); i++){ y = f(y); unsigned long long z = mr.reduce(y); q = mr.mul(q, mr.generate(x > z ? x - z : z - x)); } y = mr.reduce(y); g = binary_gcd(mr.reduce(q), n); k += m; } r <<= 1; }while(g == 1); if(g == n){ do{ ys = f(ys); unsigned long long z = mr.reduce(ys); g = binary_gcd(x > z ? x - z : z - x, n); }while(g == 1); } return g; // g == n when failure }; unsigned long long c = 1, res = n; do{ res = try_factorize(c); // c = generate_random_prime(2, n - 1); c = (c + 1) % n; }while(res == n); return res; } std::vector factorize(unsigned long long n){ if(n <= 1) return {}; unsigned long long x = rho(n); if(x == n) return {x}; auto l = factorize(x); auto r = factorize(n / x); l.insert(l.end(), r.begin(), r.end()); return l; } // {素数, 個数}の形で返す std::vector> factorize2(unsigned long long n){ auto p = factorize(n); sort(p.begin(), p.end()); std::vector> ret; for(int i : p){ if(ret.empty() || ret.back().first != i) ret.push_back({i, 1}); else ret.back().second++; } return ret; } // 素因数の集合(重複なし, ソート済)を返す std::vector prime_factor(unsigned long long n){ auto p = factorize(n); std::sort(p.begin(), p.end()); p.erase(std::unique(p.begin(), p.end()), p.end()); return p; } // 10^18以下の高度合成数 897612484786617600の約数が103680個なので全列挙して良さそう std::vector divisor(unsigned long long n){ auto p = factorize(n); std::sort(p.begin(), p.end()); std::vector> x; for(int i = 0; i < p.size(); i++){ if(!i || p[i] != p[i - 1]) x.push_back({p[i], 1}); else x.back().second++; } int sz = 1; for(auto [p_cur, cnt] : x) sz *= cnt + 1; std::vector res(sz); res[0] = 1; int r_prev = 1, r = 1; for(auto [p_cur, cnt] : x){ unsigned long long ppow = 1; for(int c = 0; c < cnt; c++){ ppow *= p_cur; for(int i = 0; i < r_prev; i++){ res[r++] = res[i] * ppow; } } r_prev = r; } return res; } int mobius_function(long long x){ auto P = rho_factorization::factorize(x); for(long long p : P) if((x / p) % p == 0) return 0; return P.size() % 2 == 0 ? 1 : -1; } unsigned long long totient_function(unsigned long long n){ unsigned long long res = n; auto prims = rho_factorization::prime_factor(n); for(auto p : prims) res -= res / p; return res; } }; // p: 素数 unsigned long long find_primitive_root(unsigned long long p){ static std::random_device seed_gen; static std::mt19937_64 engine(seed_gen()); //assert(miller_rabin_mr(p)); auto primes = rho_factorization::prime_factor(p - 1); while(true){ bool f = true; unsigned long long a = engine() % (p - 1) + 1; for(unsigned long long pk : primes){ // a ^ (p - 1) / pk ≡ 1 (mod p) -> no if(mod_pow_mr(a, (p - 1) / pk, p) == 1){ f = false; break; } } if(f) return a; } } struct bigint_prime_factor{ int sq, max_x; std::vector low; std::unordered_map high; bigint_prime_factor(int max_x): sq(sqrt(max_x)), max_x(max_x), low(sq + 1, 0){ // 篩が初期化済か assert(prime_sieve::min_factor.size() > max_x); } // O(log(x)) // x^kを掛ける void mul(int x, int k = 1){ assert(0 < x && x <= max_x); while(x > 1){ int p = prime_sieve::min_factor[x]; if(p <= sq) low[p] += k; else{ auto [itr, f] = high.emplace(x, k); if(!f) itr->second += k; } x /= p; } } // O(log(x)) // x^kで割る(割り切れない場合エラー) void div(int x, int k = 1){ assert(0 < x && x <= max_x); while(x > 1){ int p = prime_sieve::min_factor[x]; if(p <= sq){ assert(low[p] >= k); low[p] -= k; }else{ auto itr = high.find(p); assert(itr != high.end() && itr->second >= k); itr->second -= k; if(itr->second == 0) high.erase(itr); } x /= p; } } // 素数pで何回割れるか, O(1) // pが素数でない場合0 int k_divisible_prime(int p){ if(p > max_x || !prime_sieve::is_prime(p)) return 0; if(p > sq){ auto itr = high.find(p); return itr == high.end() ? 0 : itr->second; } return low[p]; } // xで何回割れるか, O(log(x)) // x = 1のときinf回割れる // @param 0 < x <= 篩の最大値 int k_divisible(int x){ assert(x > 0); static constexpr int inf = std::numeric_limits::max(); int ans = inf; for(auto [p, num] : prime_sieve::factorize(x)){ if(p > sq){ auto itr = high.find(p); if(itr == high.end()) return 0; ans = std::min(ans, itr->second / num); }else{ ans = std::min(ans, low[p] / num); } } return ans; } // rで何回割り切れるか, O(sqrt(r.max_x)以下の素数 + rの素因数の種類) int k_divisible(const bigint_prime_factor &r){ static constexpr int inf = std::numeric_limits::max(); int ans = inf; for(int i = 0; prime_sieve::primes[i] <= r.sq; i++){ int p = prime_sieve::primes[i]; if(!r.low[p]) continue; if(p <= sq){ if(low[p] < r.low[p]) return 0; ans = std::min(ans, low[p] / r.low[p]); }else{ auto itr = high.find(p); if(itr->second < r.low[p]) return 0; ans = std::min(ans, itr->second / r.low[p]); } } for(auto [p, num] : r.high){ assert(num); if(p <= sq){ if(low[p] < num) return 0; ans = std::min(ans, low[p] / num); }else{ auto itr = high.find(p); if(itr->second < num) return 0; ans = std::min(ans, itr->second / num); } } return ans; } }; // res[i] := v[i]の位数 template std::vector multiplicative_order_many(const std::vector &v){ int n = v.size(); std::vector res(n, 1); for(auto [p, q] : rho_factorization::factorize2(mint::mod() - 1)){ long long x = mint(p).pow(q).val(), y = (mint::mod() - 1) / x; for(int i = 0; i < n; i++){ long long z = x; for(int j = 0; j < q; j++){ z /= p; if(v[i].pow(y * z).val() != 1){ res[i] *= z * p; break; } } } } return res; } //rank O(1), select O(lgN), memory 2Nbit struct bitvector_memory{ static constexpr int s = 32; int n; std::vector RS{0}, table; bitvector_memory(): n(0){} bitvector_memory(const std::vector &S): n(S.size()){ int pop = 0, m = 0; for(int i = 0, t = 0; i < n; i++, t++){ if(S[i]) pop++, m += 1 << t; if(t == s - 1 || i == (n - 1)){ RS.push_back(pop); table.push_back(m); t = -1, m = 0; } } } // bit[k] bool access(int k){ assert(k < n); return (table[k / s] >> (k % s)) & 1; } // count 1, i < k int rank1(int k){ assert(k <= n); int ret = RS[k / s]; if(k % s) ret += __builtin_popcount(table[k / s] << (s - (k % s))); return ret; } // count 0, i < k int rank0(int k){ return k - rank1(k); } // kth 1, 0-indexed int select1(int k){ int L = 0, R = n + 1; while(R - L > 1){ int mid = (L + R) >> 1; if(rank1(mid) > k) R = mid; else L = mid; } return L; } // kth 0, 0-indexed int select0(int k){ int L = 0, R = n + 1; while(R - L > 1){ int mid = (L + R) >> 1; if(rank0(mid) > k) R = mid; else L = mid; } return L; } // leftmost1, k < i int succ1(int k){ return select1(rank1(k + 1)); } // leftmost0, k < i int succ0(int k){ return select0(rank0(k + 1)); } // rightmost1, i < k int pred1(int k){ int r = rank1(k); return r == 0 ? -1 : select1(r - 1); } // rightmost0, i < k int pred0(int k){ int r = rank0(k); return r == 0 ? -1 : select0(r - 1); } }; // rank O(1), select O(1), memory (2 + 32)Nbit struct bitvector_fast_select{ static constexpr int s = 32; int n; std::vector pos1, pos0; std::vector RS{0}, table; bitvector_fast_select(){} bitvector_fast_select(const std::vector &v): n(v.size()){ for(int i = 0; i < n; i++){ if(v[i]) pos1.push_back(i); else pos0.push_back(i); } int pop = 0, m = 0; for(int i = 0, t = 0; i < n; i++, t++){ if(v[i]) pop++, m += 1 << t; if(t == s - 1 || i == (n - 1)){ RS.push_back(pop); table.push_back(m); t = -1, m = 0; } } } // bit[k] bool access(int k){ assert(k < n); return (table[k / s] >> (k % s)) & 1; } // count 1, i < k int rank1(int k){ assert(k <= n); int ret = RS[k / s]; if(k % s) ret += __builtin_popcount(table[k / s] << (s - (k % s))); return ret; } // count 0, i < k int rank0(int k){ return k - rank1(k); } int select1(int k){ if(k >= pos1.size()) return n; return pos1[k]; } int select0(int k){ if(k >= pos0.size()) return n; return pos0[k]; } // leftmost1, k < i int succ1(int k){ return select1(rank1(k + 1)); } // leftmost0, k < i int succ0(int k){ return select0(rank0(k + 1)); } // rightmost1, i < k int pred1(int k){ int r = rank1(k); return r == 0 ? -1 : select1(r - 1); } // rightmost0, i < k int pred0(int k){ int r = rank0(k); return r == 0 ? -1 : select0(r - 1); } }; /* https://codeforces.com/contest/1746/my バグ // rank O(1), select O(1), memory (2 + 4 + (0 ~ 4))N bit struct bitvector{ private: constexpr static int s = 16, th_l = 255, th_m = 64, s2 = 32; int n; // rank std::vector RS{0}, table; // select struct block{ int t:2, idx:30; block(int t, int idx): t(t), idx(idx){} }; std::array, 2> B; std::array>, 2> large; std::array>, 2> medium; std::array, 2> block_idx; static uint64_t flip(uint64_t x){ constexpr static uint64_t y = 0xffffffffffffffffULL; return x ^ y; } static int find_kth_set_bit(uint64_t mask, int k){ constexpr static uint64_t m1 = 0x5555555555555555ULL; // even bits constexpr static uint64_t m2 = 0x3333333333333333ULL; // even 2-bit groups constexpr static uint64_t m4 = 0x0f0f0f0f0f0f0f0fULL; // even nibbles constexpr static uint64_t m8 = 0x00ff00ff00ff00ffULL; // even bytes int t, i = k, r = 0; uint64_t c1 = mask; uint64_t c2 = c1 - ((c1 >> 1) & m1); uint64_t c4 = ((c2 >> 2) & m2) + (c2 & m2); uint64_t c8 = ((c4 >> 4) + c4) & m4; uint64_t c16 = ((c8 >> 8) + c8) & m8; uint64_t c32 = (c16 >> 16) + c16; int c64 = (int)(((c32 >> 32) + c32) & 0x7f); t = (c32 ) & 0x3f; if (i >= t) { r += 32; i -= t; } t = (c16>> r) & 0x1f; if (i >= t) { r += 16; i -= t; } t = (c8 >> r) & 0x0f; if (i >= t) { r += 8; i -= t; } t = (c4 >> r) & 0x07; if (i >= t) { r += 4; i -= t; } t = (c2 >> r) & 0x03; if (i >= t) { r += 2; i -= t; } t = (c1 >> r) & 0x01; if (i >= t) { r += 1; } if (k >= c64) r = -1; return r; } // get[l, r], len <= 64 uint64_t get64(int l, int r){ int lb = l / s2, rb = r / s2; l %= s2, r %= s2; if(rb - lb == 2){ uint64_t ans_l = table[lb] >> l; uint32_t ans_r = table[rb] << (s2 - 1 - r); return ans_l + (table[lb + 1] << (s2 - l)) + ((uint64_t)ans_r << (r - l + s2)); }else if(rb - lb == 1){ uint64_t ans_l = table[lb] >> l; uint32_t ans_r = table[rb] << (s2 - 1 - r); return ans_l + ((uint64_t)ans_r << (r - l)); }else{ assert(lb == rb); uint32_t ans = table[lb]; ans <<= (s2 - 1 - r); ans >>= (s2 - (r - l + 1)); return ans; } } int select(int k, bool b){ int idx_l = k / s; if(idx_l >= block_idx[b].size()) return n; int ans = block_idx[b][idx_l], block_type = B[b][idx_l].t, idx_s = B[b][idx_l].idx; if(block_type == 2) ans += large[b][idx_s][k % s]; else if(block_type == 1) ans += medium[b][idx_s][k % s]; else{ int L = block_idx[b][idx_l], R = B[b][idx_l].idx; uint64_t x = b ? get64(L, R) : flip(get64(L, R)); int tmp = find_kth_set_bit(x, k % s); ans = (tmp == -1 || ans + tmp >= n ? n : ans + tmp); } return ans; } public: bitvector(){} bitvector(const std::vector &v): n(v.size()){ // rank int pop = 0, m = 0; for(int i = 0, t = 0; i < n; i++, t++){ if(v[i]) pop++, m += 1 << t; if(t == s2 - 1 || i == (n - 1)){ RS.push_back(pop); table.push_back(m); t = -1, m = 0; } } // select // 区間が[l, r]の0/1のブロックを作成 auto make_block = [&](int l, int r, bool b){ block_idx[b].push_back(l); if(r - l >= th_l){ B[b].push_back(block(2, large[b].size())); large[b].push_back(std::array()); int k = 0; for(int i = l; i <= r; i++) if(v[i] == b) large[b].back()[k++] = i - l; while(k < s) large[b].back()[k++] = r - l + 1; }else if(r - l >= th_m){ B[b].push_back(block(1, medium[b].size())); medium[b].push_back(std::array()); int k = 0; for(int i = l; i <= r; i++) if(v[i] == b) medium[b].back()[k++] = i - l; while(k < s) large[b].back()[k++] = r - l + 1;// [0, 255] }else{ B[b].push_back(block(0, r)); } }; for(int b = 0; b < 2; b++){ for(int i = 0, k = 0, l = 0; i < n; i++){ if(v[i] == b){ k++; if(k % s == 0) make_block(l, i, b), l = i + 1; } if(i == n - 1 && (block_idx[b].empty() || block_idx[b].back() != l)) make_block(l, i, b); } } } int size(){ return n; } // bit[k] bool access(int k){ assert(k < n); return (table[k / s2] >> (k % s2)) & 1; } // count 1, i < k int rank1(int k){ assert(0 <= k && k <= n); int ret = RS[k / s2]; if(k % s2) ret += __builtin_popcount(table[k / s2] << (s2 - (k % s2))); return ret; } // count 0, i < k int rank0(int k){ return k - rank1(k); } int select1(int k){ return select(k, 1); } int select0(int k){ return select(k, 0); } }; */ template struct bitvector_arbitrary_radix{ static constexpr int s = 32, sdiv = 5, smod = 31; using Z = uint32_t;//s bit int N; std::vector> B; std::vector> S; bitvector_arbitrary_radix(): N(0){} // init O(N + NR/s) bitvector_arbitrary_radix(const std::vector &v): N(v.size()){ int M = (N + s - 1) / s; std::array pop; std::array pop_small; pop.fill(0); pop_small.fill(0); B.resize(M + 1, pop); S.resize(M, pop_small); for(int i = 0, t = 0, sz = 0; i < N; i++, t++){ int x = v[i]; assert(0 <= x && x < R); pop[x]++, pop_small[x] |= (Z(1) << t); if(t == s - 1 || i == N - 1){ for(int j = 0; j < R; j++){ if(j) pop[j] += pop[j - 1], pop_small[j] |= pop_small[j - 1]; B[sz + 1][j] = pop[j] + B[sz][j]; S[sz][j] = pop_small[j]; } pop.fill(0); pop_small.fill(0); t = -1; sz++; } } } // r未満のcの数 int rank(int r, int c){ int rq = r >> sdiv, rm = r & smod; int ret = B[rq][c] - (c ? B[rq][c - 1] : 0); if(rm) ret += __builtin_popcount((S[rq][c] ^ (c ? S[rq][c - 1] : 0)) << (s - rm)); return ret; } // r未満のc未満の数 int rank_lower(int r, int c){ if(c == 0) return 0; int rq = r >> sdiv, rm = r & smod; int ret = B[rq][--c]; if(rm) ret += __builtin_popcount(S[rq][c] << (s - rm)); return ret; } }; struct wavelet_matrix{ private: using __bitvector = bitvector_memory; int n, h, inf; std::vector<__bitvector> bv; std::vector bottom_idx; void build(std::vector v){ bv.resize(h); std::vector bits(n); std::vector tmp(n), tmp_idx(n); std::iota(bottom_idx.begin(), bottom_idx.end(), 0); std::iota(tmp_idx.begin(), tmp_idx.end(), 0); std::queue> q; if(h) q.push({h - 1, 0, n}); while(!q.empty()){ auto [d, l, r] = q.front(); q.pop(); int lidx = 0, ridx = 0; for(int i = l; i < r; i++){ bool b = (v[i] >> d) & 1; bits[i] = b; if(b) tmp[ridx] = v[i], tmp_idx[ridx++] = bottom_idx[i]; else v[l + lidx] = v[i], bottom_idx[l + lidx++] = bottom_idx[i]; } for(int i = 0; i < ridx; i++){ v[l + lidx + i] = tmp[i]; bottom_idx[l + lidx + i] = tmp_idx[i]; } if(d){ int mid = l + lidx; if(mid > l) q.push({d - 1, l, mid}); if(r > mid) q.push({d - 1, mid, r}); } if(r == n) bv[d] = __bitvector(bits); } } // v[k] int __access(int k){ assert(0 <= k && k < n); int L = 0, R = n, ret = 0; for(int i = h - 1; i >= 0; i--){ int L0 = bv[i].rank0(L), R0 = bv[i].rank0(R), k0 = bv[i].rank0(k); if(bv[i].access(k)){ k += R0 - k0; L += R0 - L0; ret += 1 << i; }else{ k -= (k - L) - (k0 - L0); R = L + R0 - L0; } } return ret; } // [0, r)のcの数 int __rank(int r, int c){ if(c < 0 || c >= inf) return 0; int L = 0, R = n; for(int i = h - 1; i >= 0 && r > L; i--){ int L0 = bv[i].rank0(L), R0 = bv[i].rank0(R) - L0, r0 = bv[i].rank0(r) - L0; if((c >> i) & 1){ r += R0 - r0; L += R0; }else{ r -= (r - L) - r0; R = L + R0; } } return r - L; } // k番目のc, ない場合は-1 int __select(int k, int c){ if(c < 0 || c >= inf) return -1; int L = 0, R = n; for(int i = h - 1; i >= 0; i--){ int LR0 = bv[i].rank0(R) - bv[i].rank0(L); if((c >> i) & 1){ L += LR0; }else{ R = L + LR0; } } if(R - L <= k) return -1; return bottom_idx[L + k]; } int __range_freq(int d, int l, int r, int L, int R, int s, int t, int S, int T){ if(l == r || L == R || t <= S || T <= s) return 0; else if(s <= S && T <= t) return r - l; int L0 = bv[d].rank0(L), R0 = bv[d].rank0(R) - L0; int l0 = bv[d].rank0(l) - L0, r0 = bv[d].rank0(r) - L0; return __range_freq(d - 1, L + l0, L + r0, L, L + R0, s, t, S, (S + T) / 2) + __range_freq(d - 1, l + (R0 - l0), r + (R0 - r0), L + R0, R, s, t, (S + T) / 2, T); } std::pair __range_kth_smallest(int l, int r, int k){ if(l >= r || r - l <= k) return {-1, -1}; int L = 0, R = n, ret = 0; for(int i = h - 1; i >= 0; i--){ int L0 = bv[i].rank0(L), R0 = bv[i].rank0(R) - L0; int l0 = bv[i].rank0(l) - L0, r0 = bv[i].rank0(r) - L0; if(r0 - l0 <= k){ k -= r0 - l0; ret += 1 << i; l += R0 - l0, r += R0 - r0, L += R0; }else{ l = L + l0, r = L + r0, R = L + R0; } } return {ret, bottom_idx[l + k]}; } std::pair __range_kth_smallest_super(int d, int l, int r, int L, int R, int s, int t, int S, int T, int &k){ if(l == r || L == R || T <= s || t <= S) return {-1, -1}; else if(s <= S && T <= t){ if(r - l > k){ for(int i = d; i >= 0; i--){ int L0 = bv[i].rank0(L), R0 = bv[i].rank0(R) - L0; int l0 = bv[i].rank0(l) - L0, r0 = bv[i].rank0(r) - L0; if(r0 - l0 <= k){ k -= r0 - l0; l += R0 - l0, r += R0 - r0, L += R0; S = (S + T) / 2; }else{ l = L + l0, r = L + r0, R = L + R0; T = (S + T) / 2; } } return {S, bottom_idx[l + k]}; } k -= r - l; return {-1, -1}; } int L0 = bv[d].rank0(L), R0 = bv[d].rank0(R) - L0; int l0 = bv[d].rank0(l) - L0, r0 = bv[d].rank0(r) - L0; auto p = __range_kth_smallest_super(d - 1, L + l0, L + r0, L, L + R0, s, t, S, (S + T) / 2, k); if(p.first != -1) return p; return __range_kth_smallest_super(d - 1, l + (R0 - l0), r + (R0 - r0), L + R0, R, s, t, (S + T) / 2, T, k); } std::pair __range_kth_largest_super(int d, int l, int r, int L, int R, int s, int t, int S, int T, int &k){ if(l == r || L == R || T <= s || t <= S) return {-1, -1}; else if(s <= S && T <= t){ if(r - l > k){ for(int i = d; i >= 0; i--){ int L0 = bv[i].rank0(L), R0 = bv[i].rank0(R) - L0; int l0 = bv[i].rank0(l) - L0, r0 = bv[i].rank0(r) - L0; if((r - l) - (r0 - l0) <= k){ k -= (r - l) - (r0 - l0); l = L + l0, r = L + r0, R = L + R0; T = (S + T) / 2; }else{ l += R0 - l0, r += R0 - r0, L += R0; S = (S + T) / 2; } } return {S, bottom_idx[r - 1 - k]}; } k -= r - l; return {-1, -1}; } int L0 = bv[d].rank0(L), R0 = bv[d].rank0(R) - L0; int l0 = bv[d].rank0(l) - L0, r0 = bv[d].rank0(r) - L0; auto p = __range_kth_largest_super(d - 1, l + (R0 - l0), r + (R0 - r0), L + R0, R, s, t, (S + T) / 2, T, k); if(p.first != -1) return p; return __range_kth_largest_super(d - 1, L + l0, L + r0, L, L + R0, s, t, S, (S + T) / 2, k); } void __count_prefix(int l, int r, int x, std::vector &ret){ int L = 0, R = n; for(int i = h - 1; i >= 0; i--){ int L0 = bv[i].rank0(L), R0 = bv[i].rank0(R) - L0; int l0 = bv[i].rank0(l) - L0, r0 = bv[i].rank0(r) - L0; if((x >> i) & 1){ ret[h - 1 - i] = r0 - l0; l += R0 - l0, r += R0 - r0, L += R0; }else{ ret[h - 1 - i] = (r - l) - (r0 - l0); l = L + l0, r = L + r0, R = L + R0; } } ret[h] = r - l; } public: wavelet_matrix(): n(0){} // 値が[0, inf) wavelet_matrix(const std::vector &v, int inf): n(v.size()), inf(inf), bottom_idx(n){ assert(inf >= 0); h = 0; while((1 << h) < inf) h++; build(v); } // v[k] int access(int k){ return __access(k); } // [0, r)のcの数 int rank(int r, int c){ return __rank(r, c); } // k番目のc, ない場合は-1 int select(int k, int c){ return __select(k, c); } // k以降(k含む)のc, 無い場合は-1 int find_next(int k, int c){ if(c < 0 || c >= inf) return -1; return __select(__rank(k, c), c); } // k以前(k含む)のc, 無い場合は-1 int find_prev(int k, int c){ if(c < 0 || c >= inf) return -1; int r = __rank(k + 1, c); if(r == 0) return -1; return __select(r - 1, c); } // [l, r)の[s, t)の数 int range_freq(int l, int r, int s, int t){ assert(0 <= l && r <= n); assert(0 <= s && t <= inf); if(l >= r || s >= t) return 0; return __range_freq(h - 1, l, r, 0, n, s, t, 0, 1 << h); } // [l, r)でk番目に小さい要素とインデックス(ない場合は-1) // 値が同じ場合インデックスが小さいものを小さいとする std::pair range_kth_smallest(int l, int r, int k){ return __range_kth_smallest(l, r, k); } // [l, r)でk番目に大きい要素とインデックス(ない場合は-1) // 値が同じ場合インデックスが小さいものを小さいとする std::pair range_kth_largest(int l, int r, int k){ if(r - l <= k) return {-1, -1}; return __range_kth_smallest(l, r, r - l - 1 - k); } // [l, r)で値が[s, t)の要素でk番目に小さい // 値が同じ場合インデックスが小さいものを小さいとする std::pair range_kth_smallest_super(int l, int r, int s, int t, int k){ return __range_kth_smallest_super(h - 1, l, r, 0, n, s, t, 0, 1 << h, k); } // [l, r)で値が[s, t)の要素でk番目に大きい // 値が同じ場合インデックスが小さいものを小さいとする std::pair range_kth_largest_super(int l, int r, int s, int t, int k){ return __range_kth_largest_super(h - 1, l, r, 0, n, s, t, 0, 1 << h, k); } // 区間[l, r)のx以上の最小要素(ない場合は-1) int range_lower_bound(int l, int r, int x){ assert(0 <= l && r <= n); assert(0 <= x && x < inf); int cnt = range_freq(l, r, 0, x); return range_kth_smallest(l, r, cnt).first; } // 区間[l, r)のx以下の最大要素(ない場合は-1) int range_lower_bound_rev(int l, int r, int x){ assert(0 <= l && r <= n); assert(0 <= x && x < inf); int cnt = range_freq(l, r, 0, x + 1); if(cnt == 0) return -1; return range_kth_smallest(l, r, cnt - 1).first; } // [l, r)で値が[s, t)の要素でインデックス順にk番目, ない場合は-1 int range_select(int l, int r, int s, int t, int k){ if(range_freq(l, r, s, t) <= k) return -1; int L = l; while(r - l > 1){ int mid_r = (l + r) / 2; if(range_freq(L, mid_r, s, t) > k) r = mid_r; else l = mid_r; } return l; } // ret[i] := [l, r)でh-bitの2進数で見たときxとのlcpがちょうどiの要素の数 std::vector count_prefix(int l, int r, int x){ std::vector ret(h + 1, 0); __count_prefix(l, r, x, ret); return ret; } }; template struct compressed_wavelet_matrix{ private: std::vector rev; wavelet_matrix wm; int lb(T c){ return std::lower_bound(rev.begin(), rev.end(), c) - rev.begin(); } public: compressed_wavelet_matrix(){} compressed_wavelet_matrix(const std::vector &_v){ int n = _v.size(); std::vector v(n); std::vector> tmp(n); for(int i = 0; i < n; i++) tmp[i].first = _v[i], tmp[i].second = i; std::sort(tmp.begin(), tmp.end()); for(int i = 0; i < n; i++){ if(i == 0 || tmp[i].first != tmp[i - 1].first) rev.push_back(tmp[i].first); v[tmp[i].second] = (int)rev.size() - 1; } tmp.clear(); wm = wavelet_matrix(v, rev.size()); } // V[k] T access(int k){ return rev[wm.access(k)]; } // インデックスがr未満のcの数 int rank(int r, T c){ int idx = lb(c); if(idx == rev.size() || rev[idx] != c) return 0; return wm.rank(r, idx); } // k番目のc, 無い場合は-1 int select(int k, T c){ int idx = lb(c); if(idx == rev.size() || rev[idx] != c) return -1; return wm.select(k, idx); } // k以降(k含む)のc, 無い場合は-1 int find_next(int k, T c){ int idx = lb(c); if(idx == rev.size() || rev[idx] != c) return -1; return wm.find_next(k, idx); } // k以前(k含む)のc, 無い場合は-1 int find_prev(int k, T c){ int idx = lb(c); if(idx == rev.size() || rev[idx] != c) return -1; return wm.find_prev(k, idx); } // [l, r)の[s, t)の個数 int range_freq(int l, int r, T s, T t){ return wm.range_freq(l, r, lb(s), lb(t)); } // [l, r)でk番目に小さい要素と場所, ない場合は-1 // 値が同じ場合インデックスが小さいものを小さいとする std::pair range_kth_smallest(int l, int r, int k){ auto p = wm.range_kth_smallest(l, r, k); if(p.first == -1) return {-1, -1}; return {rev[p.first], p.second}; } // [l, r)でk番目に大きい要素と場所, ない場合は-1 // 値が同じ場合インデックスが小さいものを小さいとする std::pair range_kth_largest(int l, int r, int k){ auto p = wm.range_kth_largest(l, r, k); if(p.first == -1) return {-1, -1}; return {rev[p.first], p.second}; } // [l, r)の[s, t)の中でk番目に小さい要素と場所, ない場合は-1 // 値が同じ場合インデックスが小さいものを小さいとする std::pair range_kth_smallest_super(int l, int r, T s, T t, int k){ auto p = wm.range_kth_smallest_super(l, r, lb(s), lb(t), k); if(p.first == -1) return {-1, -1}; return {rev[p.first], p.second}; } // [l, r)の[s, t)の中でk番目に大きい要素と場所, ない場合は-1 // 値が同じ場合インデックスが小さいものを小さいとする std::pair range_kth_largest_super(int l, int r, T s, T t, int k){ auto p = wm.range_kth_largest_super(l, r, lb(s), lb(t), k); if(p.first == -1) return {-1, -1}; return {rev[p.first], p.second}; } // 区間[l, r)のx以上の最小要素(ない場合は-1) T range_lower_bound(int l, int r, T x){ int xi = lb(x); if(xi == rev.size()) return -1; int ret = wm.range_lower_bound(l, r, xi); return ret == -1 ? -1 : rev[ret]; } // 区間[l, r)のx以下の最大要素(ない場合は-1) T range_lower_bound_rev(int l, int r, T x){ int xi = std::upper_bound(rev.begin(), rev.end(), x) - rev.begin(); if(xi == 0) return -1; int ret = wm.range_lower_bound_rev(l, r, xi - 1); return ret == -1 ? -1 : rev[ret]; } // [l, r)で値が[s, t)の要素でインデックス順にk番目, ない場合は-1 int range_select(int l, int r, T s, T t, int k){ return wm.range_select(l, r, lb(s), lb(t), k); } }; template struct wavelet_matrix_arbitrary_radix{ private: static constexpr int bitlen(uint64_t x){ if(!x) return 0; return 64 - __builtin_clzll(x); } static constexpr int Rdiv = bitlen(R) - 1, Rmod = R - 1; int n, h, inf; std::vector> bv; std::vector bottom_idx; void build(const std::vector &v){ bv.resize(h); bottom_idx.resize(n); std::queue> q; if(h) q.push({h - 1, 0, n}); std::vector> cur(n); std::array>, R> next; std::vector bits(n); for(int i = 0; i < n; i++) cur[i] = {v[i], i}; while(!q.empty()){ auto [d, l, r] = q.front(); q.pop(); for(int i = l; i < r; i++){ int dir = (cur[i].first >> (d * Rdiv)) & Rmod; bits[i] = dir; next[dir].push_back({cur[i].first, cur[i].second}); } for(int i = 0, j = l; i < R; i++){ std::copy(next[i].begin(), next[i].end(), cur.begin() + j); int k = j + next[i].size(); if(d && j < k) q.push({d - 1, j, k}); j = k; next[i].clear(); } if(r == n) bv[d] = bitvector_arbitrary_radix(bits); } for(int i = 0; i < n; i++) bottom_idx[i] = cur[i].second; } int __rank(int r, int c){ int l = 0, _L = 0, _R = n; for(int d = h - 1, shift = d * Rdiv; r > l && d >= 0; d--, shift -= Rdiv){ int dir = (c >> shift) & Rmod; int Ldir = bv[d].rank_lower(_L, dir); int Ldir2 = bv[d].rank(_L, dir); int Rdir = bv[d].rank_lower(_R, dir); _L += Rdir - Ldir; _R = _L + bv[d].rank(_R, dir) - Ldir2; l = _L + bv[d].rank(l, dir) - Ldir2; r = _L + bv[d].rank(r, dir) - Ldir2; } return r - l; } int __select(int k, int c){ int _L = 0, _R = n; for(int d = h - 1, shift = d * Rdiv; d >= 0; d--, shift -= Rdiv){ int dir = (c >> shift) & Rmod; int Ldir = bv[d].rank_lower(_L, dir), Ldir2 = bv[d].rank(_L, dir), Rdir = bv[d].rank_lower(_R, dir); _L += Rdir - Ldir, _R = _L + bv[d].rank(_R, dir) - Ldir2; if(_L + k >= _R) return -1; } return bottom_idx[_L + k]; } int __range_freq(int d, int l, int r, int _L, int _R, int s, int t, int S, int T){ if(l >= r || s >= t || t <= S || T <= s) return 0; else if(s <= S && T <= t) return r - l; if(s < S) s = S; if(t > T) t = T; int shift = d * Rdiv; int sb = (s - S) >> shift, tb = (t - S) >> shift; // [0, R] if(sb == tb){ int Lsb = bv[d].rank_lower(_L, sb), Lsb2 = bv[d].rank(_L, sb), Rsb = bv[d].rank_lower(_R, sb); int Lnext = _L + Rsb - Lsb, Rnext = Lnext + bv[d].rank(_R, sb) - Lsb2; int lnext = Lnext + bv[d].rank(l, sb) - Lsb2, rnext = Lnext + bv[d].rank(r, sb) - Lsb2; return __range_freq(d - 1, lnext, rnext, Lnext, Rnext, s, t, S + (sb << shift), S + ((sb + 1) << shift)); } // f(tb) - f(sb) + [tb] - [sb] int ans = 0; ans += bv[d].rank_lower(r, tb) - bv[d].rank_lower(l, tb); ans -= bv[d].rank_lower(r, sb) - bv[d].rank_lower(l, sb); // [sb << d, s) if((sb << shift) < s){ int Lsb = bv[d].rank_lower(_L, sb), Lsb2 = bv[d].rank(_L, sb), Rsb = bv[d].rank_lower(_R, sb); int Lnext = _L + Rsb - Lsb, Rnext = Lnext + bv[d].rank(_R, sb) - Lsb2; int lnext = Lnext + bv[d].rank(l, sb) - Lsb2, rnext = Lnext + bv[d].rank(r, sb) - Lsb2; ans -= __range_freq(d - 1, lnext, rnext, Lnext, Rnext, S + (sb << shift), s, S + (sb << shift), S + ((sb + 1) << shift)); } // [tb << d, t) if((tb << shift) < t){ int Ltb = bv[d].rank_lower(_L, tb), Ltb2 = bv[d].rank(_L, tb), Rtb = bv[d].rank_lower(_R, tb); int Lnext = _L + Rtb - Ltb, Rnext = Lnext + bv[d].rank(_R, tb) - Ltb2; int lnext = Lnext + bv[d].rank(l, tb) - Ltb2, rnext = Lnext + bv[d].rank(r, tb) - Ltb2; ans += __range_freq(d - 1, lnext, rnext, Lnext, Rnext, S + (tb << shift), t, S + (tb << shift), S + ((tb + 1) << shift)); } return ans; } std::pair __range_kth_smallest(int l, int r, int k){ int _L = 0, _R = n, ans = 0; for(int d = h - 1, shift = d * Rdiv; d >= 0; d--, shift -= Rdiv){ int _l = -1, _r = R - 1; while(_r - _l > 1){ int mid = (_l + _r) >> 1; (bv[d].rank_lower(r, mid + 1) - bv[d].rank_lower(l, mid + 1) > k ? _r : _l) = mid; } int dir = _r; k -= bv[d].rank_lower(r, dir) - bv[d].rank_lower(l, dir); ans += dir << shift; int Ldir = bv[d].rank_lower(_L, dir), Ldir2 = bv[d].rank(_L, dir), Rdir = bv[d].rank_lower(_R, dir); _L += Rdir - Ldir, _R = _L + bv[d].rank(_R, dir) - Ldir2; l = _L + bv[d].rank(l, dir) - Ldir2, r = _L + bv[d].rank(r, dir) - Ldir2; } return {ans, bottom_idx[l + k]}; } public: wavelet_matrix_arbitrary_radix(){} wavelet_matrix_arbitrary_radix(const std::vector &v, int inf): n(v.size()), inf(inf){ // 高速化のためRは2べき static_assert((1 << Rdiv) == R); h = 0; long long k = 1; while(k < inf) h++, k <<= Rdiv; build(v); } // r未満のcの数 int rank(int r, int c){ assert(0 <= r && r <= n); if(c < 0 || c >= inf) return 0; return __rank(r, c); } // k番目のcのインデックス, ない場合は-1 int select(int k, int c){ assert(0 <= k); if(c < 0 || c >= inf) return -1; return __select(k, c); } // [l, r)の[s, t)の数 int range_freq(int l, int r, int s, int t){ assert(0 <= l && r <= n); assert(0 <= s && t <= inf); if(l >= r || s >= t) return 0; return __range_freq(h - 1, l, r, 0, n, s, t, 0, 1 << (h * Rdiv)); } // 区間[l, r)でk番目に小さい要素とそのインデックス, ない場合は-1 // 値が同じ場合インデックスの昇順 std::pair range_kth_smallest(int l, int r, int k){ if(r - l <= k) return {-1, -1}; return __range_kth_smallest(l, r, k); } // 区間[l, r)でk番目に大きい要素とそのインデックス, ない場合は-1 // 値が同じ場合インデックスの降順 std::pair range_kth_largest(int l, int r, int k){ if(r - l <= k) return {-1, -1}; return __range_kth_smallest(l, r, r - l - 1 - k); } // [l, r)で値が[s, t)の要素でインデックス順にk番目, ない場合は-1 int range_select(int l, int r, int s, int t, int k){ if(range_freq(l, r, s, t) <= k) return -1; int L = l; while(r - l > 1){ int mid_r = (l + r) / 2; if(range_freq(L, mid_r, s, t) > k) r = mid_r; else l = mid_r; } return l; } }; template struct compressed_wavelet_matrix_arbitrary_radix{ private: int N; wavelet_matrix_arbitrary_radix wm; std::vector rev; int lb(Val c){return std::lower_bound(rev.begin(), rev.end(), c) - rev.begin();} public: compressed_wavelet_matrix_arbitrary_radix(){} compressed_wavelet_matrix_arbitrary_radix(const std::vector &_v): N(_v.size()){ std::vector v(N); std::vector> tmp(N); for(int i = 0; i < N; i++) tmp[i] = {_v[i], i}; std::sort(tmp.begin(), tmp.end()); for(int i = 0; i < N; i++){ if(i == 0 || tmp[i].first != tmp[i - 1].first) rev.push_back(tmp[i].first); v[tmp[i].second] = (int)rev.size() - 1; } tmp.clear(); wm = wavelet_matrix_arbitrary_radix(v, rev.size()); } // val[k] Val access(int k){ return rev[wm.access(k)]; } // [0, r)のc int rank(int r, Val c){ int idx = lb(c); if(idx == rev.size() || rev[idx] != c) return 0; return wm.rank(r, idx); } // k番目のcのインデックス, ない場合は-1 int select(int k, Val c){ int idx = lb(c); if(idx == rev.size() || rev[idx] != c) return -1; return wm.select(k, idx); } // [l, r)の[s, t)の数 int range_freq(int l, int r, Val s, Val t){ return wm.range_freq(l, r, lb(s), lb(t)); } // 区間[l, r)でk番目に小さい要素とそのインデックス, ない場合は-1 // 値が同じ場合インデックスの昇順 std::pair range_kth_smallest(int l, int r, int k){ auto p = wm.range_kth_smallest(l, r, k); if(p.first == -1) return {-1, -1}; return {rev[p.first], p.second}; } // 区間[l, r)でk番目に大きい要素とそのインデックス, ない場合は-1 // 値が同じ場合インデックスの降順 std::pair range_kth_largest(int l, int r, int k){ auto p = wm.range_kth_largest(l, r, k); if(p.first == -1) return {-1, -1}; return {rev[p.first], p.second}; } // [l, r)で値が[s, t)の要素でインデックス順にk番目, ない場合は-1 int range_select(int l, int r, Val s, Val t, int k){ return wm.range_select(l, r, lb(s), lb(t), k); } }; int main(){ prime_sieve::init(200001); io_init(); int n; std::cin >> n; vector> v; range(i, 0, n){ int a, b; //a = 166320; //b = 1000000000; std::cin >> a >> b; v.push_back({a, b}); } sort(allof(v)); vector mxb(200001, 0); range(i, 0, n){ auto [a, b] = v[i]; mxb[a] = max(mxb[a], b); } int lower = 0; range(i, 1, n){ if(v[i].first == v[i - 1].first){ lower = max(lower, min(v[i].second, v[i - 1].second)); } } // x | ai // x | aj // ci := ai / g // cj := aj / g // min(floor(bi / cj), floor(bj / ci)) // gを小さく取ってしまう -> cが大きくなる -> 答えが小さくなるため問題ない // S[i] := iの倍数の集合 vector>> S(200001); range(i, 0, n){ auto [a, b] = v[i]; for(int d : prime_sieve::divisor(v[i].first)){ S[d].push_back({v[i].first / d, v[i].second}); if(a != d) lower = max(lower, min(b, mxb[d] / (a / d))); } } int ans = lower; for(int i = 200000; i > 0; i--){ if(S[i].size() <= 10){ int sz = S[i].size(); range(j, 0, sz){ auto [ci, bi] = S[i][j]; range(k, j + 1, sz){ auto [cj, bj] = S[i][k]; int x = min(bi / cj, bj / ci); ans = max(ans, x); } } }else{ sort(allof(S[i])); vector B(S[i].size()); range(j, 0, S[i].size()) B[j] = S[i][j].second; compressed_wavelet_matrix wm(B); int l = ans, r = 1000000001; while(r - l > 1){ int x = (l + r) / 2; bool f = false; range(j, 0, S[i].size()){ auto [c, b] = S[i][j]; int ridx = lower_bound(allof(S[i]), pair{b / x + 2, 0}) - S[i].begin(); // [0, r) int y = 0; if((ll)x * c <= 1000000000) y = wm.range_freq(0, ridx, 0, x * c); y -= (b >= (ll)c * x ? 1 : 0); if(y > 0){ f = true; break; } } if(f) l = x; else r = x; } ans = max(ans, l); } } std::cout << ans << '\n'; }