text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Top down segment tree supporting advanced range updates and queries
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   C: struct to combine data and lazy values
//     Required Fields:
//       Data: the data type
//       Lazy: the lazy type
//       Arg: the argument type
//     Required Functions:
//       static qdef(): returns the query default value of type Data
//       static ldef(): returns the lazy default value of type Lazy
//       static makeLazy(d, a): returns a lazy value created from the Data d
//         and Arg a
//       static breakCond(d, a): returns the condition for when to stop the
//         recusion based on Data d and Arg a
//       static tagCond(d, a): returns the condition for when to tag
//         based on Data d and Arg a
//       static merge(l, r): returns the values l of type Data merged with
//         r of type Data, must be associative
//       static applyLazy(l, r, k): returns the value r of type Lazy applied to
//         l of type Data over a segment of length k, must be associative
//       static mergeLazy(l, r): returns the values l of type Lazy merged with
//         r of type Lazy, must be associative
//     Sample Struct: supporting range modulo and range sum queries
//       struct C {
//         struct Data { int sm, mx, mn; };
//         using Lazy = int;
//         using Arg = int;
//         static Data qdef() {
//           Data ret; ret.sm = 0;
//           ret.mx = numeric_limits<int>::min();
//           ret.mn = numeric_limits<int>::max();
//           return ret;
//         }
//         static Lazy ldef() { return -1; }
//         static Lazy makeLazy(const Data &d, const Arg &a) { return a; }
//         static bool breakCond(const Data &d, const Arg &a) {
//           return d.mx < a;
//         }
//         static bool tagCond(const Data &d, const Arg &a) {
//           return d.mn == d.mx;
//         }
//         static Data merge(const Data &l, const Data &r) {
//           Data ret; ret.sm = l.sm + r.sm;
//           ret.mx = max(l.mx, r.mx); ret.mn = min(l.mn, r.mn); return ret;
//         }
//         static Data applyLazy(const Data &l, const Lazy &r, int k) {
//           Data ret; ret.sm = r * k; ret.mx = ret.mn = r; return ret;
//         }
//         static Lazy mergeLazy(const Lazy &l, const Lazy &r) { return r; }
//       };
// Constructor Arguments:
//   N: the size of the array
//   A: a vector of type C::Data
// Functions:
//   update(l, r, v): update the range [l, r] with the Arg v
//   query(l, r): queries the range [l, r] and returns the aggregate value
//   bsearchPrefix(l, r, f): returns the smallest index i in the range [l, r]
//     such that f(query(l, i)) returns true or r + 1 if none exist
//   bsearchSuffix(l, r, f): returns the largest index i in the range [l, r]
//     such that f(query(i, r)) returns true of l - 1 if none exist
// In practice, has a small constant
// Time Complexity:
//   constructor: O(N)
//   update: dependent on C
//   query, bsearchPrefix, bsearchSuffix: O(log N)
// Memory Complexity: O(N)
// Tested:
//   https://codeforces.com/contest/438/problem/D
//   https://judge.yosupo.jp/problem/range_chmin_chmax_add_range_sum
template <class C> struct SegmentTreeBeats {
  using Data = typename C::Data; using Lazy = typename C::Lazy;
  using Arg = typename C::Arg;
  struct Node {
    Data val; Lazy lz; Node() : val(C::qdef()), lz(C::ldef()) {}
  };
  int N; vector<Node> TR;
  void apply(int x, int tl, int tr, const Lazy &v) {
    TR[x].val = C::applyLazy(TR[x].val, v, tr - tl + 1);
    TR[x].lz = C::mergeLazy(TR[x].lz, v);
  }
  void propagate(int x, int tl, int tr) {
    if (TR[x].lz != C::ldef()) {
      int m = tl + (tr - tl) / 2; apply(x * 2, tl, m, TR[x].lz);
      apply(x * 2 + 1, m + 1, tr, TR[x].lz); TR[x].lz = C::ldef();
    }
  }
  void build(const vector<Data> &A, int x, int tl, int tr) {
    if (tl == tr) { TR[x].val = A[tl]; return; }
    int m = tl + (tr - tl) / 2;
    build(A, x * 2, tl, m); build(A, x * 2 + 1, m + 1, tr);
    TR[x].val = C::merge(TR[x * 2].val, TR[x * 2 + 1].val);
  }
  void update(int x, int tl, int tr, int l, int r, const Arg &v) {
    if (r < tl || tr < l || C::breakCond(TR[x].val, v)) return;
    if (l <= tl && tr <= r && C::tagCond(TR[x].val, v)) {
      Lazy lz = C::makeLazy(TR[x].val, v);
      TR[x].val = C::applyLazy(TR[x].val, lz, tr - tl + 1);
      TR[x].lz = C::mergeLazy(TR[x].lz, lz); return;
    }
    propagate(x, tl, tr); int m = tl + (tr - tl) / 2;
    update(x * 2, tl, m, l, r, v); update(x * 2 + 1, m + 1, tr, l, r, v);
    TR[x].val = C::merge(TR[x * 2].val, TR[x * 2 + 1].val);
  }
  Data query(int x, int tl, int tr, int l, int r) {
    if (r < tl || tr < l) return C::qdef();
    if (l <= tl && tr <= r) return TR[x].val;
    propagate(x, tl, tr); int m = tl + (tr - tl) / 2;
    return C::merge(query(x * 2, tl, m, l, r),
                   query(x * 2 + 1, m + 1, tr, l, r));
  }
  template <class F>
  int bsearchPrefix(int x, int tl, int tr, int l, int r, Data &agg, F f) {
    if (r < tl || tr < l) return r + 1;
    if (tl != tr) propagate(x, tl, tr);
    if (l <= tl && tr <= r) {
      Data v = C::merge(agg, TR[x].val); if (!f(v)) { agg = v; return r + 1; }
    }
    if (tl == tr) return tl;
    int m = tl + (tr - tl) / 2;
    int ret = bsearchPrefix(x * 2, tl, m, l, r, agg, f);
    return ret <= r ? ret : bsearchPrefix(x * 2 + 1, m + 1, tr, l, r, agg, f);
  }
  template <class F>
  int bsearchSuffix(int x, int tl, int tr, int l, int r, Data &agg, F f) {
    if (r < tl || tr < l) return l - 1;
    if (tl != tr) propagate(x, tl, tr);
    if (l <= tl && tr <= r) {
      Data v = C::merge(agg, TR[x].val); if (!f(v)) { agg = v; return l - 1; }
    }
    if (tl == tr) return tl;
    int m = tl + (tr - tl) / 2;
    int ret = bsearchSuffix(x * 2 + 1, m + 1, tr, l, r, agg, f);
    return l <= ret ? ret : bsearchSuffix(x * 2, tl, m, l, r, agg, f);
  }
  SegmentTreeBeats(int N) : N(N), TR(N == 0 ? 0 : 1 << __lg(N * 4 - 1)) {}
  SegmentTreeBeats(const vector<Data> &A) : SegmentTreeBeats(A.size()) {
    if (N > 0) { build(A, 1, 0, N - 1); }
  }
  void update(int l, int r, const Arg &v) { update(1, 0, N - 1, l, r, v); }
  Data query(int l, int r) { return query(1, 0, N - 1, l, r); }
  template <class F> int bsearchPrefix(int l, int r, F f) {
    Data agg = C::qdef(); return bsearchPrefix(1, 0, N - 1, l, r, agg, f);
  }
  template <class F> int bsearchSuffix(int l, int r, F f) {
    Data agg = C::qdef(); return bsearchSuffix(1, 0, N - 1, l, r, agg, f);
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Bottom up segment tree supporting range updates and range queries
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   C: struct to combine data and lazy values
//     Required Fields:
//       Data: the data type
//       Lazy: the lazy type
//     Required Functions:
//       static qdef(): returns the query default value of type Data
//       static ldef(): returns the lazy default value of type Lazy
//       static merge(l, r): returns the values l of type Data merged with
//         r of type Data, must be associative
//       static applyLazy(l, r, k): returns the value r of type Lazy applied to
//         l of type Data over a segment of length k, must be associative
//       static mergeLazy(l, r): returns the values l of type Lazy merged with
//         r of type Lazy, must be associative
//     Sample Struct: supporting range assignments and range sum queries
//       struct C {
//         using Data = int;
//         using Lazy = int;
//         static Data qdef() { return 0; }
//         static Lazy ldef() { return numeric_limits<int>::min(); }
//         static Data merge(const Data &l, const Data &r) { return l + r; }
//         static Data applyLazy(const Data &l, const Lazy &r, int k) {
//           return r * k;
//         }
//         static Lazy mergeLazy(const Lazy &l, const Lazy &r) { return r; }
//       };
// Constructor Arguments:
//   N: the size of the array
//   A: a vector of type C::Data, memory is saved if this is moved and has
//     a capacity of 2N
// Functions:
//   update(l, r, v): update the range [l, r] with the lazy value v
//   query(l, r): queries the range [l, r] and returns the aggregate value
// In practice, has a small constant, not quite as fast as fenwick trees,
//   and similar performance as top down segment trees
// Time Complexity:
//   constructor: O(N)
//   update, query: O(log N)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/lazy
//   https://mcpt.ca/problem/seq3
//   https://judge.yosupo.jp/problem/range_affine_range_sum
//   https://dmoj.ca/problem/dmopc17c4p6
template <class C> struct SegmentTreeLazyBottomUp {
  using Data = typename C::Data; using Lazy = typename C::Lazy;
  int N, lgN; vector<Data> TR; vector<Lazy> LZ;
  void apply(int i, const Lazy &v, int k) {
    TR[i] = C::applyLazy(TR[i], v, k);
    if (i < N) LZ[i] = C::mergeLazy(LZ[i], v);
  }
  void eval(int i, int k) {
    TR[i] = C::merge(TR[i * 2], TR[i * 2 + 1]);
    if (LZ[i] != C::ldef()) TR[i] = C::applyLazy(TR[i], LZ[i], k);
  }
  void propagate(int i) {
    int h = lgN + 1, k = 1 << lgN, ii = i >> h;
    for (; h > 0; ii = i >> --h, k /= 2) if (LZ[ii] != C::ldef()) {
      apply(ii * 2, LZ[ii], k); apply(ii * 2 + 1, LZ[ii], k);
      LZ[ii] = C::ldef();
    }
  }
  SegmentTreeLazyBottomUp(int N)
      : N(N), lgN(N == 0 ? 0 : __lg(N)),
        TR(N * 2, C::qdef()), LZ(N, C::ldef()) {}
  SegmentTreeLazyBottomUp(vector<Data> A)
      : N(A.size()), lgN(N == 0 ? 0 : __lg(N)), TR(move(A)), LZ(N, C::ldef()) {
    TR.resize(N * 2, C::qdef()); copy_n(TR.begin(), N, TR.begin() + N);
    for (int i = N - 1; i > 0; i--) TR[i] = C::merge(TR[i * 2], TR[i * 2 + 1]);
  }
  void update(int l, int r, const Lazy &v) {
    propagate(l += N); propagate(r += N); bool bl = 0, br = 0; int k = 1;
    for (; l <= r; l /= 2, r /= 2, k *= 2) {
      if (bl) eval(l - 1, k);
      if (br) eval(r + 1, k);
      if (l % 2) { apply(l++, v, k); bl = 1; }
      if (!(r % 2)) { apply(r--, v, k); br = 1; }
    }
    for (l--, r++; r; l /= 2, r /= 2, k *= 2) {
      if (bl) eval(l, k);
      if (br && (!bl || l != r)) eval(r, k);
    }
  }
  Data query(int l, int r) {
    propagate(l += N); propagate(r += N); Data ql = C::qdef(), qr = C::qdef();
    for (; l <= r; l /= 2, r /= 2) {
      if (l % 2) ql = C::merge(ql, TR[l++]);
      if (!(r % 2)) qr = C::merge(TR[r--], qr);
    }
    return C::merge(ql, qr);
  }
};

repo_name: wesley-a-leung/Resources

mime_type: text/x-c++


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Bottom up segment tree supporting point updates and range queries
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   C: struct to combine data and lazy values
//     Required Fields:
//       Data: the data type
//       Lazy: the lazy type
//     Required Functions:
//       static qdef(): returns the query default value of type Data
//       static merge(l, r): returns the values l of type Data merged with
//         r of type Data, must be associative
//       static applyLazy(l, r): returns the value r of type Lazy applied to
//         l of type Data, must be associative
//     Sample Struct: supporting point increments and range max queries
//       struct C {
//         using Data = int;
//         using Lazy = int;
//         static Data qdef() { return numeric_limits<int>::min(); }
//         static Data merge(const Data &l, const Data &r) {
//           return max(l, r);
//         }
//         static Data applyLazy(const Data &l, const Lazy &r) {
//           return l + r;
//         }
//       };
// Constructor Arguments:
//   N: the size of the array
//   A: a vector of type C::Data, memory is saved if this is moved and has
//     a capacity of 2N
// Functions:
//   update(i, v): updates the index i with the lazy value v
//   query(l, r): queries the range [l, r] and returns the aggregate value
// In practice, has a small constant, not quite as fast as fenwick trees,
//   but faster than top down segment trees
// Time Complexity:
//   constructor: O(N)
//   update, query: O(log N)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/ds3
//   https://dmoj.ca/problem/dmpg17g2
//   https://dmoj.ca/problem/coci17c1p5
//   https://judge.yosupo.jp/problem/point_set_range_composite
template <class C> struct SegmentTreeBottomUp {
  using Data = typename C::Data; using Lazy = typename C::Lazy;
  int N; vector<Data> TR;
  SegmentTreeBottomUp(int N) : N(N), TR(N * 2, C::qdef()) {}
  SegmentTreeBottomUp(vector<Data> A) : N(A.size()), TR(move(A))   {
    TR.resize(N * 2, C::qdef()); copy_n(TR.begin(), N, TR.begin() + N);
    for (int i = N - 1; i > 0; i--) TR[i] = C::merge(TR[i * 2], TR[i * 2 + 1]);
  }
  void update(int i, const Lazy &v) {
    for (i += N, TR[i] = C::applyLazy(TR[i], v); i /= 2;)
      TR[i] = C::merge(TR[i * 2], TR[i * 2 + 1]);
  }
  Data query(int l, int r) {
    Data ql = C::qdef(), qr = C::qdef();
    for (l += N, r += N; l <= r; l /= 2, r /= 2) {
      if (l % 2) ql = C::merge(ql, TR[l++]);
      if (!(r % 2)) qr = C::merge(TR[r--], qr);
    }
    return C::merge(ql, qr);
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Top down dynamic segment tree supporting range updates and range queries
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   IndexType: the type of the index of the array
//   LAZY: boolean to indicate whether or not range updates are enabled
//   PERSISTENT: boolean to indicate whether or not updates are persistent
//   C: struct to combine data and lazy values
//     Required Fields:
//       Data: the data type
//       Lazy: the lazy type
//     Required Functions:
//       static qdef(): returns the query default value of type Data
//       static ldef(): only required if LAZY is true, returns the lazy
//         default value of type Lazy
//       static merge(l, r): returns the values l of type Data merged with
//         r of type Data, must be associative
//       static applyLazy(l, r): only required if LAZY is false, returns the
//         value r of type Lazy applied to l of type Data, must be associative
//       static applyLazy(l, r, k): only required if LAZY is true, returns the
//         value r of type Lazy applied to l of type Data over a segment of
//         length k, must be associative
//       static mergeLazy(l, r): only required if LAZY is true, returns
//         the values l of type Lazy merged with r of type Lazy,
//         must be associative
//       static getSegmentVdef(k): returns the default value over a segment of
//         length k
//     Sample Struct: C1 supports point assignments and range sum queries,
//         where the default value of each index is 1, C2 supports
//         range assignments and range sum queries, where the default value of
//         each index is 1
//       struct C1 {
//         using Data = int;
//         using Lazy = int;
//         static Data qdef() { return 0; }
//         static Data merge(const Data &l, const Data &r) { return l + r; }
//         static Data applyLazy(const Data &l, const Lazy &r) { return r; }
//         template <class IndexType> static Data getSegmentVdef(IndexType k) {
//           return k;
//         }
//       };
//       struct C2 {
//         using Data = int;
//         using Lazy = int;
//         static Data qdef() { return 0; }
//         static Lazy ldef() { return numeric_limits<int>::min(); }
//         static Data merge(const Data &l, const Data &r) { return l + r; }
//         template <class IndexType>
//         static Data applyLazy(const Data &l, const Lazy &r, int k) {
//           return r * k;
//         }
//         static Lazy mergeLazy(const Lazy &l, const Lazy &r) { return r; }
//         template <class IndexType> static Data getSegmentVdef(IndexType k) {
//           return k;
//         }
//       };
// Constructor Arguments:
//   N: the size of the array
//   A: a vector of type C::Data
// Functions:
//   update(i, v, newRoot): updates the index i with the lazy value v
//     and creates a new version if newRoot is true
//   update(i, v, newRoot): only valid if LAZY is true, updates
//     the range [l, r] with the lazy value v and creates
//     a new version if newRoot is true
//   query(l, r, rootInd): queries the range [l, r] for the version rootInd (or
//     the latest version of rootInd is -1), and returns the aggregate value
//   bsearchPrefix(l, r, f, rootInd): returns the smallest index i in the
//     range [l, r] for the version rootInd (or the latest version of rootInd
//     is -1) such that f(query(l, i)) returns true or r + 1 if none exist
//   bsearchSuffix(l, r, f, rootInd): returns the largest index i in the
//     range [l, r] for the version rootInd (or the latest version of rootInd
//     is -1) such that f(query(i, r)) returns true of l - 1 if none exist
//   revert(rootInd): creates a new version based off of version rootInd
//   reserveNodes(k): reserves space for k nodes in the dynamic segment tree
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(1) for size constructor,
//                O(N) for iteartor and generating function constructors
//   update, query, bsearchPrefix, bsearchSuffix: O(log N) amortized unless
//     reserveNodes is called beforehand
// Memory Complexity: O(U log N) for U updates for single argument constructor,
//                    O(N + U log N) for two argument constructor
// Tested:
//   https://dmoj.ca/problem/ccc05s5 (LAZY = false, PERSISTENT = false)
//   https://dmoj.ca/problem/wc18c4s4 (LAZY = false, PERSISTENT = true,
//     bsearchPrefix, bsearchSuffix)
//   https://codeforces.com/contest/1080/problem/F
//     (LAZY = false, PERSISTENT = true)
//   https://codeforces.com/contest/915/problem/E
//     (LAZY = true, PERSISTENT = false)
//   https://www.spoj.com/problems/TTM/ (LAZY = true, PERSISTENT = true)
template <class IndexType, const bool LAZY, const bool PERSISTENT, class C>
struct DynamicSegmentTree {
#define lazy_def template <const bool _ = LAZY> typename enable_if<_>::type
#define agg_def template <const bool _ = LAZY> typename enable_if<!_>::type
  static_assert(is_integral<IndexType>::value, "IndexType must be integeral");
  using Data = typename C::Data; using Lazy = typename C::Lazy;
  template <const bool _, const int __ = 0> struct Node {
    Data val; int l, r; Node(const Data &val) : val(val), l(-1), r(-1) {}
  };
  template <const int __> struct Node<true, __> {
    Data val; Lazy lz; int l, r;
    Node(const Data &val) : val(val), lz(C::ldef()), l(-1), r(-1) {}
  };
  IndexType N; vector<int> roots; vector<Node<LAZY>> TR;
  int makeNode(int cp, IndexType tl, IndexType tr) {
    if (~cp) TR.push_back(TR[cp]);
    else TR.emplace_back(C::getSegmentVdef(tr - tl + 1));
    return int(TR.size()) - 1;
  }
  lazy_def apply(int x, IndexType tl, IndexType tr, const Lazy &v) {
    TR[x].val = C::applyLazy(TR[x].val, v, tr - tl + 1);
    TR[x].lz = C::mergeLazy(TR[x].lz, v);
  }
  agg_def apply(int x, IndexType, IndexType, const Lazy &v) {
    TR[x].val = C::applyLazy(TR[x].val, v);
  }
  lazy_def propagate(int x, IndexType tl, IndexType tr) {
    if (TR[x].lz != C::ldef()) {
      IndexType m = tl + (tr - tl) / 2; if (!~TR[x].l || PERSISTENT) {
        int nl = makeNode(TR[x].l, tl, m); TR[x].l = nl;
      }
      if (!~TR[x].r || PERSISTENT) {
        int nr = makeNode(TR[x].r, m + 1, tr); TR[x].r = nr;
      }
      apply(TR[x].l, tl, m, TR[x].lz); apply(TR[x].r, m + 1, tr, TR[x].lz);
      TR[x].lz = C::ldef();
    }
  }
  agg_def propagate(int, IndexType, IndexType) {}
  int build(const vector<Data> &A, IndexType tl, IndexType tr) {
    int x = makeNode(-1, tl, tr);
    if (tl == tr) { TR[x].val = A[tl]; return x; }
    IndexType m = tl + (tr - tl) / 2; int nl = build(A, tl, m); TR[x].l = nl;
    int nr = build(A, m + 1, tr); TR[x].r = nr;
    TR[x].val = C::merge(TR[TR[x].l].val, TR[TR[x].r].val); return x;
  }
  int update(int y, IndexType tl, IndexType tr, IndexType l, IndexType r,
             const Lazy &v, int trSz) {
    int x = !~y || (PERSISTENT && y < trSz) ? makeNode(y, tl, tr) : y;
    if (l <= tl && tr <= r) { apply(x, tl, tr, v); return x; }
    propagate(x, tl, tr); IndexType m = tl + (tr - tl) / 2;
    if (tl <= r && l <= m) {
      int nl = update(TR[x].l, tl, m, l, r, v, trSz); TR[x].l = nl;
    }
    if (m + 1 <= r && l <= tr) {
      int nr = update(TR[x].r, m + 1, tr, l, r, v, trSz); TR[x].r = nr;
    }
    TR[x].val = C::merge(~TR[x].l ? TR[TR[x].l].val
                                  : C::getSegmentVdef(m - tl + 1),
                         ~TR[x].r ? TR[TR[x].r].val
                                  : C::getSegmentVdef(tr - m));
    return x;
  }
  Data query(int x, IndexType tl, IndexType tr, IndexType l, IndexType r) {
    if (r < tl || tr < l) return C::qdef();
    if (!~x) return C::getSegmentVdef(min(r, tr) - max(l, tl) + 1);
    if (l <= tl && tr <= r) return TR[x].val;
    propagate(x, tl, tr); IndexType m = tl + (tr - tl) / 2;
    return C::merge(query(TR[x].l, tl, m, l, r),
                    query(TR[x].r, m + 1, tr, l, r));
  }
  template <class F>
  int bsearchPrefix(int y, IndexType tl, IndexType tr,
                    IndexType l, IndexType r, Data &agg, F f, IndexType &ret) {
    if (r < tl || tr < l) { ret = r + 1; return y; }
    int x = !~y ? makeNode(y, tl, tr) : y; if (tl != tr) propagate(x, tl, tr);
    if (l <= tl && tr <= r) {
      Data v = C::merge(agg, TR[x].val);
      if (!f(v)) { agg = v; ret = r + 1; return x; }
    }
    if (tl == tr) { ret = tl; return x; }
    int m = tl + (tr - tl) / 2;
    int nl = bsearchPrefix(TR[x].l, tl, m, l, r, agg, f, ret);
    TR[x].l = nl; if (ret <= r) return x;
    int nr = bsearchPrefix(TR[x].r, m + 1, tr, l, r, agg, f, ret);
    TR[x].r = nr; return x;
  }
  template <class F>
  int bsearchSuffix(int y, IndexType tl, IndexType tr,
                    IndexType l, IndexType r, Data &agg, F f, IndexType &ret) {
    if (r < tl || tr < l) { ret = l - 1; return y; }
    int x = !~y ? makeNode(y, tl, tr) : y; if (tl != tr) propagate(x, tl, tr);
    if (l <= tl && tr <= r) {
      Data v = C::merge(agg, TR[x].val);
      if (!f(v)) { agg = v; ret = l - 1; return x; }
    }
    if (tl == tr) { ret = tl; return x; }
    int m = tl + (tr - tl) / 2;
    int nr = bsearchSuffix(TR[x].r, m + 1, tr, l, r, agg, f, ret);
    TR[x].r = nr; if (l <= ret) return x;
    int nl = bsearchSuffix(TR[x].l, tl, m, l, r, agg, f, ret);
    TR[x].l = nl; return x;
  }
  DynamicSegmentTree(const vector<Data> &A) : N(A.size()) {
    if (N > 0) {
      reserveNodes(N * 2 - 1); roots.push_back(build(A, 0, N - 1));
    }
  }
  DynamicSegmentTree(IndexType N) : N(N) { roots.push_back(-1); }
  lazy_def update(IndexType l, IndexType r, const Lazy &v, bool newRoot) {
    int nr = update(roots.back(), 0, N - 1, l, r, v, TR.size());
    if (newRoot) roots.push_back(nr);
    else roots.back() = nr;
  }
  void update(IndexType i, const Lazy &v, bool newRoot) {
    int nr = update(roots.back(), 0, N - 1, i, i, v, TR.size());
    if (newRoot) roots.push_back(nr);
    else roots.back() = nr;
  }
  Data query(IndexType l, IndexType r, int rootInd = -1) {
    return query(~rootInd ? roots[rootInd] : roots.back(), 0, N - 1, l, r);
  }
  template <class F>
  IndexType bsearchPrefix(int l, int r, F f, int rootInd = -1) {
    Data agg = C::qdef(); IndexType ret = r + 1;
    int &root = ~rootInd ? roots[rootInd] : roots.back();
    root = bsearchPrefix(root, 0, N - 1, l, r, agg, f, ret);
    return ret;
  }
  template <class F>
  IndexType bsearchSuffix(int l, int r, F f, int rootInd = -1) {
    Data agg = C::qdef(); IndexType ret = l - 1;
    int &root = ~rootInd ? roots[rootInd] : roots.back();
    root = bsearchSuffix(root, 0, N - 1, l, r, agg, f, ret);
    return ret;
  }
  void revert(int rootInd) { roots.push_back(roots[rootInd]); }
  void reserveNodes(int k) { TR.reserve(k); }
#undef lazy_def
#undef agg_def
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Segment Tree supporting range updates with updates in the form of
//   adding m + b, 2m + b, 3m + b, ... to the interval [l, r],
//   and range sum queries
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of the element of the array
// Constructor Arguments:
//   N: the size of the array
//   A: a vector of type T
// Functions:
//   update(l, r, m, b): update the range [l, r] with
//     m + b, 2m + b, 3m + b, ...
//   query(l, r): returns the sum of elements in the range [l, r]
// In practice, has a small constant, but slower than FenwickTreeAffine
// Time Complexity:
//   constructor: O(N)
//   update, query: O(log N)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/acc3p4
template <class T> struct SegmentTreeAffine {
  using Pair = pair<T, T>;
  Pair add(const Pair &l, const Pair &r) {
    return Pair(l.first + r.first, l.second + r.second);
  }
  Pair mul(const Pair &l, const Pair &r) {
    return Pair(l.first * r.first, l.second * r.second);
  }
  int N; vector<Pair> TR, LZ;
  T sumTo(T k) { return k * (k + 1) / 2; }
  Pair sumBet(T l, T r) { return Pair(r - l + 1, sumTo(r) - sumTo(l - 1)); }
  void apply(int x, int tl, int tr, const Pair &v) {
    TR[x] = add(TR[x], mul(v, sumBet(tl, tr)));
    if (x < int(LZ.size())) LZ[x] = add(LZ[x], v);
  }
  void propagate(int x, int tl, int tr) {
    if (LZ[x] != Pair()) {
      int m = tl + (tr - tl) / 2; apply(x * 2, tl, m, LZ[x]);
      apply(x * 2 + 1, m + 1, tr, LZ[x]); LZ[x] = Pair();
    }
  }
  void build(const vector<T> &A, int x, int tl, int tr) {
    if (tl == tr) { TR[x] = Pair(A[tl], T()); return; }
    int m = tl + (tr - tl) / 2; build(A, x * 2, tl, m);
    build(A, x * 2 + 1, m + 1, tr); TR[x] = add(TR[x * 2], TR[x * 2 + 1]);
  }
  void update(int x, int tl, int tr, int l, int r, const Pair &v) {
    if (l <= tl && tr <= r) { apply(x, tl, tr, v); return; }
    propagate(x, tl, tr); int m = tl + (tr - tl) / 2;
    if (tl <= r && l <= m) update(x * 2, tl, m, l, r, v);
    if (m + 1 <= r && l <= tr) update(x * 2 + 1, m + 1, tr, l, r, v);
    TR[x] = add(TR[x * 2], TR[x * 2 + 1]);
  }
  Pair query(int x, int tl, int tr, int l, int r) {
    if (r < tl || tr < l) return Pair();
    if (l <= tl && tr <= r) return TR[x];
    propagate(x, tl, tr); int m = tl + (tr - tl) / 2;
    return add(query(x * 2, tl, m, l, r), query(x * 2 + 1, m + 1, tr, l, r));
  }
  SegmentTreeAffine(int N)
      : N(N), TR(N == 0 ? 0 : 1 << __lg(N * 4 - 1)),
        LZ(N == 0 ? 0 : 1 << __lg(N * 2 - 1)) {}
  SegmentTreeAffine(const vector<T> &A) : SegmentTreeAffine(A.size()) {
    if (N > 0) build(A, 1, 0, N - 1);
  }
  void update(int l, int r, T m, T b) {
    update(1, 0, N - 1, l, r, Pair((1 - l) * m + b, m));
  }
  T query(int l, int r) {
    Pair q = query(1, 0, N - 1, l, r); return q.first + q.second;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Fenwick Tree or Binary Indexed Tree supporting point updates
//   and range queries over a cumulative functor, such as sum, max, and min,
//   in any number of dimensions
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   D: the number of dimensions of the fenwick tree
//   T: the type of each element
//   Op: a struct with the cumulative operation (plus<T> by default)
//     Required Functions:
//       operator (l, r): combines the values l and r
//   Inv: a struct with the inverse cumulative operation (minus<T> by default),
//       if Op has no inverse, some dummy struct may need to be provided
//     Required Functions:
//       operator (l, r): performs the inverse operation on values l and r
// Constructor Arguments:
//   N: the size of the first dimension
//   ...args: the sizes of the second, third, fourth, etc. dimensions, with
//     the last value being the identity element (T() by default)
// Functions:
//   update(i, ...args): updates the index i in the first dimension,
//     and ...args for the next dimensions, with the last argument
//     being the value to update the index
//   queryPrefix(r, ...args): queries the cumulative value of the range [0, r]
//     in the first dimension, with the remaining arguments for the
//     next dimensions
//   query(l, r, ...args): queries the cumulative value of the range [l, r]
//     in the first dimension, and alternating arguments in args for beginning
//     and ending ranges for the next dimensions
// In practice, has a small constant, faster than segment trees and comparable
//   to the 1D version
// Time Complexity:
//   constructor: O(PI(N_i))
//   update, queryPrefix: O(PI(log(N_i)))
//   query: O(2^D * PI(log(N_i)))
// Memory Complexity: O(PI(N_i))
//   PI is the product function, N_i is the size in the ith dimension
// Tested:
//   Fuzz and Stress Tested
//   https://dmoj.ca/problem/ioi01p1
//   https://dmoj.ca/problem/gfssoc1s4
//   https://dmoj.ca/problem/ccc19s5
template <const int D, class T, class Op = plus<T>, class Inv = minus<T>>
struct FenwickTree {
  int N; vector<FenwickTree<D - 1, T, Op, Inv>> BIT; T val;
  template <class ...Args> FenwickTree(int N, Args &&...args)
      : N(N),
        BIT(N + 1, FenwickTree<D - 1, T, Op, Inv>(forward<Args>(args)...)),
        val(BIT[0].val) {}
  template <class ...Args> void update(int i, Args &&...args) {
    for (i++; i <= N; i += i & -i) BIT[i].update(forward<Args>(args)...);
  }
  template <class ...Args> T queryPrefix(int r, Args &&...args) {
    T ret = val; for (r++; r > 0; r -= r & -r)
      ret = Op()(ret, BIT[r].queryPrefix(forward<Args>(args)...));
    return ret;
  }
  template <class ...Args> T query(int l, int r, Args &&...args) {
    T ret = val; for (; l > 0; l -= l & -l)
      ret = Inv()(ret, BIT[l].query(forward<Args>(args)...));
    for (r++; r > 0; r -= r & -r)
      ret = Op()(ret, BIT[r].query(forward<Args>(args)...));
    return ret;
  }
};

template <class T, class Op, class Inv> struct FenwickTree<0, T, Op, Inv> {
  T val; FenwickTree(T qdef = T()) : val(qdef) {}
  void update(T v) { val = Op()(val, v); }
  T queryPrefix() { return val; }
  T query() { return val; }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Fenwick Tree or Binary Indexed Tree supporting range updates
//   and point queries, over a cumulative functor, such as sum, max, and min,
//   in any number of dimensions
// Indices are 0-indexed and ranges are inclusive
// Order of arguments for update differs slightly from FenwickTreeRangePoint1D
// Template Arguments:
//   D: the number of dimensions of the fenwick tree
//   T: the type of each element
//   Op: a struct with the cumulative operation (plus<T> by default)
//     Required Functions:
//       operator (l, r): combines the values l and r, must be associative and
//         commutatitve
//   Inv: a struct with the inverse cumulative operation (minus<T> by default),
//       if Op has no inverse, some dummy struct may need to be provided
//     Required Functions:
//       operator (l, r): performs the inverse operation on values l and r
// Constructor Arguments:
//   N: the size of the first dimension
//   ...args: the sizes of the second, third, fourth, etc. dimensions, with
//     the last value being the identity element (T() by default)
// Functions:
//   updateSuffix(v, l, ...args): updates the suffix starting from index l in
//     the first dimension, and ...args the next dimensions with the value a
//     being updated to op(a, v)
//   update(v, l, r, ...args): updates the range [l, r] in the first
//     dimension, and alternating arguments in args for beginning and
//     ending ranges for the next dimensions with the value a being updated to
//     op(a, v)
//   get(i, ...args): queries the value of index i in the first
//     dimension, and ...args for the next dimensions
// In practice, has a small constant, faster than segment trees and comparable
//   to the 1D version
// Time Complexity:
//   constructor: O(PI(N_i))
//   update: O(2^D * PI(log(N_i)))
//   updateSuffix, getValue: O(PI(log(N_i)))
// Memory Complexity: O(PI(N_i))
//   PI is the product function, N_i is the size in the ith dimension,
//   and D is the number of dimensions
// Tested:
//   Fuzz and Stress Tested
//   https://mcpt.ca/problem/adifferenceproblem
template <const int D, class T, class Op = plus<T>, class Inv = minus<T>>
struct FenwickTreeRangePoint {
  int N; vector<FenwickTreeRangePoint<D - 1, T, Op, Inv>> BIT; T val;
  template <class ...Args> FenwickTreeRangePoint(int N, Args &&...args)
    : N(N), BIT(N + 1, FenwickTreeRangePoint<D - 1, T, Op, Inv>(
          forward<Args>(args)...)),
      val(BIT[0].val) {}
  template <class ...Args> void updateSuffix(T v, int l, Args &&...args) {
    for (l++; l <= N; l += l & -l)
      BIT[l].updateSuffix(v, forward<Args>(args)...);
  }
  template <class ...Args> void update(T v, int l, int r, Args &&...args) {
    for (l++; l <= N; l += l & -l) BIT[l].update(v, forward<Args>(args)...);
    for (r += 2; r <= N; r += r & -r)
      BIT[r].update(Inv()(val, v), forward<Args>(args)...);
  }
  template <class ...Args> T get(int i, Args &&...args) {
    T ret = val; for (i++; i > 0; i -= i & -i)
      ret = Op()(ret, BIT[i].get(forward<Args>(args)...));
    return ret;
  }
};

template <class T, class Op, class Inv>
struct FenwickTreeRangePoint<0, T, Op, Inv> {
  T val; FenwickTreeRangePoint(T qdef = T()) : val(qdef) {}
  void updateSuffix(T v) { val = Op()(val, v); }
  void update(T v) { val = Op()(val, v); }
  T get() { return val; }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "../../PolicyBasedDataStructures.h"
#include "../../sqrt/SqrtBufferSimple.h"
#include "../../sqrt/SqrtBuffer.h"
using namespace std;
using namespace __gnu_pbds;

// A collection of sparse fenwick trees implemented in various methods
// Some are only sparse in 1 dimension, others are sparse in 2 dimensions
// Certain implementations only allow for increments and decrements of 1
// In general, offline fenwick trees are faster than sqrt implementations,
//   which are faster than the pbds implementations

// Sparse Fenwick Tree supporting point updates (with any value)
//   and range queries over a cumulative function or functor,
//   such as sum, max, and min, in 2 dimensions (sparse in 1 dimension)
// All update indices must be known beforehand
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   IndexType: the type of the index of the second dimension of the array
//   Op: a struct with the cumulative operation (plus<T> by default)
//     Required Functions:
//       operator (l, r): combines the values l and r, must be associative and
//         commutatitve
// Constructor Arguments:
//   N: the size of the first dimension of the array
//   updateInds: a vector of pairs of ints and IndexType containing the indices
//     for both dimensions for each update
//   qdef: the identity element of the operation
//   op: an instance of the Op struct
// Functions:
//   update(i, j, v): modifies the value A[i][j] to op(A[i][j], v), (i, j)
//     must be an element of updateInds
//   query(d, r): queries the cumulative value of the range [0, d] in the
//     first dimension and [0, r] in the second dimension
//   query(d, l, r, inv): queries the cumulative value of the range [0, d]
//     in the first dimension and [l, r] in the second dimension, where inv is
//     the inverse of op (minus<T>() by default)
//   query(u, d, l, r, inv): queries the cumulative value of the range [u, d]
//     in the first dimension and [l, r] in the second dimension, where inv is
//     the inverse of op (minus<T>() by default)
// In practice, has a small constant
// Time Complexity:
//   constructor: O(Q (log Q + log N) + N) for U updates
//   update, query: O(log N log U) for U updates
// Memory Complexity: O(N + U log N) for U updates
// Tested:
//   https://codeforces.com/contest/1093/problem/E
//   https://codeforces.com/contest/1523/problem/G
template <class T, class IndexType, class Op = plus<T>>
struct OfflineSemiSparseFenwickTree2D {
  static_assert(is_integral<IndexType>::value, "IndexType must be integeral");
  int N; vector<int> st, cnt; vector<IndexType> inds; vector<T> BIT;
  T qdef; Op op;
  int getInd(int i, IndexType j) {
    return upper_bound(inds.begin() + st[i], inds.begin() + st[i] + cnt[i], j)
        - inds.begin() - st[i];
  }
  OfflineSemiSparseFenwickTree2D(int N,
                                 vector<pair<int, IndexType>> updateInds,
                                 T qdef = T(), Op op = Op())
      : N(N), st(N + 1, 0), cnt(N + 1, 0), qdef(qdef), op(op) {
    sort(updateInds.begin(), updateInds.end(),
         [&] (const pair<int, IndexType> &a, const pair<int, IndexType> &b) {
      return a.second < b.second;
    });
    vector<IndexType> last(N + 1, IndexType());
    for (auto &&u : updateInds) for (int i = u.first + 1; i <= N; i += i & -i)
      if (cnt[i] == 0 || u.second != last[i]) { cnt[i]++; last[i] = u.second; }
    for (int i = 1; i <= N; i++) st[i] = st[i - 1] + cnt[i - 1];
    inds.resize(st[N] + cnt[N]); BIT.resize(st[N] + cnt[N], qdef);
    fill(cnt.begin(), cnt.end(), 0); for (auto &&u : updateInds)
      for (int i = u.first + 1; i <= N; i += i & -i)
        if (cnt[i] == 0 || u.second != inds[st[i] + cnt[i] - 1])
          inds[st[i] + cnt[i]++] = u.second;
  }
  void update(int i, IndexType j, T v) {
    for (i++; i <= N; i += i & -i)
      for (int s = st[i], c = cnt[i], y = getInd(i, j); y <= c; y += y & -y)
        BIT[s + y - 1] = op(BIT[s + y - 1], v);
  }
  T query(int d, IndexType r) {
    T ret = qdef; for (d++; d > 0; d -= d & -d)
      for (int s = st[d], y = getInd(d, r); y > 0; y -= y & -y)
        ret = op(ret, BIT[s + y - 1]);
    return ret;
  }
  template <class Inv = minus<T>>
  T query(int d, IndexType l, IndexType r, Inv inv = Inv()) {
    return inv(query(d, r), query(d, l - 1));
  }
  template <class Inv = minus<T>>
  T query(int u, int d, IndexType l, IndexType r, Inv inv = Inv()) {
    return inv(query(d, l, r), query(u - 1, l, r));
  }
};

// Sparse Fenwick Tree supporting point updates (with any value)
//   and range queries over a cumulative function or functor,
//   such as sum, max, and min, in 2 dimensions (sparse in 2 dimensions)
// All update indices must be known beforehand
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   IndexType1: the type of the index of the first dimension of the array
//   IndexType2: the type of the index of the second dimension of the array
//   Op: a struct with the cumulative operation (plus<T> by default)
//     Required Functions:
//       operator (l, r): combines the values l and r, must be associative and
//         commutatitve
// Constructor Arguments:
//   updateInds: a vector of pairs of IndexType1 and IndexType2 containing
//     the indices for both dimensions for each update
//   qdef: the identity element of the operation
//   op: an instance of the Op struct
// Functions:
//   update(i, j, v): modifies the value A[i][j] to op(A[i][j], v), (i, j)
//     must be an element of updateInds
//   query(d, r): queries the cumulative value of the range [0, d] in the
//     first dimension and [0, r] in the second dimension
//   query(d, l, r, inv): queries the cumulative value of the range [0, d]
//     in the first dimension and [l, r] in the second dimension, where inv is
//     the inverse of op (minus<T>() by default)
//   query(u, d, l, r, inv): queries the cumulative value of the range [u, d]
//     in the first dimension and [l, r] in the second dimension, where inv is
//     the inverse of op (minus<T>() by default)
// In practice, has a small constant
// Time Complexity:
//   constructor: O(U log U) for U updates
//   update, query: O((log U)^2) for U updates
// Memory Complexity: O(U log U) for U updates
// Tested:
//   https://judge.yosupo.jp/problem/point_add_rectangle_sum
//   https://codeforces.com/contest/1093/problem/E
//   https://codeforces.com/contest/1523/problem/G
template <class T, class IndexType1, class IndexType2, class Op = plus<T>>
struct OfflineSparseFenwickTree2D {
  static_assert(is_integral<IndexType1>::value,
                "IndexType1 must be integeral");
  static_assert(is_integral<IndexType2>::value,
                "IndexType2 must be integeral");
  int U; vector<int> st, cnt; vector<IndexType1> inds1;
  vector<IndexType2> inds2; vector<T> BIT; T qdef; Op op;
  int getInd1(IndexType1 i) {
    return upper_bound(inds1.begin(), inds1.end(), i) - inds1.begin();
  }
  int getInd2(int i, IndexType2 j) {
    return upper_bound(inds2.begin() + st[i],
                       inds2.begin() + st[i] + cnt[i], j)
        - inds2.begin() - st[i];
  }
  OfflineSparseFenwickTree2D(vector<pair<IndexType1, IndexType2>> updateInds,
                             T qdef = T(), Op op = Op())
      : inds1(updateInds.size()), qdef(qdef), op(op) {
    sort(updateInds.begin(), updateInds.end(),
         [&] (const pair<IndexType1, IndexType2> &a,
              const pair<IndexType1, IndexType2> &b) {
      return a.second < b.second;
    });
    for (int i = 0; i < int(updateInds.size()); i++)
      inds1[i] = updateInds[i].first;
    sort(inds1.begin(), inds1.end());
    inds1.erase(unique(inds1.begin(), inds1.end()), inds1.end());
    U = int(inds1.size()); st.assign(U + 1, 0); cnt.assign(U + 1, 0);
    vector<IndexType2> last(U + 1, IndexType2()); for (auto &&u : updateInds)
      for (int i = getInd1(u.first); i <= U; i += i & -i)
        if (cnt[i] == 0 || u.second != last[i]) {
          cnt[i]++; last[i] = u.second;
        }
    for (int i = 1; i <= U; i++) st[i] = st[i - 1] + cnt[i - 1];
    inds2.resize(st[U] + cnt[U]); BIT.resize(st[U] + cnt[U], qdef);
    fill(cnt.begin(), cnt.end(), 0); for (auto &&u : updateInds)
      for (int i = getInd1(u.first); i <= U; i += i & -i)
        if (cnt[i] == 0 || u.second != inds2[st[i] + cnt[i] - 1])
          inds2[st[i] + cnt[i]++] = u.second;
  }
  void update(IndexType1 i, IndexType2 j, T v) {
    for (int x = getInd1(i); x <= U; x += x & -x)
      for (int s = st[x], c = cnt[x], y = getInd2(x, j); y <= c; y += y & -y)
        BIT[s + y - 1] = op(BIT[s + y - 1], v);
  }
  T query(IndexType1 d, IndexType2 r) {
    T ret = qdef; for (int x = getInd1(d); x > 0; x -= x & -x)
      for (int s = st[x], y = getInd2(x, r); y > 0; y -= y & -y)
        ret = op(ret, BIT[s + y - 1]);
    return ret;
  }
  template <class Inv = minus<T>>
  T query(IndexType1 d, IndexType2 l, IndexType2 r, Inv inv = Inv()) {
    return inv(query(d, r), query(d, l - 1));
  }
  template <class Inv = minus<T>>
  T query(IndexType1 u, IndexType1 d, IndexType2 l, IndexType2 r,
          Inv inv = Inv()) {
    return inv(query(d, l, r), query(u - 1, l, r));
  }
};

// Sparse Fenwick Tree supporting point updates (with value 1 and -1)
//   and range sum queries in 2 dimensions
//   using SqrtBufferSimple (sparse in 1 dimension)
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   IndexType: the type of the index of the second dimension of the array
// Constructor Arguments:
//   N: the size of the first dimension of the array
//   SCALE: the value to scale sqrt by (for SqrtBufferSimple)
// Functions:
//   add(i, j): add 1 to the value A[i][j]
//   rem(i, j): subtract 1 from the value A[i][j]
//   query(d, r): queries the sum of the range [0, d] in the first dimension
//     and [0, r] in the second dimension
//   query(d, l, r): queries the sum of the range [0, d] in the first dimension
//     and [l, r] in the second dimension
//   query(u, d, l, r): queries the sum of the range [u, d] in the first
//     dimension and [l, r] in the second dimension
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(N)
//   add, rem: O(log N) amortized
//   query: O(log N sqrt U) amortized for U updates
// Memory Complexity: O(N + U log N) for U updates
// Tested:
//   https://codeforces.com/contest/1093/problem/E
//   https://dmoj.ca/problem/dmopc19c7p5
template <class IndexType> struct SemiSparseFenwickTree2DSimple {
  static_assert(is_integral<IndexType>::value, "IndexType must be integeral");
  int N; vector<SqrtBufferSimple<IndexType>> IN, OUT;
  SemiSparseFenwickTree2DSimple(int N, double SCALE = 1)
      : N(N), IN(N + 1, SqrtBufferSimple<IndexType>(SCALE)),
        OUT(N + 1, SqrtBufferSimple<IndexType>(SCALE)) {}
  void add(int i, IndexType j) {
    for (i++; i <= N; i += i & -i) IN[i].insert(j);
  }
  void rem(int i, IndexType j) {
    for (i++; i <= N; i += i & -i) OUT[i].insert(j);
  }
  int query(int d, IndexType r) {
    int ret = 0; for (d++; d > 0; d -= d & -d)
      ret += IN[d].aboveInd(r) - OUT[d].aboveInd(r);
    return ret;
  }
  int query(int d, IndexType l, IndexType r) {
    int ret = 0; for (d++; d > 0; d -= d & -d)
      ret += IN[d].count(l, r) - OUT[d].count(l, r);
    return ret;
  }
  int query(int u, int d, IndexType l, IndexType r) {
    return query(d, l, r) - query(u - 1, l, r);
  }
};

// Sparse Fenwick Tree supporting point updates (with any value)
//   and range sum queries in 2 dimensions using SqrtBuffer (sparse in
//   1 dimension)
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   IndexType: the type of the index of the second dimension of the array
// Constructor Arguments:
//   N: the size of the first dimension of the array
//   SCALE: the value to scale sqrt by (for SqrtBuffer)
// Functions:
//   update(i, j, v): add v to the value A[i][j]
//   query(d, r): queries the sum of the range [0, d] in the first dimension
//     and [0, r] in the second dimension
//   query(d, l, r): queries the sum of the range [0, d] in the first dimension
//     and [l, r] in the second dimension
//   query(u, d, l, r): queries the sum of the range [u, d] in the first
//     dimension and [l, r] in the second dimension
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(N)
//   update: O(log N) amortized
//   query: O(log N sqrt U) amortized for U updates
// Memory Complexity: O(N + U log N) for U updates
// Tested:
//   https://dmoj.ca/problem/apio19p3
//   https://dmoj.ca/problem/ioi01p1
template <class T, class IndexType> struct SemiSparseFenwickTree2D {
  static_assert(is_integral<IndexType>::value, "IndexType must be integeral");
  int N; vector<SqrtBuffer<T, IndexType>> BIT;
  SemiSparseFenwickTree2D(int N, double SCALE = 1)
      : N(N), BIT(N + 1, SqrtBuffer<T, IndexType>(SCALE)) {}
  void update(int i, IndexType j, T v) {
    for (i++; i <= N; i += i & -i) BIT[i].emplace(j, v);
  }
  T query(int d, IndexType r) {
    T ret = T(); for (d++; d > 0; d -= d & -d) ret += BIT[d].aboveInd(r);
    return ret;
  }
  T query(int d, IndexType l, IndexType r) {
    T ret = T(); for (d++; d > 0; d -= d & -d) ret += BIT[d].count(l, r);
    return ret;
  }
  T query(int u, int d, IndexType l, IndexType r) {
    return query(d, l, r) - query(u - 1, l, r);
  }
};

// Sparse Fenwick Tree supporting point updates (with value 1 and -1)
//   and range sum queries in 2 dimensions
//   using pbds hash_table for the first dimension,
//   and SqrtBufferSimple for the second dimension (sparse in 2 dimensions)
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   IndexType1: the type of the index of the first dimension of the array
//   IndexType2: the type of the index of the second dimension of the array
//   Container: the container to represent an unordered mapping from
//     IndexType1 to pair of SqrtBufferSimple<IndexType2>
// Constructor Arguments:
//   N: the size of the first dimension of the array
//   SCALE: the value to scale sqrt by (for SqrtBufferSimple)
// Functions:
//   add(i, j): add 1 to the value A[i][j]
//   rem(i, j): subtract 1 from the value A[i][j]
//   query(d, r): queries the sum of the range [0, d] in the first dimension
//     and [0, r] in the second dimension
//   query(d, l, r): queries the sum of the range [0, d] in the first dimension
//     and [l, r] in the second dimension
//   query(u, d, l, r): queries the sum of the range [u, d] in the first
//     dimension and [l, r] in the second dimension
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(1)
//   add, rem: O(log N) amortized
//   query: O(log N sqrt U) amortized for U updates
// Memory Complexity: O(U log N) for U updates
// Tested:
//   https://codeforces.com/contest/1093/problem/E
template <class IndexType1, class IndexType2,
          class Container = hashmap<
              IndexType1, pair<SqrtBufferSimple<IndexType2>,
                               SqrtBufferSimple<IndexType2>>>>
struct SparseFenwickTree2DSimple {
  static_assert(is_integral<IndexType1>::value,
                "IndexType1 must be integeral");
  static_assert(is_integral<IndexType2>::value,
                "IndexType2 must be integeral");
  IndexType1 N; double SCALE; Container BIT;
  SparseFenwickTree2DSimple(IndexType1 N, double SCALE = 1)
      : N(N), SCALE(SCALE) {}
  void add(IndexType1 i, IndexType2 j) {
    i++; for (IndexType1 x = i; x <= N; x += x & -x) {
      auto it = BIT.find(x); if (it == BIT.end()) {
        BIT[x] = make_pair(SqrtBufferSimple<IndexType2>(SCALE),
                           SqrtBufferSimple<IndexType2>(SCALE));
      }
      BIT[x].first.insert(j);
    }
  }
  void rem(IndexType1 i, IndexType2 j) {
    i++; for (IndexType1 x = i; x <= N; x += x & -x) {
      auto it = BIT.find(x); if (it == BIT.end()) {
        BIT[x] = make_pair(SqrtBufferSimple<IndexType2>(SCALE),
                           SqrtBufferSimple<IndexType2>(SCALE));
      }
      BIT[x].second.insert(j);
    }
  }
  int query(IndexType1 d, IndexType2 r) {
    d++; int ret = 0; for (IndexType1 x = d; x > 0; x -= x & -x) {
      auto it = BIT.find(x);
      ret += it->second.first.aboveInd(r) - it->second.second.aboveInd(r);
    }
    return ret;
  }
  int query(IndexType1 d, IndexType2 l, IndexType2 r) {
    d++; int ret = 0; for (IndexType1 x = d; x > 0; x -= x & -x) {
      auto it = BIT.find(x);
      ret += it->second.first.count(l, r) - it->second.second.count(l, r);
    }
    return ret;
  }
  int query(IndexType1 u, IndexType1 d, IndexType2 l, IndexType2 r) {
    return query(d, l, r) - query(u - 1, l, r);
  }
};

// Sparse Fenwick Tree supporting point updates (with any value)
//   and range sum queries in 2 dimensions
//   using pbds hash_table for the first dimension,
//   and SqrtBuffer for the second dimension (sparse in 2 dimensions)
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   IndexType1: the type of the index of the first dimension of the array
//   IndexType2: the type of the index of the second dimension of the array
//   Container: the container to represent an unordered mapping from
//     IndexType1 to SqrtBuffer<IndexType2, T>
// Constructor Arguments:
//   N: the size of the first dimension of the array
//   SCALE: the value to scale sqrt by (for SqrtBuffer)
// Functions:
//   update(i, j, v): add v to the value A[i][j]
//   query(d, r): queries the sum of the range [0, d] in the first dimension
//     and [0, r] in the second dimension
//   query(d, l, r): queries the sum of the range [0, d] in the first dimension
//     and [l, r] in the second dimension
//   query(u, d, l, r): queries the sum of the range [u, d] in the first
//     dimension and [l, r] in the second dimension
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(1)
//   update: O(log N) amortized
//   query: O(log N sqrt U) amortized for U updates
// Memory Complexity: O(U log N) for U update
// Tested:
//   https://dmoj.ca/problem/fallingsnowflakes
template <class T, class IndexType1, class IndexType2,
          class Container = hashmap<IndexType1, SqrtBuffer<IndexType2, T>>>
struct SparseFenwickTree2D {
  static_assert(is_integral<IndexType1>::value,
                "IndexType1 must be integeral");
  static_assert(is_integral<IndexType2>::value,
                "IndexType2 must be integeral");
  IndexType1 N; double SCALE; Container BIT;
  SparseFenwickTree2D(IndexType1 N, double SCALE = 1) : N(N), SCALE(SCALE) {}
  void update(IndexType1 i, IndexType2 j, T v) {
    i++; for (IndexType1 x = i; x <= N; x += x & -x) {
      auto it = BIT.find(x); if (it == BIT.end()) {
        (BIT[x] = SqrtBuffer<IndexType2, T>(SCALE)).emplace(j, v);
      } else it->second.emplace(j, v);
    }
  }
  T query(IndexType1 d, IndexType2 r) {
    d++; T ret = T(); for (IndexType1 x = d; x > 0; x -= x & -x) {
      auto it = BIT.find(x);
      if (it != BIT.end()) ret += it->second.aboveInd(r);
    }
    return ret;
  }
  T query(IndexType1 d, IndexType2 l, IndexType2 r) {
    d++; T ret = T(); for (IndexType1 x = d; x > 0; x -= x & -x) {
      auto it = BIT.find(x);
      if (it != BIT.end()) ret += it->second.count(l, r);
    }
    return ret;
  }
  T query(IndexType1 u, IndexType1 d, IndexType2 l, IndexType2 r) {
    return query(d, l, r) - query(u - 1, l, r);
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "FenwickTree1D.h"
#include "../../../math/LagrangePolynomialInterpolation.h"
using namespace std;

// FenwickTree supporting range updates with updates in the form of
//   adding v, 2^k v, 3^k v, ... to the range [l, r], and range sum queries
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
// Constructor Arguments:
//   N: the size of the array
// Functions:
//   update(l, r, v, k): add v, 2^k v, 3^k v, ... to the range [l, r]
//   query(r): queries the sum of the range [0, r]
//   query(l, r): queries the sum of the range [l, r]
// In practice, has a small constant, faster than segment trees
// Time Complexity:
//   constructor: O(N + K^3)
//   update: O(K (K + log N))
//   rsq: O(K log N)
// Memory Complexity: O(KN + K^3)
// Tested:
//   https://dmoj.ca/problem/tle17c3p6
template <class T> struct FenwickTreePowerSum {
  int N, K; vector<FenwickTree1D<T>> FT; vector<T> fact;
  vector<vector<T>> cof, pascal; vector<vector<vector<T>>> sumDiff;
  FenwickTreePowerSum(int N, int K)
      : N(N), K(K), FT(K + 2, FenwickTree1D<T>(N)), fact(K + 2, T(1)),
        cof(K + 1), pascal(K + 2), sumDiff(K + 1, vector<vector<T>>(K + 2)) {
    for (int i = 1; i < K + 2; i++) fact.push_back(fact.back() * T(i));
    for (int i = 0; i < K + 1; i++) {
      vector<pair<T, T>> P; P.reserve(i + 2); P.emplace_back(T(), T());
      for (int j = 1; j < i + 2; j++) {
        T p = T(1); for (int k = 1; k <= i; k++) p *= T(j);
        P.emplace_back(T(j), P.back().second + fact[i + 1] * p);
      }
      cof[i] = lagrangePolynomialInterpolation(P);
    }
    for (int k = 0; k < K + 1; k++) for (int i = 0; i < K + 2; i++)
      sumDiff[k][i].resize(K + 2 - i, T());
    for (int i = 0; i < K + 2; i++) {
      pascal[i].resize(i + 1, T(1)); for (int j = 1; j < i; j++)
        pascal[i][j] = pascal[i - 1][j - 1] + pascal[i - 1][j];
    }
    for (int k = 0; k < K + 1; k++) for (int i = 0; i < k + 2; i++) {
      T sign = T(1); for (int j = 0; j <= i; j++, sign = -sign) {
        T f = fact[K + 1] / fact[k + 1];
        sumDiff[k][i - j][j] += cof[k][i] * pascal[i][j] * f * sign;
      }
    }
  }
  void update(int l, int r, T v, int k) {
    assert(k <= K); T s = T(l - 1), len = T(r - l + 1);
    for (int i = 0; i < k + 2; i++) {
      T mul = sumDiff[k][i][k + 1 - i];
      for (int j = k - i; j >= 0; j--) mul = mul * s + sumDiff[k][i][j];
      FT[i].update(l, v * mul); FT[i].update(r + 1, -v * mul);
    }
    T mul = cof[k][k + 1];
    for (int i = k; i >= 0; i--) mul = mul * len + cof[k][i];
    mul = mul * (fact[K + 1] / fact[k + 1]); FT[0].update(r + 1, v * mul);
  }
  T query(int r) {
    T ret = FT[K + 1].query(r);
    for (int i = K; i >= 0; i--) ret = ret * r + FT[i].query(r);
    return ret / fact[K + 1];
  }
  T query(int l, int r) { return query(r) - query(l - 1); }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "FenwickTreeRangePoint.h"
using namespace std;

// Supports range additions and range sum queries in 2D where
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of the element of the grid
// Constructor Arguments:
//   N: the size of the first dimension of the array
//   M: the size of the second dimension of the array
// Functions:
//   update(u, d, l, r, v): adds v to all elements in the range [u, d] in the
//     first dimension and [l, r] in the second dimension
//   query(d, r): queries the sum of the range [0, d] in the
//     first dimension and [0, r] in the second dimension
//   query(d, l, r): queries the sum of the range [0, d] in the
//     first dimension and [l, r] in the second dimension
//   query(u, d, l, r): queries the sum of the range [u, d] in the
//     first dimension and [l, r] in the second dimension
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(NM)
//   update, query: O(log N log M)
// Memory Complexity: O(NM)
// Tested:
//   Fuzz and Stress Tested
template <class T> struct RangeAddRangeSum2D {
  int N, M; vector<FenwickTreeRangePoint<2, T>> FT;
  RangeAddRangeSum2D(int N, int M)
      : N(N), M(M), FT(4, FenwickTreeRangePoint<2, T>(N, M)) {}
  void update(int u, int d, int l, int r, const T &v) {
    FT[0].update(v, u, d, l, r);
    FT[1].update(-v * T(l - 1), u, d, l, r);
    FT[1].update(v * T(r - l + 1), u, d, r + 1, M - 1);
    FT[2].update(-v * T(u - 1), u, d, l, r);
    FT[2].update(v * T(d - u + 1), d + 1, N - 1, l, r);
    FT[3].update(v * T(u - 1) * T(l - 1), u, d, l, r);
    FT[3].update(-v * T(u - 1) * T(r - l + 1), u, d, r + 1, M - 1);
    FT[3].update(-v * T(l - 1) * T(d - u + 1), d + 1, N - 1, l, r);
    FT[3].updateSuffix(v * T(d - u + 1) * T(r - l + 1), d + 1, r + 1);
  }
  T query(int d, int r) {
    return FT[0].get(d, r) * T(d) * T(r) + FT[1].get(d, r) * T(d)
        + FT[2].get(d, r) * T(r) + FT[3].get(d, r);
  }
  T query(int d, int l, int r) { return query(d, r) - query(d, l - 1); }
  T query(int u, int d, int l, int r) {
    return query(d, l, r) - query(u - 1, l, r);
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Fenwick Tree or Binary Indexed Tree supporting point updates
//   and range queries over a cumulative function or functor,
//   such as max and min, in K dimensions, even if K is not known at
//   compile time
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   Op: a struct with the cumulative operation (plus<T> by default)
//     Required Functions:
//       operator (l, r): combines the values l and r, must be associative and
//         commutatitve
// Constructor Arguments:
//   D: a vector of length K with the size of each dimension
//   A: the flattened array, represented by a vector of type T
// Functions:
//   values(inv): returns a vector of the fenwick tree decomposed into a
//     flattened array, where inv is the inverse of op (minus<T>() by default)
//   update(ind, v): updates the indices specified by the vector ind holding
//     the index for each dimension, with the value a being modified to
//     op(a, v)
//   query(ind): queries the cumulative value of the prefix bounded by the
//     indices specified by the vector ind holding the index for each
//     dimension, where inv is the inverse of Op (minus<T> by default)
//   query(ind, inv): queries the cumulative value of the indices specified by
//     the vector ind holding a pair specifying the range for each dimension,
//     where inv is the inverse of Op (minus<T> by default)
// In practice, has a small constant, slower than FenwickTree.h
// Time Complexity:
//   constructor: O(PI(N_i))
//   update, query(ind): O(PI(log(N_i)))
//   query(ind, inv): O(2^K * PI(log(N_i)))
// Memory Complexity: O(PI(N_i))
//   PI is the product function, N_i is the size in the ith dimension
// Tested:
//   https://dmoj.ca/problem/inaho2
//   https://dmoj.ca/problem/inaho7
template <class T, class Op = plus<T>> struct FenwickTreeKD {
  int K; vector<int> D, suf; vector<T> BIT; T qdef; Op op;
  int dig(int i, int d) { return i % suf[d] / suf[d + 1]; }
  int nxt(int i) { i++; i += i & -i; return --i; }
  FenwickTreeKD(const vector<int> &D, T qdef = T(), Op op = Op())
      : K(D.size()), D(D), suf(K + 1, 1), qdef(qdef), op(op) {
    for (int i = K - 1; i >= 0; i--) suf[i] = suf[i + 1] * D[i];
    BIT.resize(suf[0], qdef);
  }
  FenwickTreeKD(const vector<int> &D, vector<T> A, T qdef = T(), Op op = Op())
      : K(D.size()), D(D), suf(K + 1, 1), BIT(move(A)), qdef(qdef), op(op) {
    for (int i = K - 1; i >= 0; i--) suf[i] = suf[i + 1] * D[i];
    assert(int(BIT.size()) == suf[0]);
    for (int d = 0; d < K; d++) for (int i = 0; i < suf[0]; i++)
      if (nxt(dig(i, d)) < D[d]) {
        int j = i + (nxt(dig(i, d)) - dig(i, d)) * suf[d + 1];
        BIT[j] = op(BIT[j], BIT[i]);
      }
  }
  template <class Inv = minus<T>> vector<T> values(Inv inv = Inv()) {
    vector<T> ret = BIT; for (int d = K - 1; d >= 0; d--)
      for (int i = suf[0] - 1; i >= 0; i--) if (nxt(dig(i, d)) < D[d]) {
        int j = i + (nxt(dig(i, d)) - dig(i, d)) * suf[d + 1];
        ret[j] = inv(ret[j], ret[i]);
      }
    return ret;
  }
  void update(const vector<int> &ind, const T &v, int d = 0, int pos = 0) {
    if (d == K) BIT[pos] = op(BIT[pos], v);
    else for (int i = ind[d] + 1; i <= D[d]; i += i & -i)
      update(ind, v, d + 1, pos + (i - 1) * suf[d + 1]);
  }
  T query(const vector<int> &ind, int d = 0, int pos = 0) {
    T ret = qdef; if (d == K) ret = op(ret, BIT[pos]);
    else for (int r = ind[d] + 1; r > 0; r -= r & r)
      ret = op(ret, query(ind, d + 1, pos + (r - 1) * suf[d + 1]));
    return ret;
  }
  template <class Inv = minus<T>>
  T query(const vector<pair<int, int>> &ind, Inv inv = Inv(),
          int d = 0, int pos = 0) {
    T ret = qdef; if (d == K) ret = op(ret, BIT[pos]);
    else {
      for (int l = ind[d].first; l > 0; l -= l & -l)
        ret = inv(ret, query(ind, inv, d + 1, pos + (l - 1) * suf[d + 1]));
      for (int r = ind[d].second + 1; r > 0; r -= r & -r)
        ret = op(ret, query(ind, inv, d + 1, pos + (r - 1) * suf[d + 1]));
    }
    return ret;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "FenwickTreeRangePoint1D.h"
using namespace std;

// Fenwick Tree supporting range updates with updates in the form of
//   adding m + b, 2m + b, 3m + b, ... to the range [l, r],
//   and range sum queries
// Indices are 0-indexed and ranges are inclusive
// numeric_limits<T>::max() * N * N must not overflow
// Template Arguments:
//   T: the type of each element
//   A: a vector of type T, memory is saved if this is moved and has
//     a capacity of N + 1
// Constructor Arguments:
//   N: the size of the array
// Functions:
//   update(l, r, m, b): add m + b, 2m + b, 3m + b, ... to the range [l, r]
//   query(r): queries the sum of the range [0, r]
//   query(l, r): queries the sum of the range [l, r]
// In practice, has a small constant, faster than SegmentTreeAffine
// Time Complexity:
//   constructor: O(N)
//   update, query: O(log N)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/acc3p4
template <class T> struct FenwickTreeAffine {
  vector<FenwickTreeRangePoint1D<T>> FT;
  FenwickTreeAffine(int N) : FT(3, FenwickTreeRangePoint1D<T>(N)) {}
  FenwickTreeAffine(vector<T> A) {
    for (auto &&a : A) a *= T(2);
    FT.reserve(3); FT.emplace_back(move(A));
    FT.emplace_back(FT[0].N); FT.emplace_back(FT[0].N);
  }
  T query(int r) {
    return (FT[2].get(r) * T(r) * T(r) + FT[1].get(r) * T(r)
            + FT[0].get(r)) / T(2);
  }
  T query(int l, int r) { return query(r) - query(l - 1); }
  void update(int l, int r, T m, T b) {
    FT[2].update(l, r, m);
    FT[1].update(l, r, m * (T(1) - T(l - 1) * T(2)) + b * T(2));
    T con1 = m * T(l - 1) * T(l - 2) + b * T(1 - l) * T(2);
    T con2 = m * T(r - l + 1) * T(r - l + 2) + b * (r - l + 1) * T(2);
    FT[0].update(l, con1); FT[0].update(r + 1, con2 - con1);
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Fenwick Tree or Binary Indexed Tree supporting range updates
//   and point queries, over a cumulative functor, such as sum, max, and min,
//   in 1 dimension
// Indices are 0-indexed and ranges are inclusive
// Order of arguments for update differs slightly from FenwickTreeRangePoint
// Template Arguments:
//   T: the type of each element
//   Op: a struct with the cumulative operation (plus<T> by default)
//     Required Functions:
//       operator (l, r): combines the values l and r, must be associative and
//         commutatitve
// Constructor Arguments:
//   N: the size of the array
//   A: a vector of type T, memory is saved if this is moved and has
//     a capacity of N + 1
//   qdef: the identity element of the operation
//   op: an instance of the Op struct
//   inv(l, r): a functor or function that performs the inverse operation of
//     op (minus<T>() by default)
// Functions:
//   values(inv): returns a vector of the fenwick tree decomposed into an
//     array, where inv is the inverse of op (minus<T>() by default)
//   update(l, v): sets A[i] to op(A[i], v) for the range [l..N - 1]
//   update(l, r, v, inv): add A[i] to op(A[i], v) for the range [l..r]
//   get(i): queries the value at index i, where inv is the inverse of
//     op (minus<T>() by default)
//   bsearch(v, cmp): returns the first index where cmp(A[i], v)
//     returns false, or N if no such index exists
// In practice, this version performs as well as the multidimensional version
// In practice, has a small constant, faster than segment trees
// Time Complexity:
//   constructor, values: O(N)
//   update, get: O(log N)
// Memory Complexity: O(N)
// Tested:
//   Fuzz and Stress Tested
//   https://mcpt.ca/problem/asquirrelproblem
//   https://codeforces.com/contest/1404/problem/C
template <class T, class Op = plus<T>> struct FenwickTreeRangePoint1D {
  int N; vector<T> BIT; Op op;
  FenwickTreeRangePoint1D(int N, T qdef = T(), Op op = Op())
      : N(N), BIT(N + 1, qdef), op(op) {}
  template <class Inv = minus<T>>
  FenwickTreeRangePoint1D(vector<T> A, T qdef = T(),
                          Op op = Op(), Inv inv = Inv())
      : N(A.size()), BIT(move(A)), op(op) {
    adjacent_difference(BIT.begin(), BIT.end(), BIT.begin(), inv);
    BIT.reserve(N + 1); BIT.insert(BIT.begin(), qdef);
    for (int i = 1; i <= N; i++) {
      int j = i + (i & -i); if (j <= N) BIT[j] = op(BIT[j], BIT[i]);
    }
  }
  template <class Inv = minus<T>> vector<T> values(Inv inv = Inv()) {
    vector<T> ret(BIT.begin() + 1, BIT.end()); for (int i = N; i >= 1; i--) {
      int j = i + (i & -i);
      if (j <= N) ret[j - 1] = inv(ret[j - 1], ret[i - 1]);
    }
    for (int i = 1; i < N; i++) ret[i] = op(ret[i - 1], ret[i]);
    return ret;
  }
  void update(int l, T v) {
    for (l++; l <= N; l += l & -l) BIT[l] = op(BIT[l], v);
  }
  template <class Inv = minus<T>>
  void update(int l, int r, T v, Inv inv = Inv()) {
    update(l, v); update(r + 1, inv(BIT[0], v));
  }
  T get(int i) {
    T ret = BIT[0]; for (i++; i > 0; i -= i & -i) ret = op(ret, BIT[i]);
    return ret;
  }
  template <class F> int bsearch(T v, F cmp) {
    T val = BIT[0]; int ind = 0; for (int j = __lg(N + 1); j >= 0; j--) {
      int i = ind + (1 << j);
      if (i <= N && cmp(op(val, BIT[i]), v)) val = op(val, BIT[ind = i]);
    }
    return ind;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "../../PolicyBasedDataStructures.h"
using namespace std;
using namespace __gnu_pbds;

// A collection of sparse fenwick trees implemented in various methods
// If bsearch is not required and online is enforced, then SqrtBufferSimple and
//   SqrtBuffer should be used instead
// If updates only involve values 1 and -1, and online is enforced,
//   then RootArray should be used instead
// In general, offline fenwick tree is faster than SqrtBufferSimple,
//   SqrtBuffer, and RootArray, which are faster than the pbds implementation

// Sparse Fenwick Tree supporting point updates (with any value)
//   and range queries over a cumulative function or functor,
//   such as sum, max, and min, in 1 dimension
// All update indices must be known beforehand
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   IndexType: the type of the index of the array
//   Op: a struct with the cumulative operation (plus<T> by default)
//     Required Functions:
//       operator (l, r): combines the values l and r, must be associative and
//         commutatitve
// Constructor Arguments:
//   N: the size of the array
//   updateInds: a vector of IndexType containing the indices for each update
//   qdef: the identity element of the operation
//   op: an instance of the Op struct
// Functions:
//   update(i, v): modifies the value A[i] to op(A[i], v), i must be an element
//     in updateInds
//   query(r): queries the cumulative value of the range [0, r]
//   query(l, r, inv): queries the cumulative value of the range [l, r], where
//     inv is the inverse of op (minus<T>() by default)
//   bsearch(v, cmp): returns the first index where cmp(op(A[0..i]), v)
//     returns false, or N if no such index exists, assumes A is sorted by cmp
// In practice, has a small constant
// Time Complexity:
//   constructor: O(Q log U) for U updates
//   update, query, bsearch: O(log U) for U updates
// Memory Complexity: O(U) for U updates
// Tested:
//   https://codeforces.com/contest/903/problem/D
//   https://dmoj.ca/problem/cco10p3
//   https://atcoder.jp/contests/dp/tasks/dp_q
template <class T, class IndexType, class Op = plus<T>>
struct OfflineSparseFenwickTree1D {
  static_assert(is_integral<IndexType>::value, "IndexType must be integeral");
  IndexType N; vector<T> BIT; vector<IndexType> inds; T qdef; Op op;
  int getInd(IndexType i) {
    return std::upper_bound(inds.begin(), inds.end(), i) - inds.begin();
  }
  OfflineSparseFenwickTree1D(IndexType N, const vector<IndexType> &updateInds,
                             T qdef = T(), Op op = Op())
      : N(N), inds(updateInds), qdef(qdef), op(op) {
    sort(inds.begin(), inds.end());
    inds.erase(unique(inds.begin(), inds.end()), inds.end());
    BIT.assign(int(inds.size()) + 1, qdef);
  }
  void update(IndexType i, T v) {
    for (int x = getInd(i); x <= int(inds.size()); x += x & -x)
      BIT[x] = op(BIT[x], v);
  }
  T query(IndexType r) {
    T ret = qdef;
    for (int x = getInd(r); x > 0; x -= x & -x) ret = op(ret, BIT[x]);
    return ret;
  }
  template <class Inv = minus<T>>
  T query(IndexType l, IndexType r, Inv inv = Inv()) {
    return inv(query(r), query(l - 1));
  }
  template <class F> IndexType bsearch(T v, F cmp) {
    T agg = qdef; int ind = 0;
    for (int j = __lg(int(inds.size()) + 1); j >= 0; j--) {
      int i = ind + (1 << j);
      if (i < int(inds.size()) && cmp(op(agg, BIT[i]), v))
        agg = op(agg, BIT[ind = i]);
    }
    return ind == int(inds.size()) ? N : inds[ind];
  }
};

// Sparse Fenwick Tree supporting point updates (with any value)
//   and range queries over a cumulative function or functor,
//   such as sum, max, and min, in 1 dimension using pbds hash_table
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   IndexType: the type of the index of the array
//   Op: a struct with the cumulative operation (plus<T> by default)
//     Required Functions:
//       operator (l, r): combines the values l and r, must be associative and
//         commutatitve
// Constructor Arguments:
//   N: the size of the array
//   qdef: the identity element of the operation
//   op: an instance of the Op struct
// Functions:
//   update(i, v): modifies the value A[i] to op(A[i], v)
//   query(r): queries the cumulative value of the range [0, r]
//   query(l, r, inv): queries the cumulative value of the range [l, r], where
//     inv is the inverse of op (minus<T>() by default)
//   bsearch(v, cmp): returns the first index where cmp(op(A[0..i]), v)
//     returns false, or N if no such index exists
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(1)
//   update, query, bsearch: O(log N) on average
// Memory Complexity: O(U log N) for U updates
// Tested:
//   https://codeforces.com/contest/903/problem/D
//   https://dmoj.ca/problem/cco10p3
//   https://atcoder.jp/contests/dp/tasks/dp_q
template <class T, class IndexType, class Container = hashmap<IndexType, T>,
          class Op = plus<T>>
struct SparseFenwickTree1D {
  static_assert(is_integral<IndexType>::value, "IndexType must be integeral");
  IndexType N; Container BIT; T qdef; Op op;
  SparseFenwickTree1D(IndexType N, T qdef = T(), Op op = Op())
      : N(N), qdef(qdef), op(op) {}
  void update(IndexType i, T v) {
    for (i++; i <= N; i += i & -i) BIT[i] = op(BIT[i], v);
  }
  T query(IndexType r) {
    T ret = qdef; for (r++; r > 0; r -= r & -r) {
      auto it = BIT.find(r); if (it != BIT.end()) ret = op(ret, it->second);
    }
    return ret;
  }
  template <class Inv = minus<T>>
  T query(IndexType l, IndexType r, Inv inv = Inv()) {
    return inv(query(r), query(l - 1));
  }
  template <class F> IndexType bsearch(T v, F cmp) {
    T agg = qdef; IndexType ind = 0;
    for (IndexType j = __lg(N + 1); j >= 0; j--) {
      IndexType i = ind + (1 << j); if (i <= N) {
        auto it = BIT.find(i); T val = it == BIT.end() ? 0 : it->second;
        if (cmp(op(agg, val), v)) { agg = op(agg, val); ind = i; }
      }
    }
    return ind;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Fenwick Tree or Binary Indexed Tree supporting point updates
//   and range queries over a cumulative function or functor,
//   such as sum, max, and min, in 1 dimension
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   Op: a struct with the cumulative operation (plus<T> by default)
//     Required Functions:
//       operator (l, r): combines the values l and r, must be associative and
//         commutatitve
// Constructor Arguments:
//   N: the size of the array
//   A: a vector of type T, memory is saved if this is moved and has
//     a capacity of N + 1
//   qdef: the identity element of the operation
//   op: an instance of the Op struct
// Functions:
//   values(inv): returns a vector of the fenwick tree decomposed into an
//     array, where inv is the inverse of op (minus<T>() by default)
//   update(i, v): modifies the value A[i] to op(A[i], v)
//   query(r): queries the cumulative value of the range [0, r]
//   query(l, r, inv): queries the cumulative value of the range [l, r], where
//     inv is the inverse of op (minus<T>() by default)
//   bsearch(v, cmp): returns the first index where cmp(op(A[0..i]), v)
//     returns false, or N if no such index exists
// In practice, has a small constant, faster than segment trees
// Time Complexity:
//   constructor, values: O(N)
//   update, query, bsearch: O(log N)
// Memory Complexity: O(N)
// Tested:
//   Fuzz and Stress Tested
//   https://dmoj.ca/problem/ds1
//   https://judge.yosupo.jp/problem/point_add_range_sum
//   https://atcoder.jp/contests/dp/tasks/dp_q
template <class T, class Op = plus<T>> struct FenwickTree1D {
  int N; vector<T> BIT; Op op;
  FenwickTree1D(int N, T qdef = T(), Op op = Op())
      : N(N), BIT(N + 1, qdef), op(op) {}
  FenwickTree1D(vector<T> A, T qdef = T(), Op op = Op())
      : N(A.size()), BIT(move(A)), op(op) {
    BIT.reserve(N + 1); BIT.insert(BIT.begin(), qdef);
    for (int i = 1; i <= N; i++) {
      int j = i + (i & -i); if (j <= N) BIT[j] = op(BIT[j], BIT[i]);
    }
  }
  template <class Inv = minus<T>> vector<T> values(Inv inv = Inv()) {
    vector<T> ret(BIT.begin() + 1, BIT.end()); for (int i = N; i >= 1; i--) {
      int j = i + (i & -i);
      if (j <= N) ret[j - 1] = inv(ret[j - 1], ret[i - 1]);
    }
    return ret;
  }
  void update(int i, T v) {
    for (i++; i <= N; i += i & -i) BIT[i] = op(BIT[i], v);
  }
  T query(int r) {
    T ret = BIT[0]; for (r++; r > 0; r -= r & -r) ret = op(ret, BIT[r]);
    return ret;
  }
  template <class Inv = minus<T>> T query(int l, int r, Inv inv = Inv()) {
    return inv(query(r), query(l - 1));
  }
  template <class F> int bsearch(T v, F cmp) {
    T agg = BIT[0]; int ind = 0; for (int j = __lg(N + 1); j >= 0; j--) {
      int i = ind + (1 << j);
      if (i <= N && cmp(op(agg, BIT[i]), v)) agg = op(agg, BIT[ind = i]);
    }
    return ind;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "FenwickTree1D.h"
using namespace std;

// Fenwick Tree or Binary Indexed Tree supporting range updates
//   and range sum queries in 1 dimension
// Indices are 0-indexed and ranges are inclusive
// numeric_limits<T>::max() * N must not overflow
// Template Arguments:
//   T: the type of each element
// Constructor Arguments:
//   N: the size of the array
//   A: a vector of type T, memory is saved if this is moved and has
//     a capacity of N + 1
// Functions:
//   update(l, r, v): add v to the range [l..r]
//   query(r): queries the sum of the range [0, r]
//   query(l, r): queries the sum of the range [l, r]
// In practice, has a small constant, faster than segment trees
// Time Complexity:
//   constructor: O(N)
//   update, query: O(log N)
// Memory Complexity: O(N)
// Tested:
//   Fuzz and Stress Tested
//   http://www.usaco.org/index.php?page=viewproblem2&cpid=973
template <class T> struct FenwickTreeRange1D {
  vector<FenwickTree1D<T>> FT;
  FenwickTreeRange1D(int N) : FT(2, FenwickTree1D<T>(N)) {}
  FenwickTreeRange1D(vector<T> A) {
    FT.reserve(2); FT.emplace_back(move(A)); FT.emplace_back(FT[0].N);
  }
  T query(int r) { return FT[1].query(r) * T(r) + FT[0].query(r); }
  T query(int l, int r) { return query(r) - query(l - 1); }
  void update(int l, int r, T v) {
    FT[1].update(l, v); FT[1].update(r + 1, -v);
    FT[0].update(l, v * T(1 - l)); FT[0].update(r + 1, v * T(r));
  }
};


text:
#pragma once
#include <bits/stdc++.h>
using namespace std;

// Fenwick Tree supporting range sum queries where all values are 0 or 1,
//   and values can be updated after they are set
// Indices are 0-indexed and ranges are inclusive
// Constructor Arguments:
//   N: the size of the array
// Functions:
//   set(i, v): sets the initial value of index i with the boolean value v
//   get(i): gets the current value of index i
//   build(): builds the fenwick tree, set should not be called afterwards,
//     does not need to be called if there are no calls to set
//   update(i, v): updates the value of index i with the boolean value v,
//     should not be called before build unless there were no calls to set
//   query(r): queries the sum of the range [0, r]
//   query(l, r): queries the sum of the range [l, r]
//   bsearch(v, cmp): returns the first index where cmp(sum(A[0..i]), v)
//     returns false, or N if no such index exists
// In practice, has a moderate constant, performs faster than a regular fenwick
//   tree and uses less memory
// Time Complexity:
//   constructor, set, get: O(1)
//   update, query, bsearch, lower_bound, upper_bound: O(log(N / 64))
//   build: O(N / 64)
// Memory Complexity: O(N / 64)
// Tested:
//   Fuzz and Stress Tested
//   https://judge.yosupo.jp/problem/predecessor_problem
struct BitFenwickTree {
  int N, M; vector<uint64_t> mask; vector<int> BIT;
  BitFenwickTree(int N) : N(N), M((N >> 6) + 1), mask(M, 0), BIT(M + 1, 0) {}
  void set(int i, bool v) {
    int j = i >> 6, k = i & 63;
    mask[j] = (mask[j] & ~(uint64_t(1) << k)) | (uint64_t(v) << k);
  }
  bool get(int i) { return (mask[i >> 6] >> (i & 63)) & 1; }
  void build() {
    for (int i = 1; i <= M; i++) {
      BIT[i] += __builtin_popcountll(mask[i - 1]); int j = i + (i & -i);
      if (j <= M) BIT[j] += BIT[i];
    }
  }
  void update(int i, bool v) {
    int j = i >> 6, k = i & 63; if (((mask[j] >> k) & 1) != v) {
      mask[j] ^= uint64_t(1) << k; int delta = v ? 1 : -1;
      for (j++; j <= M; j += j & -j) BIT[j] += delta;
    }
  }
  int query(int r) {
    r++; int j = r >> 6, k = r & 63, ret = 0;
    for (int i = j; i > 0; i -= i & -i) ret += BIT[i];
    return ret + __builtin_popcountll(mask[j] & ((uint64_t(1) << k) - 1));
  }
  int query(int l, int r) { return query(r) - query(l - 1); }
  template <class F> int bsearch(int v, F cmp) {
    int sum = 0, ind = 0; for (int j = __lg(M + 1); j >= 0; j--) {
      int i = ind + (1 << j);
      if (i <= M && cmp(sum + BIT[i], v)) sum += BIT[ind = i];
    }
    if (ind == M) return N;
    uint64_t m = mask[ind]; ind <<= 6;
    for (;; ind++, m >>= 1) if (!cmp(sum += m & 1, v)) break;
    return ind;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#if __cplusplus < 201402L
#include "../../../utils/MakeUnique.h"
#endif
using namespace std;

// Skew Heap supporting merges
// Template Arguments:
//   T: the type of each element
//   Cmp: the comparator to compare two values of type T,
//       convention is same as std::priority_queue in STL
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Functions:
//   empty(): returns whether the heap is empty
//   size(): returns the number of elements in the heap
//   top(): returns the largest element in the heap based on the comparator
//   pop(): pops the largest element in the heap based on the comparator,
//     and returns that element
//   push(key): pushes val onto the heap
//   merge(h): merges the heap h into this heap, destroys heap h
// In practice, has a moderate constant
// Time Complexity:
//   constructor, empty, size, top: O(1)
//   pop, push: O(log N) amortized
//   merge: O(log (N + M)) amortized where M is the size of heap h
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/apio16p2
template <class T, class Cmp = less<T>> struct SkewHeap {
  struct Node { T val; unique_ptr<Node> l, r; Node(const T &v) : val(v) {} };
  Cmp cmp; int cnt; unique_ptr<Node> root;
  unique_ptr<Node> merge(unique_ptr<Node> a, unique_ptr<Node> b) {
    if (!a || !b) return a ? move(a) : move(b);
    if (cmp(a->val, b->val)) a.swap(b);
    a->l.swap(a->r); a->r = merge(move(b), move(a->r)); return move(a);
  }
  SkewHeap() : cnt(0) {}
  bool empty() const { return !root; }
  int size() const { return cnt; }
  T top() const { return root->val; }
  T pop() {
    T ret = root->val; root = merge(move(root->l), move(root->r)); cnt--;
    return ret;
  }
  void push(const T &val) {
    root = merge(move(root), make_unique<Node>(val)); cnt++;
  }
  void merge(SkewHeap &h) {
    root = merge(move(root), move(h.root)); cnt += h.cnt;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "../../../utils/Random.h"
using namespace std;

// Persistent Heap supporting merges where copy assignment/constructor creates
//   a new version of the data structure
// Template Arguments:
//   T: the type of each element
//   Cmp: the comparator to compare two values of type T,
//       convention is same as std::priority_queue in STL
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Functions:
//   empty(): returns whether the current heap is empty
//   size(): returns the number of elements in the current heap
//   top(): returns the largest element in the current heap based
//     on the comparator
//   pop(): pops the largest element in the current heap based
//     on the comparator, and returns that element
//   push(key): pushes val onto the current heap
//   merge(h): merges the heap h into the current heap
// In practice, has a moderate constant
// Time Complexity:
//   constructor, empty, size, top: O(1)
//   pop, push: O(log N) expected
//   merge: O(log (N + M)) expected where M is the size of heap h
// Memory Complexity: O(N + Q log N) for Q operations
// Tested:
//   https://dmoj.ca/problem/wac4p5
template <class T, class Cmp = less<T>> struct PersistentRandomizedHeap {
  struct Node; using ptr = shared_ptr<Node>;
  struct Node {
    T val; ptr l, r;
    Node(const T &v, const ptr &l = ptr(), const ptr &r = ptr())
        : val(v), l(l), r(r) {}
  };
  Cmp cmp; int cnt; ptr root;
  ptr merge(ptr a, ptr b) {
    if (!a || !b) return a ? a : b;
    if (cmp(a->val, b->val)) a.swap(b);
    return rng() % 2 ? make_shared<Node>(a->val, a->l, merge(a->r, b))
                     : make_shared<Node>(a->val, merge(a->l, b), a->r);
  }
  PersistentRandomizedHeap() : cnt(0) {}
  bool empty() const { return !root; }
  int size() const { return cnt; }
  T top() { return root->val; }
  T pop() {
    T ret = root->val; root = merge(root->l, root->r); cnt--; return ret;
  }
  void push(const T &val) {
    root = merge(root, make_shared<Node>(val)); cnt++;
  }
  void merge(const PersistentRandomizedHeap &h) {
    root = merge(root, h.root); cnt += h.cnt;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include "../../../utils/Random.h"
using namespace std;

// Persistent Heap supporting merges and increments where
//   copy assignment/constructor creates a new version of the data structure
// Template Arguments:
//   T: the type of each element
//   Cmp: the comparator to compare two values of type T,
//       convention is same as std::priority_queue in STL
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
//   Delta: the type of the incremental element
// Functions:
//   empty(): returns whether the current heap is empty
//   size(): returns the number of elements in the current heap
//   top(): returns the largest element in the current heap based
//     on the comparator
//   pop(): pops the largest element in the current heap based
//     on the comparator, and returns that element
//   push(key): pushes val onto the current heap
//   increment(delta): increments all elements in the current heap by delta
//   merge(h): merges the heap h into the current heap
// In practice, has a moderate constant
// Time Complexity:
//   constructor, empty, size, top, increment: O(1)
//   pop, push: O(log N) expected
//   merge: O(log (N + M)) expected where M is the size of heap h
// Memory Complexity: O(N + Q log N) for Q operations
// Tested:
//   https://dmoj.ca/problem/wac4p5
template <class T, class Cmp = less<T>, class Delta = T>
struct PersistentRandomizedHeapIncremental {
  struct Node; using ptr = shared_ptr<Node>;
  struct Node {
    T val; Delta delta; ptr l, r;
    Node(const T &v, const Delta &d,
         const ptr &l = ptr(), const ptr &r = ptr())
        : val(v), delta(d), l(l), r(r) {}
  };
  Cmp cmp; Delta ddef; int cnt; ptr root;
  void propagate(ptr &a) {
    a->val += a->delta;
    if (a->l) { a->l = make_shared<Node>(*a->l); a->l->delta += a->delta; }
    if (a->r) { a->r = make_shared<Node>(*a->r); a->r->delta += a->delta; }
    a->delta = ddef;
  }
  ptr merge(ptr a, ptr b) {
    if (!a || !b) return a ? a : b;
    propagate(a); propagate(b); if (cmp(a->val, b->val)) a.swap(b);
    return rng() % 2
           ? make_shared<Node>(a->val, a->delta, a->l, merge(a->r, b))
           : make_shared<Node>(a->val, a->delta, merge(a->l, b), a->r);
  }
  PersistentRandomizedHeapIncremental(const Delta &ddef = Delta())
      : ddef(ddef), cnt(0) {}
  bool empty() const { return !root; }
  int size() const { return cnt; }
  T top() { propagate(root); return root->val; }
  T pop() {
    propagate(root); T ret = root->val;
    root = merge(root->l, root->r); cnt--; return ret;
  }
  void push(const T &val) {
    root = merge(root, make_shared<Node>(val, ddef)); cnt++;
  }
  void increment(const Delta &delta) {
    if (root) { root = make_shared<Node>(*root); root->delta += delta; }
  }
  void merge(const PersistentRandomizedHeapIncremental &h) {
    root = merge(root, h.root); cnt += h.cnt;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#if __cplusplus < 201402L
#include "../../../utils/MakeUnique.h"
#endif
using namespace std;

// Skew Heap supporting merges and increments
// Template Arguments:
//   T: the type of each element
//   Cmp: the comparator to compare two values of type T,
//       convention is same as std::priority_queue in STL
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
//   Delta: the type of the incremental element
// Functions:
//   empty(): returns whether the heap is empty
//   size(): returns the number of elements in the heap
//   top(): returns the largest element in the heap based on the comparator
//   pop(): pops the largest element in the heap based on the comparator,
//     and returns that element
//   push(key): pushes val onto the heap
//   increment(delta): increments all elements in the heap by delta
//   merge(h): merges the heap h into this heap, destroys heap h
// In practice, has a moderate constant
// Time Complexity:
//   constructor, empty, size, top, increment: O(1)
//   pop, push: O(log N) amortized
//   merge: O(log (N + M)) amortized
// Memory Complexity: O(N)
// Tested:
//   https://codeforces.com/contest/1041/problem/D
template <class T, class Cmp = less<T>, class Delta = T>
struct SkewHeapIncremental {
  Cmp cmp; Delta ddef;
  struct Node {
    T val; Delta delta; unique_ptr<Node> l, r;
    Node(const T &v, const Delta &d) : val(v), delta(d) {}
  };
  int cnt; unique_ptr<Node> root;
  void propagate(unique_ptr<Node> &a) {
    a->val += a->delta;
    if (a->l) a->l->delta += a->delta;
    if (a->r) a->r->delta += a->delta;
    a->delta = ddef;
  }
  unique_ptr<Node> merge(unique_ptr<Node> a, unique_ptr<Node> b) {
    if (!a || !b) return a ? move(a) : move(b);
    propagate(a); propagate(b); if (cmp(a->val, b->val)) a.swap(b);
    a->l.swap(a->r); a->r = merge(move(b), move(a->r)); return move(a);
  }
  SkewHeapIncremental(const Delta &ddef = Delta()) : ddef(ddef), cnt(0) {}
  bool empty() const { return !root; }
  int size() const { return cnt; }
  T top() { propagate(root); return root->val; }
  T pop() {
    propagate(root); T ret = root->val;
    root = merge(move(root->l), move(root->r)); cnt--; return ret;
  }
  void push(const T &val) {
    root = merge(move(root), make_unique<Node>(val, ddef)); cnt++;
  }
  void increment(const Delta &delta) { if (root) root->delta += delta; }
  void merge(SkewHeapIncremental &h) {
    root = merge(move(root), move(h.root)); cnt += h.cnt;
  }
};


text:
cmake_minimum_required(VERSION 3.4)
project(FeatureDetector CXX)

add_subdirectory(src/x86)
target_include_directories(FeatureDetector PUBLIC src)
if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME)
    add_executable(FeatureDetector_main src/Main.cpp)
    target_link_libraries(FeatureDetector_main PUBLIC FeatureDetector)
endif()

repo_name: Mysticial/FeatureDetector

mime_type: text/xml


text: Program('bin/FeatureDetector', ['src/Main.cpp', 'src/x86/cpu_x86.cpp'])


text:
Microsoft Visual Studio Solution File, Format Version 12.00
# Visual Studio 2013
VisualStudioVersion = 12.0.31101.0
MinimumVisualStudioVersion = 10.0.40219.1
Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "FeatureDetector", "FeatureDetector\FeatureDetector.vcxproj", "{5AF4062D-386F-48AA-B824-2C10A4B1BA97}"
EndProject
Global
	GlobalSection(SolutionConfigurationPlatforms) = preSolution
		Debug|Win32 = Debug|Win32
		Release|Win32 = Release|Win32
	EndGlobalSection
	GlobalSection(ProjectConfigurationPlatforms) = postSolution
		{5AF4062D-386F-48AA-B824-2C10A4B1BA97}.Debug|Win32.ActiveCfg = Debug|Win32
		{5AF4062D-386F-48AA-B824-2C10A4B1BA97}.Debug|Win32.Build.0 = Debug|Win32
		{5AF4062D-386F-48AA-B824-2C10A4B1BA97}.Release|Win32.ActiveCfg = Release|Win32
		{5AF4062D-386F-48AA-B824-2C10A4B1BA97}.Release|Win32.Build.0 = Release|Win32
	EndGlobalSection
	GlobalSection(SolutionProperties) = preSolution
		HideSolutionNode = FALSE
	EndGlobalSection
EndGlobal


text:
# FeatureDetector

A simple app that detects what hardware features are supported by the host CPU and operating system.

-----

Sample output for a Haswell processor:

````
CPU Vendor String: GenuineIntel

CPU Vendor:
    AMD         = No
    Intel       = Yes

OS Features:
    64-bit      = Yes
    OS AVX      = Yes
    OS AVX512   = No

Hardware Features:
    MMX         = Yes
    x64         = Yes
    ABM         = Yes
    RDRAND      = Yes
    BMI1        = Yes
    BMI2        = Yes
    ADX         = No
    MPX         = No
    PREFETCHWT1 = No

SIMD: 128-bit
    SSE         = Yes
    SSE2        = Yes
    SSE3        = Yes
    SSSE3       = Yes
    SSE4a       = No
    SSE4.1      = Yes
    SSE4.2      = Yes
    AES-NI      = Yes
    SHA         = No

SIMD: 256-bit
    AVX         = Yes
    XOP         = No
    FMA3        = Yes
    FMA4        = No
    AVX2        = Yes

SIMD: 512-bit
    AVX512-F    = No
    AVX512-CD   = No
    AVX512-PF   = No
    AVX512-ER   = No
    AVX512-VL   = No
    AVX512-BW   = No
    AVX512-DQ   = No
    AVX512-IFMA = No
    AVX512-VBMI = No

Summary:
    Safe to use AVX:     Yes
    Safe to use AVX512:  No


text:
/* Main.cpp
 *
 * Author           : Alexander J. Yee
 * Date Created     : 04/17/2015
 * Last Modified    : 04/17/2015
 *
 */

#include <iostream>
using std::cout;
using std::endl;

#include "x86/cpu_x86.h"
using namespace FeatureDetector;

int main(){

    cout << "CPU Vendor String: " << cpu_x86::get_vendor_string() << endl;
    cout << endl;

    cpu_x86::print_host();

#if _WIN32
    system("pause");
#endif
}


text:
/* cpu_x86_Linux.ipp
 *
 * Author           : Alexander J. Yee
 * Date Created     : 04/12/2014
 * Last Modified    : 04/12/2014
 *
 */

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//  Dependencies
#include <cpuid.h>
#include "cpu_x86.h"
namespace FeatureDetector{
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void cpu_x86::cpuid(int32_t out[4], int32_t eax, int32_t ecx){
    __cpuid_count(eax, ecx, out[0], out[1], out[2], out[3]);
}
uint64_t xgetbv(unsigned int index){
    uint32_t eax, edx;
    __asm__ __volatile__("xgetbv" : "=a"(eax), "=d"(edx) : "c"(index));
    return ((uint64_t)edx << 32) | eax;
}
#define _XCR_XFEATURE_ENABLED_MASK  0
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//  Detect 64-bit
bool cpu_x86::detect_OS_x64(){
    //  We only support x64 on Linux.
    return true;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
}


text:
/* cpu_x86.h
 *
 * Author           : Alexander J. Yee
 * Date Created     : 04/12/2014
 * Last Modified    : 04/12/2014
 *
 */

#pragma once
#ifndef _cpu_x86_H
#define _cpu_x86_H
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//  Dependencies
#include <stdint.h>
#include <string>
namespace FeatureDetector{
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
struct cpu_x86{
    //  Vendor
    bool Vendor_AMD;
    bool Vendor_Intel;

    //  OS Features
    bool OS_x64;
    bool OS_AVX;
    bool OS_AVX512;

    //  Misc.
    bool HW_MMX;
    bool HW_x64;
    bool HW_ABM;
    bool HW_RDRAND;
    bool HW_RDSEED;
    bool HW_BMI1;
    bool HW_BMI2;
    bool HW_ADX;
    bool HW_MPX;
    bool HW_PREFETCHW;
    bool HW_PREFETCHWT1;
    bool HW_RDPID;

    //  SIMD: 128-bit
    bool HW_SSE;
    bool HW_SSE2;
    bool HW_SSE3;
    bool HW_SSSE3;
    bool HW_SSE41;
    bool HW_SSE42;
    bool HW_SSE4a;
    bool HW_AES;
    bool HW_SHA;

    //  SIMD: 256-bit
    bool HW_AVX;
    bool HW_XOP;
    bool HW_FMA3;
    bool HW_FMA4;
    bool HW_AVX2;

    //  SIMD: 512-bit
    bool HW_AVX512_F;
    bool HW_AVX512_CD;

    //  Knights Landing
    bool HW_AVX512_PF;
    bool HW_AVX512_ER;

    //  Skylake Purley
    bool HW_AVX512_VL;
    bool HW_AVX512_BW;
    bool HW_AVX512_DQ;

    //  Cannon Lake
    bool HW_AVX512_IFMA;
    bool HW_AVX512_VBMI;

    //  Knights Mill
    bool HW_AVX512_VPOPCNTDQ;
    bool HW_AVX512_4VNNIW;
    bool HW_AVX512_4FMAPS;

    //  Cascade Lake
    bool HW_AVX512_VNNI;

    //  Cooper Lake
    bool HW_AVX512_BF16;

    //  Ice Lake
    bool HW_AVX512_VBMI2;
    bool HW_GFNI;
    bool HW_VAES;
    bool HW_AVX512_VPCLMUL;
    bool HW_AVX512_BITALG;

public:
    cpu_x86();
    void detect_host();

    void print() const;
    static void print_host();

    static void cpuid(int32_t out[4], int32_t eax, int32_t ecx);
    static std::string get_vendor_string();

private:
    static void print(const char* label, bool yes);

    static bool detect_OS_x64();
    static bool detect_OS_AVX();
    static bool detect_OS_AVX512();
};
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
}
#endif