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

// Palindromic Tree Node backed by a map
// Template Arguments:
//   _T: the type of the element in the string/array for the Palindromic Tree
// Constructor Arguments:
//   len: the length of the longest palindromic substring represented
//     by this node
// Fields:
//   T: the data type of each element in the string/array
//   len: the length of the longest palindromic substring this
//     node represents
//   link: the index of the suffix link of this node (the node with the
//     longest palindromic suffix of this node)
//   qlink: the index of the quick link of this node (the node with the
//     longest palindromic suffix which has a different preceding character
//     as the link node)
// Functions:
//   getEdge(a): returns the index of the other node on the
//     incident edge with the element a or 1 if no such edge exists
//   setEdge(a, n): sets the other node on the incident edge
//     with the element a to n
// Time Complexity:
//   constructor: O(1)
//   getEdge, setEdge: O(log E) where E is the number of edges
//     incident to this node
// Memory Complexity: O(E) where E is the number of edges incident to this node
// Tested:
//   https://dmoj.ca/problem/apio14p1
template <class _T> struct PalTreeMapNode {
  using T = _T; int len, link, qlink; map<T, int> to;
  PalTreeMapNode(int len) : len(len), link(1), qlink(1) {}
  int getEdge(const T &a) const {
    auto it = to.find(a); return it == to.end() ? 1 : it->second;
  }
  void setEdge(const T &a, int n) { to[a] = n; }
};

// Palindromic Tree Node  backed by an array
// Template Arguments:
//   _T: the type of the element in the string/array for the Palindromic Tree
//   ALPHABET_SIZE: the size of the alphabet
//   OFFSET: the offset for the start of the alphabet
// Constructor Arguments:
//   len: the length of the longest palindromic substring represented
//     by this node
// Fields:
//   T: the data type of each element in the string/array
//   len: the length of the longest palindromic substring this
//     node represents
//   link: the index of the suffix link of this node (the node with the
//     longest palindromic suffix of this node)
//   qlink: the index of the quick link of this node (the node with the
//     longest palindromic suffix which has a different preceding character
//     as the link node)
// Functions:
//   getEdge(a): returns the index of the other node on the
//     incident edge with the element a or 1 if no such edge exists
//   setEdge(a, n): sets the other node on the incident edge
//     with the element a to n
// Time Complexity:
//   constructor: O(ALPHABET_SIZE)
//   getEdge, setEdge: O(1)
// Memory Complexity: O(ALPHABET_SIZE)
// Tested:
//   https://dmoj.ca/problem/mmcc15p3
//   https://dmoj.ca/problem/apio14p1
template <class _T, const int ALPHABET_SIZE, const _T OFFSET>
struct PalTreeArrayNode {
  using T = _T; int len, link, qlink; array<int, ALPHABET_SIZE> to;
  PalTreeArrayNode(int len) : len(len), link(1), qlink(1) { to.fill(1); }
  int getEdge(const T &a) const { return to[a - OFFSET]; }
  void setEdge(const T &a, int n) { to[a - OFFSET] = n; }
};

// Palindromic Tree with two roots at 0 (with length -1) and 1 (with length 0)
// Template Arguments:
//   Node: a node class
//     Required Fields:
//       T: the data type of each element in the string/array
//       len: the length of the longest palindromic substring this
//         node represents
//       link: the index of the suffix link of this node (the node with the
//         longest palindromic suffix of this node)
//       qlink: the index of the quick link of this node (the node with the
//         longest palindromic suffix which has a different preceding character
//         as the link node)
//     Required Functions:
//       constructor(len): initializes the node with a length of len with
//         link, qlink, and edges pointing to 1 by default
//       getEdge(a): returns the index of the other node on the
//         incident edge with the element a
//       setEdge(a, n): sets the other node on the incident edge
//         with the element a to n
// Constructor Arguments:
//   def: the default value of type T, cannot be in the string/array
// Fields:
//   S: the current string with the default character at the front and
//     between words
//   TR: the vector of all nodes in the tree
//   last: a vector of integers with the last node after each addition or
//     terminate call
// Functions:
//   add(a): adds the element a to the tree
//   terminate(): terminates the current word, returning to the root (which
//     allows for another word to be added with add(a))
//   undo(): undoes the last added element, or a terminate operation
// In practice, has a small constant
// Time Complexity:
//   constructor: time complexity of node constructor
//   add: O(log N) + time complexity of node constructor
//        + time complexity of getEdge in Node
//   terminate: O(1)
//   undo: time complexity of setEdge in Node
// Memory Complexity: O(N) * memory complexity of node, after N calls to add
// Tested:
//   https://dmoj.ca/problem/mmcc15p3
//   https://dmoj.ca/problem/apio14p1
template <class Node> struct PalindromicTree {
  using T = typename Node::T;
  T def; vector<T> S; vector<Node> TR; vector<int> last, modified;
  PalindromicTree(const T &def)
      : def(def), S(1, def), TR(vector<Node>{Node(-1), Node(0)}), last(1, 1) {
    TR[1].link = TR[1].qlink = 0;
  }
  int getLink(int x, int i) {
    while (S[i - 1 - TR[x].len] != S[i])
      x = S[i - 1 - TR[TR[x].link].len] == S[i] ? TR[x].link : TR[x].qlink;
    return x;
  }
  void add(const T &a) {
    int i = S.size(); S.push_back(a); int p = getLink(last.back(), i);
    modified.push_back(-1); if (TR[p].getEdge(a) == 1) {
      int u = TR.size(); TR.emplace_back(TR[p].len + 2);
      TR[u].link = TR[getLink(TR[p].link, i)].getEdge(a);
      T b = S[i - TR[TR[u].link].len], c = S[i - TR[TR[TR[u].link].link].len];
      TR[u].qlink = b == c ? TR[TR[u].link].qlink : TR[u].link;
      TR[modified.back() = p].setEdge(a, u);
    }
    last.push_back(TR[p].getEdge(a));
  }
  void terminate() {
    S.push_back(def); last.push_back(1); modified.push_back(-1);
  }
  void undo() {
    if (modified.back() != -1) {
      TR[modified.back()].setEdge(S.back(), 1); TR.pop_back();
    }
    S.pop_back(); last.pop_back(); modified.pop_back();
  }
};


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

// Modified from https://github.com/cheran-senthil/PyRival/blob/master/pyrival/strings/suffix_array.py,
//   which has an Apache 2.0 license
// Suffix Array Induced Sort to sort suffixes of an array in
//   lexicographical order
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   _T: the type of each element in the array
// Constructor Arguments:
//   S: a vector of type _T
// Fields:
//   T: the type of the character/element in the array
//   N: the length of the array
//   rnk: a vector of the ranks of the suffixes (rnk[i] is the rank of the
//     suffix starting from index i)
//   ind: a vector of the indices in the original array of the suffixes
//     sorted in lexicographical order (ind[i] is the index in original array
//     of the ith lexicographically smallest suffix)
//   LCP: a vector of the longest common prefixes between the suffixes when
//     sorted in lexicographical order (LCP[i] is the longest common prefix of
//     the ith and (i + 1)th lexicographically smallest suffix, with LCP[N - 1]
//     being 0)
// In practice, has a moderate constant, usually faster than SuffixArray
// Time Complexity:
//   constructor: O(N + K) where K is the range of the array
// Memory Complexity: O(N + K)
// Tested:
//   Fuzz and Stress Tested
//   https://judge.yosupo.jp/problem/suffixarray
//   https://judge.yosupo.jp/problem/number_of_substrings
//   https://dmoj.ca/problem/coci06c5p6
//   https://dmoj.ca/problem/ccc20s3
template <class _T> struct SAISSuffixArray {
  using T = _T;
  static vector<int> SAIS(const vector<int> &S) {
    if (S.empty()) return vector<int>();
    int N = S.size(), K = *max_element(S.begin(), S.end()) + 1;
    vector<bool> isL(N + 1, false); vector<int> ind(N), cnt(K + 1, 0), lms;
    lms.reserve(N); for (auto &&a : S) cnt[a + 1]++;
    partial_sum(cnt.begin(), cnt.end(), cnt.begin());
    isL[N - 1] = true; for (int i = N - 2; i >= 0; i--)
      isL[i] = S[i] == S[i + 1] ? isL[i + 1] : S[i] > S[i + 1];
    auto c = [&] (int i) { return i > 0 && !isL[i] && isL[i - 1]; };
    for (int i = 1; i < N; i++) if (c(i)) lms.push_back(i);
    auto IS = [&] {
      vector<int> tmp(cnt.begin() + 1, cnt.end());
      fill(ind.begin(), ind.end(), -1);
      for (int i = int(lms.size()) - 1; i >= 0; i--)
        ind[--tmp[S[lms[i]]]] = lms[i];
      tmp = vector<int>(cnt.begin(), cnt.end() - 1);
      for (int i = -1, j; i < N; i++)
        if ((j = i < 0 ? N - 1 : ind[i] - 1) >= 0 && isL[j])
          ind[tmp[S[j]]++] = j;
      tmp = vector<int>(cnt.begin() + 1, cnt.end());
      for (int i = N - 1, j; i >= 0; i--)
        if ((j = ind[i] - 1) >= 0 && !isL[j]) ind[--tmp[S[j]]] = j;
    };
    if (int(lms.size()) > 1) {
      IS(); vector<int> tmp(ind.begin(), ind.end());
      for (int i = 0, j = 0, a = -1; i < N; i++) if (c(tmp[i])) {
        for (int b = tmp[i]; a >= 0 && S[a] == S[b];) {
          a++; b++; if (c(a) || c(b)) { j -= int(c(a) && c(b)); break; }
        }
        ind[a = tmp[i]] = ++j;
      }
      tmp = vector<int>(lms.size());
      for (int i = 0; i < int(tmp.size()); i++) tmp[i] = ind[lms[i]];
      tmp = SAIS(tmp);
      for (int i = 0; i < int(tmp.size()); i++) tmp[i] = lms[tmp[i]];
      lms = tmp;
    }
    IS(); return ind;
  }
  static vector<int> init(const vector<T> &S) {
    if (S.empty()) return vector<int>();
    T offset = *min_element(S.begin(), S.end()); vector<int> A(S.size());
    for (int i = 0; i < int(S.size()); i++) A[i] = S[i] - offset;
    return A;
  }
  int N; vector<int> ind, rnk, LCP;
  SAISSuffixArray(const vector<T> &S) : N(S.size()) {
    ind = SAIS(init(S)); rnk.resize(N); LCP.resize(N);
    for (int i = 0; i < N; i++) rnk[ind[i]] = i;
    for (int i = 0, k = 0; i < N; i++) {
      if (rnk[i] == N - 1) { LCP[rnk[i]] = k = 0; continue; }
      int j = ind[rnk[i] + 1];
      while (i + k < N && j + k < N && S[i + k] == S[j + k]) k++;
      if ((LCP[rnk[i]] = k) > 0) k--;
    }
  }
};

repo_name: wesley-a-leung/Resources

mime_type: text/x-c++


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

// Factorizes a array S into Lyndon words w1w2w3... where w1, w2, w3, ... are
//   in non-increasing order
// Guaranteed to exist and is unique
// A Lyndon word is lexicographically smaller than all of its nontrivial
//   suffixes
// Template Arguments:
//   T: the type of each element in the array
// Function Arguments:
//   S: a vector of type T
// Return Value: a vector of integers representing the lengths of each factor,
//   where the sum of the values in the vector is equal to
//   the length of the array
// In practice, has a small constant
// Time Complexity: O(N)
// Memory Complexity: O(N)
// Tested:
//   https://cses.fi/problemset/task/1110/
template <class T> vector<int> lyndonFactorization(const vector<T> &S) {
  int N = S.size(); vector<int> ret; for (int i = 0; i < N;) {
    int j = i + 1, k = i;
    for (; j < N && S[k] <= S[j]; j++) k = S[k] < S[j] ? i : k + 1;
    for (; i <= k; i += j - k) ret.push_back(j - k);
  }
  return ret;
}


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

// Computes the longest common prefix of two suffixes of an array
//   be using the Fischer Heun Structure over the LCP array generated by the
//   Suffix Array
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   SuffixArray: a generic suffix array to be used (should be either
//       SuffixArray or SAISSuffixArray)
//     Required Fields:
//       N: the length of the array
//       rnk: a vector of the ranks of the suffixes (rnk[i] is the rank of the
//         suffix starting from index i)
//       ind: a vector of the indices in the original array of the suffixes
//         sorted in lexicographical order (ind[i] is the index in original
//         array of the ith lexicographically smallest suffix)
//       LCP: a vector of the longest common prefixes between the suffixes when
//         sorted in lexicographical order (LCP[i] is the longest common prefix
//         of the ith and (i + 1)th lexicographically smallest suffix, with
//         LCP[N - 1] being 0)
//     Required Functions:
//       constructor(S): construts a suffix array from the array S
// Constructor Arguments:
//   S: a vector of type SuffixArray::T
// Fields:
//   SA: the associated suffix array constructed from the array
// Functions:
//   lcpRnk(i, j): computes the longest common prefix of the ith and jth
//     lexicographically least suffixes
//   lcp(i, j): computes the longest common prefix of the suffixes starting
//     from index i and j
// In practice, the constructor has a very small constant, lcp has a
//   moderate constant, still faster than using segment trees
// Time Complexity:
//   constructor: time complexity of SuffixArray constructor
//   lcpRnk, lcp: O(1)
// Memory Complexity: memory complexity of SuffixArray
// Tested:
//   Fuzz and Stress Tested
//   https://dmoj.ca/problem/ccc20s3
template <class SuffixArray> struct LongestCommonPrefix {
  SuffixArray SA; FischerHeunStructure<int, greater<int>> FHS;
  LongestCommonPrefix(const vector<typename SuffixArray::T> &S)
    : SA(S), FHS(SA.LCP) {}
  int lcpRnk(int i, int j) {
    if (i > j) swap(i, j);
    return i == j ? int(SA.N) - SA.ind[j] : FHS.query(i, j - 1);
  }
  int lcp(int i, int j) { return lcpRnk(SA.rnk[i], SA.rnk[j]); }
};


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

// Suffix Array using Sadakane's algorithm to sort suffixes of an array in
//   lexicographical order
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   _T: the type of each element in the array
// Constructor Arguments:
//   S: a vector of type _T
// Fields:
//   T: the type of the character/element in the array
//   N: the length of the array
//   rnk: a vector of the ranks of the suffixes (rnk[i] is the rank of the
//     suffix starting from index i)
//   ind: a vector of the indices in the original array of the suffixes
//     sorted in lexicographical order (ind[i] is the index in original array
//     of the ith lexicographically smallest suffix)
//   LCP: a vector of the longest common prefixes between the suffixes when
//     sorted in lexicographical order (LCP[i] is the longest common prefix of
//     the ith and (i + 1)th lexicographically smallest suffix, with LCP[N - 1]
//     being 0)
// In practice, has a very small constant, usually slower than
//   SAISSuffixArray.h
// Time Complexity:
//   constructor: O(N (log N)^2)
// Memory Complexity: O(N)
// Tested:
//   Fuzz and Stress Tested
//   https://judge.yosupo.jp/problem/suffixarray
//   https://judge.yosupo.jp/problem/number_of_substrings
//   https://dmoj.ca/problem/coci06c5p6
//   https://dmoj.ca/problem/ccc20s3
template <class _T> struct SuffixArray {
  using T = _T; int N; vector<int> ind, rnk, LCP;
  SuffixArray(const vector<T> &S)
      : N(S.size()), ind(N + 1), rnk(N + 1), LCP(N + 1) {
    vector<int> &tmp = LCP; iota(ind.begin(), ind.end(), 0);
    sort(ind.begin(), ind.begin() + N, [&] (int a, int b) {
      return S[a] < S[b];
    });
    rnk[ind[N]] = -1; for (int i = 0; i < N; i++) {
      rnk[ind[i]] = i > 0 && S[ind[i]] == S[ind[i - 1]] ? rnk[ind[i - 1]] : i;
    }
    for (int h = 1; h < N; h += h) for (int l = 0, r = 1; r <= N; r++) {
      if (rnk[ind[r - 1]] != rnk[ind[r]] && l + 1 < r) {
        sort(ind.begin() + l, ind.begin() + r, [&] (int a, int b) {
          return rnk[h + a] < rnk[h + b];
        });
        tmp[l] = l; for (int j = l + 1; j < r; j++)
          tmp[j] = rnk[h + ind[j - 1]] < rnk[h + ind[j]] ? j : tmp[j - 1];
        for (l++; l < r; l++) rnk[ind[l]] = tmp[l];
      } else if (rnk[ind[r - 1]] != rnk[ind[r]]) l++;
    }
    ind.pop_back(); rnk.pop_back(); tmp.pop_back();
    for (int i = 0, k = 0; i < N; i++) {
      if (rnk[i] == N - 1) { LCP[rnk[i]] = k = 0; continue; }
      int j = ind[rnk[i] + 1];
      while (i + k < N && j + k < N && S[i + k] == S[j + k]) k++;
      if ((LCP[rnk[i]] = k) > 0) k--;
    }
  }
};


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

// Finds the longest common subarray of two arrays
// Indices are 0-indexed
// Template Arguments:
//   SuffixArray: a generic suffix array to be used (should be either
//       SuffixArray or SAISSuffixArray)
//     Required Fields:
//       N: the length of the array
//       rnk: a vector of the ranks of the suffixes (rnk[i] is the rank of the
//         suffix starting from index i)
//       ind: a vector of the indices in the original array of the suffixes
//         sorted in lexicographical order (ind[i] is the index in original
//         array of the ith lexicographically smallest suffix)
//       LCP: a vector of the longest common prefixes between the suffixes when
//         sorted in lexicographical order (LCP[i] is the longest common prefix
//         of the ith and (i + 1)th lexicographically smallest suffix, with
//         LCP[N - 1] being 0)
//     Required Functions:
//       constructor(S): construts a suffix array from the array S
// Constructor Arguments:
//   A: the first vector of type SuffixArray::T
//   B: the second vector of type SuffixArray::T
//   sep: a separator character that does not appear in either array
// Fields:
//   SA: the associated suffix array constructed from the array A,
//     concatenated with the separator, conatenated with the array B
//   lcs: the longest common subarray of A and B
// In practice, the constructor has a very small constant
// Time Complexity:
//   constructor: time complexity of SuffixArray for an array of
//     length N + M + 1
// Memory Complexity: memory complexity of SuffixArray for an array of
//   length N + M + 1
// Tested:
//   https://www.spoj.com/problems/LCS/
template <class SuffixArray> struct LongestCommonSubstringOfTwo {
  using T = typename SuffixArray::T; SuffixArray SA; vector<T> lcs;
  vector<T> init(vector<T> A, const vector<T> &B, T sep) {
    A.reserve(A.size() + B.size() + 1); A.push_back(sep);
    A.insert(A.end(), B.begin(), B.end()); return A;
  }
  LongestCommonSubstringOfTwo(const vector<T> &A, const vector<T> &B, T sep)
      : SA(init(A, B, sep)) {
    pair<int, int> mx(-1, -1); int N = A.size();
    for (int i = 0; i < N + int(B.size()); i++)
      if (SA.ind[i] != N && SA.ind[i + 1] != N
          && (SA.ind[i] < N) != (SA.ind[i + 1] < N))
        mx = max(mx, make_pair(SA.LCP[i], SA.ind[i]));
    if (mx.first != -1) {
      lcs.reserve(mx.first); for (int i = 0; i < mx.first; i++) {
        int j = mx.second + i; lcs.push_back(j < N ? A[j] : B[j - 1 - N]);
      }
    }
  }
};


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

// Computes the hash of a 2D array to allow for easy computation of
//   submatrix hashes
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   HTYPE: the type of the hash
//   HASHES: the number of hashes to compute
// Construtor Arguments:
//   A: the elements in the 2D array
//   offset: the amount to subtract from each element, a - offset + 1 should be
//     positive
//   ROWBASE: array<HTYPE, HASHES> that contains the base for the rows for
//     each hash, must be less than the repective mod and positive (preferably
//     in the range [mod / 4, mod / 4 * 3])
//   COLBASE: array<HTYPE, HASHES> that contains the base for the columns for
//     each hash, must be less than the repective mod and positive (preferably
//     in the range [mod / 4, mod / 4 * 3])
//   MOD: array<HTYPE, HASHES> that contains the mods for each hash, should
//     be large (products should be least (NM)^2) and prime
// Fields:
//   N: the number of rows in the 2D array
//   M: the number of columns in the 2D array
//   ROWBASE: array<HTYPE, HASHES> with the bases for the rows for each hash
//   COLBASE: array<HTYPE, HASHES> with the bases for the columns for each hash
//   MOD: array<HTYPE, HASHES> with the mods for each hash
//   H: a vector of vectors of array<HTYPE, HASHES> of size N + 1 by M + 1 with
//     the suffix hashes starting at each cell
//   ROWPOW: a vector of array<HTYPE, HASHES> of size N + 1 with the ith
//     representing row_base^i modulo mod
//   COLPOW: a vector of array<HTYPE, HASHES> of size M + 1 with the jth
//     representing col_base^j modulo mod
// Functions:
//   initBase(MOD): returns an array<HTYPE, HASHES> with randomly initialized
//     values for the bases for each mod in the range [mod / 4, mod / 4 * 3]
//   getHash(u, d, l, r): returns a type HASH with the hashes in the submatrix
//     [u, d] x [l, r]
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(NM HASHES)
//   getHash: O(HASHES)
// Memory Complexity: O(N HASHES)
template <class T, class HTYPE, const int HASHES> struct Hashing2D {
  using arr = array<HTYPE, HASHES>; int N, M; arr ROWBASE, COLBASE, MOD;
  vector<vector<arr>> H; vector<arr> ROWPOW, COLPOW;
  struct HASH : public arr { int rows, cols; };
  Hashing2D(const vector<vector<T>> &A, T offset, const arr &ROWBASE,
            const arr &COLBASE, const arr &MOD)
      : N(A.size()), M(N == 0 ? 0 : A[0].size()), ROWBASE(ROWBASE),
        COLBASE(COLBASE), MOD(MOD), H(N + 1, vector<arr>(M + 1)),
        ROWPOW(N + 1), COLPOW(M + 1) {
    ROWPOW[0].fill(HTYPE(1)); COLPOW[0].fill(HTYPE(1));
    for (int i = 0; i <= N; i++) H[i][M].fill(HTYPE());
    for (int j = 0; j <= M; j++) H[N][j].fill(HTYPE());
    for (int i = 1; i <= N; i++) for (int h = 0; h < HASHES; h++)
      ROWPOW[i][h] = mulMod(ROWPOW[i - 1][h], ROWBASE[h], MOD[h]);
    for (int j = 1; j <= M; j++) for (int h = 0; h < HASHES; h++)
      COLPOW[j][h] = mulMod(COLPOW[j - 1][h], COLBASE[h], MOD[h]);
    for (int i = N - 1; i >= 0; i--) for (int j = M - 1; j >= 0; j--)
      for (int h = 0; h < HASHES; h++) {
        H[i][j][h] = posMod(HTYPE(A[i][j] - offset) + HTYPE(1), MOD[h]);
        H[i][j][h] = addMod(H[i][j][h],
                            mulMod(H[i + 1][j][h], ROWBASE[h], MOD[h]),
                            MOD[h]);
        H[i][j][h] = addMod(H[i][j][h],
                            mulMod(H[i][j + 1][h], COLBASE[h], MOD[h]),
                            MOD[h]);
        H[i][j][h] = subMod(H[i][j][h],
                            mulMod(mulMod(H[i + 1][j + 1][h], ROWBASE[h],
                                          MOD[h]),
                                   COLBASE[h], MOD[h]),
                            MOD[h]);
      }
  }
  static arr initBase(const arr &MOD) {
    arr base;
    for (int h = 0; h < HASHES; h++) base[h] = uniform_int_distribution<HTYPE>(
        MOD[h] / 4, MOD[h] / 4 * 3)(rng64);
    return base;
  }
  Hashing2D(const vector<vector<T>> &A, T offset, const arr &MOD)
      : Hashing2D(A, offset, initBase(MOD), initBase(MOD), MOD) {}
  HASH getHash(int u, int d, int l, int r) {
    HASH ret; ret.rows = ++d - u; ret.cols = ++r - l;
    for (int h = 0; h < HASHES; h++) {
      ret[h] = H[u][l][h];
      ret[h] = subMod(ret[h], mulMod(H[d][l][h], ROWPOW[d - u][h], MOD[h]),
                      MOD[h]);
      ret[h] = subMod(ret[h], mulMod(H[u][r][h], COLPOW[r - l][h], MOD[h]),
                      MOD[h]);
      ret[h] = addMod(ret[h],
                      mulMod(mulMod(H[d][r][h], ROWPOW[d - u][h], MOD[h]),
                             COLPOW[r - l][h], MOD[h]),
                      MOD[h]);
    }
    return ret;
  }
};


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

// Trie Node backed by a map
// Template Arguments:
//   _T: the type of the element in the string/array for the Trie
// Fields:
//   T: the data type of each element in the string/array
//   link: a link to the parent in the trie
// Functions:
//   hasEdge(a): returns whether there exists an edge to an adjacent node
//     with the element a
//   getEdge(a): returns the index of the other node on the
//     incident edge with the element a
//   setEdge(a, n): sets the other node on the incident edge
//     with the element a to n
// Time Complexity:
//   constructor: O(1)
//   hasEdge, getEdge, setEdge: O(log E) where E is the number of edges
//     incident to this node
// Memory Complexity: O(E) where E is the number of edges incident to this node
// Tested:
//   https://csacademy.com/contest/round-77/task/expected-lcp/
template <class _T> struct TrieMapNode {
  using T = _T; int link; map<T, int> to;
  TrieMapNode() : link(-1) {}
  bool hasEdge(const T &a) const { return to.count(a); }
  int getEdge(const T &a) const { return to.at(a); }
  void setEdge(const T &a, int n) { to[a] = n; }
};

// Trie Node backed by an array
// Template Arguments:
//   _T: the type of the element in the string/array for the Trie
//   ALPHABET_SIZE: the size of the alphabet
//   OFFSET: the offset for the start of the alphabet
// Fields:
//   T: the data type of each element in the string/array
//   link: a link to the parent in the trie
// Functions:
//   hasEdge(a): returns whether there exists an edge to an adjacent node
//     with the element a
//   getEdge(a): returns the index of the other node on the
//     incident edge with the element a
//   setEdge(a, n): sets the other node on the incident edge
//     with the element a to n
// Time Complexity:
//   constructor: O(ALPHABET_SIZE)
//   hasEdge, getEdge, setEdge: O(1)
// Memory Complexity: O(ALPHABET_SIZE)
// Tested:
//   https://csacademy.com/contest/round-77/task/expected-lcp/
template <class _T, const int ALPHABET_SIZE, const int OFFSET>
struct TrieArrayNode {
  using T = _T; int link; array<int, ALPHABET_SIZE> to;
  TrieArrayNode() : link(-1) { to.fill(-1); }
  bool hasEdge(const T &a) const { return to[a - OFFSET] != -1; }
  int getEdge(const T &a) const { return to[a - OFFSET]; }
  void setEdge(const T &a, int n) { to[a - OFFSET] = n; }
};

// Trie with the root node at 0
// Template Arguments:
//   Node: a node class
//     Required Fields:
//       T: the data type of each element in the string/array
//       link: a link to the parent in the trie
//     Required Functions:
//       constructor(): initializes the node with a with no links or edges
//       hasEdge(a): returns whether there exists an edge to an adjacent node
//         with the element a
//       getEdge(a): returns the index of the other node on the
//         incident edge with the element a
//       setEdge(a, n): sets the other node on the incident edge
//         with the element a to n
// Fields:
//   TR: the vector of all nodes in the automaton
//   last: the last node representing the current string
// Functions:
//   add(a): adds the element a to the trie
//   terminate(): terminates the current word, returning to the root (which
//     allows for another word to be added with add(a))
// Time Complexity:
//   constructor: time complexity of node constructor
//   add: O(1) + time complexity of node constructor
//        + time complexity of getEdge in Node
//   terminate: O(1)
// Memory Complexity: O(N) * memory complexity of node, after N calls to add
// Tested:
//   https://csacademy.com/contest/round-77/task/expected-lcp/
template <class Node> struct Trie {
  using T = typename Node::T; vector<Node> TR; int last;
  Trie() : TR(1, Node()), last(0) {}
  void add(const T &a) {
    if (TR[last].hasEdge(a)) last = TR[last].getEdge(a);
    else {
      int u = TR.size(); TR.emplace_back();
      TR[TR[u].link = last].setEdge(a, u); last = u;
    }
  }
  void terminate() { last = 0; }
};


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

// Suffix Automaton Node backed by a map
// Template Arguments:
//   _T: the type of the element in the string/array for the Suffix Automaton
// Constructor Arguments:
//   len: the length of the longest substring represented by this node
// Fields:
//   T: the data type of each element in the string/array
//   len: the length of the longest substring this node represents
//   link: the index of the suffix link of this node (the node with the
//     longest suffix of this node)
// Functions:
//   hasEdge(a): returns whether there exists an edge to an adjacent node
//     with the element a
//   getEdge(a): returns the index of the other node on the
//     incident edge with the element a
//   setEdge(a, n): sets the other node on the incident edge
//     with the element a to n
// Time Complexity:
//   constructor: O(1)
//   hasEdge, getEdge, setEdge: O(log E) where E is the number of edges
//     incident to this node
// Memory Complexity: O(E) where E is the number of edges incident to this node
// Tested:
//   https://www.spoj.com/problems/LCS/
//   https://open.kattis.com/problems/stringmultimatching
template <class _T> struct SAMMapNode {
  using T = _T; int len, link; map<T, int> to;
  SAMMapNode(int len) : len(len), link(-1) {}
  bool hasEdge(const T &a) const { return to.count(a); }
  int getEdge(const T &a) const { return to.at(a); }
  void setEdge(const T &a, int n) { to[a] = n; }
};

// Suffix Automaton Node backed by an array
// Template Arguments:
//   _T: the type of the element in the string/array for the Suffix Automaton
//   ALPHABET_SIZE: the size of the alphabet
//   OFFSET: the offset for the start of the alphabet
// Constructor Arguments:
//   len: the length of the longest substring represented by this node
// Fields:
//   T: the data type of each element in the string/array
//   len: the length of the longest substring this node represents
//   link: the index of the suffix link of this node (the node with the
//     longest suffix of this node)
// Functions:
//   hasEdge(a): returns whether there exists an edge to an adjacent node
//     with the element a
//   getEdge(a): returns the index of the other node on the
//     incident edge with the element a
//   setEdge(a, n): sets the other node on the incident edge
//     with the element a to n
// Time Complexity:
//   constructor: O(ALPHABET_SIZE)
//   hasEdge, getEdge, setEdge: O(1)
// Memory Complexity: O(ALPHABET_SIZE)
// Tested:
//   https://www.spoj.com/problems/LCS/
//   https://dmoj.ca/problem/coci14c5p6
//   https://dmoj.ca/problem/coci11c5p6
template <class _T, const int ALPHABET_SIZE, const _T OFFSET>
struct SAMArrayNode {
  using T = _T; int len, link; array<int, ALPHABET_SIZE> to;
  SAMArrayNode(int len) : len(len), link(-1) { to.fill(-1); }
  bool hasEdge(const T &a) const { return to[a - OFFSET] != -1; }
  int getEdge(const T &a) const { return to[a - OFFSET]; }
  void setEdge(const T &a, int n) { to[a - OFFSET] = n; }
};

// Suffix Automaton with the root node at 0
// Each distinct path from the root is a substring of the set of words inserted
// Template Arguments:
//   Node: a node class
//     Required Fields:
//       T: the data type of each element in the string/array
//       len: the length of the longest substring this node represents
//       link: the index of the suffix link of this node (the node with the
//         longest suffix of this node)
//     Required Functions:
//       constructor(len): initializes the node with a length of len with no
//         links or edges
//       hasEdge(a): returns whether there exists an edge to an adjacent node
//         with the element a
//       getEdge(a): returns the index of the other node on the
//         incident edge with the element a
//       setEdge(a, n): sets the other node on the incident edge
//         with the element a to n
// Fields:
//   TR: the vector of all nodes in the automaton
//   last: the last node representing the current string
// Functions:
//   add(a): adds the element a to the automaton
//   terminate(): terminates the current word, returning to the root (which
//     allows for another word to be added with add(a))
// In practice, has a moderate constant
// Time Complexity:
//   constructor: time complexity of node constructor
//   add: O(1) amortized + time complexity of node constructor
//        + time complexity of getEdge in Node
//   terminate: O(1)
// Memory Complexity: O(N) * memory complexity of node, after N calls to add
// Tested:
//   https://www.spoj.com/problems/LCS/
//   https://dmoj.ca/problem/coci14c5p6
//   https://open.kattis.com/problems/stringmultimatching
//   https://dmoj.ca/problem/coci11c5p6
template <class Node> struct SuffixAutomaton {
  using T = typename Node::T; vector<Node> TR; int last;
  SuffixAutomaton() : TR(1, Node(0)), last(0) {}
  void add(const T &a) {
    int u = -1; if (!TR[last].hasEdge(a)) {
      u = TR.size(); TR.emplace_back(TR[last].len + 1);
      for (; last != -1 && !TR[last].hasEdge(a); last = TR[last].link)
        TR[last].setEdge(a, u);
      if (last == -1) { TR[last = u].link = 0; return; }
    }
    int p = TR[last].getEdge(a); if (TR[p].len == TR[last].len + 1) {
      (u == -1 ? last : TR[last = u].link) = p; return;
    }
    int q = TR.size(); TR.push_back(TR[p]); TR[q].len = TR[last].len + 1;
    TR[p].link = q; if (u != -1) TR[u].link = q;
    while (last != -1 && TR[last].hasEdge(a) && TR[last].getEdge(a) == p) {
      TR[last].setEdge(a, q); last = TR[last].link;
    }
    last = u == -1 ? q : u;
  }
  void terminate() { last = 0; }
};


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

// Finds the minimum lexicographical rotation of an array
// Template Arguments:
//   T: the type of each element in the array
// Function Arguments:
//   S: a vector of type T
// Return Value: the index of the first element of the minimum
//   lexicographical rotation of the array
// In practice, has a small constant
// Time Complexity: O(N)
// Memory Complexity: O(1)
// Tested:
//   https://cses.fi/problemset/task/1110/
template <class T> int minRotation(const vector<T> &S) {
  int N = S.size(), i = 0;
  auto ind = [&] (int i) { return i < N ? i : i - N; };
  for (int j = 0; j < N; j++) for (int k = 0; k < N; k++) {
    auto &a = S[ind(i + k)], &b = S[ind(j + k)];
    if (i + k == j || a < b) { j += max(0, k - 1); break; }
    if (a > b) { i = j; break; }
  }
  return i;
}


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

// Computes the hash of an array to allow for easy computation of
//   subarray hashes
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   HTYPE: the type of the hash
//   HASHES: the number of hashes to compute
// Construtor Arguments:
//   S: the elements in the array
//   offset: the amount to subtract from each element, a - offset + 1 should be
//     positive
//   BASE: array<HTYPE, HASHES> that contains the base for each hash, must
//     be less than the repective mod and positive (preferably in the range
//     [mod / 4, mod / 4 * 3])
//   MOD: array<HTYPE, HASHES> that contains the mods for each hash, should
//     be large (products should be least N^2) and prime
// Fields:
//   N: the length of the array
//   BASE: array<HTYPE, HASHES> with the bases for each hash
//   MOD: array<HTYPE, HASHES> with the mods for each hash
//   H: a vector of array<HTYPE, HASHES> of size N + 1 with the suffix hashes
//     starting at each index
//   POW: a vector of array<HTYPE, HASHES> of size N + 1 with the ith
//     representing base^i modulo mod
// Functions:
//   initBase(MOD): returns an array<HTYPE, HASHES> with randomly initialized
//     values for the bases for each mod in the range [mod / 4, mod / 4 * 3]
//   getHash(l, r): returns a type HASH with the hashes in the subarray [l, r]
//   merge(h1, h2): returns a type HASH of 2 hashes merged together as if
//     the subarray represented by h2 was immediately after h1
//   concat(l1, r1, l2, r2): returns a type HASH of
//     the subarray [l1, r1] contactenated with the subarray [l2, r2]
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(N HASHES)
//   getHash, merge, concat: O(HASHES)
// Memory Complexity: O(N HASHES)
// Tested:
//   https://dmoj.ca/problem/globexcup19j5hard
template <class T, class HTYPE, const int HASHES> struct Hashing {
  using arr = array<HTYPE, HASHES>; int N; arr BASE, MOD; vector<arr> H, POW;
  struct HASH : public arr { int len; };
  Hashing(const vector<T> &S, T offset, const arr &BASE, const arr &MOD)
      : N(S.size()), BASE(BASE), MOD(MOD), H(N + 1), POW(N + 1) {
    POW[0].fill(HTYPE(1)); H[N].fill(HTYPE());
    for (int i = 1; i <= N; i++) for (int h = 0; h < HASHES; h++)
      POW[i][h] = mulMod(POW[i - 1][h], BASE[h], MOD[h]);
    for (int i = N - 1; i >= 0; i--) for (int h = 0; h < HASHES; h++)
      H[i][h] = addMod(mulMod(H[i + 1][h], BASE[h], MOD[h]),
                       posMod(HTYPE(S[i] - offset) + HTYPE(1), MOD[h]),
                       MOD[h]);
  }
  static arr initBase(const arr &MOD) {
    arr base;
    for (int h = 0; h < HASHES; h++) base[h] = uniform_int_distribution<HTYPE>(
        MOD[h] / 4, MOD[h] / 4 * 3)(rng64);
    return base;
  }
  Hashing(const vector<T> &S, T offset, const arr &MOD)
      : Hashing(S, offset, initBase(MOD), MOD) {}
  HASH getHash(int l, int r) {
    HASH ret; ret.len = ++r - l; for (int h = 0; h < HASHES; h++)
      ret[h] = subMod(H[l][h], mulMod(H[r][h], POW[r - l][h], MOD[h]), MOD[h]);
    return ret;
  }
  HASH merge(const HASH &h1, const HASH &h2) {
    HASH ret; ret.len = h1.len + h2.len; for (int h = 0; h < HASHES; h++)
      ret[h] = addMod(h1[h], mulMod(h2[h], POW[h1.len][h], MOD[h]), MOD[h]);
    return ret;
  }
  HASH concat(int l1, int r1, int l2, int r2) {
    return merge(getHash(l1, r1), getHash(l2, r2));
  }
};


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

// Finds the longest common subarray of multiple arrays
// Indices are 0-indexed
// Template Arguments:
//   SuffixArray: a generic suffix array to be used (should be either
//       SuffixArray or SAISSuffixArray)
//     Required Fields:
//       N: the length of the array
//       rnk: a vector of the ranks of the suffixes (rnk[i] is the rank of the
//         suffix starting from index i)
//       ind: a vector of the indices in the original array of the suffixes
//         sorted in lexicographical order (ind[i] is the index in original
//         array of the ith lexicographically smallest suffix)
//       LCP: a vector of the longest common prefixes between the suffixes when
//         sorted in lexicographical order (LCP[i] is the longest common prefix
//         of the ith and (i + 1)th lexicographically smallest suffix, with
//         LCP[N - 1] being 0)
//     Required Functions:
//       constructor(S): construts a suffix array from the array S
// Constructor Arguments:
//   A: a vector of the arrays
//   sep: a vector of max(0, A.size() - 1) different separator character
//     that do not appear in any of the arrays
// Fields:
//   SA: the associated suffix array constructed from the arrays
//     concatenated together, separated by the separator character
//   lcs: the longest common subarray of all arrays
// In practice, the constructor has a very small constant
// Time Complexity:
//   constructor: time complexity of SuffixArray for the length of the
//     concatenated array
// Memory Complexity: memory complexity of SuffixArray for the length of the
//   concatenated array
// Tested:
//   https://www.spoj.com/problems/LCS2/
//   https://open.kattis.com/problems/longestcommonsubstring
template <class SuffixArray> struct LongestCommonSubstring {
  using T = typename SuffixArray::T;
  vector<int> ind; vector<T> lcs; vector<T> S; SuffixArray SA;
  vector<T> init(const vector<vector<T>> &A, const vector<T> &sep) {
    int len = max(0, int(A.size()) - 1); for (auto &&a : A) len += a.size();
    vector<T> S; S.reserve(len); ind.resize(len, -1);
    for (int i = 0; i < int(A.size()); i++) {
      if (i > 0) S.push_back(sep[i - 1]);
      fill(ind.begin() + S.size(), ind.begin() + S.size() + A[i].size(), i);
      S.insert(S.end(), A[i].begin(), A[i].end());
    }
    return S;
  }
  LongestCommonSubstring(const vector<vector<T>> &A, const vector<T> &sep)
      : ind(), S(init(A, sep)), SA(S) {
    int K = A.size(); if (K == 0) return;
    if (K == 1) { lcs = A[0]; return; }
    vector<int> freq(K, 0), dq(SA.N); pair<int, int> mx(-1, -1);
    for (int l = 0, r = 0, front = 0, back = 0, cnt = 0; r < SA.N; r++) {
      if (ind[SA.ind[r]] != -1) {
        cnt += freq[ind[SA.ind[r]]]++ == 0; if (cnt == K) {
          while (true) {
            int i = ind[SA.ind[l]]; if (i != -1 && freq[i] == 1) break;
            l++; if (i != -1) freq[i]--;
          }
          while (dq[front] < l) front++;
          mx = max(mx, make_pair(SA.LCP[dq[front]], SA.ind[dq[front]]));
        }
      }
      while (front < back && SA.LCP[dq[back - 1]] >= SA.LCP[r]) back--;
      dq[back++] = r;
    }
    if (mx.first != -1) {
      lcs.reserve(mx.first);
      for (int i = 0; i < mx.first; i++) lcs.push_back(S[mx.second + i]);
    }
  }
};


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

// Computes the longest palindromic subarray centered at each half index
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element in the array
// Constructor Arguments:
//   S: a vector of type T
// Fields:
//   N: the length of the array
//   p: the length of the longest palindromic substring/subarray centered
//     at each half index (including paddings at each end)
// Functions:
//   lps(): returns a pair containing the starting index and length (in the
//     original string) of the longest palindromic substring/subarray, picking
//     the earliest string if there are multiple of the same length
//   lps(i): returns a pair containing the starting index and length (in the
//     original string) of the longest palindromic substring/subarray centered
//     at index i/2.0
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(N)
//   lps(): O(N)
//   lps(i): O(1)
// Memory Complexity: O(S)
// Tested:
//   https://judge.yosupo.jp/problem/enumerate_palindromes
//   https://www.spoj.com/problems/LPS/
template <class T> struct ManacherPalindrome {
  int N; vector<int> p;
  ManacherPalindrome(const vector<T> &S) : N(S.size()), p(N * 2 + 1, 0) {
    vector<T> SS(N * 2 + 1, T());
    for (int i = 0; i < N; i++) SS[i * 2 + 1] = S[i];
    for (int i = 0, cen = 0, mxr = 0; i < N * 2 + 1; i++) {
      if (mxr > i) p[i] = min(mxr - i, p[cen * 2 - i]);
      int l = i - p[i], r = i + p[i];
      for (; l > 0 && r < N * 2 && SS[l - 1] == SS[r + 1]; l--, r++) p[i]++;
      if (r > mxr) { cen = i; mxr = r; }
    }
  }
  pair<int, int> lps() {
    int len = 0, cen = 0;
    for (int i = 0; i < N * 2 + 1; i++) if (p[i] > len) len = p[cen = i];
    return make_pair((cen - len) / 2, len);
  }
  pair<int, int> lps(int i) {
    int len = p[i + 1], cen = i + 1; return make_pair((cen - len) / 2, len);
  }
};


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

// Two players alternate taking coins from either end of a row. If both players
//   play optimally, what is the value of player 1's total minus
//   player 2's total?
// The standard dp recurrence is
//   dp[l][r] = l > r ? 0 : max(A[l] - dp[l + 1][r], A[r] - dp[l][r - 1])
// This implementation is generalized to support any number of rows, possibly
//   with a fixed endpoint on the right
// The value of a single player can be recovered by solving the system of
//   equations: P1 + P2 = A, P1 - P2 = B
// Template Arguments:
//   T: the type of the value of each coin
// Function Arguments:
//   A: a vector of vectors of type T, with each vector representing a row
//   fixed: a vector of booleans of the same length of A representing whether
//     the right endpoint of that row is fixed
// Return Value: player 1's total minus player 2's total
// In practice, has a small constant
// Time Complexity: O(N log N) for N total coins or O(N) if there is a
//                  single row
// Memory Complexity: O(N) for N total coins
// Tested:
//   https://dmoj.ca/problem/ahardergame
//   https://atcoder.jp/contests/dp/tasks/dp_l
template <class T>
T termity(const vector<vector<T>> &A, const vector<bool> &fixed) {
  int N = 0; T ret = T(); for (auto &&a : A) N += a.size();
  int mult = N % 2 == 0 ? 1 : -1; vector<T> B; B.reserve(N);
  for (int i = 0; i < int(A.size()); i++) {
    int front = B.size(); for (auto &&a : A[i]) {
      B.push_back(a); while (int(B.size()) - front >= 3
          && B[B.size() - 3] <= B[B.size() - 2]
          && B[B.size() - 2] >= B[B.size() - 1]) {
        B[B.size() - 3] -= B[B.size() - 2] - B[B.size() - 1];
        B.pop_back(); B.pop_back();
      }
    }
    if (fixed[i]) while (int(B.size()) - front >= 2
        && B[B.size() - 2] <= B[B.size() - 1]) {
      ret += mult * (B[B.size() - 2] - B[B.size() - 1]);
      B.pop_back(); B.pop_back();
    }
  }
  if (int(A.size()) == 1) {
    for (int l = 0, r = int(B.size()) - 1, m = 1; l <= r; m *= -1) {
      if (B[l] >= B[r]) ret += m * B[l++];
      else ret += m * B[r--];
    }
  } else {
    sort(B.begin(), B.end(), greater<T>());
    for (int i = 0, m = 1; i < int(B.size()); i++, m *= -1) ret += m * B[i];
  }
  return ret;
}


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

// Solves the maximum disjoint intervals problem
// Given a set of intervals in the form [L, R], find the maximum number of
//   disjoint intervals
// Maximum number of disjoint intervals is equivalent to the minimum number of
//   points to cover each interval (with the points being the right endpoints
//   of the disjoint intervals)
// Range is modified in-place
// Template Arguments:
//   T: the type of the endpoints of the intervals
//   Cmp: the comparator to compare two points
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Function Arguments:
//   A: a reference to a vector of pairs with the first element being the
//     inclusive left bound of the interval and the second element being the
//     inclusive right bound of the interval
//   cmp: an instance of the Cmp struct
// Return Value: a reference to the modified vector
// In practice, has a very small constant
// Time Complexity: O(N log N)
// Memory Complexity: O(1)
// Tested:
//   https://codeforces.com/contest/1141/problem/F2
//   https://oj.uz/problem/view/COCI21_planine
template <class T, class Cmp = less<T>>
vector<pair<T, T>> &maxDisjointIntervals(vector<pair<T, T>> &A,
                                         Cmp cmp = Cmp()) {
  sort(A.begin(), A.end(), [&] (const pair<T, T> &a, const pair<T, T> &b) {
    return cmp(a.second, b.second);
  });
  int i = 0; for (int l = 0, r = 0, N = A.size(); l < N; l = r, i++) {
    A[i] = A[l]; for (r = l + 1; r < N && !cmp(A[i].second, A[r].first); r++);
  }
  A.erase(A.begin() + i, A.end()); return A;
}


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

// Solves the minimum interval cover problem
// Given a set of intervals in the form [L, R], find the minimum number of
//   intervals to cover a target interval
// Range is modified in-place
// Template Arguments
//   T: the type of the endpoints of the intervals
//   Cmp: the comparator to compare two points
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
//   T: the type of the points
// Function Arguments:
//   A: a reference to a vector of pairs with the first element being the
//     inclusive left bound of the interval and the second element being the
//     inclusive right bound of the interval
//   target: the target interval to cover with the first element being the
//     inclusive left bound of the target interval and the second element
//     being the inclusive right bound of the interval
//   cmp: an instance of the Cmp struct
// Return Value: a reference to the modified vector
// In practice, has a very small constant
// Time Complexity: O(N log N)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/intervalcover
template <class T, class Cmp = less<T>>
vector<pair<T, T>> &minIntervalCover(vector<pair<T, T>> &A,
                                     pair<T, T> target, Cmp cmp = Cmp()) {
  sort(A.begin(), A.end(), [&] (const pair<T, T> &a, const pair<T, T> &b) {
    return cmp(a.first, b.first);
  });
  bool first = true; int i = 0, N = A.size();
  for (int j = 0; j < N && (first || cmp(target.first, target.second));) {
    if (cmp(target.first, A[j].first)) { A.clear(); return A; }
    else {
      for (A[i] = A[j]; j < N && !cmp(target.first, A[j].first); j++)
        if (cmp(A[i].second, A[j].second)) A[i] = A[j];
      target.first = A[i++].second; first = false;
    }
  }
  if (first || cmp(target.first, target.second)) i = 0;
  A.erase(A.begin() + i, A.end()); return A;
}


text:
#pragma once
#include <bits/stdc++.h>

// Unsyncs the C++ and C std output streams and unties cin to
//   speed up standard input
void unSyncUntie() {
  std::ios::sync_with_stdio(0); std::cin.tie(0);
}

// Functions for fast IO
namespace IO {
#define INTERACTIVE_INPUT 0
  template <class T> struct is_iterator {
    template <class U, typename std::enable_if<
        !std::is_convertible<U, const char *>::value, int>::type = 0>
    static constexpr auto has_indirection(int)
        ->decltype(*std::declval<U>(), bool()) {
      return true;
    }
    template <class>
    static constexpr bool has_indirection(long) { return false; }
    static constexpr bool value = has_indirection<T>(0);
  };
  constexpr const int _bufSize = 1 << 16, _maxNumLength = 128;
  char _inputBuf[_bufSize + 1], *_inputPtr = _inputBuf, _sign, _c, _last = -1;
  char *_tempInputBuf = nullptr, _tempOutputBuf[_maxNumLength];
  char _outputBuf[_bufSize], _numBuf[_maxNumLength], _fill = ' ';
  int _cur, _tempOutputPtr = 0, _outputPtr = 0, _cnt, _numPtr = 0;
  int _width = 0, _precision = 9;
  const char *_delimiter = " "; unsigned long long _precisionBase = 1000000000;
  FILE *_input = stdin, *_output = stdout, *_error = stderr;
#if INTERACTIVE_INPUT
  char _getchar() { return _last = getchar(); }
#else
  char _getchar() {
    if (!*_inputPtr)
      _inputBuf[fread(_inputPtr = _inputBuf, 1, _bufSize, _input)] = '\0';
    return _last = *_inputPtr++;
  }
#endif
  char _getcharskipr() {
    while (_getchar() == '\r');
    return _last;
  }
  template <class I> void _readSigned(I &x) {
    while ((x = _getchar()) <= ' ');
    if ((_sign = x == '-')) x = _getchar();
    for (x -= '0'; (_c = _getchar()) >= '0'; x = x * 10 + _c - '0');
    if (_sign) x = -x;
  }
  template <class UI> void _readUnsigned(UI &x) {
    while ((x = _getchar()) <= ' ');
    for (x -= '0'; (_c = _getchar()) >= '0'; x = x * 10 + _c - '0');
  }
  template <class F> void _readFloatingPoint(F &x) {
    for (x = 0; (_c = _getchar()) <= ' ';);
    if ((_sign = _c == '-')) _c = _getchar();
    if (_c >= '0')
      for (x = _c - '0'; (_c = _getchar()) >= '0'; x = x * 10 + _c - '0');
    if (_c == '.') {
      F _div = 1.0;
      for (; (_c = _getchar()) >= '0'; x += (_c - '0') / (_div *= 10));
    }
    if (_sign) x = -x;
  }
  void setLength(int x) {
    if (_tempInputBuf) delete[](_tempInputBuf);
    _tempInputBuf = new char[x + 1];
  }
  template <class I> typename std::enable_if<std::is_integral<I>::value
      && std::is_signed<I>::value>::type read(I &x) { _readSigned(x); }
  template <class UI> typename std::enable_if<std::is_integral<UI>::value
      && std::is_unsigned<UI>::value>::type read(UI &x) {
    _readUnsigned(x);
  }
#if __SIZEOF_INT128__
  void read(__int128_t &x) { _readSigned(x); }
  void read(__uint128_t &x) { _readUnsigned(x); }
#endif
  template <class F>
  typename std::enable_if<std::is_floating_point<F>::value>::type read(F &x) {
    _readFloatingPoint(x);
  }
  void read(char &x) { while ((x = _getchar()) <= ' '); }
  void read(char *x) {
    _cur = 0; do { _c = _getchar(); } while (_c <= ' ');
    do { x[_cur++] = _c; } while ((_c = _getchar()) > ' ');
    x[_cur] = '\0';
  }
  void readln(char *x) {
    if (_last == '\r') _getcharskipr();
    for (_cur = 0; (_c = _getcharskipr()) != '\n' && _c; x[_cur++] = _c);
    x[_cur] = '\0';
  }
  void read(std::string &x) {
    if (!_tempInputBuf) assert(0);
    read(_tempInputBuf); x = std::string(_tempInputBuf, _cur);
  }
  void readln(std::string &x) {
    if (!_tempInputBuf) assert(0);
    readln(_tempInputBuf); x = std::string(_tempInputBuf, _cur);
  }
  template <class T> typename std::enable_if<
      is_iterator<typename T::iterator>::value>::type read(T &x);
  template <class T1, class T2>
  void read(std::pair<T1, T2> &x) { read(x.first); read(x.second); }
  template <class T> void read(std::complex<T> &x) {
    T _re, _im; read(_re); read(_im); x.real(_re); x.imag(_im);
  }
  template <class T, class ...Ts> void read(T &x, Ts &&...xs);
  template <class It>
  typename std::enable_if<is_iterator<It>::value>::type read(It st, It en) {
    for (It _i = st; _i != en; _i++) read(*_i);
  }
  template <class It, class ...Ts>
  typename std::enable_if<is_iterator<It>::value>::type read(
      It st, It en, Ts &&...xs) {
    read(st, en); read(std::forward<Ts>(xs)...);
  }
  template <class T> typename std::enable_if<
      is_iterator<typename T::iterator>::value>::type read(T &x) {
    for (auto &&_i : x) read(_i);
  }
  template <class T, class...Ts> void read(T &x, Ts &&...xs) {
    read(x); read(std::forward<Ts>(xs)...);
  }
  void setInput(FILE *file) { *(_inputPtr = _inputBuf) = 0; _input = file; }
  void setInput(const std::string &s) {
    *(_inputPtr = _inputBuf) = 0; _input = fopen(s.c_str(), "r");
  }
  int _flushBuf() {
    fwrite(_outputBuf, 1, _outputPtr, _output); return _outputPtr = 0;
  }
  void flush() { _flushBuf(); fflush(_output); }
  int _putchar(char x) {
    _outputBuf[_outputPtr == _bufSize ? _flushBuf() : _outputPtr] = x;
    return _outputPtr++;
  }
  void _writeTempBuf(char x) { _tempOutputBuf[_tempOutputPtr++] = x; }
  void _writeOutput() {
    for (int _i = 0; _i < _tempOutputPtr; _putchar(_tempOutputBuf[_i++]));
    _tempOutputPtr = 0;
  }
  void _fillBuf(int x) { for (int _i = 0; _i < x; _i++) _putchar(_fill); }
  void _flushNumBuf() { for (; _numPtr; _writeTempBuf(_numBuf[--_numPtr])); }
  void _flushTempBuf() {
    int _tempLen = _tempOutputPtr; _fillBuf(_width - _tempLen);
    _writeOutput(); _fillBuf(-_width - _tempLen);
  }
  void setPrecision(int x) {
    _precision = x; _precisionBase = 1;
    for (int _i = 0; _i < x; _i++, _precisionBase *= 10);
  }
  void setWidth(int x) { _width = x; }
  void setFill(char x) { _fill = x; }
  void setDelimiter(const char *x) { _delimiter = x; }
  void setDelimiter(const std::string &x) { _delimiter = x.c_str(); }
  void writeDelimiter() {
    for (const char *_p = _delimiter; *_p; _putchar(*_p++));
  }
  template <class T> void _writeNum(const T &x, int digits) {
    _cnt = 0;
    for (T _y = x; _y; _y /= 10, _cnt++) _numBuf[_numPtr++] = '0' + _y % 10;
    for (; _cnt < digits; _cnt++) _numBuf[_numPtr++] = '0';
    _flushNumBuf();
  }
  template <class F> void _writeFloatingPoint(const F &x) {
    unsigned long long _I = x, _F = (x - _I) * _precisionBase + F(0.5);
    if (_F >= _precisionBase) { _I++; _F = 0; }
    _writeNum(_I, 1); _writeTempBuf('.'); _writeNum(_F, _precision);
  }
  void write(const bool &x) {
    if (x) _writeTempBuf('1');
    else _writeTempBuf('0');
    _flushTempBuf();
  }
  void write(const char &x) { _writeTempBuf(x); _flushTempBuf(); }
  void write(const char *x) {
    int _slen = strlen(x); _fillBuf(_width - _slen);
    for (const char *_p = x; *_p; _putchar(*_p++));
    _fillBuf(-_width - _slen);
  }
  void write(const std::string &x) {
    _fillBuf(_width - int(x.length()));
    for (int _i = 0; _i < int(x.length()); _putchar(x[_i++]));
    _fillBuf(-_width - int(x.length()));
  }
  template <class I>
  typename std::enable_if<std::is_integral<I>::value
      && std::is_signed<I>::value>::type write(const I &x) {
    using UI = typename std::make_unsigned<I>::type;
    if (x < 0) { _writeTempBuf('-'); _writeNum(UI(-x), 1); }
    else { _writeNum(UI(x), 1); }
    _flushTempBuf();
  }
  template <class UI>
  typename std::enable_if<std::is_integral<UI>::value
      && std::is_unsigned<UI>::value >::type write(const UI &x) {
    _writeNum(x, 1); _flushTempBuf();
  }
  template <class F>
  typename std::enable_if<std::is_floating_point<F>::value>::type write(
      const F &x) {
    if (std::isnan(x)) write("NaN");
    else if (std::isinf(x)) write("Inf");
    else if (x < 0) { _writeTempBuf('-'); _writeFloatingPoint(-x); }
    else _writeFloatingPoint(x);
    _flushTempBuf();
  }
#if __SIZEOF_INT128__
  void write(const __int128_t &x) {
    if (x < 0) { _writeTempBuf('-'); _writeNum(__uint128_t(-x), 1); }
    else _writeNum(__uint128_t(x), 1);
    _flushTempBuf();
  }
  void write(const __uint128_t &x) { _writeNum(x, 1); _flushTempBuf(); }
#endif
  template <class T> typename std::enable_if<
      is_iterator<typename T::iterator>::value>::type write(const T &x);
  template <class T1, class T2> void write(const std::pair<T1, T2> &x) {
    write(x.first); writeDelimiter(); write(x.second);
  }
  template <class T> void write(const std::complex<T> &x) {
    write(x.real()); writeDelimiter(); write(x.imag());
  }
  template <class T, class ...Ts> void write(const T &x, Ts &&...xs);
  template <class It>
  typename std::enable_if<is_iterator<It>::value>::type write(It st, It en) {
    bool _first = 1; for (It _i = st; _i != en; _i++) {
      if (_first) _first = 0;
      else writeDelimiter();
      write(*_i);
    }
  }
  template <class It, class ...Ts>
  typename std::enable_if<is_iterator<It>::value>::type write(
      It st, It en, Ts &&...xs) {
    write(st, en); writeDelimiter(); write(std::forward<Ts>(xs)...);
  }
  template <class T> typename std::enable_if<
      is_iterator<typename T::iterator>::value>::type write(const T &x) {
    bool _first = 1; for (auto &&_i : x) {
      if (_first) _first = 0;
      else writeDelimiter();
      write(_i);
    }
  }
  template <class T, class ...Ts> void write(const T &x, Ts &&...xs) {
    write(x); writeDelimiter(); write(std::forward<Ts>(xs)...);
  }
  void writeln() { _putchar('\n'); }
  template <class ...Ts>
  void writeln(Ts &&...xs) { write(std::forward<Ts>(xs)...); _putchar('\n'); }
  struct IOManager {
    ~IOManager() { flush(); if (_tempInputBuf) delete[](_tempInputBuf); }
  };
  std::unique_ptr<IOManager> iomanager = std::unique_ptr<IOManager>(
      new IOManager());
  void setOutput(FILE *file) { flush(); _output = file; }
  void setOutput(const std::string &s) {
    flush(); _output = fopen(s.c_str(), "w");
  }
  template <class ...Ts> void debug(Ts &&...xs) {
    FILE *_temp = _output; setOutput(_error);
    write(std::forward<Ts>(xs)...); setOutput(_temp);
  }
  template <class...Ts> void debugln(Ts &&...xs) {
    FILE *_temp = _output; setOutput(_error);
    writeln(std::forward<Ts>(xs)...); setOutput(_temp);
  }
  void setError(FILE *file) { flush(); _error = file; }
  void setError(const std::string &s) {
    flush(); _error = fopen(s.c_str(), "w");
  }
#undef INTERACTIVE_INPUT
}


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

// With C++ 11 not having std::make_unique defined, a simple implementation
//   of make_unique is provided
// Template Arguments:
//   T: the type of the object to construct
//   ...Args: variadic arguments of the types of the arguments being passed to
//     the object being constructor
// Return Value:
//   a unique_ptr of type T
// Tested:
//   https://dmoj.ca/problem/set
template <class T, class ...Args> unique_ptr<T> _make_unique(Args &&...args) {
  return unique_ptr<T>(new T(forward<Args>(args)...));
}

#define make_unique _make_unique


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

// Functions for __int128_t
__int128_t abs(__int128_t a) { return a >= 0 ? a : -a; }
#define FUN(name) \
  long double name(__int128_t a) { return name((long double)(a)); }
FUN(sqrt) FUN(sin) FUN(cos) FUN(tan) FUN(asin) FUN(acos) FUN(atan)
#undef FUN
long double atan2(__int128_t y, __int128_t x) {
  return atan2((long double)(y), (long double)(x));
}

// 128-bit integers IO operations
__uint128_t stoui128(const string &s) {
  __uint128_t x = 0;
  for (int i = 0; i < int(s.length()); i++) x = x * 10 + (s[i] - '0');
  return x;
}

__int128_t stoi128(const string &s) {
  __int128_t x = 0; int i = 0; bool neg = false;
  if (s[0] == '-') { neg = true; i++; }
  for (; i < int(s.length()); i++) x = x * 10 + (s[i] - '0');
  if (neg) x *= -1;
  return x;
}

istream &operator >> (istream &stream, __uint128_t &x) {
  string s; stream >> s; x = stoui128(s); return stream;
}

istream &operator >> (istream &stream, __int128_t &x) {
  string s; stream >> s; x = stoi128(s); return stream;
}

string to_string(__uint128_t x) {
  if (x == 0) return string("0");
  string s; for (; x > 0; x /= 10) s.push_back('0' + x % 10);
  reverse(s.begin(), s.end()); return s;
}

string to_string(__int128_t x) {
  string s = to_string(__uint128_t(x >= 0 ? x : -x));
  if (x < 0) s.insert(s.begin(), '-');
  return s;
}

ostream &operator << (ostream &stream, __uint128_t x) {
  return stream << to_string(x);
}

ostream &operator << (ostream &stream, __int128_t x) {
  return stream << to_string(x);
}


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

// Static allocator to improve the speed of memory access
// Tested:
//   https://dmoj.ca/problem/ioi12p3
constexpr const int MB = 200;
char buf[MB << 20]; size_t buf_ind = sizeof(buf);

// Overloading new and delete operators
void *operator new (size_t n) { return (void *) &buf[buf_ind -= n]; }
void *operator new[] (size_t n) { return (void *) &buf[buf_ind -= n]; }
void operator delete (void *) {}
void operator delete[] (void *) {}
void operator delete (void *, size_t) {}
void operator delete[] (void *, size_t) {}

// Allocator class to be used with stl data structures
template <class T> struct StaticAllocator {
  typedef T value_type; StaticAllocator() {}
  template <class U> StaticAllocator(const U &) {}
  T *allocate(size_t n) {
    (buf_ind -= n * sizeof(T)) &= 0 - alignof(T); return (T*) (buf + buf_ind);
  }
  void deallocate(T *, size_t) {}
};

// 32-bit pointer
#define CMP(op, body) bool operator op (small_ptr p) const { return body; }
template <class T> struct small_ptr {
  unsigned ind;
  small_ptr(T *p = 0) : ind(p ? unsigned((char *) p - buf) : 0) {}
  T &operator * () const { return *(T *) (buf + ind); }
  T *operator -> () const { return &**this; }
  T &operator [] (int a) const { return (&**this)[a]; }
  explicit operator bool () const { return ind; }
  bool operator < (small_ptr p) const { return ind < p.ind; }
  CMP(<=, !(p < *this)) CMP(>, p < *this) CMP(>=, !(*this < p))
  CMP(==, !(*this < p) && !(p < *this)) CMP(!=, *this < p || p < *this)
};
#undef CMP


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

// Functions for epsilon comparison of floating points
using T = long double;
constexpr const T EPS = 1e-9;
bool lt(T a, T b) { return a + EPS < b; }
bool le(T a, T b) { return !lt(b, a); }
bool gt(T a, T b) { return lt(b, a); }
bool ge(T a, T b) { return !lt(a, b); }
bool eq(T a, T b) { return !lt(a, b) && !lt(b, a); }
bool ne(T a, T b) { return lt(a, b) || lt(b, a); }
int sgn(T a) { return lt(a, 0) ? -1 : lt(0, a) ? 1 : 0; }
struct eps_lt { bool operator () (T a, T b) const { return lt(a, b); } };
struct eps_le { bool operator () (T a, T b) const { return !lt(b, a); } };
struct eps_gt { bool operator () (T a, T b) const { return lt(b, a); } };
struct eps_ge { bool operator () (T a, T b) const { return !lt(a, b); } };
struct eps_eq {
  bool operator () (T a, T b) const { return !lt(a, b) && !lt(b, a); }
};
struct eps_ne {
  bool operator () (T a, T b) const { return lt(a, b) || lt(b, a); }
};


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

// 32-bit and 64-bit random number generators using a time-based seed sequence
// Tested:
//   https://dmoj.ca/problem/set
//   https://atcoder.jp/contests/agc026/tasks/agc026_c
//   https://dmoj.ca/problem/ds4
//   https://judge.yosupo.jp/problem/associative_array
seed_seq seq {
  uint64_t(chrono::duration_cast<chrono::nanoseconds>(
      chrono::steady_clock::now().time_since_epoch()).count()),
  uint64_t(__builtin_ia32_rdtsc()),
  uint64_t(uintptr_t(make_unique<char>().get()))
};

mt19937 rng(seq); mt19937_64 rng64(seq);


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

// Combine struct used for maximum non empty subarray
//   for Segment Trees, Dynamic Range Operations, Link Cut Trees, etc
// Template Arguments:
//   T: the type of the element
// Functions:
//   makeData(v): returns a MaxSubarraySumCombine<T>::Data initialized with v
// Time Complexity:
//   all functions: O(1)
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://www.spoj.com/problems/GSS1/
//   https://dmoj.ca/problem/dmpg17g2
//   https://mcpt.ca/problem/seq3
//   https://dmoj.ca/problem/acc1p1
//   https://dmoj.ca/problem/noi05p2
template <class T> struct MaxSubarraySumCombine {
  struct Data { T pre, suf, sum, maxSum; };
  using Lazy = T;
  static Data makeData(const T &v) {
    Data ret; ret.pre = ret.suf = ret.sum = ret.maxSum = v; return ret;
  }
  static Data qdef() {
    Data ret = makeData(numeric_limits<T>::lowest()); ret.sum = T();
    return ret;
  }
  static Lazy ldef() { return numeric_limits<T>::lowest(); }
  static Data merge(const Data &l, const Data &r) {
    if (l.maxSum == numeric_limits<T>::lowest()) return r;
    if (r.maxSum == numeric_limits<T>::lowest()) return l;
    Data ret; ret.pre = max(l.pre, l.sum + r.pre);
    ret.suf = max(l.suf + r.sum, r.suf); ret.sum = l.sum + r.sum;
    ret.maxSum = max(max(l.maxSum, r.maxSum), l.suf + r.pre); return ret;
  }
  template <class IndexType>
  static Data applyLazy(const Data &, const Lazy &r, IndexType k) {
    Data ret; ret.pre = ret.suf = ret.maxSum = max(Lazy(r * k), r);
    ret.sum = r * k; return ret;
  }
  static Lazy mergeLazy(const Lazy &, const Lazy &r) { return r; }
  static void revData(Data &v) { swap(v.pre, v.suf); }
  template <class IndexType> static Data getSegmentVdef(IndexType) {
    return makeData(IndexType());
  }
};


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

// Maintains the aggregate value of a sequence of elements over an associative
//   operation where elements can be pushed to the back and popped from the
//   front of the queue
// Template Arguments:
//   T: the type of each element
//   Op: a struct with the operation (can also be of type
//       std::function<T(T, T)>); in practice, custom struct is faster than
///      std::function
//     Required Functions:
//       operator (l, r): merges the values l and r, must be associative
// Constructor Arguments:
//   qdef: the query default value
//   op: an instance of the Op struct
// Functions:
//   reserve(v): reserves space for N elements
//   push(v): pushes the value v into the queue
//   getAgg(): returns the aggregate value of the elements in the queue
//     aggregated in the order they were pushed
//   pop(): pops from the front of the queue
// In practice, has a moderate constant
// Time Complexity:
//   constructor, getAgg: O(1)
//   reserve: O(N)
//   push: O(1) if reserved is called beforehand, O(1) amortized otherwise
//   pop: O(1) amortized
// Memory Complexity: O(N)
// Tested:
//   https://judge.yosupo.jp/problem/queue_operate_all_composite
template <class T, class Op> struct SWAG {
  vector<T> q; T qdef, backAgg; int front, mid; Op op;
  SWAG(const T &qdef, Op op = Op())
      : qdef(qdef), backAgg(qdef), front(0), mid(0), op(op) {}
  void reserve(int N) { q.reserve(N); }
  void push(const T &v) { q.push_back(v); backAgg = op(backAgg, v); }
  T getAgg() const { return front == mid ? backAgg : op(q[front], backAgg); }
  void pop() {
    if (front++ < mid) return;
    for (int i = int(q.size()) - 2; i >= mid; i--) q[i] = op(q[i], q[i + 1]);
    mid = q.size(); backAgg = qdef;
  }
};


text:
#pragma once
#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/priority_queue.hpp>
#include <ext/pb_ds/tree_policy.hpp>
#include "RandomizedHash.h"
using namespace std;
using namespace __gnu_pbds;

// Policy-based data structures for ordered and unordered sets and maps,
//   as well as priority queues

// Unordered hashmap with randomized hash
// API is similar to set and map in C++ 03,
//   except there is no count function (use member function find instead)
// https://gcc.gnu.org/onlinedocs/libstdc++/ext/pb_ds/gp_hash_table.html
//   and https://gcc.gnu.org/onlinedocs/libstdc++/ext/pb_ds/container_base.html
//   contain most of the public methods
// In practice, has a moderate constant, faster than unordered_set
//   and unordered_map (at the expense of using more memory), but slower than
//   a sorted vector + binary search
// Time Complexity: O(1) on average, O(N) worst case
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/set
//   https://atcoder.jp/contests/agc026/tasks/agc026_c
//   https://judge.yosupo.jp/problem/associative_array
template <class K, class H = rand_hash, class ...Ts>
using hashset = gp_hash_table<K, null_type, H, Ts ...>;
template <class K, class V, class H = rand_hash, class ...Ts>
using hashmap = gp_hash_table<K, V, H, Ts ...>;

// Ordered treeset and treemap
// API is similar to set and map in C++ 03
// There are additional functions to find the kth element with find_by_order(k)
//   and find the 0-indexed rank of a key with order_of_key(key)
// https://gcc.gnu.org/onlinedocs/libstdc++/ext/pb_ds/tree.html
//   https://gcc.gnu.org/onlinedocs/libstdc++/ext/pb_ds/basic_tree.html
//   and https://gcc.gnu.org/onlinedocs/libstdc++/ext/pb_ds/container_base.html
//   contain most of the public methods
// In practice, has a moderate constant, slower that set and map as well as
//   custom written balanced binary search trees
// Time Complexity: O(log N)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/ds4
//   https://codeforces.com/contest/1093/problem/E
//   http://www.usaco.org/index.php?page=viewproblem2&cpid=898
template <class K, class C = less<K>, class ...Ts>
using treeset = tree<K, null_type, C, rb_tree_tag,
                     tree_order_statistics_node_update, Ts ...>;
template <class K, class V, class C = less<K>, class ...Ts>
using treemap = tree<K, V, C, rb_tree_tag,
                     tree_order_statistics_node_update, Ts ...>;

// Priority Queue
// API is similar to std::priority_queue in C++ 03
// There are additional functions to modify an existing key with
//   PQ.modify(iter, key) where ptr is an iterator returned from push,
//   and merge heap B into A with A.join(B)
// https://gcc.gnu.org/onlinedocs/libstdc++/ext/pb_ds/priority_queue.html
//   contains most of the public methods
// There are 5 different types of tags to specify the underlying data structure
//   of the heap, which can be seen at
//   https://gcc.gnu.org/onlinedocs/libstdc++/ext/pb_ds/pq_design.html
// https://gcc.gnu.org/onlinedocs/libstdc++/ext/pb_ds/pq_performance_tests.html
//   contains a summary of the time complexities of each tag
// On average, pairing_heap_tag performs the best out of all tags
// In practice, has a moderate constant, slower than std::priority_queue
//   join is only slightly faster than small to large merging when
//   amortization applies
// Time Complexity (for pairing_heap_tag):
//   push, join: O(1)
//   pop: O(log N)
//   modify, erase: O(log N) amortized, O(N) worst case
// Memory Complexity (for pairing_heap_tag): O(N)
// Tested:
//   https://dmoj.ca/problem/apio16p2
template <class T, class C = less<T>, class Tag = pairing_heap_tag>
using pbdsheap = __gnu_pbds::priority_queue<T, C, Tag>;


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

// Maintains the difference between adjacent elements in a sorted multiset
// Template Arguments:
//   T: the type of the elements being stored
// Fields:
//   vals: a mapping from the values in the multiset to its frequency
//   diffs: a mapping from the difference of each pair of adjacent elements in
//     the multiset to its frequency
// Functions:
//   insert(v): inserts the value v into the multiset
//   erase(v); erases the value of v from the multiset assuming v exists
// In practice, has a moderate constant
// Time Complexity:
//   insert, erase: O(log N)
// Memory Complexity: O(N)
// Tested:
//   https://codeforces.com/contest/1418/problem/D
//   https://dmoj.ca/problem/coci19c3p3
//   https://dmoj.ca/problem/year2018p6
template <class T> struct SetDifferenceMaintenance {
  map<T, int> vals, diffs;
  void insert(T v) {
    auto nxt = vals.lower_bound(v); if (nxt != vals.end() && nxt->first == v) {
      ++nxt->second; ++diffs[T()]; return;
    }
    if (nxt != vals.begin()) {
      auto prv = prev(nxt); diffs[v - prv->first]++; if (nxt != vals.end()) {
        auto it = diffs.find(nxt->first - prv->first);
        if (--it->second == 0) diffs.erase(it);
      }
    }
    if (nxt != vals.end()) ++diffs[nxt->first - v];
    vals.emplace_hint(nxt, v, 1);
  }
  void erase(T v) {
    auto cur = vals.find(v); if (cur->second >= 2) {
      --cur->second; auto it = diffs.find(T());
      if (--it->second == 0) diffs.erase(it);
      return;
    }
    auto nxt = vals.erase(cur); if (nxt != vals.end()) {
      auto it = diffs.find(nxt->first - v);
      if (--it->second == 0) diffs.erase(it);
    }
    if (nxt != vals.begin()) {
      auto prv = prev(nxt); auto it = diffs.find(v - prv->first);
      if (--it->second == 0) diffs.erase(it);
      if (nxt != vals.end()) ++diffs[nxt->first - prv->first];
    }
  }
};


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

// Sparse Table supporting associative and idempotent range queries on
//   a static array
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   Op: a struct with the operation (can also be of type
//       std::function<T(T, T)>); in practice, custom struct is faster than
///      std::function
//     Required Functions:
//       operator (l, r): combines the values l and r, must be
//         associative and idempotent
// Constructor Arguments:
//   A: a vector of type T
//   op: an instance of the Op struct
// Functions:
//   query(l, r): returns the aggregate value of the elements in
//     the range [l, r]
// In practice, the constructor has a small constant,
//   query has a moderate constant, but still faster than segment trees,
//   slightly faster than Fischer Heun Structure, and performs similarly to
//   disjoint sparse tables
// Time Complexity:
//   constructor: O(N log N)
//   query: O(1)
// Memory Complexity: O(N log N)
// Tested:
//   Fuzz and Stress Tested
//   https://dmoj.ca/problem/ncco3d2p1
//   https://www.spoj.com/problems/RMQSQ/
//   https://judge.yosupo.jp/problem/staticrmq
template <class T, class Op> struct SparseTable {
  int N; vector<vector<T>> ST; Op op;
  SparseTable(vector<T> A, Op op = Op())
      : N(A.size()), ST(N == 0 ? 0 : __lg(N) + 1, move(A)), op(op) {
    for (int i = 0; i < int(ST.size()) - 1; i++) {
      ST[i + 1] = ST[0]; for (int j = 0; j < N; j++)
        ST[i + 1][j] = op(ST[i][j], ST[i][min(j + (1 << i), N - 1)]);
    }
  }
  T query(int l, int r) {
    int i = __lg(r - l + 1); return op(ST[i][l], ST[i][r - (1 << i) + 1]);
  }
};


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

// Helper struct to compare two pairs
template <class Cmp> struct PairCmp {
  Cmp cmp; PairCmp(Cmp cmp = Cmp()) : cmp(cmp) {}
  template <class T>
  bool operator () (const pair<T, T> &a, const pair<T, T> &b) const {
    if (cmp(a.first, b.first)) return true;
    if (cmp(b.first, a.first)) return false;
    return cmp(a.second, b.second);
  }
};

// Struct representing no operation
struct NoOp {
  template <class T> void operator () (const T &, const T &) const {}
};

// Adding and removing half-open intervals from a set such that the combined
//   intervals are disjoint
// Template Arguments:
//   T: the type of the endpoints of the intervals
//   Cmp: the comparator to compare two points
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
//   Add: a struct with a callback function (can also be of type
//       std::function<T(T, T)>); in practice, custom struct is faster than
//       std::function
//     Required Functions:
//       operator (l, r): adds the interval [l, r) to the set of
//         disjoint intervals
//   Rem: a struct with a callback function (can also be of type
//       std::function<T(T, T)>); in practice, custom struct is faster than
//       std::function
//     Required Functions:
//       operator (l, r): removes the interval [l, r) from the set of
//         disjoint intervals
// Constructor Arguments:
//   cmp: an instance of Cmp
//   add: an instance of Add
//   rem: an instance of Rem
// Functions:
//   addInterval(L, R): adds an interval [L, R) to the set
//   removeInterval(L, R): removes the interval [L, R) from the set
// In practice, has a moderate constant
// Time Complexity:
//   addInterval, removeInterval: O(log N) amortized
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/art6
//   http://www.usaco.org/index.php?page=viewproblem2&cpid=973
template <class T, class Cmp = less<T>, class Add = NoOp, class Rem = NoOp>
struct IntervalUnion : public set<pair<T, T>, PairCmp<Cmp>> {
  Cmp cmp; Add add; Rem rem;
  IntervalUnion(Cmp cmp = Cmp(), Add add = Add(), Rem rem = Rem())
      : set<pair<T, T>, PairCmp<Cmp>>(PairCmp<Cmp>(cmp)),
        cmp(cmp), add(add), rem(rem) {}
  typename set<pair<T, T>, PairCmp<Cmp>>::iterator addInterval(T L, T R) {
    if (!cmp(L, R) && !cmp(R, L)) return this->end();
    auto it = this->lower_bound(make_pair(L, R)), before = it;
    while (it != this->end() && !cmp(R, it->first)) {
      R = max(R, it->second, cmp); rem(it->first, it->second);
      before = it = this->erase(it);
    }
    if (it != this->begin() && !cmp((--it)->second, L)) {
      L = min(L, it->first, cmp); R = max(R, it->second, cmp);
      rem(it->first, it->second); this->erase(it);
    }
    add(L, R); return this->emplace_hint(before, L, R);
  }
  void removeInterval(T L, T R) {
    if (!cmp(L, R) && !cmp(R, L)) return;
    auto it = addInterval(L, R); auto r2 = it->second;
    if (!cmp(it->first, L) && !cmp(L, it->first)) {
      rem(it->first, it->second); this->erase(it);
    } else (T &) it->second = L;
    if (cmp(R, r2) || cmp(r2, R)) { add(R, r2); this->emplace(R, r2); }
  }
};

// Given a set of intervals (by the PairCmp struct), combine
//   them into disjoint half-open intervals of the form [L, R)
// Range is modified in-place
// Template Arguments:
//   T: the type of the endpoints of the intervals
//   Cmp: the comparator to compare two points
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Function Arguments:
//   A: a reference to a vector of pairs with the first element being the
//     inclusive left bound of the interval and the second element being the
//     exclusive right bound of the interval
//   cmp: an instance of the Cmp struct
// Return Value: a reference to the modified vector
// In practice, has a small constant
// Time Complexity: O(N log N)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/art6
//   https://open.kattis.com/problems/drawingcircles
template <class T, class Cmp = less<T>>
vector<pair<T, T>> &intervalUnion(vector<pair<T, T>> &A, Cmp cmp = Cmp()) {
  sort(A.begin(), A.end(), [&] (const pair<T, T> &a, const pair<T, T> &b) {
    if (cmp(a.first, b.first)) return true;
    if (cmp(b.first, a.first)) return false;
    return cmp(a.second, b.second);
  });
  int i = 0; for (int l = 0, r = 0, N = A.size(); l < N; l = r, i++) {
    A[i] = A[l]; for (r = l + 1; r < N && !cmp(A[i].second, A[r].first); r++)
      A[i].second = max(A[i].second, A[r].second, cmp);
  }
  A.erase(A.begin() + i, A.end()); return A;
}


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

// Fischer Heun Structure in 2 dimensions supporting range maximum queries on
//   a static 2D array
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   Cmp: the comparator to compare two values,
//       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
//   mask_t: the type to store a bitmask, should have around log(N) bits
// Constructor Arguments:
//   A: a 2D vector of elements of type T
//   cmp: an instance of the Cmp struct
// Functions:
//   query(u, d, l, r): returns the maximum element (based on the comparator)
//     in the range [u, d] in the first dimension and [l, r] in
//     the second dimension
// In practice, the constructor has a moderate constant and is significantly
//   faster than sparse table's constructor, query has a moderate constant and
//   is slightly slower than sparse table's query
// Time Complexity:
//   constructor: O(N log (N) (M + M / B log (M / B))), where B is
//     the number of bits in mask_t
//   queryInd, query: O(1) assuming bitshift for mask_t is O(1)
// Memory Complexity: O(N log (N) (M + M / B log (M / B))) assuming mask_t is
//   O(1) memory
// Tested:
//   Fuzz and Stress Tested
//   http://www.usaco.org/index.php?page=viewproblem2&cpid=972
template <class T, class Cmp = less<T>, class mask_t = uint32_t>
struct FischerHeunStructure2D {
  int N, M, lgN;
  vector<vector<FischerHeunStructure<T, Cmp, mask_t>>> FHS; Cmp cmp;
  T cmpVal(const T &a, const T &b) { return cmp(a, b) ? b : a; }
  FischerHeunStructure2D(vector<vector<T>> A, Cmp cmp = Cmp())
      : N(A.size()), M(N == 0 ? 0 : A[0].size()),
        lgN(N == 0 ? 0 : __lg(N) + 1), FHS(lgN), cmp(cmp) {
    for (int i = 0; i < lgN; i++) {
      FHS[i].reserve(N);
      for (int j = 0; j < N; j++) FHS[i].emplace_back(A[j], cmp);
      if (i < lgN - 1) for (int j = 0; j < N; j++) for (int k = 0; k < M; k++)
        A[j][k] = cmpVal(A[j][k], A[min(j + (1 << i), N - 1)][k]);
    }
  }
  T query(int u, int d, int l, int r) {
    int i = __lg(d - u + 1);
    return cmpVal(FHS[i][u].query(l, r), FHS[i][d - (1 << i) + 1].query(l, r));
  }
};


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

// A map that supports insertions of a key value pair (k, v) and
//   queries for the maximum value (based on VCmp) of v for
//   all keys less than or equal to (based on KCmp) k
// Template Arguments:
//   K: the type of the key
//   V: the type of the value
//   KCmp: the comparator for the key, convention is the same as
//       std::map in STL
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
//   VCmp: the comparator for the value, convention is the same as
//       std::priority_queue in STL
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Constructor Arguments:
//   kcmp: an instance of KCmp
//   vcmp: an instance of VCmp
//   NEG_INF: a value for negative infinity
// Functions:
//   add(k, v): adds the key value pair (k, v) and maintains the pareto
//     optimums
//   queryPrefix(k): returns the maximum value (based on VCmp) of v for
//     all keys less than or equal to (based on KCmp) k, or NEG_INF if no such
//     keys exist
//   queryProperPrefix(k): returns the maximum value (based on VCmp) of v for
//     all keys less than or equal to (based on KCmp) k, or NEG_INF if no such
//     keys exist
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(1)
//   add: O(log N) amortized
//   queryPrefix, queryProperPrefix: O(log N)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/dpq
template <class K, class V, class KCmp = less<K>, class VCmp = less<V>>
struct ParetoMap : public map<K, V, KCmp> {
  using M = map<K, V, KCmp>; VCmp vcmp; V NEG_INF;
  ParetoMap(KCmp kcmp = KCmp(), VCmp vcmp = VCmp(),
            V NEG_INF = numeric_limits<V>::lowest())
      : M(kcmp), vcmp(vcmp), NEG_INF(NEG_INF) {}
  void add(K k, V v) {
    auto it = M::lower_bound(k);
    if (it != M::begin() && !VCmp()(prev(it)->second, v)) return;
    while (it != M::end() && !VCmp()(v, it->second)) it = M::erase(it);
    M::emplace_hint(it, k, v);
  }
  V queryPrefix(K k) const {
    auto it = M::upper_bound(k);
    return it == M::begin() ? NEG_INF : prev(it)->second;
  }
  V queryProperPrefix(K k) const {
    auto it = M::lower_bound(k);
    return it == M::begin() ? NEG_INF : prev(it)->second;
  }
};


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

// Wavelet Matrix supporting rank and select operations for a subarray
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of the element of the array
//   Cmp: the comparator to compare two elements
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Constructor Arguments:
//   A: a vector of type T of the values in the array
//   cmp: an instance of the Cmp struct
// Functions:
//   rank(l, r, k): returns the number of elements less than k (using the
//     comparator) in the range [l, r]
//   count(l, r, lo, hi) returns the number of elements not less than lo and
//     not greater than hi (using the comparator) in the range [l, r]
//   select(l, r, k): returns the kth element (0-indexed) sorted by the
//     comparator if the range [l, r] was sorted
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(N log N)
//   rank, count, select: O(log N)
// Memory Complexity: O(N + (N log N) / 64)
// Tested:
//   https://www.spoj.com/problems/KQUERY/ (rank/count)
//   https://www.spoj.com/problems/KQUERYO/ (rank/count)
//   https://codeforces.com/contest/1284/problem/D (rank/count)
//   https://www.spoj.com/problems/MKTHNUM/ (select)
//   https://judge.yosupo.jp/problem/range_kth_smallest (select)
template <class T, class Cmp = less<T>> struct WaveletMatrix {
#define clt [&] (const T &a, const T &b) { return cmp(a, b); }
#define cle [&] (const T &a, const T &b) { return !cmp(b, a); }
  int N, H; Cmp cmp; vector<int> mid; vector<BitPrefixSumArray> B; vector<T> S;
  WaveletMatrix(vector<T> A, Cmp cmp = Cmp())
      : N(A.size()), H(N == 0 ? 0 : __lg(N) + 1), cmp(cmp), mid(H),
        B(H, BitPrefixSumArray(N)), S(move(A)) {
    vector<T> temp = S; sort(S.begin(), S.end(), cmp); vector<int> C(N);
    for (int i = 0; i < N; i++)
      C[i] = lower_bound(S.begin(), S.end(), temp[i], cmp) - S.begin();
    for (int h = H - 1; h >= 0; h--) {
      int ph = 1 << h; for (int i = 0; i < N; i++) B[h].set(i, C[i] <= ph - 1);
      mid[h] = stable_partition(C.begin(), C.end(), [&] (int v) {
                                  return v <= ph - 1;
                                }) - C.begin();
      B[h].build(); for (int i = mid[h]; i < N; i++) C[i] -= ph;
    }
  }
  template <class F> int cnt(int l, int r, const T &v, F f) {
    int ret = 0; for (int cur = 0, h = H - 1; h >= 0; h--) {
      int ph = 1 << h, ql = B[h].query(l - 1), qr = B[h].query(r);
      if (cur + ph - 1 >= N || f(v, S[cur + ph - 1])) { l = ql; r = qr - 1; }
      else { cur += ph; ret += qr - ql; l += mid[h] - ql; r += mid[h] - qr; }
    }
    return ret;
  }
  int rank(int l, int r, const T &v) { return cnt(l, r, v, cle); }
  int count(int l, int r, const T &lo, const T &hi) {
    return cnt(l, r, hi, clt) - cnt(l, r, lo, cle);
  }
  T select(int l, int r, int k) {
    int cur = 0; for (int h = H - 1; h >= 0; h--) {
      int ql = B[h].query(l - 1), qr = B[h].query(r);
      if (k < qr - ql) { l = ql; r = qr - 1; }
      else { cur += 1 << h; k -= qr - ql; l += mid[h] - ql; r += mid[h] - qr; }
    }
    return S[cur];
  }
#undef clt
#undef cle
};


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

// 2D Sparse Table supporting associative and idempotent range queries on
//   a static 2D array
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   Op: a struct with the operation (can also be of type
//       std::function<T(T, T)>); in practice, custom struct is faster than
///      std::function
//     Required Functions:
//       operator (l, r): combines the values l and r, must be
//         associative and idempotent
// Constructor Arguments:
//   A: a 2D vector of elements of type T
//   op: an instance of the Op struct
// Functions:
//   query(u, d, l, r): returns the aggregate value of the elements in
//     the range [u, d] in the first dimension and [l, r] in
//     the second dimension
// In practice, the constructor has a small constant,
//   query has a moderate constant, but still faster than segment trees
// Time Complexity:
//   constructor: O(NM log N log M)
//   query: O(1)
// Memory Complexity: O(NM log N log M)
// Tested:
//   Fuzz and Stress Tested
//   https://dmoj.ca/problem/2drmq
//   http://www.usaco.org/index.php?page=viewproblem2&cpid=972
template <class T, class Op> struct SparseTable2D {
  int N, M, lgN, lgM; vector<vector<vector<vector<T>>>> ST; Op op;
  SparseTable2D(vector<vector<T>> A, Op op = Op())
      : N(A.size()), M(N == 0 ? 0 : A[0].size()),
        lgN(N == 0 ? 0 : __lg(N) + 1), lgM(M == 0 ? 0 : __lg(M) + 1),
        ST(lgN, vector<vector<vector<T>>>(lgM, move(A))), op(op) {
    for (int ir = 0; ir < lgN; ir++) {
      for (int ic = 0; ic < lgM - 1; ic++)
        for (int jr = 0; jr < N; jr++) for (int jc = 0; jc < M; jc++)
          ST[ir][ic + 1][jr][jc] = op(
              ST[ir][ic][jr][jc], ST[ir][ic][jr][min(jc + (1 << ic), M - 1)]);
      if (ir < lgN - 1) for (int jr = 0; jr < N; jr++)
        for (int jc = 0; jc < M; jc++) ST[ir + 1][0][jr][jc] = op(
            ST[ir][0][jr][jc], ST[ir][0][min(jr + (1 << ir), N - 1)][jc]);
    }
  }
  T query(int u, int d, int l, int r) {
    int ir = __lg(d - u + 1), ic = __lg(r - l + 1);
    return op(op(ST[ir][ic][u][l],
                 ST[ir][ic][d - (1 << ir) + 1][r - (1 << ic) + 1]),
              op(ST[ir][ic][d - (1 << ir) + 1][l],
                 ST[ir][ic][u][r - (1 << ic) + 1]));
  }
};


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

// Disjoint Sparse Table supporting associative range queries on a static array
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   Op: a struct with the operation (can also be of type
//       std::function<T(T, T)>); in practice, custom struct is faster than
///      std::function
//     Required Functions:
//       operator (l, r): combines the values l and r, must be associative
// Constructor Arguments:
//   A: a vector of type T
//   op: an instance of the Op struct
// Functions:
//   query(l, r): returns the aggregate value of the elements in
//     the range [l, r]
// In practice, the constructor has a small constant,
//   query has a moderate constant, but still faster than segment trees,
//   slightly faster than Fischer Heun Structure, and performs similarly to
//   sparse tables
// Time Complexity:
//   constructor: O(N log N)
//   query: O(1)
// Memory Complexity: O(N log N)
// Tested:
//   Fuzz and Stress Tested
//   https://www.spoj.com/problems/GSS1/
//   https://www.spoj.com/problems/GSS5/
//   https://judge.yosupo.jp/problem/staticrmq
template <class T, class Op> struct DisjointSparseTable {
  int N; vector<vector<T>> ST; Op op;
  DisjointSparseTable(vector<T> A, Op op = Op())
      : N(A.size()), ST(N == 0 ? 0 : __lg(N * 2 - 1) + 1, move(A)), op(op) {
    for (int i = 1; i < int(ST.size()); i++) {
      ST[i] = ST[0]; for (int j = 0, len = 1 << i; j < N; j += len) {
        int k = min(j + len / 2, N) - 1; T agg = ST[i][k] = ST[0][k];
        for (k--; k >= j; k--) agg = ST[i][k] = op(ST[0][k], agg);
        k = j + len / 2; if (k < N) agg = ST[i][k] = ST[0][k];
        for (k++; k < min(j + len, N); k++) agg = ST[i][k] = op(agg, ST[0][k]);
      }
    }
  }
  T query(int l, int r) {
    if (l == r) return ST[0][l];
    int i = __lg(l ^ r) + 1; return op(ST[i][l], ST[i][r]);
  }
};


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

// Prefix Sum Array supporting range sum queries where all values are 0 or 1,
//   and all values are set before all queries
// Indices are 0-indexed and ranges are inclusive
// Constructor Arguments:
//   N: the size of the array
// Functions:
//   set(i, v): sets the index i with the boolean value v
//   get(i): gets the value of index i
//   build(): builds the prefix sum array, must be called before query can
//     be called after a set operation
//   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 moderate constant, performs faster than a regular prefix
//   sum array and uses less memory
// Time Complexity:
//   constructor, set, get, query: O(1)
//   build: O(N / 64)
// Memory Complexity: O(N / 64)
// Tested:
//   Fuzz and Stress Tested
struct BitPrefixSumArray {
  int M; vector<uint64_t> mask; vector<int> pre;
  BitPrefixSumArray(int N) : M((N >> 6) + 1), mask(M, 0), pre(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 j = 0; j < M; j++)
      pre[j + 1] = pre[j] + __builtin_popcountll(mask[j]);
  }
  int query(int r) {
    r++; int j = r >> 6, k = r & 63;
    return pre[j] + __builtin_popcountll(mask[j] & ((uint64_t(1) << k) - 1));
  }
  int query(int l, int r) { return query(r) - query(l - 1); }
};


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

// Supports online queries and FIFO operations on a multiset of elements,
//   where the underlying data structure only supports LIFO operations
// Template Arguments:
//   S: struct to maintain a multiset of elements
//   Required Fields:
//     T: the type of each element
//   Required Functions:
//     constructor(...args): takes any number of arguments (arguments are
//       passed from constructor of FIFOWithLIFO)
//     push(v): pushes the value v into the multiset
//     pop(): pops the most recent element pushed to the multiset
//     query(...args): queries the multiset
//   Sample Struct: supports queries for graph connectivity
//     struct S {
//       using T = pair<int, int>;
//       UnionFindUndo uf;
//       S(int V) : uf(V) {}
//       void push(const T &e) { uf.join(e.first, e.second); }
//       void pop() { uf.undo(); }
//       bool query(int v, int w) { return uf.connected(v, w); }
//     };
// Constructor Arguments:
//   ...args: arguments to be forwarded to the constructor of S
// Functions:
//   push(v): pushes the value v into the multiset
//   pop(): pops the least recent element pushed into the multiset
//   query(...args): inheirited from S, queries the multiset
// In practice, has a small constant
// Time Complexity:
//   push: O(P) where P is the time complexity of S.push
//   pop: O(P log N) amortized where P is the time complexity of
//     S.push and S.pop, with S.push and S.pop being called at most O(log N)
//     times in total for each element, over all operations
//   query: O(T) where T is the time complexity of S.query
// Memory Complexity: O(N)
// Tested:
//   https://codeforces.com/contest/1386/problem/C
template <class S> struct FIFOWithLIFO : public S {
  using T = typename S::T; vector<pair<T, bool>> stk; vector<T> A, B; int cntA;
  template <class ...Args>
  FIFOWithLIFO(Args &&...args) : S(forward<Args>(args)...), cntA(0) {}
  void push(const T &v, bool a = false) { stk.emplace_back(v, a); S::push(v); }
  void pop() {
    if (cntA == 0) {
      reverse(stk.begin(), stk.end()); cntA = stk.size();
      for (int i = 0; i < cntA; i++) S::pop();
      for (auto &&s : stk) { S::push(s.first); s.second = true; }
    }
    for (; !stk.back().second; stk.pop_back(), S::pop())
      B.push_back(stk.back().first);
    int m = cntA & -cntA; for (int i = 0; i < m; i++, stk.pop_back(), S::pop())
      A.push_back(stk.back().first);
    for (; !B.empty(); B.pop_back()) push(B.back());
    for (; !A.empty(); A.pop_back()) push(A.back(), true);
    stk.pop_back(); S::pop(); cntA--;
  }
};


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

const size_t RANDOM = uniform_int_distribution<size_t>(
    0, numeric_limits<size_t>::max())(rng64);

// Randomized hash for types with std::hash defined and for pairs
// Tested:
//   https://dmoj.ca/problem/set
//   https://atcoder.jp/contests/agc026/tasks/agc026_c
//   https://judge.yosupo.jp/problem/associative_array
struct rand_hash {
  static uint64_t splitmix64(uint64_t x) {
    x += 0x9e3779b97f4a7c15; x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
    x = (x ^ (x >> 27)) * 0x94d049bb133111eb; return x ^ (x >> 31);
  }
  template <class T> size_t operator () (const T &x) const {
    return splitmix64(hash<T>()(x) + RANDOM);
  }
  template <class T1, class T2>
  size_t operator () (const pair<T1, T2> &p) const {
    return 31 * operator ()(p.first) + operator ()(p.second);
  }
};

// unordered_set and unordered_map with a randomized hash
// In practice, has a moderate constant, but still slower than
//   pbds::gp_hash_table
// Tested:
//   https://dmoj.ca/problem/set
//   https://atcoder.jp/contests/agc026/tasks/agc026_c
//   https://judge.yosupo.jp/problem/associative_array
template <class K, class H = rand_hash, class ...Ts>
using uset = unordered_set<K, H, Ts ...>;
template <class K, class V, class H = rand_hash, class ...Ts>
using umap = unordered_map<K, V, H, Ts ...>;