・テンプレート
#ifndef ONLINE_JUDGE
#define _GLIBCXX_DEBUG
#endif
#include <bits/stdc++.h>
using namespace std;

#define pass (void)0
#define INF (1<<30)-1
#define INFLL (1LL<<60)-1
#define rep(i, n) for (int i = 0; i < (int)(n); i++)
#define repr(i, n) for (int i = (int)(n) - 1; i >= 0; i--)
#define rep2(i, a, b) for (int i = (int)(a); i < (int)(b); i++)
#define repr2(i, a, b) for (int i = (int)(b) - 1; i >= (int)(a); i--)
#define all(x) (x).begin(), (x).end()
#define rall(x) (x).rbegin(), (x).rend()
#define sz(x) ((int)(x).size())
#define YesNo(cond) cout << ((cond) ? "Yes\n" : "No\n")
#define YESNO(cond) cout << ((cond) ? "YES\n" : "NO\n")

using ll = long long;
using pii = pair<int,int>;
using pll = pair<ll,ll>;
using vi = vector<int>;
using vl = vector<ll>;
using vvi = vector<vi>;
using vvl = vector<vl>;

template <typename T> void print(const T& value) { cout << value << "\n"; }
template <typename T> void print(const vector<T>& vec) { for (auto& v : vec) cout << v << " "; cout << "\n"; }
template <typename T> void input(vector<T>& vec) { for (auto& v : vec) cin >> v; };
template <typename T> bool chmin(T& a, const T& b) { if (a > b) { a = b; return true; } return false; }
template <typename T> bool chmax(T& a, const T& b) { if (a < b) { a = b; return true; } return false; }

#include <atcoder/all>
using namespace atcoder;
using mint = modint998244353;
using vm = vector<mint>;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cout << fixed << setprecision(10);
}

・遅延セグ木の関数
// 区間加算・区間最小値取得
ll op(ll a, ll b){ return min(a, b); }
ll e(){ return INFLL; }
ll mapping(ll f, ll x){ return f+x; }
ll composition(ll f, ll g){ return f+g; }
ll id(){ return 0; }

// 区間加算・区間最大値取得
ll op(ll a, ll b){ return max(a, b); }
ll e(){ return -INFLL; }
ll mapping(ll f, ll x){ return f+x; }
ll composition(ll f, ll g){ return f+g; }
ll id(){ return 0; }

// 区間加算・区間和取得
struct S{
    ll value;
    int size;
};
S op(S a, S b){ return {a.value+b.value, a.size+b.size}; }
S e(){ return {0, 0}; }
S mapping(ll f, S x){ return {x.value + f*x.size, x.size}; }
ll composition(ll f, ll g){ return f+g; }
ll id(){ return 0; }

// 区間変更・区間最小値取得
ll op(ll a, ll b){ return min(a, b); }
ll e(){ return INFLL; }
ll mapping(ll f, ll x){ return (f == INFLL ? x : f); }
ll composition(ll f, ll g){ return (f == INFLL ? g : f); }
ll id(){ return INFLL; }

// 区間変更・区間最大値取得
ll op(ll a, ll b){ return max(a, b); }
ll e(){ return -INFLL; }
ll mapping(ll f, ll x){ return (f == -INFLL ? x : f); }
ll composition(ll f, ll g){ return (f == -INFLL ? g : f); }
ll id(){ return -INFLL; }

// 区間変更・区間和取得
struct S{
    ll value;
    int size;
};
S op(S a, S b){ return {a.value+b.value, a.size+b.size}; }
S e(){ return {0, 0}; }
S mapping(ll f, S x){
    if(f != INFLL) x.value = f*x.size;
    return x;
}
ll composition(ll f, ll g){ return (f == INFLL ? g : f); }
ll id(){ return INFLL; }

・動的セグメントツリー
template <class S, S (*op)(S, S), S (*e)()> class dynamic_segtree {
    public:
        dynamic_segtree(size_t n) : n(n), root(nullptr) {}
    
        void set(size_t p, S x) {
            assert(p < n);
            set(root, 0, n, p, x);
        }
    
        S get(size_t p) const {
            assert(p < n);
            return get(root, 0, n, p);
        }
    
        S prod(size_t l, size_t r) const {
            assert(l <= r && r <= n);
            return prod(root, 0, n, l, r);
        }
    
        S all_prod() const { return root ? root->product : e(); }
    
        void reset(size_t l, size_t r) {
            assert(l <= r && r <= n);
            return reset(root, 0, n, l, r);
        }
    
        template <bool (*f)(S)> size_t max_right(size_t l) const {
            return max_right(l, [](S x) { return f(x); });
        }
    
        template <class F> size_t max_right(size_t l, const F& f) const {
            assert(l <= n);
            S product = e();
            assert(f(product));
            return max_right(root, 0, n, l, f, product);
        }
    
        template <bool (*f)(S)> size_t min_left(size_t r) const {
            return min_left(r, [](S x) { return f(x); });
        }
    
        template <class F> size_t min_left(size_t r, const F& f) const {
            assert(r <= n);
            S product = e();
            assert(f(product));
            return min_left(root, 0, n, r, f, product);
        }
    
    private:
        struct node;
        using node_ptr = std::unique_ptr<node>;
    
        struct node {
            size_t index;
            S value, product;
            node_ptr left, right;
    
            node(size_t index, S value)
                : index(index),
                  value(value),
                  product(value),
                  left(nullptr),
                  right(nullptr) {}
    
            void update() {
                product = op(op(left ? left->product : e(), value),
                             right ? right->product : e());
            }
        };
    
        const size_t n;
        node_ptr root;
    
        void set(node_ptr& t, size_t a, size_t b, size_t p, S x) const {
            if (!t) {
                t = std::make_unique<node>(p, x);
                return;
            }
            if (t->index == p) {
                t->value = x;
                t->update();
                return;
            }
            size_t c = (a + b) >> 1;
            if (p < c) {
                if (t->index < p) std::swap(t->index, p), std::swap(t->value, x);
                set(t->left, a, c, p, x);
            } else {
                if (p < t->index) std::swap(p, t->index), std::swap(x, t->value);
                set(t->right, c, b, p, x);
            }
            t->update();
        }
    
        S get(const node_ptr& t, size_t a, size_t b, size_t p) const {
            if (!t) return e();
            if (t->index == p) return t->value;
            size_t c = (a + b) >> 1;
            if (p < c) return get(t->left, a, c, p);
            else return get(t->right, c, b, p);
        }
    
        S prod(const node_ptr& t, size_t a, size_t b, size_t l, size_t r) const {
            if (!t || b <= l || r <= a) return e();
            if (l <= a && b <= r) return t->product;
            size_t c = (a + b) >> 1;
            S result = prod(t->left, a, c, l, r);
            if (l <= t->index && t->index < r) result = op(result, t->value);
            return op(result, prod(t->right, c, b, l, r));
        }
    
        void reset(node_ptr& t, size_t a, size_t b, size_t l, size_t r) const {
            if (!t || b <= l || r <= a) return;
            if (l <= a && b <= r) {
                t.reset();
                return;
            }
            size_t c = (a + b) >> 1;
            reset(t->left, a, c, l, r);
            reset(t->right, c, b, l, r);
            t->update();
        }
    
        template <class F>
        size_t max_right(const node_ptr& t, size_t a, size_t b,
                         size_t l, const F& f, S& product) const {
            if (!t || b <= l) return n;
            if (f(op(product, t->product))) {
                product = op(product, t->product);
                return n;
            }
            size_t c = (a + b) >> 1;
            size_t result = max_right(t->left, a, c, l, f, product);
            if (result != n) return result;
            if (l <= t->index) {
                product = op(product, t->value);
                if (!f(product)) return t->index;
            }
            return max_right(t->right, c, b, l, f, product);
        }
    
        template <class F>
        size_t min_left(const node_ptr& t, size_t a, size_t b,
                        size_t r, const F& f, S& product) const {
            if (!t || r <= a) return 0;
            if (f(op(t->product, product))) {
                product = op(t->product, product);
                return 0;
            }
            size_t c = (a + b) >> 1;
            size_t result = min_left(t->right, c, b, r, f, product);
            if (result != 0) return result;
            if (t->index < r) {
                product = op(t->value, product);
                if (!f(product)) return t->index + 1;
            }
            return min_left(t->left, a, c, r, f, product);
        }
};

・IntervalSet
// 区間集合 [l, r] をまとめて管理するクラス
// add/remove は変化した要素数を返す
// add2/remove2 は実際に追加・削除された区間のリストを返す
struct IntervalSet {
    set<pair<long long, long long>> s;

    long long add(long long l, long long r) {
        long long cnt = r - l + 1;
        auto it = s.lower_bound({l, -1});
        if (it != s.begin()) {
            auto pit = prev(it);
            if (pit->second + 1 >= l) {
                l = pit->first;
                r = max(r, pit->second);
                cnt -= min(r, pit->second) - l + 1;
                s.erase(pit);
            }
        }
        while (it != s.end() && it->first <= r + 1) {
            cnt -= min(r, it->second) - it->first + 1;
            r = max(r, it->second);
            it = s.erase(it);
        }
        s.insert({l, r});
        return cnt;
    }

    long long remove(long long l, long long r) {
        long long cnt = 0;
        auto it = s.lower_bound({l, -1});
        if (it != s.begin()) {
            auto pit = prev(it);
            if (pit->second >= l) {
                long long L = pit->first, R = pit->second;
                cnt += min(r, R) - l + 1;
                s.erase(pit);
                if (L <= l - 1) s.insert({L, l - 1});
                if (r < R) s.insert({r + 1, R});
            }
        }
        while (it != s.end() && it->first <= r) {
            long long L = it->first, R = it->second;
            cnt += min(R, r) - L + 1;
            it = s.erase(it);
            if (r < R) {
                s.insert({r + 1, R});
                break;
            }
        }
        return cnt;
    }

    vector<pair<long long,long long>> add2(long long l, long long r) {
        vector<pair<long long,long long>> ans;
        auto it = s.lower_bound({l, -1});
        long long pl = l;

        if (it != s.begin()) {
            auto pit = prev(it);
            if (pit->second + 1 >= l) {
                long long L = pit->first, R = pit->second;
                l = L;
                r = max(r, R);
                pl = R + 1;
                s.erase(pit);
            }
        }
        while (it != s.end() && it->first <= r + 1) {
            long long L = it->first, R = it->second;
            if (pl < L) ans.push_back({pl, L - 1});
            pl = R + 1;
            r = max(r, R);
            it = s.erase(it);
        }
        s.insert({l, r});
        if (pl <= r) ans.push_back({pl, r});
        return ans;
    }

    vector<pair<long long,long long>> remove2(long long l, long long r) {
        vector<pair<long long,long long>> ans;
        auto it = s.lower_bound({l, -1});
        if (it != s.begin()) {
            auto pit = prev(it);
            if (pit->second >= l) {
                long long L = pit->first, R = pit->second;
                s.erase(pit);
                if (L <= l - 1) s.insert({L, l - 1});
                ans.push_back({l, min(R, r)});
                if (r < R) s.insert({r + 1, R});
            }
        }
        while (it != s.end() && it->first <= r) {
            long long L = it->first, R = it->second;
            ans.push_back({L, min(R, r)});
            it = s.erase(it);
            if (r < R) {
                s.insert({r + 1, R});
                break;
            }
        }
        return ans;
    }
};

・UnionFind
// Implement (union by size) + (path compression)
// Reference:
// Zvi Galil and Giuseppe F. Italiano,
// Data structures and algorithms for disjoint set union problems
struct dsu {
  public:
    dsu() : _n(0) {}
    dsu(int n) : _n(n), parent_or_size(n, -1) {}

    int merge(int a, int b) {
        assert(0 <= a && a < _n);
        assert(0 <= b && b < _n);
        int x = leader(a), y = leader(b);
        if (x == y) return x;
        if (-parent_or_size[x] < -parent_or_size[y]) std::swap(x, y);
        parent_or_size[x] += parent_or_size[y];
        parent_or_size[y] = x;
        return x;
    }

    bool same(int a, int b) {
        assert(0 <= a && a < _n);
        assert(0 <= b && b < _n);
        return leader(a) == leader(b);
    }

    int leader(int a) {
        assert(0 <= a && a < _n);
        if (parent_or_size[a] < 0) return a;
        return parent_or_size[a] = leader(parent_or_size[a]);
    }

    int size(int a) {
        assert(0 <= a && a < _n);
        return -parent_or_size[leader(a)];
    }

    std::vector<std::vector<int>> groups() {
        std::vector<int> leader_buf(_n), group_size(_n);
        for (int i = 0; i < _n; i++) {
            leader_buf[i] = leader(i);
            group_size[leader_buf[i]]++;
        }
        std::vector<std::vector<int>> result(_n);
        for (int i = 0; i < _n; i++) {
            result[i].reserve(group_size[i]);
        }
        for (int i = 0; i < _n; i++) {
            result[leader_buf[i]].push_back(i);
        }
        result.erase(
            std::remove_if(result.begin(), result.end(),
                           [&](const std::vector<int>& v) { return v.empty(); }),
            result.end());
        return result;
    }

  private:
    int _n;
    // root node: -1 * component size
    // otherwise: parent
    std::vector<int> parent_or_size;
};

・素数判定
bool isPrime(long long n) {
    if (n <= 1) return false;
    for (long long i = 2; i <= n / i; i++) {
        if (n % i == 0) return false;
    }
    return true;
}

・エラトステネスのふるい
vector<int> eratosthenes(int n) {
    vector<bool> sieve(n + 1, true);

    for (int i = 0; i * 1LL * i <= n; i++) {
        if (i < 2) {
            sieve[i] = false;
        } else if (sieve[i]) {
            for (int j = 2; j <= n / i; j++) {
                sieve[i * j] = false;
            }
        }
    }

    vector<int> ans;
    for (int i = 0; i <= n; i++) {
        if (sieve[i]) ans.push_back(i);
    }
    return ans;
}

・10進数→n進数
string base_n(long long num_10, long long n) {
    if (n <= 1) return "";  // 1進数以下は定義しない前提

    string res = "";
    while (num_10 > 0) {
        long long d = num_10 % n;
        if (d >= 10) return "";
        res += char('0' + d);
        num_10 /= n;
    }

    if (res.empty()) return 0;  // num_10 = 0 の場合

    reverse(res.begin(), res.end());
    return res;  // 数字として返す（先頭ゼロなし）
}

・素因数分解
vector<pair<long long, long long>> factorization(long long n) {
    vector<pair<long long, long long>> res;
    long long temp = n;

    // 上限は ceil(sqrt(n))
    for (long long i = 2; i * i <= temp; i++) {
        if (temp % i == 0) {
            long long cnt = 0;
            while (temp % i == 0) {
                temp /= i;
                cnt++;
            }
            res.push_back({i, cnt});
        }
    }

    if (temp != 1) {
        res.push_back({temp, 1});
    }

    if (res.empty()) {
        res.push_back({n, 1});
    }

    return res;
}

・約数列挙
vector<long long> divisor(long long n) {
    vector<long long> ans;
    for (long long i = 1; i * i <= n; i++) {
        if (n % i == 0) {
            ans.push_back(i);
            if (i * i != n) {
                ans.push_back(n / i);
            }
        }
    }
    return ans;
}

・2点間の距離
double distance(double x1, double y1, double x2, double y2) {
    return sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
}
long long distance2(long long x1, long long y1, long long x2, long long y2) {
    return (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);
}

・トリミング
vector<string> trimming(const vector<string>& A) {
    int H = A.size();
    int W = A[0].size();

    int minH = H, maxH = -1;
    int minW = W, maxW = -1;

    // # の bounding box を探す
    for (int i = 0; i < H; i++) {
        for (int j = 0; j < W; j++) {
            if (A[i][j] == '#') {
                minH = min(minH, i);
                maxH = max(maxH, i);
                minW = min(minW, j);
                maxW = max(maxW, j);
            }
        }
    }

    vector<string> res;

    // # が1つも無い場合 → 空の結果を返す（Python版準拠）
    if (maxH == -1) return res;

    for (int i = minH; i <= maxH; i++) {
        res.push_back(A[i].substr(minW, maxW - minW + 1));
    }

    return res;
}

・回転
template <typename T>
vector<vector<T>> rotateR(const vector<vector<T>>& A) {
    int H = A.size();
    int W = A[0].size();
    vector<vector<T>> res(W, vector<T>(H));
    for (int i = 0; i < H; i++) {
        for (int j = 0; j < W; j++) {
            res[j][H - 1 - i] = A[i][j];  // 時計回り
        }
    }
    return res;
}

template <typename T>
vector<vector<T>> rotateL(const vector<vector<T>>& A) {
    int H = A.size();
    int W = A[0].size();
    vector<vector<T>> res(W, vector<T>(H));
    for (int i = 0; i < H; i++) {
        for (int j = 0; j < W; j++) {
            res[W - 1 - j][i] = A[i][j];  // 反時計回り
        }
    }
    return res;
}

・セグメント木
namespace internal {

// @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;
}

// @param n `1 <= n`
// @return minimum non-negative `x` s.t. `(n & (1 << x)) != 0`
int bsf(unsigned int n) {
#ifdef _MSC_VER
    unsigned long index;
    _BitScanForward(&index, n);
    return index;
#else
    return __builtin_ctz(n);
#endif
}

}  // namespace internal

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) {
        assert(0 <= p && p < _n);
        p += size;
        d[p] = x;
        for (int i = 1; i <= log; i++) update(p >> i);
    }

    S get(int p) {
        assert(0 <= p && p < _n);
        return d[p + size];
    }

    S prod(int l, int 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) {
        return max_right(l, [](S x) { return f(x); });
    }
    template <class F> int max_right(int l, F f) {
        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) {
        return min_left(r, [](S x) { return f(x); });
    }
    template <class F> int min_left(int r, F f) {
        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]); }
};

・遅延セグメント木
namespace internal {

// @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;
}

// @param n `1 <= n`
// @return minimum non-negative `x` s.t. `(n & (1 << x)) != 0`
int bsf(unsigned int n) {
#ifdef _MSC_VER
    unsigned long index;
    _BitScanForward(&index, n);
    return index;
#else
    return __builtin_ctz(n);
#endif
}

}  // namespace internal
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();
    }
};

・強連結成分分解
namespace internal {

template <class E> struct csr {
    std::vector<int> start;
    std::vector<E> elist;
    csr(int n, const std::vector<std::pair<int, E>>& edges)
        : start(n + 1), elist(edges.size()) {
        for (auto e : edges) {
            start[e.first + 1]++;
        }
        for (int i = 1; i <= n; i++) {
            start[i] += start[i - 1];
        }
        auto counter = start;
        for (auto e : edges) {
            elist[counter[e.first]++] = e.second;
        }
    }
};

// Reference:
// R. Tarjan,
// Depth-First Search and Linear Graph Algorithms
struct scc_graph {
  public:
    scc_graph(int n) : _n(n) {}

    int num_vertices() { return _n; }

    void add_edge(int from, int to) { edges.push_back({from, {to}}); }

    // @return pair of (# of scc, scc id)
    std::pair<int, std::vector<int>> scc_ids() {
        auto g = csr<edge>(_n, edges);
        int now_ord = 0, group_num = 0;
        std::vector<int> visited, low(_n), ord(_n, -1), ids(_n);
        visited.reserve(_n);
        auto dfs = [&](auto self, int v) -> void {
            low[v] = ord[v] = now_ord++;
            visited.push_back(v);
            for (int 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) {
                    int u = visited.back();
                    visited.pop_back();
                    ord[u] = _n;
                    ids[u] = group_num;
                    if (u == v) break;
                }
                group_num++;
            }
        };
        for (int 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};
    }

    std::vector<std::vector<int>> scc() {
        auto ids = scc_ids();
        int group_num = ids.first;
        std::vector<int> counts(group_num);
        for (auto x : ids.second) counts[x]++;
        std::vector<std::vector<int>> groups(ids.first);
        for (int i = 0; i < group_num; i++) {
            groups[i].reserve(counts[i]);
        }
        for (int i = 0; i < _n; i++) {
            groups[ids.second[i]].push_back(i);
        }
        return groups;
    }

  private:
    int _n;
    struct edge {
        int to;
    };
    std::vector<std::pair<int, edge>> edges;
};

}  // namespace internal

struct scc_graph {
  public:
    scc_graph() : internal(0) {}
    scc_graph(int n) : internal(n) {}

    void add_edge(int from, int to) {
        int n = internal.num_vertices();
        assert(0 <= from && from < n);
        assert(0 <= to && to < n);
        internal.add_edge(from, to);
    }

    std::vector<std::vector<int>> scc() { return internal.scc(); }

  private:
    internal::scc_graph internal;
};

・繰り返し二乗法
long long pow(long long x, long long n) {
    long long ret = 1;
    while (n > 0) {
        if (n & 1) ret = ret * x % MOD;  // n の最下位bitが 1 ならば x^(2^i) をかける
        x = x * x % MOD;
        n >>= 1;  // n を1bit 左にずらす
    }
    return ret;
}

・拡張ユークリッドの互除法
pair<long long, long long> extended_euclid(long long a, long long b) {
    if (b == 0) {
        return {1, 0};
    } else {
        auto p = extended_euclid(b, a % b);
        long long xd = p.first;
        long long yd = p.second;
        return {yd, xd - (a / b) * yd};
    }
}

・4方向の差分
vi dh = {-1, 0, 1, 0};
vi dw = {0, 1, 0, -1};

・範囲内の判定
0<=nh && nh<H && 0<=nw && nw<W

・最大流
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++; }
};

}  // namespace internal

template <class Cap> struct mf_graph {
  public:
    mf_graph() : _n(0) {}
    mf_graph(int n) : _n(n), g(n) {}

    int add_edge(int from, int to, Cap cap) {
        assert(0 <= from && from < _n);
        assert(0 <= to && to < _n);
        assert(0 <= cap);
        int m = int(pos.size());
        pos.push_back({from, int(g[from].size())});
        int from_id = int(g[from].size());
        int to_id = int(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 {
        int from, to;
        Cap cap, flow;
    };

    edge get_edge(int i) {
        int m = int(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() {
        int m = int(pos.size());
        std::vector<edge> result;
        for (int i = 0; i < m; i++) {
            result.push_back(get_edge(i));
        }
        return result;
    }
    void change_edge(int i, Cap new_cap, Cap new_flow) {
        int m = int(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(int s, int t) {
        return flow(s, t, std::numeric_limits<Cap>::max());
    }
    Cap flow(int s, int t, Cap flow_limit) {
        assert(0 <= s && s < _n);
        assert(0 <= t && t < _n);
        assert(s != t);

        std::vector<int> level(_n), iter(_n);
        internal::simple_queue<int> que;

        auto bfs = [&]() {
            std::fill(level.begin(), level.end(), -1);
            level[s] = 0;
            que.clear();
            que.push(s);
            while (!que.empty()) {
                int 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, int v, Cap up) {
            if (v == s) return up;
            Cap res = 0;
            int level_v = level[v];
            for (int& i = iter[v]; i < int(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(int s) {
        std::vector<bool> visited(_n);
        internal::simple_queue<int> que;
        que.push(s);
        while (!que.empty()) {
            int 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:
    int _n;
    struct _edge {
        int to, rev;
        Cap cap;
    };
    std::vector<std::pair<int, int>> pos;
    std::vector<std::vector<_edge>> g;
};

・最小費用流
template <class Cap, class Cost> struct mcf_graph {
  public:
    mcf_graph() {}
    mcf_graph(int n) : _n(n), g(n) {}

    int add_edge(int from, int to, Cap cap, Cost cost) {
        assert(0 <= from && from < _n);
        assert(0 <= to && to < _n);
        int m = int(pos.size());
        pos.push_back({from, int(g[from].size())});
        int from_id = int(g[from].size());
        int to_id = int(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 {
        int from, to;
        Cap cap, flow;
        Cost cost;
    };

    edge get_edge(int i) {
        int m = int(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,
        };
    }
    std::vector<edge> edges() {
        int m = int(pos.size());
        std::vector<edge> result(m);
        for (int i = 0; i < m; i++) {
            result[i] = get_edge(i);
        }
        return result;
    }

    std::pair<Cap, Cost> flow(int s, int t) {
        return flow(s, t, std::numeric_limits<Cap>::max());
    }
    std::pair<Cap, Cost> flow(int s, int t, Cap flow_limit) {
        return slope(s, t, flow_limit).back();
    }
    std::vector<std::pair<Cap, Cost>> slope(int s, int t) {
        return slope(s, t, std::numeric_limits<Cap>::max());
    }
    std::vector<std::pair<Cap, Cost>> slope(int s, int t, Cap flow_limit) {
        assert(0 <= s && s < _n);
        assert(0 <= t && t < _n);
        assert(s != t);
        // variants (C = maxcost):
        // -(n-1)C <= dual[s] <= dual[i] <= dual[t] = 0
        // reduced cost (= e.cost + dual[e.from] - dual[e.to]) >= 0 for all edge
        std::vector<Cost> dual(_n, 0), dist(_n);
        std::vector<int> pv(_n), pe(_n);
        std::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;
                int 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()) {
                int v = que.top().to;
                que.pop();
                if (vis[v]) continue;
                vis[v] = true;
                if (v == t) break;
                // dist[v] = shortest(s, v) + dual[s] - dual[v]
                // dist[v] >= 0 (all reduced cost are positive)
                // dist[v] <= (n-1)C
                for (int i = 0; i < int(g[v].size()); i++) {
                    auto e = g[v][i];
                    if (vis[e.to] || !e.cap) continue;
                    // |-dual[e.to] + dual[v]| <= (n-1)C
                    // cost <= C - -(n-1)C + 0 = nC
                    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;
            }

            for (int v = 0; v < _n; v++) {
                if (!vis[v]) continue;
                // dual[v] = dual[v] - dist[t] + dist[v]
                //         = dual[v] - (shortest(s, t) + dual[s] - dual[t]) + (shortest(s, v) + dual[s] - dual[v])
                //         = - shortest(s, t) + dual[t] + shortest(s, v)
                //         = shortest(s, v) - shortest(s, t) >= 0 - (n-1)C
                dual[v] -= dist[t] - dist[v];
            }
            return true;
        };
        Cap flow = 0;
        Cost cost = 0, prev_cost_per_flow = -1;
        std::vector<std::pair<Cap, Cost>> result;
        result.push_back({flow, cost});
        while (flow < flow_limit) {
            if (!dual_ref()) break;
            Cap c = flow_limit - flow;
            for (int v = t; v != s; v = pv[v]) {
                c = std::min(c, g[pv[v]][pe[v]].cap);
            }
            for (int 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:
    int _n;

    struct _edge {
        int to, rev;
        Cap cap;
        Cost cost;
    };

    std::vector<std::pair<int, int>> pos;
    std::vector<std::vector<_edge>> g;
};

・BIT
namespace internal {

#ifndef _MSC_VER
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;

#else

template <class T> using is_integral = typename std::is_integral<T>;

template <class T>
using is_signed_int =
    typename std::conditional<is_integral<T>::value && std::is_signed<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,
                              std::true_type,
                              std::false_type>::type;

template <class T>
using to_unsigned = typename std::conditional<is_signed_int<T>::value,
                                              std::make_unsigned<T>,
                                              std::common_type<T>>::type;

#endif

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;

}  // namespace internal

// Reference: https://en.wikipedia.org/wiki/Fenwick_tree
template <class T> struct fenwick_tree {
    using U = internal::to_unsigned_t<T>;

  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] += U(x);
            p += p & -p;
        }
    }

    T sum(int l, int r) {
        assert(0 <= l && l <= r && r <= _n);
        return sum(r) - sum(l);
    }

  private:
    int _n;
    std::vector<U> data;

    U sum(int r) {
        U s = 0;
        while (r > 0) {
            s += data[r - 1];
            r -= r & -r;
        }
        return s;
    }
};

・UnionFind（値を保持）
template <typename T>
struct UnionFindValue {
    int n;
    vector<int> parent_size;  // 負ならサイズ、正なら親
    vector<T> A;              // 各リーダーが持つ値

    UnionFindValue(int n, const vector<T>& A_) 
        : n(n), parent_size(n, -1), A(A_) {}

    int leader(int a) {
        if (parent_size[a] < 0) return a;
        return parent_size[a] = leader(parent_size[a]);
    }

    void merge(int a, int b) {
        a = leader(a);
        b = leader(b);
        if (a == b) return;

        T c = A[a];
        
        // union by size
        if (-parent_size[a] < -parent_size[b]) swap(a, b);

        // a が大きい側
        parent_size[a] += parent_size[b];
        parent_size[b] = a;

        A[a] = c;
    }

    // A[x] を取得
    T operator[](int x) {
        return A[leader(x)];
    }

    // 値を上書き
    void update(int x, T v) {
        A[leader(x)] = v;
    }

    bool same(int a, int b) {
        return leader(a) == leader(b);
    }

    int size(int x) {
        return -parent_size[leader(x)];
    }

    vector<vector<int>> groups() {
        vector<vector<int>> res(n);
        for (int i = 0; i < n; i++) {
            res[leader(i)].push_back(i);
        }
        vector<vector<int>> ans;
        for (auto &g : res) {
            if (!g.empty()) ans.push_back(g);
        }
        return ans;
    }
};

・UnionFind（値を演算）
template <typename T>
struct UnionFindOp {
    int n;
    vector<int> parent_size;  
    vector<T> A;              
    function<T(T,T)> op;      

    UnionFindOp(int n, const vector<T>& A_, function<T(T,T)> op_) 
        : n(n), parent_size(n, -1), A(A_), op(op_) {}

    int leader(int a) {
        if (parent_size[a] < 0) return a;
        return parent_size[a] = leader(parent_size[a]);
    }

    void merge(int a, int b) {
        a = leader(a);
        b = leader(b);
        if (a == b) return;

        T l = A[a];
        T r = A[b];

        // union by size
        if (-parent_size[a] < -parent_size[b]) swap(a, b);

        parent_size[a] += parent_size[b];
        parent_size[b] = a;

        // 値を演算してマージ
        A[a] = op(l, r);
    }

    T operator[](int x) {
        return A[leader(x)];
    }

    void update(int x, T v) {
        A[leader(x)] = v;
    }

    bool same(int a, int b) {
        return leader(a) == leader(b);
    }

    int size(int x) {
        return -parent_size[leader(x)];
    }

    vector<vector<int>> groups() {
        vector<vector<int>> res(n);
        for (int i = 0; i < n; i++) {
            res[leader(i)].push_back(i);
        }
        vector<vector<int>> ans;
        for (auto &g : res) {
            if (!g.empty()) ans.push_back(g);
        }
        return ans;
    }
};

・セグメント木（区間演算、一点取得）
template <typename T>
struct SegTreeApplyRange {
    int n;
    int num;
    vector<T> tree;
    function<T(T, T)> op;
    T ide;

    SegTreeApplyRange(const vector<T>& init_val, function<T(T,T)> op_, T ide_)
        : n(init_val.size()), op(op_), ide(ide_) {
        
        num = 1;
        while (num < n) num <<= 1;

        tree.assign(2 * num, ide);

        for (int i = 0; i < n; i++) {
            tree[num + i] = init_val[i];
        }
    }

    // apply op(x, tree[i]) on range [l, r)
    void update(int l, int r, T x) {
        l += num;
        r += num;
        while (l < r) {
            if (l & 1) {
                tree[l] = op(x, tree[l]);
                l++;
            }
            if (r & 1) {
                r--;
                tree[r] = op(x, tree[r]);
            }
            l >>= 1;
            r >>= 1;
        }
    }

    // get point k (propagating all ops from root to leaf)
    T query(int k) {
        k += num;
        T ans = op(ide, tree[k]);
        while (k > 1) {
            k >>= 1;
            ans = op(ans, tree[k]);
        }
        return ans;
    }
};

・セグメント木（区間更新、一点取得）
template <typename T>
struct SegTreeRangeAssign {
    int n;
    int num;
    vector<T> tree;
    vector<int> rank;
    int cnt;

    SegTreeRangeAssign(const vector<T>& init_val)
        : n(init_val.size()), cnt(0)
    {
        num = 1;
        while (num < n) num <<= 1;

        tree.assign(2 * num, T());
        rank.assign(2 * num, -1);

        for (int i = 0; i < n; i++) {
            tree[num + i] = init_val[i];
            rank[num + i] = cnt;
        }
        cnt++;
    }

    // assign x to range [l, r)
    void update(int l, int r, T x) {
        l += num;
        r += num;
        while (l < r) {
            if (l & 1) {
                tree[l] = x;
                rank[l] = cnt;
                l++;
            }
            if (r & 1) {
                r--;
                tree[r] = x;
                rank[r] = cnt;
            }
            l >>= 1;
            r >>= 1;
        }
        cnt++;
    }

    // get point k (latest assignment on path)
    T query(int k) {
        k += num;
        T ans = tree[k];
        int mx = rank[k];

        while (k > 1) {
            k >>= 1;
            if (mx < rank[k]) {
                ans = tree[k];
                mx = rank[k];
            }
        }
        return ans;
    }
};

・二項係数（手動MOD）
struct Comb {
    long long MOD;
    int N;
    vector<long long> fact, fact_inv;

    Comb(int N, long long MOD) : N(N), MOD(MOD), fact(N+1), fact_inv(N+1) {
        fact[0] = 1;
        for (int i = 1; i <= N; i++) fact[i] = fact[i-1] * i % MOD;

        // fermat: inv = fact[N]^(MOD-2)
        auto modpow = [&](long long a, long long e) {
            long long r = 1;
            while (e) {
                if (e & 1) r = r * a % MOD;
                a = a * a % MOD;
                e >>= 1;
            }
            return r;
        };

        fact_inv[N] = modpow(fact[N], MOD-2);
        for (int i = N; i >= 1; i--) {
            fact_inv[i-1] = fact_inv[i] * i % MOD;
        }
    }

    long long C(long long n, long long k) const {
        if (n < 0 || k < 0 || n < k) return 0;
        return fact[n] * fact_inv[k] % MOD * fact_inv[n-k] % MOD;
    }

    long long P(long long n, long long k) const {
        if (n < 0 || k < 0 || n < k) return 0;
        return fact[n] * fact_inv[n-k] % MOD;
    }

    long long H(long long n, long long k) const {
        if (n < 0 || k < 0) return 0;
        if (n == 0 && k == 0) return 1;
        return C(n + k - 1, k);
    }
};

・二項係数（mint）
template <class mint>
struct CombMint {
    int N;
    vector<mint> fact, fact_inv;

    CombMint(int N) : N(N), fact(N+1), fact_inv(N+1) {
        fact[0] = 1;
        for (int i = 1; i <= N; i++) fact[i] = fact[i-1] * i;

        fact_inv[N] = fact[N].inv();
        for (int i = N; i >= 1; i--) {
            fact_inv[i-1] = fact_inv[i] * i;
        }
    }

    mint C(long long n, long long k) const {
        if (n < 0 || k < 0 || n < k) return mint(0);
        return fact[n] * fact_inv[k] * fact_inv[n-k];
    }

    mint P(long long n, long long k) const {
        if (n < 0 || k < 0 || n < k) return mint(0);
        return fact[n] * fact_inv[n-k];
    }

    mint H(long long n, long long k) const {
        if (n < 0 || k < 0) return mint(0);
        if (n == 0 && k == 0) return mint(1);
        return C(n + k - 1, k);
    }
};

・全方位木DP
vector<vector<int>> G;
vector<vector<long long>> dp, L, R;
vector<bool> visited;

void LRset(int n){
    int m = G[n].size();
    for(int i = 1; i < m; i++){
        L[n][i] = L[n][i-1] + dp[n][i-1];
    }
    for(int i = m-2; i >= 0; i--){
        R[n][i] = R[n][i+1] + dp[n][i+1];
    }
}

void rerooting(int n, int pre, int no){
    visited[n] = true;

    int preno = 0;
    for(int i = 0; i < (int)G[n].size(); i++){
        if(G[n][i] == pre){
            preno = i;
            break;
        }
    }

    dp[n][preno] = L[pre][no] + R[pre][no] + 1;
    LRset(n);

    for(int i = 0; i < (int)G[n].size(); i++){
        int nxt = G[n][i];
        if(!visited[nxt]){
            rerooting(nxt, n, i);
        }
    }
}

dp.assign(N, {});
L.assign(N, {});
R.assign(N, {});
for(int i = 0; i < N; i++){
    int m = G[i].size();
    dp[i].assign(m, 0LL);
    L[i].assign(m, 0LL);
    R[i].assign(m, 0LL);
}

// root = 0 の初期 L,R
LRset(0);

visited.assign(N, false);
visited[0] = true;

for(int i = 0; i < (int)G[0].size(); i++){
    int nxt = G[0][i];
    if(!visited[nxt]){
        rerooting(nxt, 0, i);
    }
}

・Z-algorithm
vector<int> z_algo(const string &S) {
    int N = (int)S.size();
    vector<int> A(N);
    A[0] = N;

    int i = 1, j = 0;
    while (i < N) {
        while (i + j < N && S[j] == S[i + j]) j++;
        if (j == 0) {
            i++;
            continue;
        }
        A[i] = j;
        int k = 1;
        while (i + k < N && k < j && A[k] < j - k) {
            A[i + k] = A[k];
            k++;
        }
        i += k;
        j -= k;
    }
    return A;
}

・RollingHash
struct RollingHash {
    vector<long long> hash;
    vector<long long> powv;

    RollingHash(const string &S, bool type) {
        hash.assign(S.size() + 1, 0);
        powv.assign(S.size() + 1, 1);

        int a = (type == false ? 'a' - 1 : 'A' - 1);

        for (int i = 0; i < (int)S.size(); i++) {
            hash[i + 1] = (hash[i] * 26 + (S[i] - a)) % MOD;
            powv[i + 1] = (powv[i] * 26) % MOD;
        }
    }

    // [l, r) のハッシュ
    long long query(int l, int r) const {
        long long res = hash[r] - (hash[l] * powv[r - l]) % MOD;
        if (res < 0) res += MOD;  // マイナス回避
        return res;
    }
};

・最小共通祖先
struct LCA {
    int V;
    int LOG;
    vector<vector<int>> parent;
    vector<int> depth;

    LCA(const vector<vector<int>>& G, int root = 0) {
        V = (int)G.size();
        LOG = 1;
        while ((1 << LOG) < V) LOG++;

        parent.assign(LOG, vector<int>(V, -1));
        depth.assign(V, -1);

        bfs(G, root);

        for (int i = 0; i + 1 < LOG; i++) {
            for (int v = 0; v < V; v++) {
                if (parent[i][v] != -1) {
                    parent[i + 1][v] = parent[i][ parent[i][v] ];
                }
            }
        }
    }

    void bfs(const vector<vector<int>>& G, int root) {
        queue<int> q;
        depth[root] = 0;
        q.push(root);

        while (!q.empty()) {
            int u = q.front(); q.pop();
            for (int v : G[u]) {
                if (depth[v] == -1) {
                    depth[v] = depth[u] + 1;
                    parent[0][v] = u;
                    q.push(v);
                }
            }
        }
    }

    int lca(int a, int b) const {
        if (depth[a] < depth[b]) swap(a, b);

        int diff = depth[a] - depth[b];
        for (int i = 0; i < LOG; i++) {
            if (diff & (1 << i)) {
                a = parent[i][a];
            }
        }

        if (a == b) return a;

        for (int i = LOG - 1; i >= 0; i--) {
            if (parent[i][a] != parent[i][b]) {
                a = parent[i][a];
                b = parent[i][b];
            }
        }
        return parent[0][a];
    }

    int dist(int a, int b) const {
        int c = lca(a, b);
        return depth[a] + depth[b] - 2 * depth[c];
    }

    bool is_ancestor(int u, int v) const {
        return lca(u, v) == u;
    }

    int kth_ancestor(int v, int k) const {
        for (int i = 0; i < LOG; i++) {
            if (k & (1 << i)) {
                v = parent[i][v];
                if (v == -1) break;
            }
        }
        return v;
    }

    int jump(int u, int v) const {
        if (u == v) return u;
        int c = lca(u, v);
        if (c == u) {
            return kth_ancestor(v, dist(u, v) - 1);
        } else {
            return parent[0][u];
        }
    }

    int jump_k(int u, int v, int k) const {
        int d = dist(u, v);
        if (d <= k) return v;

        int c = lca(u, v);
        int du = depth[u] - depth[c];

        if (k <= du) {
            return kth_ancestor(u, k);
        } else {
            return kth_ancestor(v, d - k);
        }
    }

    bool on_path(int u, int v, int x) const {
        return dist(u, x) + dist(x, v) == dist(u, v);
    }

    pair<int,int> create_path(int u, int v, int x) const {
        if (on_path(u, v, x)) return {u, v};
        if (on_path(u, x, v)) return {u, x};
        if (on_path(v, x, u)) return {v, x};
        return {-1, -1};
    }
};

・行列の積
vector<vector<mint>> matrix_mul(const vector<vector<mint>> &a,
                                const vector<vector<mint>> &b) {
    int H1 = a.size();
    int W1 = a[0].size();
    int W2 = b[0].size();

    vector<vector<mint>> ans(H1, vector<mint>(W2, 0));

    for (int i = 0; i < H1; i++) {
        for (int k = 0; k < W1; k++) {
            mint aik = a[i][k];
            for (int j = 0; j < W2; j++) {
                ans[i][j] += aik * b[k][j];
            }
        }
    }
    return ans;
}

・エラトステネスのふるいを用いた高速素因数分解と高速約数列挙
struct Eratosthenes {
    int n;
    vector<bool> isPrime;
    vector<int> minfactor;

    Eratosthenes(int n) : n(n) {
        isPrime.assign(n + 1, true);
        minfactor.assign(n + 1, -1);

        if (n >= 0) isPrime[0] = false;
        if (n >= 1) {
            isPrime[1] = false;
            minfactor[1] = 1;
        }

        for (int i = 2; i <= n; i++) {
            if (isPrime[i]) {
                minfactor[i] = i;
                for (long long j = 1LL * i * 2; j <= n; j += i) {
                    isPrime[j] = false;
                    if (minfactor[j] == -1) minfactor[j] = i;
                }
            }
        }
    }

    vector<pair<int,int>> factorize(long long x) const {
        vector<pair<int,int>> res;
        while (x > 1) {
            int p = minfactor[x];
            int cnt = 0;
            while (minfactor[x] == p) {
                x /= p;
                cnt++;
            }
            res.emplace_back(p, cnt);
        }
        return res;
    }

    vector<long long> divisor(long long x) const {
        vector<long long> ans = {1};
        auto pf = factorize(x);

        for (auto [p, c] : pf) {
            int L = ans.size();
            for (int i = 0; i < L; i++) {
                long long v = 1;
                for (int t = 0; t < c; t++) {
                    v *= p;
                    ans.push_back(ans[i] * v);
                }
            }
        }
        return ans;
    }
};

・高速ゼータ変換、高速メビウス変換
vector<bool> Eratosthenes(int n) {
    vector<bool> isPrime(n + 1, true);
    if (n >= 0) isPrime[0] = false;
    if (n >= 1) isPrime[1] = false;
    for (int i = 2; i <= n; i++) {
        if (isPrime[i]) {
            for (long long j = 1LL * i * 2; j <= n; j += i) {
                isPrime[j] = false;
            }
        }
    }
    return isPrime;
}

void fast_zeta(vector<long long> &F) {
    int N = F.size();
    vector<bool> isPrime = Eratosthenes(N - 1);

    for (int i = 2; i < N; i++) {
        if (!isPrime[i]) continue;
        for (int j = (N - 1) / i; j >= 1; j--) {
            F[j] += F[j * i];
        }
    }
}

void fast_mobius(vector<long long> &F) {
    int N = F.size();
    vector<bool> isPrime = Eratosthenes(N - 1);

    for (int i = 2; i < N; i++) {
        if (!isPrime[i]) continue;
        for (int j = 1; j <= (N - 1) / i; j++) {
            F[j] -= F[j * i];
        }
    }
}

・メビウス関数値の列挙
struct Eratosthenes {
    int n;
    vector<bool> isPrime;
    vector<int> mobius;

    Eratosthenes(int n) : n(n) {
        isPrime.assign(n + 1, true);
        mobius.assign(n + 1, 1);

        if (n >= 0) isPrime[0] = false;
        if (n >= 1) isPrime[1] = false;

        for (int i = 2; i <= n; i++) {
            if (isPrime[i]) {
                mobius[i] = -1;
                for (int j = 2 * i; j <= n; j += i) {
                    isPrime[j] = false;
                    if ((j / i) % i == 0) {
                        mobius[j] = 0;
                    } else {
                        mobius[j] *= -1;
                    }
                }
            }
        }
    }
};

・マージソートツリー
struct MergeSortTree {
    int n, num;
    vector<vector<long long>> tree;
    vector<vector<long long>> cum;

    MergeSortTree(const vector<long long> &init_val) {
        n = init_val.size();
        num = 1;
        while (num < n) num <<= 1;
        tree.assign(num * 2, {});
        cum.assign(num * 2, {0});

        // 葉に初期値
        for (int i = 0; i < n; i++) {
            tree[num + i].push_back(init_val[i]);
        }

        // 内部ノードのマージ
        for (int i = num - 1; i > 0; i--) {
            int L = 2 * i, R = 2 * i + 1;
            tree[i].resize(tree[L].size() + tree[R].size());
            merge(tree[L].begin(), tree[L].end(), tree[R].begin(), tree[R].end(), tree[i].begin());
        }

        // 累積和
        for (int i = 1; i < num * 2; i++) {
            cum[i].resize(tree[i].size() + 1);
            for (size_t j = 0; j < tree[i].size(); j++) {
                cum[i][j + 1] = cum[i][j] + tree[i][j];
            }
        }
    }

    // A[i] <= X の個数
    long long query1(int l, int r, long long X) {
        long long ans = 0;
        l += num; r += num;
        while (l < r) {
            if (l & 1) ans += upper_bound(tree[l].begin(), tree[l].end(), X) - tree[l].begin(), l++;
            if (r & 1) r--, ans += upper_bound(tree[r].begin(), tree[r].end(), X) - tree[r].begin();
            l >>= 1; r >>= 1;
        }
        return ans;
    }

    // A[i] >= X の個数
    long long query2(int l, int r, long long X) {
        long long ans = 0;
        l += num; r += num;
        while (l < r) {
            if (l & 1) ans += tree[l].size() - (lower_bound(tree[l].begin(), tree[l].end(), X) - tree[l].begin()), l++;
            if (r & 1) r--, ans += tree[r].size() - (lower_bound(tree[r].begin(), tree[r].end(), X) - tree[r].begin());
            l >>= 1; r >>= 1;
        }
        return ans;
    }

    // leftX <= A[i] <= rightX の個数
    long long query3(int l, int r, long long leftX, long long rightX) {
        long long ans = 0;
        l += num; r += num;
        while (l < r) {
            if (l & 1) {
                int left = lower_bound(tree[l].begin(), tree[l].end(), leftX) - tree[l].begin();
                int right = upper_bound(tree[l].begin(), tree[l].end(), rightX) - tree[l].begin();
                ans += right - left;
                l++;
            }
            if (r & 1) {
                r--;
                int left = lower_bound(tree[r].begin(), tree[r].end(), leftX) - tree[r].begin();
                int right = upper_bound(tree[r].begin(), tree[r].end(), rightX) - tree[r].begin();
                ans += right - left;
            }
            l >>= 1; r >>= 1;
        }
        return ans;
    }

    // A[i] <= X の総和
    long long query4(int l, int r, long long X) {
        long long ans = 0;
        l += num; r += num;
        while (l < r) {
            if (l & 1) {
                int b = upper_bound(tree[l].begin(), tree[l].end(), X) - tree[l].begin();
                ans += cum[l][b];
                l++;
            }
            if (r & 1) {
                r--;
                int b = upper_bound(tree[r].begin(), tree[r].end(), X) - tree[r].begin();
                ans += cum[r][b];
            }
            l >>= 1; r >>= 1;
        }
        return ans;
    }

    // A[i] >= X の総和
    long long query5(int l, int r, long long X) {
        long long ans = 0;
        l += num; r += num;
        while (l < r) {
            if (l & 1) {
                int b = lower_bound(tree[l].begin(), tree[l].end(), X) - tree[l].begin();
                ans += cum[l].back() - cum[l][b];
                l++;
            }
            if (r & 1) {
                r--;
                int b = lower_bound(tree[r].begin(), tree[r].end(), X) - tree[r].begin();
                ans += cum[r].back() - cum[r][b];
            }
            l >>= 1; r >>= 1;
        }
        return ans;
    }

    // leftX <= A[i] <= rightX の総和
    long long query6(int l, int r, long long leftX, long long rightX) {
        long long ans = 0;
        l += num; r += num;
        while (l < r) {
            if (l & 1) {
                int left = lower_bound(tree[l].begin(), tree[l].end(), leftX) - tree[l].begin();
                int right = upper_bound(tree[l].begin(), tree[l].end(), rightX) - tree[l].begin();
                ans += cum[l][right] - cum[l][left];
                l++;
            }
            if (r & 1) {
                r--;
                int left = lower_bound(tree[r].begin(), tree[r].end(), leftX) - tree[r].begin();
                int right = upper_bound(tree[r].begin(), tree[r].end(), rightX) - tree[r].begin();
                ans += cum[r][right] - cum[r][left];
            }
            l >>= 1; r >>= 1;
        }
        return ans;
    }
};

・LowLink
struct LowLink {
    int N, M;
    vector<vector<pair<int,int>>> G;
    vector<int> ord, low;
    vector<bool> visited;
    int cnt;
    vector<int> articulation_points;
    vector<int> bridge;

    LowLink(int N) : N(N), M(0), G(N), cnt(0) {}

    void add_edge(int u, int v) {
        G[u].emplace_back(v, M);
        G[v].emplace_back(u, M);
        M++;
    }

    void dfs(int n, int parent) {
        visited[n] = true;
        ord[n] = low[n] = cnt++;
        bool is_ap = false;
        int child = 0;

        for (auto [v, idx] : G[n]) {
            if (!visited[v]) {
                child++;
                dfs(v, idx);
                low[n] = min(low[n], low[v]);
                if (parent != -1 && ord[n] <= low[v]) is_ap = true;
                if (ord[n] < low[v]) bridge.push_back(idx);
            } else if (idx != parent) {
                low[n] = min(low[n], ord[v]);
            }
        }

        if (parent == -1 && child >= 2) is_ap = true;
        if (is_ap) articulation_points.push_back(n);
    }

    pair<vector<int>, vector<int>> find() {
        ord.assign(N, -1);
        low.assign(N, -1);
        visited.assign(N, false);
        articulation_points.clear();
        bridge.clear();
        cnt = 0;

        for (int i = 0; i < N; i++) {
            if (!visited[i]) dfs(i, -1);
        }

        sort(articulation_points.begin(), articulation_points.end());
        sort(bridge.begin(), bridge.end());
        return {articulation_points, bridge};
    }
};

・Trie木
struct Trie {
    int N, base;
    int nex;
    vector<vector<int>> tree;
    vector<int> arrive, last;

    // type=0 -> 'a' 基準, type=1 -> 'A' 基準
    Trie(int N, int type) : N(N), tree(N+1, vector<int>(26, -1)), 
                            arrive(N+1, 0), last(N+1, 0), nex(1) {
        base = (type == 0 ? 'a' : 'A');
    }

    void insert(const string &S) {
        int now = 0;
        arrive[now]++;
        for (char s : S) {
            int c = s - base;
            if (tree[now][c] == -1) {
                tree[now][c] = nex;
                now = nex++;
            } else {
                now = tree[now][c];
            }
            arrive[now]++;
        }
        last[now]++;
    }

    int search(const string &S) const {
        int now = 0;
        for (char s : S) {
            int c = s - base;
            if (tree[now][c] == -1) return 0;
            now = tree[now][c];
        }
        return last[now];
    }

    int start_with(const string &S) const {
        int now = 0;
        for (char s : S) {
            int c = s - base;
            if (tree[now][c] == -1) return 0;
            now = tree[now][c];
        }
        return arrive[now];
    }

    int count() const {
        return arrive[0];
    }
};

・座標圧縮
template<typename T>
vector<T> compress(vector<T> A) {
    sort(A.begin(), A.end());
    A.erase(unique(A.begin(), A.end()), A.end());
    return A;
}

・最小共通祖先（重み付き）
struct LCA {
    int V;
    int LOG;
    vector<vector<int>> parent;
    vector<int> depth;
    vector<long long> depth2; // 重み付き距離

    LCA(const vector<vector<pair<int, long long>>>& G, int root = 0) {
        V = G.size();
        LOG = 1;
        while ((1 << LOG) < V) LOG++;

        parent.assign(LOG, vector<int>(V, -1));
        depth.assign(V, -1);
        depth2.assign(V, 0);

        bfs(G, root);

        for (int i = 0; i + 1 < LOG; i++) {
            for (int v = 0; v < V; v++) {
                if (parent[i][v] != -1) {
                    parent[i + 1][v] = parent[i][parent[i][v]];
                }
            }
        }
    }

    void bfs(const vector<vector<pair<int, long long>>>& G, int root) {
        queue<int> q;
        depth[root] = 0;
        depth2[root] = 0;
        q.push(root);

        while (!q.empty()) {
            int u = q.front(); q.pop();
            for (auto [v, w] : G[u]) {
                if (depth[v] == -1) {
                    depth[v] = depth[u] + 1;
                    depth2[v] = depth2[u] + w;
                    parent[0][v] = u;
                    q.push(v);
                }
            }
        }
    }

    int lca(int a, int b) {
        if (depth[a] < depth[b]) swap(a, b);

        int diff = depth[a] - depth[b];
        for (int i = 0; i < LOG; i++) {
            if (diff & (1 << i)) {
                a = parent[i][a];
            }
        }

        if (a == b) return a;

        for (int i = LOG - 1; i >= 0; i--) {
            if (parent[i][a] != parent[i][b]) {
                a = parent[i][a];
                b = parent[i][b];
            }
        }
        return parent[0][a];
    }

    int dist(int a, int b) {
        int c = lca(a, b);
        return depth[a] + depth[b] - 2 * depth[c];
    }

    long long dist2(int a, int b) {
        int c = lca(a, b);
        return depth2[a] + depth2[b] - 2 * depth2[c];
    }

    bool is_ancestor(int u, int v) {
        return lca(u, v) == u;
    }

    int kth_ancestor(int v, int k) {
        for (int i = 0; i < LOG; i++) {
            if (k & (1 << i)) {
                v = parent[i][v];
                if (v == -1) break;
            }
        }
        return v;
    }

    int jump(int u, int v) {
        if (u == v) return u;
        int c = lca(u, v);
        if (c == u) {
            return kth_ancestor(v, dist(u, v) - 1);
        } else {
            return parent[0][u];
        }
    }

    int jump_k(int u, int v, int k) {
        int d = dist(u, v);
        if (d <= k) return v;

        int c = lca(u, v);
        int du = depth[u] - depth[c];

        if (k <= du) {
            return kth_ancestor(u, k);
        } else {
            return kth_ancestor(v, d - k);
        }
    }

    bool on_path(int u, int v, int x) {
        return dist(u, x) + dist(x, v) == dist(u, v);
    }

    pair<int,int> create_path(int u, int v, int x) {
        if (on_path(u, v, x)) return {u, v};
        if (on_path(u, x, v)) return {u, x};
        if (on_path(v, x, u)) return {v, x};
        return {-1, -1};
    }
};

・DSU on Tree
vector<vector<int>> G;
vector<int> sz;
vector<vector<int>> A;

void get_size(int n, int p) {
    for (int v : G[n]) {
        if (v == p) continue;
        get_size(v, n);
        sz[n] += sz[v];
    }
}

void dfs(int n, int p, bool keep) {
    int MAX = -1;
    int bigv = -1;

    for (int v : G[n]) {
        if (v == p) continue;
        if (sz[v] > MAX) {
            MAX = sz[v];
            bigv = v;
        }
    }

    for (int v : G[n]) {
        if (v == p || v == bigv) continue;
        dfs(v, n, false);
    }

    if (bigv != -1) {
        dfs(bigv, n, true);
        swap(A[n], A[bigv]);
    }

    for (int v : G[n]) {
        if (v == p || v == bigv) continue;
        for (int x : A[v]) {
            A[n].push_back(x);
        }
    }

    A[n].push_back(n);

    if (!keep) {
        for (int v : A[n]) {
        }
    }
}

・Mo's algorithm
template<class T, class AddL, class AddR, class RemoveL, class RemoveR>
struct Mo {
    int N;
    int route;
    AddL add_L;
    AddR add_R;
    RemoveL remove_L;
    RemoveR remove_R;

    Mo(int N, AddL add_L, AddR add_R, RemoveL remove_L, RemoveR remove_R, int route = -1)
        : N(N), add_L(add_L), add_R(add_R), remove_L(remove_L), remove_R(remove_R) {
        if (route == -1) this->route = sqrt(N) + 1;
        else this->route = route;
    }

    vector<T> solve(const vector<pair<int,int>>& query) {
        int Q = query.size();
        vector<vector<tuple<int,int,int>>> bucket(N / route + 1);

        for (int i = 0; i < Q; i++) {
            auto [l, r] = query[i];
            bucket[l / route].emplace_back(i, l, r);
        }

        vector<T> ans(Q, -1);

        T now = 0;
        int left = 0;
        int right = -1;

        for (int idx = 0; idx < (int)bucket.size(); idx++) {
            auto &q = bucket[idx];

            sort(q.begin(), q.end(), [&](auto &a, auto &b) {
                if (idx % 2) return get<2>(a) > get<2>(b);
                else return get<2>(a) < get<2>(b);
            });

            for (auto [i, l, r] : q) {
                while (l < left) {
                    left--;
                    now = add_L(now, left);
                }
                while (right < r) {
                    right++;
                    now = add_R(now, right);
                }
                while (left < l) {
                    now = remove_L(now, left);
                    left++;
                }
                while (r < right) {
                    now = remove_R(now, right);
                    right--;
                }
                ans[i] = now;
            }
        }

        return ans;
    }
};

auto addL = [&](long long now, int idx) {
    // 処理
    return now;
};

auto addR = [&](long long now, int idx) {
    return now;
};

auto removeL = [&](long long now, int idx) {
    return now;
};

auto removeR = [&](long long now, int idx) {
    return now;
};

・最大値（最小値）Top2
struct Top2 {
    pair<ll, int> best1, best2;

    Top2() {
        best1 = {-INFLL, -1};
        best2 = {-INFLL, -1};
    }

    void add(ll val, int genre) {
        auto [b1v, b1g] = best1;
        auto [b2v, b2g] = best2;

        if (genre == b1g) {
            if (val > b1v) {
                best1 = {val, genre};
            }
        } else {
            if (val > b1v) {
                best2 = best1;
                best1 = {val, genre};
            } else if (val > b2v) {
                best2 = {val, genre};
            }
        }
    }

    ll get_exclude(int genre) const {
        if (best1.second != genre) return best1.first;
        return best2.first;
    }
};

・IntervalMap
struct IntervalMap {
    struct Node {
        ll l, r, x;
        bool operator<(const Node& other) const {
            return l < other.l;
        }
    };
    set<Node> s;

    void add_interval(long long l, long long r, long long x) {
    }

    void remove_interval(long long l, long long r, long long x) {
    }

    // 先に[l, r]をremoveで消してから呼び出す
    void add(long long l, long long r, long long x) {
        s.insert({l, r, x});
        add_interval(l, r, x);
    }

    long long remove(long long l, long long r) {
        long long cnt = 0;
        auto it = s.lower_bound({l, -1, -1});
        if (it != s.begin()) {
            auto pit = prev(it);
            if (pit->r >= l) {
                long long L = pit->l, R = pit->r, X = pit->x;
                cnt += min(r, R) - l + 1;
                remove_interval(L, R, X);
                s.erase(pit);
                if (L <= l - 1) {
                    s.insert({L, l - 1, X});
                    add_interval(L, l - 1, X);
                }
                if (r < R) {
                    s.insert({r + 1, R, X});
                    add_interval(r + 1, R, X);
                }
            }
        }
        while (it != s.end() && it->l <= r) {
            long long L = it->l, R = it->r, X = it->x;
            cnt += min(R, r) - L + 1;
            remove_interval(L, R, X);
            it = s.erase(it);
            if (r < R) {
                s.insert({r + 1, R, X});
                add_interval(r + 1, R, X);
                break;
            }
        }
        return cnt;
    }

    vector<tuple<long long,long long,long long>> remove2(long long l, long long r) {
        vector<tuple<long long,long long,long long>> ans;
        auto it = s.lower_bound({l, -1, -1});
        if (it != s.begin()) {
            auto pit = prev(it);
            if (pit->r >= l) {
                long long L = pit->l, R = pit->r, X = pit->x;
                remove_interval(L, R, X);
                s.erase(pit);
                if (L <= l - 1) {
                    s.insert({L, l - 1, X});
                    add_interval(L, l - 1, X);
                }
                ans.push_back({l, min(R, r), X});
                if (r < R) {
                    s.insert({r + 1, R, X});
                    add_interval(r + 1, R, X);
                }
            }
        }
        while (it != s.end() && it->l <= r) {
            long long L = it->l, R = it->r, X = it->x;
            ans.push_back({L, min(R, r), X});
            remove_interval(L, R, X);
            it = s.erase(it);
            if (r < R) {
                s.insert({r + 1, R, X});
                add_interval(r + 1, R, X);
                break;
            }
        }
        return ans;
    }
};

・SortedMultiset
// https://github.com/tatyam-prime/SortedSet/blob/main/SortedMultiset.py
template<class T>
class SortedMultiset {
private:
    static constexpr int BUCKET_RATIO = 16;
    static constexpr int SPLIT_RATIO = 24;

    vector<vector<T>> a;
    int sz = 0;

    void build(const vector<T>& v) {
        int n = v.size();
        sz = n;
        if (n == 0) {
            a.clear();
            return;
        }

        int num_bucket = (int)ceil(sqrt((double)n / BUCKET_RATIO));
        a.resize(num_bucket);

        for (int i = 0; i < num_bucket; i++) {
            int l = n * i / num_bucket;
            int r = n * (i + 1) / num_bucket;
            a[i] = vector<T>(v.begin() + l, v.begin() + r);
        }
    }

    pair<int,int> position(const T& x) const {
        // return {bucket index, position in bucket}
        for (int i = 0; i < (int)a.size(); i++) {
            if (x <= a[i].back()) {
                return {i, lower_bound(a[i].begin(), a[i].end(), x) - a[i].begin()};
            }
        }

        int i = (int)a.size() - 1;
        return {i, (int)a[i].size()};
    }

    T pop_internal(int b, int i) {
        T ans = a[b][i];
        a[b].erase(a[b].begin() + i);
        sz--;

        if (a[b].empty()) {
            a.erase(a.begin() + b);
        }

        return ans;
    }

public:
    SortedMultiset() {}

    SortedMultiset(vector<T> v) {
        sort(v.begin(), v.end());
        build(v);
    }

    int size() const {
        return sz;
    }

    bool empty() const {
        return sz == 0;
    }

    bool contains(const T& x) const {
        if (sz == 0) return false;

        auto [b, i] = position(x);

        return i < (int)a[b].size() && a[b][i] == x;
    }

    int count(const T& x) const {
        return index_right(x) - index(x);
    }

    void add(const T& x) {
        if (sz == 0) {
            a = {{x}};
            sz = 1;
            return;
        }

        auto [b, i] = position(x);

        a[b].insert(a[b].begin() + i, x);
        sz++;

        if ((int)a[b].size() > (int)a.size() * SPLIT_RATIO) {
            int mid = a[b].size() >> 1;

            vector<T> left(a[b].begin(), a[b].begin() + mid);
            vector<T> right(a[b].begin() + mid, a[b].end());

            a.erase(a.begin() + b);
            a.insert(a.begin() + b, right);
            a.insert(a.begin() + b, left);
        }
    }

    bool discard(const T& x) {
        if (sz == 0) return false;

        auto [b, i] = position(x);

        if (i == (int)a[b].size() || a[b][i] != x) {
            return false;
        }

        pop_internal(b, i);
        return true;
    }

    optional<T> lt(const T& x) const {
        for (int i = (int)a.size() - 1; i >= 0; i--) {
            if (a[i][0] < x) {
                auto it = lower_bound(a[i].begin(), a[i].end(), x);
                --it;
                return *it;
            }
        }
        return nullopt;
    }

    optional<T> le(const T& x) const {
        for (int i = (int)a.size() - 1; i >= 0; i--) {
            if (a[i][0] <= x) {
                auto it = upper_bound(a[i].begin(), a[i].end(), x);
                --it;
                return *it;
            }
        }
        return nullopt;
    }

    optional<T> gt(const T& x) const {
        for (int i = 0; i < (int)a.size(); i++) {
            if (a[i].back() > x) {
                auto it = upper_bound(a[i].begin(), a[i].end(), x);
                return *it;
            }
        }
        return nullopt;
    }

    optional<T> ge(const T& x) const {
        for (int i = 0; i < (int)a.size(); i++) {
            if (a[i].back() >= x) {
                auto it = lower_bound(a[i].begin(), a[i].end(), x);
                return *it;
            }
        }
        return nullopt;
    }

    T operator[](int i) const {
        if (i < 0) i += sz;

        for (auto &bucket : a) {
            if (i < (int)bucket.size()) {
                return bucket[i];
            }
            i -= bucket.size();
        }

        throw out_of_range("index out of range");
    }

    T pop(int i = -1) {
        if (i < 0) i += sz;

        for (int b = 0; b < (int)a.size(); b++) {
            if (i < (int)a[b].size()) {
                return pop_internal(b, i);
            }
            i -= a[b].size();
        }

        throw out_of_range("index out of range");
    }

    int index(const T& x) const {
        int ans = 0;

        for (auto &bucket : a) {
            if (bucket.back() >= x) {
                return ans + (lower_bound(bucket.begin(), bucket.end(), x) - bucket.begin());
            }
            ans += bucket.size();
        }

        return ans;
    }

    int index_right(const T& x) const {
        int ans = 0;

        for (auto &bucket : a) {
            if (bucket.back() > x) {
                return ans + (upper_bound(bucket.begin(), bucket.end(), x) - bucket.begin());
            }
            ans += bucket.size();
        }

        return ans;
    }

    vector<T> to_vector() const {
        vector<T> res;

        for (auto &bucket : a) {
            for (auto &x : bucket) {
                res.push_back(x);
            }
        }

        return res;
    }
};

・最小共通祖先（Sparse Table）
struct LCA {
    vector<vector<int>> G;
    int N;
    vector<int> A;
    vector<int> IN;
    vector<int> depth;

    int L;
    vector<int> log;
    vector<vector<int>> sparse;

    LCA(const vector<vector<int>>& G, int root = 0)
        : G(G), N((int)G.size()),
          IN(N), depth(N)
    {
        dfs(root, -1, 0);

        L = (int)A.size();

        log.assign(L + 1, 0);
        for (int i = 2; i <= L; i++) {
            log[i] = log[i >> 1] + 1;
        }

        sparse.resize(log[L] + 1);
        sparse[0] = A;

        int b = 1;
        for (int i = 0; i < log[L]; i++) {
            sparse[i + 1].resize((int)sparse[i].size() - b);

            for (int j = 0; j + b < (int)sparse[i].size(); j++) {
                int p = sparse[i][j];
                int q = sparse[i][j + b];
                sparse[i + 1][j] =
                    (depth[p] <= depth[q] ? p : q);
            }

            b <<= 1;
        }
    }

    void dfs(int n, int p, int d) {
        depth[n] = d;
        IN[n] = (int)A.size();
        A.push_back(n);

        for (int v : G[n]) {
            if (v == p) continue;
            dfs(v, n, d + 1);
            A.push_back(n);
        }
    }

    int lca(int u, int v) {
        int l = IN[u];
        int r = IN[v];

        if (l > r) {
            swap(l, r);
        }

        int k = log[r - l + 1];

        int a = sparse[k][l];
        int b = sparse[k][r - (1 << k) + 1];

        return depth[a] <= depth[b] ? a : b;
    }

    int dist(int a, int b) {
        return depth[a] + depth[b] - depth[lca(a, b)]*2;
    }
};

・Segment Tree Beats
// https://smijake3.hatenablog.com/entry/2019/04/28/021457
class SegmentTreeBeats {
    int n, n0;
    vl max_v, smax_v, max_c;
    vl min_v, smin_v, min_c;
    vl change, add;
    vl SUM;

    void update_node_max(int k, ll x) {
        SUM[k] += (x-max_v[k])*max_c[k];

        if (min_v[k] == max_v[k]) {
            min_v[k] = x;
        }
        else if (smin_v[k] == max_v[k]) {
            smin_v[k] = x;
        }

        max_v[k] = x;
    }

    void update_node_min(int k, ll x) {
        SUM[k] += (x-min_v[k])*min_c[k];

        if (max_v[k] == min_v[k]) {
            max_v[k] = x;
        }
        else if (smax_v[k] == min_v[k]) {
            smax_v[k] = x;
        }

        min_v[k] = x;
    }

    void update_node_change(int k, ll x, int size) {
        max_v[k] = min_v[k] = x;
        smax_v[k] = -INFLL;
        smin_v[k] = INFLL;
        max_c[k] = min_c[k] = size;
        SUM[k] = x*size;

        change[k] = x;
        add[k] = 0;
    }

    void update_node_add(int k, ll x, int size) {
        SUM[k] += x*size;

        if (max_v[k] != -INFLL) max_v[k] += x;
        if (smax_v[k] != -INFLL) smax_v[k] += x;
        if (min_v[k] != INFLL) min_v[k] += x;
        if (smin_v[k] != INFLL) smin_v[k] += x;

        add[k] += x;
    }

    void push(int k, int l, int r) {
        int m = (l+r)/2;

        if (change[k] != INFLL) {
            update_node_change(2*k+1, change[k], m-l);
            update_node_change(2*k+2, change[k], r-m);
            change[k] = INFLL;
        }

        if (add[k] != 0) {
            update_node_add(2*k+1, add[k], m-l);
            update_node_add(2*k+2, add[k], r-m);
            add[k] = 0;
        }

        if (max_v[k] < max_v[2*k+1]) {
            update_node_max(2*k+1, max_v[k]);
        }
        if (max_v[k] < max_v[2*k+2]) {
            update_node_max(2*k+2, max_v[k]);
        }

        if (min_v[k] > min_v[2*k+1]) {
            update_node_min(2*k+1, min_v[k]);
        }
        if (min_v[k] > min_v[2*k+2]) {
            update_node_min(2*k+2, min_v[k]);
        }
    }

    void update(int k) {
        SUM[k] = SUM[2*k+1]+SUM[2*k+2];

        if (max_v[2*k+1] < max_v[2*k+2]) {
            max_v[k] = max_v[2*k+2];
            max_c[k] = max_c[2*k+2];
            smax_v[k] = max(max_v[2*k+1], smax_v[2*k+2]);
        }
        else if (max_v[2*k+1] > max_v[2*k+2]) {
            max_v[k] = max_v[2*k+1];
            max_c[k] = max_c[2*k+1];
            smax_v[k] = max(smax_v[2*k+1], max_v[2*k+2]);
        }
        else {
            max_v[k] = max_v[2*k+1];
            max_c[k] = max_c[2*k+1]+max_c[2*k+2];
            smax_v[k] = max(smax_v[2*k+1], smax_v[2*k+2]);
        }

        if (min_v[2*k+1] < min_v[2*k+2]) {
            min_v[k] = min_v[2*k+1];
            min_c[k] = min_c[2*k+1];
            smin_v[k] = min(smin_v[2*k+1], min_v[2*k+2]);
        }
        else if (min_v[2*k+1] > min_v[2*k+2]) {
            min_v[k] = min_v[2*k+2];
            min_c[k] = min_c[2*k+2];
            smin_v[k] = min(smin_v[2*k+1], min_v[2*k+2]);
        }
        else {
            min_v[k] = min_v[2*k+1];
            min_c[k] = min_c[2*k+1] + min_c[2*k+2];
            smin_v[k] = min(smin_v[2*k+1], smin_v[2*k+2]);
        }
    }

    void _update_min(ll x, int a, int b, int k, int l, int r) {
        if (b <= l || r <= a || max_v[k] <= x) {
            return;
        }
        if (a <= l && r <= b && smax_v[k] < x) {
            update_node_max(k, x);
            return;
        }

        push(k, l, r);
        _update_min(x, a, b, 2*k+1, l, (l+r)/2);
        _update_min(x, a, b, 2*k+2, (l+r)/2, r);
        update(k);
    }

    void _update_max(ll x, int a, int b, int k, int l, int r) {
        if (b <= l || r <= a || min_v[k] >= x) {
            return;
        }
        if (a <= l && r <= b && smin_v[k] > x) {
            update_node_min(k, x);
            return;
        }

        push(k, l, r);
        _update_max(x, a, b, 2*k+1, l, (l+r)/2);
        _update_max(x, a, b, 2*k+2, (l+r)/2, r);
        update(k);
    }

    void _update_change(ll x, int a, int b, int k, int l, int r) {
        if (b <= l || r <= a) {
            return;
        }

        if (a <= l && r <= b) {
            update_node_change(k, x, r-l);
            return;
        }

        push(k, l, r);
        _update_change(x, a, b, 2*k+1, l, (l+r)/2);
        _update_change(x, a, b, 2*k+2, (l+r)/2, r);
        update(k);
    }

    void _update_add(ll x, int a, int b, int k, int l, int r) {
        if (b <= l || r <= a) {
            return;
        }

        if (a <= l && r <= b) {
            update_node_add(k, x, r-l);
            return;
        }

        push(k, l, r);
        _update_add(x, a, b, 2*k+1, l, (l+r)/2);
        _update_add(x, a, b, 2*k+2, (l+r)/2, r);
        update(k);
    }

    ll _query_max(int a, int b, int k, int l, int r) {
        if (b <= l || r <= a) {
            return -INFLL;
        }
        if (a <= l && r <= b) {
            return max_v[k];
        }

        push(k, l, r);
        ll lv = _query_max(a, b, 2*k+1, l, (l+r)/2);
        ll rv = _query_max(a, b, 2*k+2, (l+r)/2, r);
        return max(lv, rv);
    }

    ll _query_min(int a, int b, int k, int l, int r) {
        if (b <= l || r <= a) {
            return INFLL;
        }
        if (a <= l && r <= b) {
            return min_v[k];
        }

        push(k, l, r);
        ll lv = _query_min(a, b, 2*k+1, l, (l+r)/2);
        ll rv = _query_min(a, b, 2*k+2, (l+r)/2, r);
        return min(lv, rv);
    }

    ll _query_sum(int a, int b, int k, int l, int r) {
        if (b <= l || r <= a) {
            return 0;
        }
        if (a <= l && r <= b) {
            return SUM[k];
        }

        push(k, l, r);
        ll lv = _query_sum(a, b, 2*k+1, l, (l+r)/2);
        ll rv = _query_sum(a, b, 2*k+2, (l+r)/2, r);
        return lv + rv;
    }

public:
    SegmentTreeBeats(int n) : SegmentTreeBeats(vl(n, 0)) {}

    SegmentTreeBeats(vl A) : n(A.size()) {
        n0 = 1;
        while (n0 < n) n0 <<= 1;

        int size = n0 << 1;

        max_v.assign(size, -INFLL);
        smax_v.assign(size, -INFLL);
        max_c.assign(size, 0);
        min_v.assign(size, INFLL);
        smin_v.assign(size, INFLL);
        min_c.assign(size, 0);
        change.assign(size, INFLL);
        add.assign(size, 0);
        SUM.assign(size, 0);

        for (int i = 0; i < n; i++) {
            max_v[n0-1+i] = A[i];
            smax_v[n0-1+i] = -INFLL;
            max_c[n0-1+i] = 1;
            min_v[n0-1+i] = A[i];
            smin_v[n0-1+i] = INFLL;
            min_c[n0-1+i] = 1;
            SUM[n0-1+i] = A[i];
        }

        for (int i = n; i < n0; i++) {
            max_v[n0-1+i] = smax_v[n0-1+i] = -INFLL;
            max_c[n0-1+i] = 0;
            min_v[n0-1+i] = smin_v[n0-1+i] = INFLL;
            min_c[n0-1+i] = 0;
            SUM[n0-1+i] = 0;
        }

        for (int i = n0-2; i >= 0; i--) {
            update(i);
        }
    }

    void update_min(int a, int b, ll x) {
        _update_min(x, a, b, 0, 0, n0);
    }

    void update_max(int a, int b, ll x) {
        _update_max(x, a, b, 0, 0, n0);
    }

    void update_change(int a, int b, ll x) {
        _update_change(x, a, b, 0, 0, n0);
    }

    void update_add(int a, int b, ll x) {
        _update_add(x, a, b, 0, 0, n0);
    }

    ll query_max(int a, int b) {
        return _query_max(a, b, 0, 0, n0);
    }

    ll query_min(int a, int b) {
        return _query_min(a, b, 0, 0, n0);
    }

    ll query_sum(int a, int b) {
        return _query_sum(a, b, 0, 0, n0);
    }
};

・オイラーツアー（3種類）
// 頂点・部分木用
struct EulerTourVertexSubtree {
    int N;
    vector<int> IN, OUT, REV;

    EulerTourVertexSubtree(const vector<vector<int>>& G, int root = 0) {
        N = G.size();
        IN.assign(N, -1);
        OUT.assign(N, -1);
        REV.assign(N, -1);
        dfs(G, root, -1, 0);
    }

    int dfs(const vector<vector<int>>& G, int n, int pre, int time) {
        IN[n] = time;
        REV[time] = n;
        for (int v : G[n]) {
            if (v == pre) continue;
            time = dfs(G, v, n, time + 1);
        }
        OUT[n] = time;
        return time;
    }

    int index(int n) {
        return IN[n];
    }

    pair<int, int> interval(int n) {
        return {IN[n], OUT[n]};
    }

    int index_rev(int n) {
        return REV[n];
    }
};

// 辺用
struct EulerTourEdge {
    int N;
    vector<int> IN, OUT, REV;

    EulerTourEdge(const vector<vector<pair<int, int>>>& G, int root = 0) {
        N = G.size();
        IN.assign(N - 1, -1);
        OUT.assign(N - 1, -1);
        REV.assign(2 * (N - 1), -1);
        dfs(G, root, -1, 0);
    }

    int dfs(const vector<vector<pair<int, int>>>& G, int n, int pre, int time) {
        for (auto [v, eidx] : G[n]) {
            if (v == pre) continue;
            IN[eidx] = time;
            REV[time] = eidx;
            time++;
            time = dfs(G, v, n, time);
            OUT[eidx] = time;
            REV[time] = eidx;
            time++;
        }
        return time;
    }

    int index_in(int n) {
        return IN[n];
    }

    int index_out(int n) {
        return OUT[n];
    }

    pair<int, int> interval(int n) {
        return {IN[n], OUT[n]};
    }

    int index_rev(int n) {
        return REV[n];
    }
};

// 頂点・行き帰り用
struct EulerTourVertex {
    int N;
    vector<int> IN, OUT, REV;

    EulerTourVertex(const vector<vector<int>>& G, int root = 0) {
        N = G.size();
        IN.assign(N, -1);
        OUT.assign(N, -1);
        REV.assign(2 * N, -1);
        dfs(G, root, -1, 0);
    }

    int dfs(const vector<vector<int>>& G, int n, int pre, int time) {
        IN[n] = time;
        REV[time] = n;
        time++;
        for (int v : G[n]) {
            if (v == pre) continue;
            time = dfs(G, v, n, time);
        }
        OUT[n] = time;
        REV[time] = n;
        time++;
        return time;
    }

    int index_in(int n) {
        return IN[n];
    }

    int index_out(int n) {
        return OUT[n];
    }

    pair<int, int> interval(int n) {
        return {IN[n], OUT[n]};
    }

    int index_rev(int n) {
        return REV[n];
    }
};

・BIT2D
struct BIT2D {
    enum class Mode { GE, LE, BETWEEN };

    int n;
    vector<ll> a;
    vector<vector<ll>> xs, candidates;
    vector<atcoder::fenwick_tree<int>> bit_cnt;
    vector<atcoder::fenwick_tree<ll>> bit_sum;
    bool built = false;

    BIT2D(const vector<ll>& a)
        : n((int)a.size()), a(a), xs(n + 1), candidates(n) {
        for (int p = 0; p < n; p++) {
            reserve_update(p, a[p]);
        }
    }

    // build前に、位置pに代入する可能性のある値xを登録
    void reserve_update(int p, ll x) {
        assert(!built && 0 <= p && p < n);
        candidates[p].push_back(x);

        for (int i = p + 1; i <= n; i += i & -i) {
            xs[i].push_back(x);
        }
    }

    // 構築
    void build() {
        assert(!built);

        for (auto& v : candidates) {
            sort(v.begin(), v.end());
            v.erase(unique(v.begin(), v.end()), v.end());
        }

        bit_cnt.resize(n + 1);
        bit_sum.resize(n + 1);

        for (int i = 1; i <= n; i++) {
            auto& v = xs[i];
            sort(v.begin(), v.end());
            v.erase(unique(v.begin(), v.end()), v.end());

            bit_cnt[i] = atcoder::fenwick_tree<int>((int)v.size());
            bit_sum[i] = atcoder::fenwick_tree<ll>((int)v.size());
        }

        built = true;
        for (int p = 0; p < n; p++) {
            add_value(p, a[p], 1);
        }
    }

    // a[p] = x
    void set(int p, ll x) {
        assert(built && 0 <= p && p < n);
        if (a[p] == x) return;

        // 更新値は位置pに事前登録されている必要がある
        assert(binary_search(
            candidates[p].begin(), candidates[p].end(), x
        ));

        add_value(p, a[p], -1);
        add_value(p, x, 1);
        a[p] = x;
    }

    // [l,r) 内の値 >= x の (個数, 合計)
    pair<ll, ll> query_ge(int l, int r, ll x) {
        return range_query(l, r, x, 0, Mode::GE);
    }

    // [l,r) 内の値 <= x の (個数, 合計)
    pair<ll, ll> query_le(int l, int r, ll x) {
        return range_query(l, r, 0, x, Mode::LE);
    }

    // [l,r) 内の x <= 値 <= y の (個数, 合計)
    pair<ll, ll> query_range(int l, int r, ll x, ll y) {
        assert(built && 0 <= l && l <= r && r <= n);
        if (x > y) return {0, 0};
        return range_query(l, r, x, y, Mode::BETWEEN);
    }

    // 個数のみ
    ll count_ge(int l, int r, ll x) {
        return query_ge(l, r, x).first;
    }

    ll count_le(int l, int r, ll x) {
        return query_le(l, r, x).first;
    }

    ll count_range(int l, int r, ll x, ll y) {
        return query_range(l, r, x, y).first;
    }

    // 合計のみ
    ll sum_ge(int l, int r, ll x) {
        return query_ge(l, r, x).second;
    }

    ll sum_le(int l, int r, ll x) {
        return query_le(l, r, x).second;
    }

    ll sum_range(int l, int r, ll x, ll y) {
        return query_range(l, r, x, y).second;
    }

    // a[p] を取得
    ll get(int p) const {
        assert(0 <= p && p < n);
        return a[p];
    }

private:
    // 位置pの値xをd個追加（d = +1 または -1）
    void add_value(int p, ll x, int d) {
        for (int i = p + 1; i <= n; i += i & -i) {
            const auto& v = xs[i];
            int j = lower_bound(v.begin(), v.end(), x) - v.begin();

            bit_cnt[i].add(j, d);
            bit_sum[i].add(j, (ll)d * x);
        }
    }

    // [0,r) における条件を満たす値の (個数, 合計)
    pair<ll, ll> prefix_query(int r, ll x, ll y, Mode mode) {
        ll cnt = 0, sum = 0;

        for (int i = r; i > 0; i -= i & -i) {
            const auto& v = xs[i];
            int L = 0, R = (int)v.size();

            if (mode != Mode::LE) {
                L = lower_bound(v.begin(), v.end(), x) - v.begin();
            }
            if (mode != Mode::GE) {
                R = upper_bound(v.begin(), v.end(), y) - v.begin();
            }

            cnt += bit_cnt[i].sum(L, R);
            sum += bit_sum[i].sum(L, R);
        }

        return {cnt, sum};
    }

    // [l,r) における条件を満たす値の (個数, 合計)
    pair<ll, ll> range_query(int l, int r, ll x, ll y, Mode mode) {
        assert(built && 0 <= l && l <= r && r <= n);

        auto R = prefix_query(r, x, y, mode);
        auto L = prefix_query(l, x, y, mode);

        return {R.first - L.first, R.second - L.second};
    }
};

・BITSortedMultiset
struct BITSortedMultiset {
    int N;
    vector<ll> A;
    vector<SortedMultiset<ll>> bit;

    BITSortedMultiset(const vector<ll>& A)
        : N(A.size()), A(A), bit(N+1) {
        for (int p = 0; p < N; p++) {
            for (int i = p+1; i <= N; i += i&-i) {
                bit[i].add(A[p]);
            }
        }
    }

    // A[p] = x
    void set(int p, ll x) {
        assert(0 <= p && p < N);
        if (A[p] == x) return;

        for (int i = p+1; i <= N; i += i&-i) {
            bit[i].discard(A[p]);
            bit[i].add(x);
        }

        A[p] = x;
    }

    // [0,r) 内の x以上の個数
    ll prefix_count_ge(int r, ll x) const {
        assert(0 <= r && r <= N);
        ll ans = 0;

        for (int i = r; i > 0; i -= i&-i) {
            ans += bit[i].size() - bit[i].index(x);
        }

        return ans;
    }

    // [0,r) 内の x以下の個数
    ll prefix_count_le(int r, ll x) const {
        assert(0 <= r && r <= N);
        ll ans = 0;

        for (int i = r; i > 0; i -= i&-i) {
            ans += bit[i].index_right(x);
        }

        return ans;
    }

    // [0,r) 内の x以上y以下の個数
    ll prefix_count_range(int r, ll x, ll y) const {
        assert(0 <= r && r <= N);
        if (x > y) return 0;

        ll ans = 0;

        for (int i = r; i > 0; i -= i&-i) {
            ans += bit[i].index_right(y) - bit[i].index(x);
        }

        return ans;
    }

    // [l,r) 内の x以上の個数
    ll count_ge(int l, int r, ll x) const {
        assert(0 <= l && l <= r && r <= N);
        return prefix_count_ge(r, x) - prefix_count_ge(l, x);
    }

    // [l,r) 内の x以下の個数
    ll count_le(int l, int r, ll x) const {
        assert(0 <= l && l <= r && r <= N);
        return prefix_count_le(r, x) - prefix_count_le(l, x);
    }

    // [l,r) 内の x以上y以下の個数
    ll count_range(int l, int r, ll x, ll y) const {
        assert(0 <= l && l <= r && r <= N);
        if (x > y) return 0;
        return prefix_count_range(r, x, y)
             - prefix_count_range(l, x, y);
    }

    ll get(int p) const {
        assert(0 <= p && p < N);
        return A[p];
    }
};