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

// Supports adding lines in the form f(x) = mx + b, finding
//   the maximum value of f(x) at any given x, and removing the last line added
// Template Arguments:
//   T: the type of the slope (m) and intercept (b) of the line, as well as
//     the type of the function argument (x), must be able to store m * b
//   Cmp: the comparator to compare two f(x) 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
// Constructor Arguments:
//   cmp: an instance of the Cmp struct
// Functions:
//   addLine(m, b): adds a line in the form f(x) = mx + b to the set of lines,
//     lines must be added in the order of slope sorted by Cmp
//   undo(): removed the last line added
//   getMax(x): finds the maximum x value (based on the comparator) for all
//     inserted lines
//   size(): returns the number of lines in the convex hull
//   reserve(N): reserves space for N lines in the convex hull
// In practice, has a moderate constant
// Time Complexity:
//   addLine: O(log N) if reserve is called beforehand,
//            O(log N) amortized otherwise
//   getMax: O(log N)
//   undo, size: O(1)
//   reserve: O(N)
// Memory Complexity: O(N) for N lines in the convex hull
// Tested:
//   https://oj.uz/problem/view/CEOI09_harbingers
template <class T, class Cmp = less<T>> struct ConvexHullTrickUndo {
  struct Line {
    T m, b; Line(T m = T(), T b = T()) : m(m), b(b) {}
    T eval(T x) const { return m * x + b; }
  };
  vector<pair<int, Line>> history; vector<Line> L; Cmp cmp; int back;
  ConvexHullTrickUndo(Cmp cmp = Cmp()) : cmp(cmp), back(0) {}
  int size() const { return back; }
  void addLine(T m, T b) {
    int i = back; if (i >= 1)
      i = bsearch<LAST>(1, i + 1, [&] (int j) {
        return cmp(m, L[j - 1].m) || cmp(L[j - 1].m, m) || cmp(b, L[j - 1].b);
      });
    if (i >= 2)
      i = bsearch<LAST>(2, i + 1, [&] (int j) {
        T c1 = (L[j - 1].m - m) * (L[j - 2].b - b);
        T c2 = (L[j - 1].b - b) * (L[j - 2].m - m);
        return c1 > c2;
      });
    if (i == int(L.size())) L.emplace_back();
    history.emplace_back(back, L[i]); L[i] = Line(m, b); back = i + 1;
  }
  void undo() { tie(back, L[back - 1]) = history.back(); history.pop_back(); }
  T getMax(T x) const {
    return L[bsearch<FIRST>(0, back - 1, [&] (int i) {
      return cmp(L[i + 1].eval(x), L[i].eval(x));
    })].eval(x);
  }
  void reserve(int N) { L.reserve(N); history.reserve(N); }
};


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

// Functions for a 2D point
// * operator between 2 points is complex number multiplication
// / operator between 2 points is complex number division
#define OP(op, U, a, x, y) pt operator op (U a) const { return pt(x, y); } \
  pt &operator op##= (U a) { return *this = *this op a; }
#define CMP(op, body) bool operator op (pt p) const { return body; }
struct pt {
  T x, y; constexpr pt(T x = 0, T y = 0) : x(x), y(y) {}
  pt operator + () const { return *this; }
  pt operator - () const { return pt(-x, -y); }
  OP(+, pt, p, x + p.x, y + p.y) OP(-, pt, p, x - p.x, y - p.y)
  OP(*, T, a, x * a, y * a) OP(/, T, a, x / a, y / a)
  friend pt operator * (T a, pt p) { return pt(a * p.x, a * p.y); }
  bool operator < (pt p) const { return eq(x, p.x) ? lt(y, p.y) : lt(x, p.x); }
  CMP(<=, !(p < *this)) CMP(>, p < *this) CMP(>=, !(*this < p))
  CMP(==, !(*this < p) && !(p < *this)) CMP(!=, *this < p || p < *this)
  OP(*, pt, p, x * p.x - y * p.y, y * p.x + x * p.y)
  OP(/, pt, p, (x * p.x + y * p.y) / (p.x * p.x + p.y * p.y),
               (y * p.x - x * p.y) / (p.x * p.x + p.y * p.y))
};
#undef OP
#undef CMP
istream &operator >> (istream &stream, pt &p) { return stream >> p.x >> p.y; }
ostream &operator << (ostream &stream, pt p) {
  return stream << p.x << ' ' << p.y;
}
pt conj(pt a) { return pt(a.x, -a.y); }
T dot(pt a, pt b) { return a.x * b.x + a.y * b.y; }
T cross(pt a, pt b) { return a.x * b.y - a.y * b.x; }
T norm(pt a) { return dot(a, a); }
T abs(pt a) { return sqrt(norm(a)); }
T arg(pt a) { return atan2(a.y, a.x); }
pt polar(T r, T theta) { return r * pt(cos(theta), sin(theta)); }
T distSq(pt a, pt b) { return norm(b - a); }
T dist(pt a, pt b) { return abs(b - a); }
T ang(pt a, pt b) { return arg(b - a); }
// sign of ang, area2, ccw: 1 if counterclockwise, 0 if collinear,
//   -1 if clockwise
T ang(pt a, pt b, pt c) {
  a -= b; c -= b; return arg(pt(dot(c, a), cross(c, a)));
}
// twice the signed area of triangle a, b, c
T area2(pt a, pt b, pt c) { return cross(b - a, c - a); }
int ccw(pt a, pt b, pt c) { return sgn(area2(a, b, c)); }
// a rotated theta radians around p
pt rot(pt a, pt p, T theta) { return (a - p) * pt(polar(T(1), theta)) + p; }
// rotated 90 degrees ccw
pt perp(pt a) { return pt(-a.y, a.x); }

repo_name: wesley-a-leung/Resources

mime_type: text/x-c++


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

// Binary search over a range for the first or last occurrence of
//   a return value for a boolean function
// Template Arguments:
//   ISFIRST: boolean of whether or not the first of last occurrence is being
//     searched for
//   T: the type of the range to search over, must be integral
//   F: the type of the function that is being searched over
// Function Argyments:
//   lo: the inclusive lower bound
//   hi: the exclusive upper bound
//   f: the function to search over
// Return Value:
//   If ISFIRST is true:
//     Returns the first value in the range [lo, hi) where f(x) is true
//       if no value in [lo, hi) satisfies f(x), then it returns hi
//       assumes that all values where f(x) is true are greater than all values
//       where f(x) is false
//   IF ISFIRST is false:
//     Returns the last value in the range [lo, hi) where f(x) is true
//       if no value in [lo, hi) satisfies f(x), then it returns lo - 1
//       assumes that all values where f(x) is true are less than all values
//       where f(x) is false
// In practice, has a small constant
// Time Complexity: O(log (hi - lo)) * (cost to compute f(x))
// Memory Complexity: O(1)
// Tested:
//   https://mcpt.ca/problem/lcc18c5s3
//   https://dmoj.ca/problem/apio19p3
//   https://dmoj.ca/problem/pib20p2
const bool FIRST = true, LAST = false;
template <const bool ISFIRST, class T, class F> T bsearch(T lo, T hi, F f) {
  static_assert(is_integral<T>::value, "T must be integral");
  hi--; while (lo <= hi) {
    T mid = lo + (hi - lo) / 2;
    if (f(mid) == ISFIRST) hi = mid - 1;
    else lo = mid + 1;
  }
  return ISFIRST ? lo : hi;
}


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

// Ternary search for the maximum of a function
// Template Arguments:
//   T: the type of the range to search over, must be a floating point type
//   F: the type of the function that is being searched over
//   Cmp: the comparator to compare two f(x) 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
// Function Arguments:
//   lo: the inclusive lower bound
//   hi: the inclusive upper bound
//   f: the function to search over
//   iters: the number of iterations to run the ternary search
//   cmp: an instance of the Cmp class
// Return Value: the value x in the range [lo, hi] such that f(x) is
//   maximum based on the comparator
// In practice, has a small constant
// Time Complexity: O(iters) * (cost to compute f(x))
// Memory Complexity: O(1)
// Tested:
//   https://codeforces.com/contest/578/problem/C
template <class T, class F, class Cmp = less<T>>
T tsearch(T lo, T hi, F f, int iters, Cmp cmp = less<T>()) {
  static_assert(is_floating_point<T>::value,
                "T must be a floating point type");
  for (int it = 0; it < iters; it++) {
    T m1 = lo + (hi - lo) / 3, m2 = hi - (hi - lo) / 3;
    if (cmp(f(m1), f(m2))) lo = m1;
    else hi = m2;
  }
  return lo + (hi - lo) / 2;
}


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

// Binary search a boolean function over a floating point range
// Template Arguments:
//   T: the type of the range to search over, must be a floating point type
//   F: the type of the function that is being searched over
// Function Arguments:
//   lo: the inclusive lower bound
//   hi: the inclusive upper bound
//   f: the function to search over
//   iters: the number of iterations to run the ternary search
// Return Value: the value x in the range [lo, hi] such that f(x - EPS) is
//   false, while f(x + EPS) is true for some small value EPS
// In practice, has a small constant
// Time Complexity: O(iters) * (cost to compute f(x))
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/globexcup19s2
template <class T, class F> T bsearchFloat(T lo, T hi, F f, int iters) {
  static_assert(is_floating_point<T>::value,
                "T must be a floating point type");
  for (int it = 0; it < iters; it++) {
    T mid = lo + (hi - lo) / 2;
    if (f(mid)) hi = mid;
    else lo = mid;
  }
  return lo + (hi - lo) / 2;
}


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

// Computes a longest increasing subsequence in an array
// Template Arguments:
//   T: the type of each element in the array
// Function Arguments:
//   A: a vector of type T
// Return Value: a vector of indices of the element in the longest increasing
//   subsequence of the array, guaranteed that the values of those indices form
//   the lexicographically least subsequence
// In practice, has a small constant, slower than
//   LongestIncreasingSubsequence.h
// Time Complexity: O(N log N)
// Memory Complexity: O(N)
// Tested:
//   Fuzz and Stress Tested
//   https://dmoj.ca/problem/lis
//   https://open.kattis.com/problems/longincsubseq
template <class T>
vector<int> longestIncreasingSubsequenceFenwick(const vector<T> &A) {
  vector<T> tmp = A; int N = A.size(), mx = 0; sort(tmp.begin(), tmp.end());
  tmp.erase(unique(tmp.begin(), tmp.end()), tmp.end());
  vector<int> ret(N), ind(N + 1, -1); for (int i = 0; i < N; i++)
    ret[i] = lower_bound(tmp.begin(), tmp.end(), A[i]) - tmp.begin();
  struct Max { int operator () (int a, int b) { return max(a, b); } };
  FenwickTree1D<int, Max> FT(N, 0); for (int i = 0; i < N; i++) {
    int len = FT.query(ret[i] - 1) + 1; if (ind[len] == -1) ind[mx = len] = i;
    else if (!(A[ind[len]] < A[i])) ind[len] = i;
    FT.update(ret[i], len); ret[i] = ind[len - 1];
  }
  for (int i = N, k = ind[mx], p; k != -1; k = p) { p = ret[k]; ret[--i] = k; }
  ret.erase(ret.begin(), ret.end() - mx); return ret;
}


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

// Sum over subsets (or supersets) (similar to FST)
// Transforms an array a into an array a' such that a'[i] is equal to the
//   aggregate of all a[j] such that i | j == i if TYPE == SUBSET
//   or i & j == i if TYPE == SUPERSET
// Template Arguments:
//   TYPE: SUBSET or SUPERSET
//   T: the type of each element
// Function Arguments:
//   a: a reference to the vector of type T to transform
//   op: the op being performed (this can also be an inverse operation)
// Time Complexity: O(N log N) where N = size(a)
// Memory Complexity: O(1)
// Tested:
//   https://judge.yosupo.jp/problem/bitwise_and_convolution
//   https://csacademy.com/contest/round-53/task/maxor/
const bool SUBSET = true, SUPERSET = false;
template <const bool TYPE, class T, class Op> void sos(vector<T> &a, Op op) {
  int N = a.size(); assert(!(N & (N - 1)));
  for (int m = 1; m < N; m <<= 1) for (int mask = 0; mask < N; mask++)
    if (bool(mask & m) == TYPE) a[mask] = op(a[mask], a[mask ^ m]);
}


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

// Solves the 0-1 knapsack problem (each item can appear either 0 or 1 times)
// Template Arguments:
//   V: the value type
// Function Arguments:
//   A: a vetor of pairs, with the first element of each pair being the weight
//     of type int, and the second being the value of type V
//   M: type maximum weight the knapsack can hold
//   NEG_INF: a value of type V for negative infinity
// Return Value: a vector dp of size M + 1 with dp[i] being the maximum value
//   that a knapsack with weights summing to exactly i has (or NEG_INF if
//   that sum is not possible)
// In practice, has a very small constant
// Time Complexity: O(NM)
// Memory Complexity: O(M)
// Tested:
//   https://atcoder.jp/contests/dp/tasks/dp_d
template <class V>
vector<V> zeroOneKnapsack(const vector<pair<int, V>> &A, int M,
                          V NEG_INF = numeric_limits<V>::lowest()) {
  vector<V> dp(M + 1, NEG_INF); dp[0] = V(); for (auto &&a : A)
    for (int j = M; j >= a.first; j--) if (dp[j - a.first] > NEG_INF)
      dp[j] = max(dp[j], dp[j - a.first] + a.second);
  return dp;
}

// Solves the dual of the 0-1 knapsack problem (each item can appear
//   either 0 or 1 times)
// Template Arguments:
//   W: the weight type
// Function Arguments:
//   A: a vetor of pairs, with the first element of each pair being the weight
//     of type W, and the second being the value of type int
//   K: type maximum value the knapsack can hold
//   INF: a value of type W for infinity
// Return Value: a vector dp of size K + 1 with dp[i] being the minimum weight
//   that a knapsack with values summing to exactly i has (or INF if that sum
//   is not possible)
// In practice, has a very small constant
// Time Complexity: O(NK)
// Memory Complexity: O(K)
// Tested:
//   https://atcoder.jp/contests/dp/tasks/dp_e
template <class W> vector<W> zeroOneKnapsackDual(
    const vector<pair<W, int>> &A, int K, W INF = numeric_limits<W>::max()) {
  vector<W> dp(K + 1, INF); dp[0] = W(); for (auto &&a : A)
    for (int j = K; j >= a.second; j--) if (dp[j - a.second] < INF)
      dp[j] = min(dp[j], dp[j - a.second] + a.first);
  return dp;
}


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

// Solves the multisubset sum problem with the occurrence of each sum (elements
//   can be chosen multiple times)
// Template Arguments:
//   U: the type of the count
//   T: the type of each element in the array
// Function Arguments:
//   A: a vector of type T
//   M: the maximum sum to count the number of multisubsets
// Return Value: a vector dp of type U with length M + 1 with dp[i] equal to
//   the number of multisubsets that have a sum of i
// In practice, has a very small constant
// Time Complexity: O(NM)
// Memory Complexity: O(M)
// Tested:
//   https://cses.fi/problemset/task/1636
template <class U, class T>
vector<U> multisubsetSumCount(const vector<T> &A, int M) {
  vector<U> dp(M + 1, U()); dp[0] = U(1); for (auto &&a : A)
    for (int j = a; j <= M; j++) dp[j] += dp[j - a];
  return dp;
}

// Solves the multisubset sum problem (elements can be chosen multiple times)
// Template Arguments:
//   U: the type representing the boolean
//   T: the type of each element in the array
// Function Arguments:
//   A: a vector of type T
//   M: the maximum sum to determine if a subset exists
// Return Value: a boolean vector dp with length M + 1 with dp[i] equal to
//   1 if a multisubset has a sum of i, 0 otherwise
// In practice, has a very small constant, faster if U is int, but less memory
//   is used if U is bool (as long as the compiler optimizes vector<booL>);
//   char gives mixed results depending on whether vector<char> is compiled
//   as a string (better if it is not a string)
// Time Complexity: O(NM)
// Memory Complexity: O(M / 64)
// Tested:
//   Stress Tested
//   https://mcpt.ca/problem/acoinproblem
template <class U = bool, class T>
vector<U> multisubsetSum(const vector<T> &A, int M) {
  vector<U> dp(M + 1, U(0)); dp[0] = U(1); for (auto &&a : A)
    for (int j = a; j <= M; j++) dp[j] = dp[j] || dp[j - a];
  return dp;
}


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

// Solves the subset sum problem with the occurrence of each sum (elements
//   can be chosen 0 or 1 times)
// Template Arguments:
//   U: the type of the count
//   T: the type of each element in the array
// Function Arguments:
//   A: a vector of type T
//   M: the maximum sum to count the number of subsets
// Return Value: a vector dp of type U with length M + 1 with dp[i] equal to
//   the number of subsets that have a sum of i
// In practice, has a very small constant
// Time Complexity: O(NM)
// Memory Complexity: O(M)
// Tested:
//   Fuzz Tested
template <class U, class T>
vector<U> subsetSumCount(const vector<T> &A, int M) {
  vector<U> dp(M + 1, U()); dp[0] = U(1); for (auto &&a : A)
    for (int j = M; j >= a; j--) dp[j] += dp[j - a];
  return dp;
}

// Solves the subset sum problem (elements can be chosen 0 or 1 times)
// Template Arguments:
//   M: the maximum sum to determine if a subset exists
//   T: the type of each element in the array
// Function Arguments:
//   A: a vector of type T
//   en: an iterator pointing to after the last element in the array
// Return Value: a bitset dp with length M + 1 with dp[i] equal to
//   whether a subset has a sum of i
// In practice, has a very small constant
// Time Complexity: O(NM / 64)
// Memory Complexity: O(M / 64)
// Tested:
//   Fuzz Tested
//   https://cses.fi/problemset/task/1745/
template <const int M, class T> bitset<M + 1> subsetSum(const vector<T> &A) {
  bitset<M + 1> dp; dp[0] = true; for (auto &&a : A) dp |= dp << a;
  return dp;
}


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

// Computes the maximum sum of a non empty, non consecutive subsequence
// Template Arguments:
//   T: the type of each element in the array
// Function Arguments:
//   A: a vector of type T
//   NEG_INF: a value for negative infinity
// Return Value: the maximum sum of a non empty, non consecutive subsequence
// In practice, has a small constant
// Time Complexity: O(N)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/seq0
template <class T> T maxNonConsecutiveSum(
    const vector<T> &A, T NEG_INF = numeric_limits<T>::lowest()) {
  T incl = NEG_INF, excl = NEG_INF; for (auto &&a : A) {
    T v = a; if (excl > NEG_INF) v = max(v, v + excl);
    excl = max(incl, excl); incl = v;
  }
  return max(incl, excl);
}


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

// Solves the bounded knapsack problem (each item can appear up to a specified
//   number of times)
// Template Arguments:
//   V: the value type
// Function Arguments:
//   A: a vector of tuples, with the first element of each tuple being the
//     weight of type int, the second being the value of type V, and the third
//     being the frequency of type int
//   M: type maximum weight the knapsack can hold
//   NEG_INF: a value of type V for negative infinity
// Return Value: a vector dp of size M + 1 with dp[i] being the maximum value
//   that a knapsack with weights summing to exactly i has (or NEG_INF if
//   that sum is not possible)
// In practice, has a small constant
// Time Complexity: O(NM)
// Memory Complexity: O(M)
// Tested:
//   https://dmoj.ca/problem/knapsack
template <class V>
vector<V> boundedKnapsack(const vector<tuple<int, V, int>> &A, int M,
                          V NEG_INF = numeric_limits<V>::lowest()) {
  vector<V> dp(M + 1, NEG_INF), q(M + 1, V()), dq(M + 1, V()); dp[0] = V();
  for (auto &&a : A) {
    int w = get<0>(a), f = get<2>(a); V v = get<1>(a);
    if (w <= M) for (int s = 0; s < w; s++) {
      V alpha = V(); int ql = 0, qr = 0, dql = 0, dqr = 0;
      for (int j = s; j <= M; j += w) {
        alpha += v; V a = dp[j] <= NEG_INF ? NEG_INF : dp[j] - alpha;
        while (dql < dqr && dq[dqr - 1] < a) dqr--;
        q[qr++] = dq[dqr++] = a;
        while (qr - ql > f + 1) if (q[ql++] == dq[dql]) dql++;
        dp[j] = dq[dql] <= NEG_INF ? NEG_INF : dq[dql] + alpha;
      }
    }
  }
  return dp;
}

// Solves the dual of the bounded knapsack problem (each item can appear up to
//   a specified number of times)
// Template Arguments:
//   W: the weight type
// Function Arguments:
//   A: a vector of tuples, with the first element of each tuple being the
//     weight of type W, the second being the value of type int, and the third
//     being the frequency of type int
//   K: type maximum value the knapsack can hold
//   INF: a value of type W for infinity
// Return Value: a vector dp of size K + 1 with dp[i] being the minimum weight
//   that a knapsack with values summing to exactly i has (or INF if that sum
//   is not possible)
// In practice, has a small constant
// Time Complexity: O(NK)
// Memory Complexity: O(K)
// Tested:
//   https://dmoj.ca/problem/knapsack (Subtask 1)
template <class W>
vector<W> boundedKnapsackDual(const vector<tuple<W, int, int>> &A, int K,
                              W INF = numeric_limits<W>::max()) {
  vector<W> dp(K + 1, INF), q(K + 1, W()), dq(K + 1, W()); dp[0] = W();
  for (auto &&a : A) {
    int v = get<1>(a), f = get<2>(a); W w = get<0>(a);
    if (v <= K) for (int s = 0; s < v; s++) {
      W alpha = W(); int ql = 0, qr = 0, dql = 0, dqr = 0;
      for (int j = s; j <= K; j += v) {
        alpha += w; W a = dp[j] >= INF ? INF : dp[j] - alpha;
        while (dql < dqr && dq[dqr - 1] > a) dqr--;
        q[qr++] = dq[dqr++] = a;
        while (qr - ql > f + 1) if (q[ql++] == dq[dql]) dql++;
        dp[j] = dq[dql] >= INF ? INF : dq[dql] + alpha;
      }
    }
  }
  return dp;
}


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

// Computes the longest common subsequence between 2 arrays
// Template Arguments:
//   T: the type of the elements in the arrays
// Function Arguments:
//   A: the first array
//   B: the second array
// Return Value: the length of the longest common subsequence between the
//   2 arrays
// In practice, has a small constant
// Time Complexity: O(NM) where N = en1 - st1 and M = en2 - st2
// Memory Complexity: O(M)
// Tested:
//   https://dmoj.ca/problem/lcs
template <class T>
int longestCommonSubsequence(const vector<T> &A, const vector<T> &B) {
  int N = A.size(), M = B.size();
  vector<vector<int>> dp(2, vector<int>(M + 1, 0));
  for (int i = 1; i <= N; i++) {
    int cur = i % 2, prv = 1 - cur; for (int j = 1; j <= M; j++)
      dp[cur][j] = A[i - 1] == B[j - 1] ? dp[prv][j - 1] + 1
                                        : max(dp[cur][j - 1], dp[prv][j]);
  }
  return dp[N % 2][M];
}


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

// Knuth's Dynamic Programming Optimization
// Must satisfy dp[l][r] = max(dp[l][m] + dp[m][r] + cost(l, m, r))
//   for l <= m <= r and
//   max(l + 1, opt[l][r - 1]) <= opt[l][r] <= min(opt[l + 1][r], r - 1),
//   where opt[l][r] is the optimal value of m for dp[l][r]
// Template Arguments:
//   T: the type of the return value of the cost function
//   F: the type of the function used to compute the cost function
//   Cmp: the comparator to compare two T 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
// Function Arguments:
//   N: the size of the array, must be positive
//   f(l, m, r): computes the cost function for the range [l, r) with the
//     midpoint m
//   cmp: an instance of the Cmp struct
// Return Value: a dp array of size N x (N + 1) with dp[l][r] equal to the
//   maximum dp value (based on the comparator) of the range [l, r)
// In practice, has a small constant
// Time Complexity: O(N^2)
// Memory Complexity: O(N^2)
// Tested:
//   https://atcoder.jp/contests/dp/tasks/dp_n
template <class T, class F, class Cmp = less<T>>
vector<vector<T>> knuth(int N, F f, Cmp cmp = Cmp()) {
  assert(N > 0); vector<vector<T>> dp(N, vector<T>(N + 1, T()));
  vector<vector<int>> opt(N, vector<int>(N + 1));
  for (int l = N - 1; l >= 0; l--) for (int r = l; r <= N; r++) {
    if (r - l <= 1) { opt[l][r] = l; continue; }
    int st = max(l + 1, opt[l][r - 1]), en = min(opt[l + 1][r], r - 1);
    for (int m = st; m <= en; m++) {
      T cost = dp[l][m] + dp[m][r] + f(l, m, r);
      if (m == st || cmp(dp[l][r], cost)) { dp[l][r] = cost; opt[l][r] = m; }
    }
  }
  return dp;
}


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

// Solves the unbounded knapsack problem (each item can appear any number
//   of times)
// Template Arguments:
//   V: the value type
// Function Arguments:
//   A: a vetor of pairs, with the first element of each pair being the weight
//     of type int, and the second being the value of type V
//   M: type maximum weight the knapsack can hold
//   NEG_INF: a value of type V for negative infinity
// Return Value: a vector dp of size M + 1 with dp[i] being the maximum value
//   that a knapsack with weights summing to exactly i has (or NEG_INF if
//   that sum is not possible)
// In practice, has a very small constant
// Time Complexity: O(NM)
// Memory Complexity: O(M)
// Tested:
//   https://www.spoj.com/problems/DBALLZ/
template <class V>
vector<V> unboundedKnapsack(const vector<pair<int, V>> &A, int M,
                            V NEG_INF = numeric_limits<V>::lowest()) {
  vector<V> dp(M + 1, NEG_INF); dp[0] = V(); for (auto &&a : A)
    for (int j = a.first; j <= M; j++) if (dp[j - a.first] > NEG_INF)
      dp[j] = max(dp[j], dp[j - a.first] + a.second);
  return dp;
}

// Solves the dual of the unbounded knapsack problem (each item can appear
//   any number of times)
// Template Arguments:
//   W: the weight type
// Function Arguments:
//   A: a vetor of pairs, with the first element of each pair being the weight
//     of type W, and the second being the value of type int
//   K: type maximum value the knapsack can hold
//   INF: a value of type W for infinity
// Return Value: a vector dp of size K + 1 with dp[i] being the minimum weight
//   that a knapsack with values summing to exactly i has (or INF if that sum
//   is not possible)
// In practice, has a very small constant
// Time Complexity: O(NK)
// Memory Complexity: O(K)
// Tested:
//   https://www.spoj.com/problems/PIGBANK/
template <class W> vector<W> unboundedKnapsackDual(
    const vector<pair<W, int>> &A, int K, W INF = numeric_limits<W>::max()) {
  vector<W> dp(K + 1, INF); dp[0] = W(); for (auto &&a : A)
  for (int j = a.second; j <= K; j++) if (dp[j - a.second] < INF)
    dp[j] = min(dp[j], dp[j - a.second] + a.first);
  return dp;
}


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

// Computes a longest increasing subsequence in an array
// Template Arguments:
//   T: the type of each element in the array
// Function Arguments:
//   A: the array
// Return Value: a vector of indices of the element in a longest increasing
//   subsequence of the array, the indices is guaranteed to be the
//   lexicographically largest vector of indices
// In practice, has a very small constant, faster than the Fenwick Tree version
// Time Complexity: O(N log N)
// Memory Complexity: O(N)
// Tested:
//   Fuzz and Stress Tested
//   https://dmoj.ca/problem/lis
//   https://open.kattis.com/problems/longincsubseq
template <class T>
vector<int> longestIncreasingSubsequence(const vector<T> &A) {
  int N = A.size(), k = -1; vector<int> ret(N);
  vector<pair<T, int>> dp; dp.reserve(N); for (int i = 0; i < N; i++) {
    int j = lower_bound(dp.begin(), dp.end(), make_pair(A[i], -1))
        - dp.begin();
    ret[i] = j == 0 ? -1 : dp[j - 1].second;
    if (j >= int(dp.size())) dp.emplace_back(A[i], i);
    else dp[j] = make_pair(A[i], i);
    if (j == int(dp.size()) - 1) k = i;
  }
  for (int i = N, p; k != -1; k = p) { p = ret[k]; ret[--i] = k; }
  ret.erase(ret.begin(), ret.end() - dp.size()); return ret;
}


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

// Divide and Conquer Dynamic Programming Optimization
// Must satisfy dp[k][i] = max(dp[k - 1][j] + cost(k, j, i))
//   for 0 <= j <= i <= N and opt[k][i] <= opt[k][i + 1] where opt[k][i] is
//   the optimal value of j for dp[k][i]
// Template Arguments:
//   T: the type of the return value of the cost function
//   F: the type of the function used to compute the cost function
//   Cmp: the comparator to compare two T 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
// Function Arguments:
//   K: the number of iterations, must be positive
//   N: the size of the array, must be positive
//   f(k, j, i): computes the cost function for the kth iteration (1-indexed)
//     for the transition from j to i (0 <= j <= i <= N)
//   cmp: an instance of the Cmp struct
// Return Value: a vector of length K + 1 with the value of dp[k][N]
//   for 1 <= k <= N with dp[0][N] = 0
// In practice, has a small constant
// Time Complexity: O(KN log N)
// Memory Complexity: O(N)
// Tested:
//   https://codeforces.com/contest/321/problem/E
//   https://www.acmicpc.net/problem/7052
template <class T, class F, class Cmp = less<T>>
vector<T> divideAndConquer(int K, int N, F f, Cmp cmp = Cmp()) {
  assert(N > 0 && K > 0); vector<vector<T>> dp(2, vector<T>(N + 1, T()));
  for (int i = 0; i <= N; i++) dp[1][i] = f(1, 0, i);
  vector<T> ret(K + 1, T()); ret[1] = dp[1][N]; vector<int> mid(N + 1);
  vector<pair<int, int>> bnd(N + 1, make_pair(0, N)), opt = bnd;
  for (int i = 0, b = 1; i <= N; i++) {
    int l, r; tie(l, r) = bnd[i]; int m = mid[i] = l + (r - l) / 2;
    if (l <= m - 1) bnd[b++] = make_pair(l, m - 1);
    if (m + 1 <= r) bnd[b++] = make_pair(m + 1, r);
  }
  for (int k = 2; k <= K; k++) {
    int cur = k % 2, prv = 1 - cur; for (int i = 0, b = 1; i <= N; i++) {
      int l, r, m = mid[i], optl, optr;
      tie(l, r) = bnd[i]; tie(optl, optr) = opt[i]; int optm = optl;
      for (int j = optl, en = min(m, optr); j <= en; j++) {
        T cost = dp[prv][j] + f(k, j, m);
        if (j == optl || cmp(dp[cur][m], cost)) {
          dp[cur][m] = cost; optm = j;
        }
      }
      if (l <= m - 1) opt[b++] = make_pair(optl, optm);
      if (m + 1 <= r) opt[b++] = make_pair(optm, optr);
    }
    ret[k] = dp[cur][N];
  }
  return ret;
}


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

// Computes the maximum sum of a non empty subarray
// Template Arguments:
//   T: the type of each element in the array
// Function Arguments:
//   A: a vector of type T
//   NEG_INF: a value for negative infinity
// Return Value: a tuple with the maximum sum of a non empty subarray, and the
//   inclusive left and right indices of the subarray
// In practice, has a small constant
// Time Complexity: O(N)
// Memory Complexity: O(1)
// Tested:
//   https://mcpt.ca/problem/subarraymaximization
template <class T> tuple<T, int, int> maxSubarraySum(
    const vector<T> &A, T NEG_INF = numeric_limits<T>::lowest()) {
  tuple<T, int, int> ret = make_tuple(NEG_INF, 0, -1); T sm = T();
  for (int i = 0, l = 0; i < int(A.size()); i++) {
    sm += A[i]; if (get<0>(ret) < sm) ret = make_tuple(sm, l, i);
    if (sm < T()) { sm = T(); l = i + 1; }
  }
  return ret;
}


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

// Computes the minimum edit distance between 2 arrays, allowing for custom
//   penalties for copying, replacement, insertion, and deletion
// Template Arguments:
//   T: the type of each element in the arrays
//   U: the type of the penalties
// Function Arguments:
//   A: the first array of type T
//   B: the second array of type T
// Return Value: the minimum cost to edit the first array to become the
//   second array
// In practice, has a small constant
// Time Complexity: O(NM)
// Memory Complexity: O(M)
// Tested:
//   https://www.spoj.com/problems/EDIST/
template <class T, class U = int>
U minEditDistance(const vector<T> &A, const vector<T> &B, U cpyPen = U(),
    U repPen = U(1), U insPen = U(1), U delPen = U(1)) {
  int N = A.size(), M = B.size();
  vector<vector<U>> dp(2, vector<U>(M + 1, U()));
  for (int i = 1; i <= M; i++) dp[0][i] = dp[0][i - 1] + insPen;
  for (int i = 1; i <= N; i++) {
    int cur = i % 2, prv = 1 - cur; dp[cur][0] = dp[prv][0] + delPen;
    for (int j = 1; j <= M; j++) {
      dp[cur][j] = min(min(dp[prv][j - 1] + repPen, dp[cur][j - 1] + insPen),
                       dp[prv][j] + delPen);
      if (A[i - 1] == B[j - 1])
        dp[cur][j] = min(dp[cur][j], dp[prv][j - 1] + cpyPen);
    }
  }
  return dp[N % 2][M];
}


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

// Counts inversions and sorts an array using a Fenwick Tree with a comparator
// Template Arguments:
//   It: the type of the iterator
//   Cmp: the comparator to compare two values
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Function Arguments:
//   st: an iterator pointing to the first element in the array
//   en: an iterator pointing to after the last element in the array
//   cmp: an instance of the Cmp struct
// Return Value: the number of inversions in the array
// In practice, has a moderate constant, slower than the merge sort version
// Time Complexity: O(N log N)
// Memory Complexity: O(N)
// Tested:
//   https://www.spoj.com/problems/INVCNT/
//   https://codeforces.com/problemsets/acmsguru/problem/99999/180
template <class It, class Cmp = less<typename iterator_traits<It>::value_type>>
long long countInversionsFenwick(It st, It en, Cmp cmp = Cmp()) {
  vector<typename iterator_traits<It>::value_type> tmp(st, en);
  int N = en - st; long long ret = 0; sort(tmp.begin(), tmp.end(), cmp);
  tmp.erase(unique(tmp.begin(), tmp.end()), tmp.end());
  FenwickTree1D<int> FT(tmp.size()); for (int i = N - 1; i >= 0; i--) {
    int c = lower_bound(tmp.begin(), tmp.end(), st[i], cmp) - tmp.begin();
    ret += FT.query(c - 1); FT.update(c, 1);
  }
  return ret;
}


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

// Sorts an array using shell sort with a comparator
// Shell sort is not stable
// Template Arguments:
//   It: the type of the iterator
//   Cmp: the comparator to compare two values
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Function Arguments:
//   st: an iterator pointing to the first element in the array
//   en: an iterator pointing to after the last element in the array
//   cmp: an instance of the Cmp struct
// In practice, has a very small constant, slower than std::sort by
//   a factor of 2 at N = 1e7
// Time Complexity: O(N^(4/3))
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/bf1hard
template <class It, class Cmp = less<typename iterator_traits<It>::value_type>>
void shellSort(It st, It en, Cmp cmp = Cmp()) {
  int n = en - st, h = 1; while (h < n * 4 / 9) h = h * 9 / 4 + 1;
  for (; h >= 1; h = h * 4 / 9) for (int i = h; i < n; i++)
    for (int j = i; j >= h && cmp(st[j], st[j - h]); j -= h)
      swap(st[j], st[j - h]);
}


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

// Helper function for inversion counting
template <class InIt, class OutIt, class Cmp>
long long countInversions(InIt st1, OutIt st2, int N, Cmp cmp, int len) {
  long long ret = 0; for (int lo = 0; lo < N; lo += len + len) {
    int mid = min(lo + len, N), hi = min(mid + len, N);
    for (int i = lo, j = mid, k = lo; k < hi; k++) {
      if (i >= mid) st2[k] = st1[j++];
      else if (j >= hi) { st2[k] = st1[i++]; ret += j - mid; }
      else if (cmp(st1[j], st1[i])) { st2[k] = st1[j++]; ret += mid - i; }
      else { st2[k] = st1[i++]; ret += j - mid; }
    }
  }
  return ret;
}

// Counts inversions and sorts an array using merge sort with a comparator
// Template Arguments:
//   It: the type of the iterator
//   Cmp: the comparator to compare two values
//     Required Functions:
//       operator (a, b): returns true if and only if a compares less than b
// Function Arguments:
//   st: an iterator pointing to the first element in the array
//   en: an iterator pointing to after the last element in the array
//   cmp: an instance of the Cmp struct
// Return Value: the number of inversions in the array
// In practice, has a small constant, faster than the Fenwick Tree version
// Time Complexity: O(N log N)
// Memory Complexity: O(N)
// Tested:
//   https://www.spoj.com/problems/INVCNT/
//   https://codeforces.com/problemsets/acmsguru/problem/99999/180
template <class It, class Cmp = less<typename iterator_traits<It>::value_type>>
long long countInversions(It st, It en, Cmp cmp = Cmp()) {
  vector<typename iterator_traits<It>::value_type> tmp(st, en);
  int N = en - st; auto aux = tmp.begin(); long long ret = 0; bool flag = true;
  for (int len = 1; len < N; len *= 2, flag = !flag)
    ret += flag ? countInversions(aux, st, N, cmp, len)
                : countInversions(st, aux, N, cmp, len);
  if (flag) copy_n(aux, N, st);
  return ret / 2;
}


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

// Mo's algorithm to answer offline ranges queries over
//   a static array A of length N where elements can be only added to a
//   multiset-like object
// Template Arguments:
//   S: struct to maintain a multiset of elements
//   Required Fields:
//     T: the type of each element
//     R: the type of the return value for each query
//     Q: the query object that contains information for each query
//       Required Fields:
//         l: the left endpoint of the query range
//         r: the right endpoint of the query range
//   Required Functions:
//     constructor(...args): takes any number of arguments (arguments are
//       passed from constructor of Mo)
//     add(v): adds the value v to the multiset
//     save(): saves the current state of the multiset
//     rollback(): rollbacks the multiset to the saved state
//     reset(): resets the multiset to its initial state
//     query(q): returns the answer of type R for the query q of type Q to
//       the multiset, which is guaranteed to contain all elements in the range
//       [q.l, q.r] and only those elements
//   Sample Struct: supporting queries for 0-1 knapsack problem with a size of
//       q.w over a subarray of items with a weight and value
//     struct S {
//       using T = pair<int, int>; using R = int;
//       struct Q { int l, r, w; };
//       int MAXW; vector<int> dp, saved;
//       S(int MAXW) : MAXW(MAXW), dp(MAXW + 1) {}
//       void add(const T &v) {
//         for (int j = MAXW; j >= v.first; j--)
//           dp[j] = max(dp[j], dp[j - v.first] + v.second);
//       }
//       void save() { saved = dp; }
//       void rollback() { dp = saved; }
//       void reset() { fill(dp.begin(), dp.end(), 0); }
//       R query(const Q &q) const { return dp[q.w]; }
//     };
// Constructor Arguments:
//   A: a vector of type S::T of the values in the array
//   queries: a vector of type S::Q representing the queries
//   SCALE: the value to scale sqrt by
//   ...args: arguments to pass to the constructor of S
// Fields:
//   ans: a vector of integers with the answer for each query
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(C + K ((log K + A sqrt N) + R + T) + D sqrt N)
//     for K queries where C is the time complexity of S's constructor,
//     A is the time complexity of S.add, T is the time complexity of S.query,
//     R is the time complexity of S.save and S.rollback,
//     and D is the time complexity of S.reset
// Memory Complexity: O(K + M) for K queries where M is the memory complexity
//   of S
// Tested:
//   https://www.spoj.com/problems/DQUERY/
//   https://dmoj.ca/problem/vpex1p5
template <class S> struct MoRollback {
  using T = typename S::T; using R = typename S::R; using Q = typename S::Q;
  struct Query {
    Q q; int i, b; Query(const Q &q, int i, int b) : q(q), i(i), b(b) {}
    pair<int, int> getPair() const { return make_pair(b, q.r); }
    bool operator < (const Query &o) const { return getPair() < o.getPair(); }
  };
  vector<R> ans;
  template <class ...Args> MoRollback(
      const vector<T> &A, const vector<Q> &queries, double SCALE = 1,
      Args &&...args) {
    int N = A.size(), K = queries.size(), bsz = max(1.0, sqrt(N) * SCALE);
    vector<Query> q; q.reserve(K); S s(forward<Args>(args)...);
    int l = 0, r = -1, last = -1; for (int i = 0; i < K; i++)
      q.emplace_back(queries[i], i, queries[i].l / bsz);
    sort(q.begin(), q.end()); for (auto &&qi : q) {
      if (qi.b != last) {
        l = min(N, ((last = qi.b) + 1) * bsz); r = l - 1; s.reset(); s.save();
      }
      if (qi.q.r < l) { for (int i = qi.q.l; i <= qi.q.r; i++) s.add(A[i]); }
      else {
        while (r < qi.q.r) s.add(A[++r]);
        s.save(); for (int i = l - 1; i >= qi.q.l; i--) s.add(A[i]);
      }
      R res = s.query(qi.q); if (ans.empty()) ans.resize(K, res);
      ans[qi.i] = res; s.rollback();
    }
  }
};


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

// Supports online queries for the number of distinct elements in the
//   range [l, r] for a static array A of length N
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
// Constructor Arguments:
//   A: a vector of type T of the values in the array
// Functions:
//   query(l, r): returns the number of distinct values in the range [l, r]
// In practice, has a moderate constant, faster than mo, but slower than the
//   offline version
// Time Complexity:
//   constructor: O(N log N)
//   query: O(log N)
// Memory Complexity: O(N + (N log N) / 64)
// Tested:
//   https://www.acmicpc.net/problem/14898
template <class T> struct CountDistinct {
  WaveletMatrix<int> wm;
  vector<int> init(const vector<T> &A) {
    vector<T> temp = A; sort(temp.begin(), temp.end());
    temp.erase(unique(temp.begin(), temp.end()), temp.end());
    vector<int> last(temp.size(), -1), ret(A.size());
    for (int i = 0; i < int(A.size()); i++) {
      int c = lower_bound(temp.begin(), temp.end(), A[i]) - temp.begin();
      ret[i] = last[c]; last[c] = i;
    }
    return ret;
  }
  CountDistinct(const vector<T> &A) : wm(init(A)) {}
  int query(int l, int r) { return wm.rank(l, r, l); }
};


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

// Mo's algorithm to answer offline ranges queries over
//   a static array A of length N where elements can be added to and removed
//   from a multiset-like object
// Template Arguments:
//   S: struct to maintain a multiset of the elements in the subarray [l, r]
//   Required Fields:
//     T: the type of each element
//     R: the type of the return value for each query
//     Q: the query object that contains information for each query
//       Required Fields:
//         l: the left endpoint of the query range
//         r: the right endpoint of the query range
//   Required Functions:
//     constructor(...args): takes any number of arguments (arguments are
//       passed from constructor of Mo)
//     addLeft(v): adds the value v to the multiset by expanding the subarray
//       [l, r] to [l - 1, r]
//     addRight(v): adds the value v to the multiset by expanding the subarray
//       [l, r] to [l, r + 1]
//     removeLeft(v): removes the value v from the multiset assuming
//       it exists by contracting the subarray [l, r] to [l + 1, r]
//     removeRight(v): removes the value v from the multiset assuming
//       it exists by contracting the subarray [l, r] to [l, r - 1]
//     query(q): returns the answer of type R for the query q of type Q to
//       the multiset, which is guaranteed to contain all elements in the range
//       [q.l, q.r] and only those elements
//   Sample Struct: supporting queries for the number of inversions
//       in an integer array with a small maximum value
//     struct S {
//       using T = int; using R = long long;
//       struct Q { int l, r; };
//       int U; FenwickTree1D<int> ft; long long inversions;
//       S(const vector<T> &A)
//           : U(*max_element(A.begin(), A.end()) + 1), ft(U),
//             inversions(0) {}
//       void addLeft(const T &v) {
//         inversions += ft.query(0, v - 1); ft.update(v, 1);
//       }
//       void addRight(const T &v) {
//         inversions += ft.query(v + 1, U - 1); ft.update(v, 1);
//       }
//       void removeLeft(const T &v) {
//         inversions -= ft.query(0, v - 1); ft.update(v, -1);
//       }
//       void removeRight(const T &v) {
//         inversions -= ft.query(v + 1, U - 1); ft.update(v, -1);
//       }
//       R query(const Q &q) const { return inversions; }
//     };
// Constructor Arguments:
//   A: a vector of type S::T of the values in the array
//   queries: a vector of type S::Q representing the queries
//   SCALE: the value to scale sqrt by
//   ...args: arguments to pass to the constructor of S
// Fields:
//   ans: a vector of integers with the answer for each query
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(C + K ((log K + U sqrt N) + T))
//     for K queries where C is the time complexity of S's constructor,
//     U is the time complexity of S.addLeft, S.addRight,
//     S.removeLeft, S.removeRight,
//     and T is the time complexity of S.query
// Memory Complexity: O(K + M) for K queries where M is the memory complexity
//   of S
// Tested:
//   https://www.spoj.com/problems/DQUERY/
//   https://atcoder.jp/contests/abc174/tasks/abc174_f
//   https://judge.yosupo.jp/problem/static_range_inversions_query
template <class S> struct Mo {
  using T = typename S::T; using R = typename S::R; using Q = typename S::Q;
  struct Query {
    Q q; int i, b; Query(const Q &q, int i, int b) : q(q), i(i), b(b) {}
    pair<int, int> getPair() const { return make_pair(b, b % 2 ? -q.r : q.r); }
    bool operator < (const Query &o) const { return getPair() < o.getPair(); }
  };
  vector<R> ans;
  template <class ...Args> Mo(const vector<T> &A, const vector<Q> &queries,
                              double SCALE = 2, Args &&...args) {
    int K = queries.size(), bsz = max(1.0, sqrt(A.size()) * SCALE);
    vector<Query> q; q.reserve(K); S s(forward<Args>(args)...);
    for (int i = 0; i < K; i++)
      q.emplace_back(queries[i], i, queries[i].l / bsz);
    sort(q.begin(), q.end()); int l = 0, r = l - 1; for (auto &&qi : q) {
      while (l > qi.q.l) s.addLeft(A[--l]);
      while (r < qi.q.r) s.addRight(A[++r]);
      while (l < qi.q.l) s.removeLeft(A[l++]);
      while (r > qi.q.r) s.removeRight(A[r--]);
      R res = s.query(qi.q); if (ans.empty()) ans.resize(K, res);
      ans[qi.i] = res;
    }
  }
};


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

// Uses divide and conquer to answer offline ranges queries over
//   a multiset where elements can be added or removed at anytime, but the
//   underlying data structure is only able to add and delete elements in LIFO
//   order, as well as save and restore its state
// Template Arguments:
//   S: struct to maintain a multiset of elements
//   Required Fields:
//     T: the type of each element
//     R: the type of the return value for each query
//     Q: the query object that contains information for each query
//   Required Functions:
//     constructor(...args): takes any number of arguments (arguments are
//       passed from LIFOSetDivAndConq::solveQueries)
//     push(v): adds the value v to the multiset
//     pop(): pops the most recent element pushed to the multiset
//     query(q): returns the answer of type R for the query q of type Q to
//       the multiset
//   Sample Struct: supports queries for the sum of element in the same
//       connected component as a vertex in a graph, where edges are
//       added and removed
//     struct S {
//       using T = pair<int, int>; using R = long long;
//       struct Q { int v; };
//       WeightedUnionFindUndo<R, plus<R>> uf;
//       S(const vector<R> &W) : uf(W) {}
//       void push(const T &e) { uf.join(e.first, e.second); }
//       void pop() { uf.undo(); }
//       R query(const Q &q) { return uf.getWeight(q.v); }
//     };
// Fields:
//   ans: a vector of integers with the answer for each query
// Functions:
//   addElement(v): adds v to the multiset
//   removeElement(v): removes v from the multiset
//   addQuery(q): adds a query of type S::Q
//   solveQueries(...args): solves all queries asked so far, with
//     ...args being passed to the constructor of S
// In practice, has a small constant
// Time Complexity:
//   addElement, removeElement: O(1) amortized
//   solveQueries: O(C + P E log E + KT)
//     for K queries and E total elements where C is the time complexity
//     of S's constructor, P is the time complexity of S.push and S.pop,
//     T is the time complexity of S.query, with S.push and S.pop being called
//     at most O(log K) times in total for each element, over all operations
// Memory Complexity: O(K + E + M) for K queries and E total elements, where
//   M is the memory complexity of S
// Tested:
//   https://judge.yosupo.jp/problem/dynamic_graph_vertex_add_component_sum
template <class S> struct LIFOSetDivAndConq {
  using T = typename S::T; using R = typename S::R; using Q = typename S::Q;
  vector<T> add, rem; vector<Q> queries; vector<R> ans; vector<char> type;
  void addElement(const T &v) { add.push_back(v); type.push_back(1); }
  void removeElement(const T &v) { rem.push_back(v); type.push_back(-1); }
  void addQuery(const Q &q) { queries.push_back(q); type.push_back(0); }
  template <class ...Args> void solveQueries(Args &&...args) {
    int E = type.size(); ans.clear(); ans.reserve(queries.size());
    vector<T> A = add; sort(A.begin(), A.end());
    vector<pair<int, int>> events; events.reserve(E);
    vector<int> last(A.size(), INT_MAX); S s(forward<Args>(args)...);
    for (int i = 0, addInd = 0, remInd = 0, queryInd = 0; i < E; i++) {
      if (type[i] == 1) {
        int j = lower_bound(A.begin(), A.end(), add[addInd++]) - A.begin();
        events.emplace_back(j, last[j]); last[j] = i;
      } else if (type[i] == -1) {
        int j = lower_bound(A.begin(), A.end(), rem[remInd++]) - A.begin();
        events.emplace_back(j, last[j]); int temp = events[last[j]].second;
        events[last[j]].second = i; last[j] = temp;
      } else events.emplace_back(queryInd++, i);
    }
    function<void(int, int)> dc = [&] (int l, int r) {
      if (l == r) {
        if (events[l].second == l)
          ans.push_back(s.query(queries[events[l].first]));
        return;
      }
      int m = l + (r - l) / 2; int cnt = 0; for (int i = m + 1; i <= r; i++)
        if (events[i].second < l) { s.push(A[events[i].first]); cnt++; }
      dc(l, m); for (int i = 0; i < cnt; i++) s.pop();
      cnt = 0; for (int i = l; i <= m; i++)
        if (events[i].second > r) { s.push(A[events[i].first]); cnt++; }
      dc(m + 1, r); for (int i = 0; i < cnt; i++) s.pop();
    };
    if (E > 0) dc(0, E - 1);
  }
};


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

// Supports online queries for the number of distinct elements in the
//   range [l, r] for an array A of length N where the sequence of all update
//   indices and values are known beforehand
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   C: a container representing a mapping from type T to a set of integers
// Constructor Arguments:
//   A: a vector of type T of the initial values in the array
//   updates: a vector of pairs of type int and T with the first element
//     representing the index being updated, and the second element being the
//     the value it is updated to
//   ...args: arguments to pass to the constructor of an instance of type C
// Functions:
//   processNextUpdate(): processes the next update (starting with the intial
//     array and proceeding in the order given by the updates vector)
//   query(l, r): returns the number of distinct values in the range [l, r] in
//     the current array
// In practice, has a small constant, slgithly faster than the online version
// Time Complexity:
//   constructor: O((N + U) log (N + U)) for U total updates
//   query: O(log(N + U) log((N + U) / 64)) for U total updates
// Memory Complexity: O(N + ((N + U) log (N + U)) / 64) for U total updates
// Tested:
//   https://dmoj.ca/problem/mnyc17p6
template <class T, class C = map<T, set<int>>>
struct CountDistinctOfflineUpdates {
  struct R {
    using Data = int; using Lazy = int;
    static Data qdef() { return 0; }
    static Data merge(const Data &l, const Data &r) { return l + r; }
    static Data applyLazy(const Data &l, const Lazy &r) { return l + r; }
    BitFenwickTree FT;
    R(const vector<Data> &A) : FT(A.size()) {
      for (int i = 0; i < int(A.size()); i++) FT.set(i, A[i]);
      FT.build();
    }
    void update(int i, const Lazy &val) { FT.update(i, val); }
    Data query(int l, int r) { return FT.query(l, r); }
  };
  C M; int N, cur, ucnt; vector<int> en; vector<pair<int, int>> wupds;
  WaveletMatrixAggregationOfflineUpdates<int, R> wm;
  void add(int i, const T &v, vector<pair<int, int>> &wus) {
    set<int> &mv = M[v]; auto it = mv.upper_bound(i);
    int j = it == mv.begin() ? -1 : *prev(it); wus.emplace_back(i, j);
    if (it != mv.end()) wus.emplace_back(*it, i);
    mv.insert(it, i);
  }
  void rem(int i, const T &v, vector<pair<int, int>> &wus) {
    set<int> &mv = M[v]; auto it = mv.erase(mv.find(i));
    int j = it == mv.begin() ? -1 : *prev(it); wus.emplace_back(i, N);
    if (it != mv.end()) wus.emplace_back(*it, j);
  }
  vector<pair<int, int>> init1(vector<T> A,
                              const vector<pair<int, T>> &updates) {
    vector<pair<int, int>> ret;
    for (int i = 0; i < int(A.size()); i++) add(i, A[i], ret);
    cur = ret.size(); for (int i = 0; i < int(updates.size()); i++) {
      int j = updates[i].first; rem(j, A[j], ret);
      add(j, A[j] = updates[i].second, ret); en[i] = ret.size();
    }
    return ret;
  }
  WaveletMatrixAggregationOfflineUpdates<int, R> init2(const vector<T> &A) {
    int N = A.size(); vector<vector<int>> vals(N);
    for (auto &&wu : wupds) vals[wu.first].push_back(wu.second);
    return WaveletMatrixAggregationOfflineUpdates<int, R>(vector<int>(N, N),
                                                          vector<int>(N, 1),
                                                          vals);
  }
  template <class ...Args>
  CountDistinctOfflineUpdates(const vector<T> &A,
                              const vector<pair<int, T>> &updates,
                              Args &&...args)
      : M(forward<Args>(args)...), N(A.size()), cur(0), ucnt(0),
        en(updates.size(), 0), wupds(init1(A, updates)), wm(init2(A)) {
    for (int k = 0; k < cur; k++) wm.updateToNextKey(wupds[k].first);
  }
  void processNextUpdate() {
    for (; cur < en[ucnt]; cur++) wm.updateToNextKey(wupds[cur].first);
    ucnt++;
  }
  int query(int l, int r) { return wm.query(l, r, l - 1); }
};


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

// Supports offline queries for the number of distinct elements in the
//   range [l, r] for a static array A of length N
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
// Constructor Arguments:
//   A: a vector of type T of the values in the array
//   queries: a vector of pairs containing the inclusive endpoints of
//     the queries
// Fields:
//   ans: a vector of integers with the answer for each query
// In practice, has a small constant, faster than Mo and the online version
// Time Complexity:
//   constructor: O(N log N + Q (log Q + log N))
// Memory Complexity: O(N + Q)
// Tested:
//   https://www.spoj.com/problems/DQUERY/
//   https://atcoder.jp/contests/abc174/tasks/abc174_f
template <class T> struct CountDistinctOffline {
  vector<int> ans;
  CountDistinctOffline(const vector<T> &A,
                       const vector<pair<int, int>> &queries)
      : ans(queries.size()) {
    int N = A.size(); vector<T> temp = A; sort(temp.begin(), temp.end());
    temp.erase(unique(temp.begin(), temp.end()), temp.end());
    vector<int> last(temp.size(), -1);
    vector<tuple<int, int, int>> q; int i = 0; BitFenwickTree ft(N);
    for (auto &&qi : queries) q.emplace_back(qi.first, qi.second, i++);
    sort(q.rbegin(), q.rend()); i = N - 1; for (auto &&qi : q) {
      int l, r, ind; tie(l, r, ind) = qi; for (; i >= l; i--) {
        int c = lower_bound(temp.begin(), temp.end(), A[i]) - temp.begin();
        if (last[c] != -1) ft.update(last[c], 0);
        ft.update(last[c] = i, 1);
      }
      ans[ind] = ft.query(r);
    }
  }
};


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

// Supports online queries for the number of distinct elements in the
//   range [l, r] for an array A of length N with point updates
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element
//   C: a container representing a mapping from type T to a set of integers
// Constructor Arguments:
//   A: a vector of type T of the initial values in the array
//   SCALE: the value to scale sqrt by
//   ...args: arguments to pass to the constructor of an instance of type C
// Functions:
//   update(i, v): updates index i with the value v
//   query(l, r): returns the number of distinct values in the range [l, r]
// In practice, has a very small constant, slightly slower than the
//   offline version
// Time Complexity:
//   constructor: O(N log N)
//   update: O(log N)
//   query: O(log N sqrt(N + U)) for U updates
// Memory Complexity: O((N + U) log N) for U updates
// Tested:
//   https://dmoj.ca/problem/mnyc17p6
template <class T, class C = map<T, set<int>>> struct CountDistinctUpdates {
  vector<T> A; C M; SemiSparseFenwickTree2DSimple<int> ft;
  void add(int i, const T &v) {
    set<int> &mv = M[v]; auto it = mv.upper_bound(++i);
    int j = it == mv.begin() ? 0 : *prev(it);
    ft.add(j, i); if (it != mv.end()) { ft.rem(j, *it); ft.add(i, *it); }
    mv.insert(it, i);
  }
  void rem(int i, const T &v) {
    set<int> &mv = M[v]; auto it = mv.erase(mv.find(++i));
    int j = it == mv.begin() ? 0 : *prev(it);
    ft.rem(j, i); if (it != mv.end()) { ft.rem(i, *it); ft.add(j, *it); }
  }
  void update(int i, const T &v) { rem(i, A[i]); add(i, A[i] = v); }
  int query(int l, int r) { return ft.query(l, l + 1, r + 1); }
  template <class ...Args>
  CountDistinctUpdates(const vector<T> &A, double SCALE = 1, Args &&...args)
      : A(A), M(forward<Args>(args)...), ft(A.size() + 1, SCALE) {
    for (int i = 0; i < int(A.size()); i++) add(i, A[i]);
  }
};


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

// Fast Subset Transform (similar to SumOverSubsets)
// Template Arguments:
//   TYPE: the type of susbet transform (OR, AND, or XOR)
//   T: the type of each element
// Function Arguments:
//   a: a reference to the vector of type T to transform
//   inv: a boolean indicating whether the inverse transform should be
//     performed or not
// In practice, has a small constant
// Time Complexity: O(N log N) where N = size(a)
// Memory Complexity: O(1)
// Tested:
//   https://judge.yosupo.jp/problem/subset_convolution
//   https://judge.yosupo.jp/problem/bitwise_and_convolution
//   https://judge.yosupo.jp/problem/bitwise_xor_convolution
//   https://csacademy.com/contest/round-53/task/maxor/
const int OR = 0, AND = 1, XOR = 2;
template <const int TYPE, class T> void fst(vector<T> &a, bool inv) {
  int N = a.size(); assert(!(N & (N - 1)));
  for (int len = 1; len < N; len <<= 1) for (int i = 0; i < N; i += len << 1)
    for (int j = 0; j < len; j++) {
      T &u = a[i + j], &v = a[len + i + j];
      if (TYPE == OR) tie(u, v) = inv ? make_pair(v, u - v)
                                      : make_pair(u + v, u);
      else if (TYPE == AND) tie(u, v) = inv ? make_pair(v - u, u)
                                            : make_pair(v, u + v);
      else if (TYPE == XOR) tie(u, v) = make_pair(u + v, u - v);
      else assert(false);
    }
  if (TYPE == XOR && inv) for (auto &&ai : a) ai /= T(N);
}


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

// Computes n!
// Template Arguments:
//   T: the type of n
// Function Arguments:
//   n: the value to find the factorial of, must be non negative
// Return Value: n!
// In practice has a very small constant
// Time Complexity: O(n)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> T factorial(T n) {
  T ret = 1; for (T i = 2; i <= n; i++) ret *= i;
  return ret;
}

// Computes n! modulo mod
// Template Arguments:
//   T: the type of n and mod
// Function Arguments:
//   n: the value to find the factorial of, must be non negative
//   mod: the modulo, must be not less than 2
// Return Value: n! modulo mod
// In practice has a very small constant
// Time Complexity: O(n)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> T factorialMod(T n, T mod) {
  T ret = 1; for (T i = 2; i <= n; i++) ret = mulMod(ret, i, mod);
  return ret;
}

// Computes n permute k (n! / (n - k)!)
// Template Arguments:
//   T: the type of n and k
// Function Arguments:
//   n: the total number of elements, must be non negative
//   k: the number of elements to permute, must be in the range [0, n]
// Return Value: n permute k
// In practice has a very small constant
// Time Complexity: O(min(k, n - k))
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> T permute(T n, T k) {
  T ret = 1; for (T i = 0; i < k; i++) ret = ret * (n - i);
  return ret;
}

// Computes n permute k (n! / (n - k)!) modulo mod
// Template Arguments:
//   T: the type of n, k, and mod
// Function Arguments:
//   n: the total number of elements, must be non negative
//   k: the number of elements to permute
//   mod: the modulo, must be not less than 2
// Return Value: n permute k modulo mod
// In practice has a very small constant
// Time Complexity: O(min(k, n - k))
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> T permuteMod(T n, T k, T mod) {
  T ret = 1; for (T i = 0; i < k; i++) ret = mulMod(ret, n - i, mod);
  return ret;
}

// Computes n choose k (n! / (k! (n - k)!))
// Template Arguments:
//   T: the type of n and k
// Function Arguments:
//   n: the total number of elements, must be non negative
//   k: the number of elements to choose, must be in the range [0, n]
// Return Value: n choose k
// In practice has a very small constant
// Time Complexity: O(min(k, n - k))
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> T choose(T n, T k) {
  if (k > n - k) k = n - k;
  T ret = 1; for (T i = 0; i < k; i++) ret = ret * (n - i) / (i + 1);
  return ret;
}

// Computes n choose k (n! / (k! (n - k)!)) modulo a prime p
// Template Arguments:
//   T: the type of n, k, and p
// Function Arguments:
//   n: the total number of elements, must be non negative
//   k: the number of elements to choose, must be in the range [0, n]
//   p: the prime modulo, must be greater than n
// Return Value: n permute k modulo p
// In practice has a very small constant
// Time Complexity: O(min(k, n - k) + log(p))
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> T chooseModPrime(T n, T k, T p) {
  if (k > n - k) k = n - k;
  T num = 1, den = 1; for (T i = 0; i < k; i++) {
    num = mulMod(num, n - i, p); den = mulMod(den, i + 1, p);
  }
  return divModPrime(num, den, p);
}

// Computes n multichoose k ((n + k - 1)! / (k! (n - 1)!))
// Template Arguments:
//   T: the type of n and k
// Function Arguments:
//   n: the total number of elements, must be positive
//   k: the number of elements to choose with repetitions, must be non negative
// Return Value: n choose k
// In practice has a very small constant
// Time Complexity: O(min(k, n - 1))
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> T multiChoose(T n, T k) { return choose(n + k - 1, k); }

// Computes  n multichoose k ((n + k - 1)! / (k! (n - 1)!)) modulo a prime p
// Template Arguments:
//   T: the type of n, k, and p
// Function Arguments:
//   n: the total number of elements, must be positive
//   k: the number of elements to choose with repetitions, must be non negative
//   p: the prime modulo
// Return Value: n permute k modulo p
// In practice has a very small constant
// Time Complexity: O(min(k, n - 1) + log(p))
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> T multiChooseModPrime(T n, T k, T p) {
  return chooseModPrime(n + k - 1, k, p);
}

// Combinatorics struct to compute factorials, permutations, and combinations
// Template Arguments:
//   T: the type of the value of N!
// Constructor Arguments:
//   N: the maximum value to compute factorials of
// Functions:
//   factorial(n): returns n! for 0 <= n <= N
//   permute(n, k): returns n permute k (n! / (n - k)!) for 0 <= k <= n <= N
//   choose(n, k): returns n permute k (n! / (k! (n - k)!))
//     for 0 <= k <= n <= N
//   multiChoose(n, k): returns n permute k ((n + k - 1)! / (k! (n - 1)!))
//     for n >= 1, k >= 0, n + k - 1 <= N
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(N)
//   factorial, permute, choose, multiChoose: O(1)
// Memory Complexity: O(N)
// Tested:
//   Fuzz Tested
template <class T> struct Combinatorics {
  vector<T> fact; Combinatorics(int N) : fact(N + 1, 1) {
    for (int i = 1; i <= N; i++) fact[i] = fact[i - 1] * i;
  }
  T factorial(int n) { return fact[n]; }
  T permute(int n, int k) { return fact[n] / fact[n - k]; }
  T choose(int n, int k) { return fact[n] / fact[k] / fact[n - k]; }
  T multiChoose(int n, int k) { return choose(n + k - 1, k); }
};

// Combinatorics struct to compute factorials, permutations, and combinations
//   modulo a prime P
// Template Arguments:
//   T: the type of the value of N! and P
// Constructor Arguments:
//   N: the maximum value to compute factorials of
//   P: the prime modulo P, must be greater than N
// Functions:
//   factorial(n): returns n! modulo P for 0 <= n <= N
//   invFactorial(n): returns multiplicative inverse of n! modulo P
//     for 0 <= n <= N
//   permute(n, k): returns n permute k (n! / (n - k)!) modulo P
//     for 0 <= k <= n <= N
//   choose(n, k): returns n permute k (n! / (k! (n - k)!)) modulo P
//     for 0 <= k <= n <= N
//   multiChoose(n, k): returns n permute k ((n + k - 1)! / (k! (n - 1)!))
//     modulo P for n >= 1, k >= 0, n + k - 1 <= N
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(N + log P)
//   factorial, invFactorial, permute, choose, multiChoose: O(1)
// Memory Complexity: O(N)
// Tested:
//   Fuzz Tested
template <class T> struct CombinatoricsModPrime {
  vector<T> fact, invFact; T P;
  CombinatoricsModPrime(int N, T P) : fact(N + 1, 1), invFact(N + 1), P(P) {
    for (int i = 1; i <= N; i++) fact[i] = mulMod(fact[i - 1], T(i), P);
    invFact[N] = mulInvModPrime(fact[N], P); for (int i = N - 1; i >= 0; i--)
      invFact[i] = mulMod(invFact[i + 1], T(i + 1), P);
  }
  T factorial(int n) { return fact[n]; }
  T invFactorial(int n) { return invFact[n]; }
  T permute(int n, int k) { return mulMod(fact[n], invFact[n - k], P); }
  T choose(int n, int k) {
    return mulMod(mulMod(fact[n], invFact[k], P), invFact[n - k], P);
  }
  T multiChoose(int n, int k) { return choose(n + k - 1, k); }
};

// Computes a row of Pascal's triangle
// Template Arguments:
//   T: the type of the values in the row of Pascal's triangle
// Function Arguments:
//   N: the row of Pascal's triangle to compute
// Return Value: a vector of length N + 1 with the values in the Nth row of
//   Pascal's triangle
// In practice, has a very small constant
// Time Complexity: O(N)
// Memory Complexity: O(N)
// Tested:
//   Fuzz Tested
template <class T> vector<T> pascalsRow(int N) {
  vector<T> C(N + 1, 1);
  for (int j = 1; j <= N; j++) C[j] = C[j - 1] * (N - j + 1) / j;
  return C;
}

// Computes a row of Pascal's triangle modulo a prime p
// Template Arguments:
//   T: the type of the values in the row of Pascal's triangle
// Function Arguments:
//   N: the row of Pascal's triangle to compute
//   p: the prime modulo, must be greater than N
// Return Value: a vector of length N + 1 with the values in the Nth row of
//   Pascal's triangle modulo a prime p
// In practice, has a small constant
// Time Complexity: O(N log N)
// Memory Complexity: O(N)
template <class T> vector<T> pascalsRowModPrime(int N, T p) {
  vector<T> C(N + 1, 1); for (int j = 1; j <= N; j++)
    C[j] = divModPrime(mulMod(C[j - 1], T(N - j + 1), p), T(j), p);
  return C;
}

// Computes Pascal's triangle
// Template Arguments:
//   T: the type of the values in Pascal's triangle
// Function Arguments:
//   N: the maximum row of Pascal's triangle to compute
// Return Value: a vector of vectors length N + 1 with the ith vector having
//   length i + 1 with the values in the ith row of Pascal's triangle
// In practice, has a very small constant
// Time Complexity: O(N^2)
// Memory Complexity: O(N^2)
// Tested:
//   Fuzz Tested
template <class T> vector<vector<T>> pascalsTriangle(int N) {
  vector<vector<T>> C(N + 1); for (int i = 0; i <= N; i++) {
    C[i] = vector<T>(i + 1, 1); for (int j = 1; j < i; j++)
      C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
  }
  return C;
}

// Computes Pascal's triangle modulo mod
// Template Arguments:
//   T: the type of the values in Pascal's triangle and mod
// Function Arguments:
//   N: the maximum row of Pascal's triangle to compute
//   mod: the modulo, must be not less than 2
// Return Value: a vector of vectors length N + 1 with the ith vector having
//   length i + 1 with the values in the ith row of Pascal's triangle
//   modulo mod
// In practice, has a very small constant
// Time Complexity: O(N^2)
// Memory Complexity: O(N^2)
// Tested:
//   Fuzz Tested
template <class T> vector<vector<T>> pascalsTriangleMod(int N, T mod) {
  vector<vector<T>> C(N + 1); for (int i = 0; i <= N; i++) {
    C[i] = vector<T>(i + 1, 1); for (int j = 1; j < i; j++)
      C[i][j] = addMod(C[i - 1][j - 1], C[i - 1][j], mod);
  }
  return C;
}

// Computes the sum from 1 to n inclusive
// Template Arguments:
//   T: the type of n
// Function Arguments:
//   n: the maximum bound, must be non negative
// Return Value: the sum from 1 to n inclusive
// In practice has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> T sumTo(T n) { return n * (n + 1) / 2; }

// Computes the sum from a to b inclusive
// Template Arguments:
//   T: the type of n
// Function Arguments:
//   a: the lower bound, must be non negative
//   b: the upper bound, must be non negative and not less than a
// Return Value: the sum from a to b inclusive
// In practice has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> T sumBetween(T a, T b) { return sumTo(b) - sumTo(a - 1); }

// Computes the nth term of an arithmetic sequence with starting value a1 and
//   common difference d
// Template Arguments:
//   T: the type of a1 and d
//   U: the type of n
// Function Arguments:
//   a1: the starting value
//   d: the common difference
//   n: the index of the term to find, must be positive
// Return Value: the nth term of the arithmetic sequence
// In practice has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T, class U>
T arithSeq(T a1, T d, U n) { return a1 + d * (n - 1); }

// Computes the sum of an arithmetic sequence with starting value a1 and
//   common difference d from term 1 to n
// Template Arguments:
//   T: the type of a1 and d
//   U: the type of n
// Function Arguments:
//   a1: the starting value
//   d: the common difference
//   n: the maximum index of the term to sum, must be positive
// Return Value: the sum of the arithmetic sequence from term 1 to n
// In practice has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
//   Fuzz Tested
template <class T, class U>
T arithSeries(T a1, T d, U n) { return n * (a1 + arithSeq(a1, d, n)) / 2; }

// Computes the nth term of an geometric sequence with starting value a1 and
//   common ratio r
// Template Arguments:
//   T: the type of a1 and r
//   U: the type of n
// Function Arguments:
//   a1: the starting value
//   r: the common ratio
//   n: the index of the term to find, must be positive
// Return Value: the nth term of the geometric sequence
// In practice has a small constant
// Time Complexity: O(log n)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T, class U>
T geoSeq(T a1, T r, U n) { return a1 * pow2(r, n - 1); }

// Computes the sum of a geometric sequence with starting value a1 and
//   common ratio r from term 1 to n
// Template Arguments:
//   T: the type of a1 and r
//   U: the type of n
// Function Arguments:
//   a1: the starting value
//   r: the common ratio
//   n: the maximum index of the term to sum, must be positive
// Return Value: the sum of the geometric sequence from term 1 to n
// In practice has a small constant
// Time Complexity: O(log n)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T, class U> T geoSeries(T a1, T r, U n) {
  return r == 1 ? a1 * n : a1 * (T(1) - pow2(r, n)) / (T(1) - r);
}

// Computes the sum of a infinite geometric sequence with starting value a1 and
//   common ratio r where -1 < r < 1
// Template Arguments:
//   T: the type of a1 and r
// Function Arguments:
//   a1: the starting value
//   r: the common ratio
// Return Value: the sum of the infinite geometric sequence
// In practice has a small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
template <class T>
T infGeoSeries(T a1, T r) { assert(-1 < r && r < 1); return a1 / (T(1) - r); }

// Computes the sum of floor((a * i + b) / m) for 0 <= i < n
// Equivalent to the number of integer points (x, y) where 0 <= x < n and
//   0 <= y <= (a * i + b) / m
// Template Arguments:
//   T: an integral type
// Function Arguments:
//   n: the upper bound for x
//   a: the linear factor's numerator
//   b: the constant factor's numerator
//   m: the common denominator
// Return Value: the sum of floor((a * i + b) / m) for 0 <= i < n
// In practice, has a small constant
// Time Complexity: O(log n)
// Memory Complexity: O(1)
// Tested:
//   https://judge.yosupo.jp/problem/sum_of_floor_of_linear
template <class T> T sumFloorLinear(T n, T a, T b, T m) {
  static_assert(is_integral<T>::value, "T must be integral");
  if (a == 0) return (b / m) * n;
  if (a >= m || b >= m) return ((a / m) * (n - 1) + 2 * (b / m)) * n / 2
                               + sumFloorLinear(n, a % m, b % m, m);
  return sumFloorLinear((a * n + b) / m, m, (a * n + b) % m, a);
}


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

// Data structure representing a fraction
// Template Arguments:
//   T: the type of the numerator and denominator, must be integral
// Constructor Arguments:
//   num: the numerator
//   den: the denominator, should be non negative
// Fields:
//   num: the numerator
//   den: the denominator
// Functions:
//   reduce(): reduces the fraction to lowest terms with a positive denominator
//   +, +=, -, -=, *, *=, /, /=: standard arithmetic operators on fractions
//   <, <=, >, >=, ==, !=: comparison operators
//   >>, <<: input and output operators
// Time Complexity:
//   constructor, +, +=, -, -=, *, *=, /, /=, >>, <<: O(1)
//   reduce, <, <=, >, >=, ==, !=: O(log(max(num, den)))
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/bfs19p4
template <class T> struct Fraction {
  static_assert(is_integral<T>::value, "T must be an integral type");
  using F = Fraction<T>; T num, den;
  constexpr Fraction(T num = T(), T den = T(1)) : num(num), den(den) {}
  F reduce() const {
    T g = gcd(num, den);
    return den >= 0 ? F(num / g, den / g) : F(-num / g, -den / g);
  }
  F operator + (const F &f) const { return F(*this) += f; }
  F &operator += (const F &f) {
    num = (num * f.den + f.num * den); den = den * f.den; return *this;
  }
  F operator - (const F &f) const { return F(*this) -= f; }
  F &operator -= (const F &f) {
    num = (num * f.den - f.num * den); den = den * f.den; return *this;
  }
  F operator + () const { return *this; }
  F operator - () const { return F(-num, den); }
  F operator * (const F &f) const { return F(*this) *= f; }
  F &operator *= (const F &f) {
    num = num * f.num; den = den * f.den; return *this;
  }
  F operator / (const F &f) const { return F(*this) /= f; }
  F &operator /= (const F &f) {
    T t_num = num * f.den, t_den = den * f.num; num = t_num; den = t_den;
    return *this;
  }
  bool operator < (const F &f) const {
    F a = den >= 0 ? *this : F(-num, -den);
    F b = f.den >= 0 ? f : F(-f.num, -f.den);
    return a.num * b.den < b.num * a.den;
  }
  bool operator <= (const F &f) const { return !(f < *this); }
  bool operator > (const F &f) const { return f < *this; }
  bool operator >= (const F &f) const { return !(*this < f); }
  bool operator == (const F &f) const { return !(*this < f) && !(f < *this); }
  bool operator != (const F &f) const { return *this < f || f < *this; }
  friend istream &operator >> (istream &stream, F &f) {
    return stream >> f.num >> f.den;
  }
  friend ostream &operator << (ostream &stream, const F &f) {
    return stream << f.num << '/' << f.den;
  }
};

#define ftype long double
template <class T> Fraction<T> abs(Fraction<T> a) { return a >= 0 ? a : -a; }
#define FUN(name) \
  template <class T> ftype name(Fraction<T> a) { \
    return name((ftype)(a.num) / (ftype)(a.den)); \
  }
FUN(sqrt) FUN(sin) FUN(cos) FUN(tan) FUN(asin) FUN(acos) FUN(atan)
#undef FUN
template <class T> ftype atan2(Fraction<T> y, Fraction<T> x) {
  return atan2((ftype)(y.num) / (ftype)(y.den),
               (ftype)(x.num) / (ftype)(x.den));
}
#undef ftype


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

// Karatsuba Multiplication
// Template Arguments:
//   T: the type of each element
//   ItA: the type of the iterator pointing to the first multiplyPolynomial
//   ItB: the type of the iterator pointing to the second polynomial
//   ItRes: the type of the iterator pointing to the result polynomial
// Function Arguments:
//   n: the length of the polynomial, must be a power of 2
//   a: a pointer to the first polynomial of length n, a[i] stores
//     the coefficient of x^i
//   b: a pointer to the second polynomial of length n, b[i] stores
//     the coefficient of x^i
//   res: a pointer to the result polynomial of length 2n which stores the
//     result of a multipled by b, res[i] stores the coefficient of x^i
// In practice, has a small constant, slower than FFT and NTT
// Time Complexity: O(n^log_2(3))
// Memory Complexity: O(n)
// Tested:
//   https://open.kattis.com/problems/polymul2
//   https://dmoj.ca/problem/atimesb
const int KARATSUBA_CUTOFF = 16;
template <class T, class ItA, class ItB, class ItRes>
void karatsuba(int n, ItA a, ItB b, ItRes res) {
  if (n <= KARATSUBA_CUTOFF) {
    fill(res, res + n * 2, 0); for (int i = 0; i < n; i++)
      for (int j = 0; j < n; j++) res[i + j] += a[i] * b[j];
    return;
  }
  assert(!(n & (n - 1))); int k = n / 2; vector<T> tmp(n, T()), c(n, T());
  auto atmp = tmp.begin(), btmp = atmp + k; for (int i = 0; i < k; i++) {
    atmp[i] = a[i] + a[i + k]; btmp[i] = b[i] + b[i + k];
  }
  karatsuba<T>(k, atmp, btmp, c.begin()); karatsuba<T>(k, a, b, res);
  karatsuba<T>(k, a + k, b + k, res + n); for (int i = 0; i < k; i++) {
    T t = res[i + k]; res[i + k] += c[i] - res[i] - res[i + k * 2];
    res[i + k * 2] += c[i + k] - t - res[i + k * 3];
  }
}

// Polynomial Multiplication
// Template Arguments:
//   T: the type of each element
// Functions Arguments:
//   a: a vector of type T representing the first polynomial, a[i] stores
//     the coefficient of x^i
//   b: a vector of type T representing the second polynomial, b[i] stores
//     the coefficient of x^i
// Return Value: a vector of type T representing the polynomial a times b with
//   no trailing zeros
// In practice, has a small constant
// Time Complexity: O(n^log_2(3)) where n = max(size(a), size(b))
// Memory Complexity: O(n)
// Tested:
//   https://open.kattis.com/problems/polymul2
template <class T> vector<T> mulPolyKaratsuba(vector<T> a, vector<T> b) {
  int n = max(a.size(), b.size()); while (n & (n - 1)) n++;
  a.resize(n, T()); b.resize(n, T()); vector<T> res(n * 2, T());
  karatsuba<T>(n, a.begin(), b.begin(), res.begin());
  while (int(res.size()) > 1 && res.back() == T()) res.pop_back();
  return res;
}

// Integer Multiplication
// Template Arguments:
//   T: the type of each digit
// Functions Arguments:
//   a: a vector of type T representing the first integer, a[i] stores
//     the digit of BASE^i
//   b: a vector of type T representing the second polynomial, b[i] stores
//     the digit of BASE^i
//   BASE: the base of the integer
// Return Value: a vector of type T representing the integer a times b with
//   no trailing zeros
// In practice, has a small constant
// Time Complexity: O(n^log_2(3)) where n = max(size(a), size(b))
// Memory Complexity: O(n)
// Tested:
//   https://dmoj.ca/problem/atimesb
template <class T>
vector<T> mulIntKaratsuba(vector<T> a, vector<T> b, T BASE = T(10)) {
  vector<T> res = mulPolyKaratsuba(move(a), move(b)); T carry = T();
  res.reserve(res.size() + 1); for (int i = 0; i < int(res.size()); i++) {
    T cur = res[i] + carry; res[i] = T(cur % BASE); carry = T(cur / BASE);
  }
  if (carry > T()) res.push_back(carry);
  return res;
}


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

// Computes the floor of the sqrt of an integer
// Template Arguments:
//   T: the type of x, must be integral
// Function Arguments:
//   x: the value to compute the sqrt of, must be non-negative
// Return Value: returns the floor of the sqrt of x
// In practice, has a moderate constant, slower than binary search, but this
//   does not overflow
// Time Complexity: O(log x)
// Memory Complexity: O(1)
// Tested:
//   https://codeforces.com/gym/100753/problem/F
template <class T> T isqrt(T x) {
  static_assert(is_integral<T>::value, "T must be integral");
  assert(x >= 0); if (x <= 1) return x;
  T y = x, z = 0; while (y > z) z = x / (y = (y + z) / 2);
  return y;
}


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

// Computes the nth tetration of a modulo mod
// Assumes 0^0 = 1
// Template Arguments:
//   T: the type of base and mod
//   U: the type of the number of iterations
// Function Arguments:
//   a: the base
//   n: the number of iterations
//   mod: the modulo
// Return Value: the nth tetration of a modulo mod,
//   a value in the range [0, mod)
// In practice, has a small constant
// Time Complexity: O(sqrt(mod) log(mod))
// Memory Complexity: O(log(mod))
// Tested:
//   https://judge.yosupo.jp/problem/tetration_mod
template <class T, class U> T tetraMod(T a, U n, T mod) {
  if (mod == T(1)) return T(0);
  if (a == T(0)) return (T(1) - T(n % U(2))) % mod;
  if (n == U(0)) return T(1) % mod;
  if (n == U(1)) return a % mod;
  if (n == U(2)) return powMod(a % mod, a, mod);
  T p = phi(mod), ret = tetraMod(a, n - 1, p); if (ret == 0) ret += p;
  return powMod(a % mod, ret, mod);
}

// Helper function returning a * b % mod + mod if a * b >= mod,
//   a * b % mod otherwise
template <class T> T mulMod2(T a, T b, T mod) {
  return (a *= b) >= mod ? a % mod + mod : a;
}

// Helper function returning base^pow % mod + mod if base^pow >= mod,
//   base^pow % mod otherwise
template <class T, class U> T powMod2(T base, U pow, T mod) {
  if (base >= mod) base = base % mod + mod;
  T x = 1; while (true) {
    if (pow % 2 == 1) x = mulMod2(x, base, mod);
    if ((pow /= 2) == 0) break;
    base = mulMod2(base, base, mod);
  }
  return x;
}

// Helper function
template <class T> T powSeqModRec(const vector<T> &A, int i, T mod) {
  if (i == int(A.size()) - 1) return A[i];
  T p = phi(mod), res = powSeqModRec(A, i + 1, p);
  return res == T(0) ? T(1) : powMod2(A[i], res, mod);
}

// Computes A[0]^A[1]^...^A[N - 1] modulo mod
// Assumes 0^0 = 1
// Template Arguments:
//   T: the type of each element in the sequence
// Function Arguments:
//   A: the sequence
//   mod: the modulo
// Return Value: A[0]^A[1]^...^A[N - 1] modulo mod,
//   a value in the range [0, mod)
// In practice, has a small constant
// Time Complexity: O(N + sqrt(mod) log(mod))
// Memory Complexity: O(N + log(mod))
// Tested:
//   https://dmoj.ca/problem/dmopc20c3p5
template <class T> T powSeqMod(const vector<T> &A, T mod) {
  return powSeqModRec(A, 0, mod) % mod;
}

// Finds an integer x such that (x^2 = a) modulo a prime p
// If x is a solution, then (p - x) % p is the other solution
// Template Arguments:
//   T: the type of a and p
// Function Arguments:
//   a: the argument, must be in the range [0, p)
//   p: the prime modulo
// Return Value: the square root of a modulo p in the range [0, p), or -1 if
//   it does not exist,
// In practice, has a small constant
// Time Complexity: O(log(p)^2)
// Memory Complexity: O(1)
// Tested:
//   https://judge.yosupo.jp/problem/sqrt_mod
template <class T> T sqrtMod(T a, T p) {
  if (a == T(0)) return 0;
  if (powMod(a, (p - 1) / 2, p) != T(1)) return -1;
  if (p % 4 == 3) return powMod(a, (p + 1) / 4, p);
  T s = p - 1, n = 2; int r = 0; for (; s % 2 == 0; s /= 2) r++;
  while (powMod(n % p, (p - 1) / 2, p) != p - 1) n++;
  T x = powMod(a, (s + 1) / 2, p), b = powMod(a, s, p);
  T g = powMod(n % p, s, p); while (true) {
    T t = b, m = T(0); for (; m < r && t != T(1); m++) t = mulMod(t, t, p);
    if (m == T(0)) return x;
    T gs = g; for (int i = 0; i < r - m - 1; i++) gs = mulMod(gs, gs, p);
    g = mulMod(gs, gs, p); x = mulMod(x, gs, p); b = mulMod(b, g, p); r = m;
  }
  return -1;
}

// Finds the smallest non negative integer c such that
//   (a^c = b) modulo mod
// Assumes 0^0 = 1
// Template Arguments:
//   T: the type of a, b, and mod
// Function Arguments:
//   a: the base, must be in the range [0, mod)
//   b: the argument, must be in the range [0, mod)
//   mod: the modulo
// Return Value: the discrete log of b in base a modulo mod, or -1 if
//   it does not exist
// In practice, has a small constant
// Time Complexity: O(sqrt(mod) log(mod))
// Memory Complexity: O(sqrt(mod))
// Tested:
//   https://judge.yosupo.jp/problem/discrete_logarithm_mod
template <class T> T discreteLogMod(T a, T b, T mod) {
  if (b == T(1) % mod) return 0;
  T n = max(T(0), T(sqrt(mod)) - 1), e = T(1), f = T(1), j = T(1);
  vector<pair<T, T>> A; while (n * n <= mod) n++;
  while (j <= n && (e = f = mulMod(e, a, mod)) != b)
    A.emplace_back(mulMod(e, b, mod), j++);
  if (e == b) return j;
  sort(A.begin(), A.end()); if (gcd(mod, e) == gcd(mod, b)) {
    for (T i = 2; i <= n + 1; i++) {
      e = mulMod(e, f, mod);
      auto it = lower_bound(A.begin(), A.end(), make_pair(e, T(0)));
      if (it != A.end() && it->first == e) return n * i - it->second;
    }
  }
  return -1;
}


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

// Given a set of N points, P[i] = (X[i], Y[i]), compute a polynomial Q of
//   degree N - 1, such that Q passes through them:
//   Q(x) = A[0] * x^0 + ... + A[N - 1] * x^(N - 1).
// For real numbers, pick X[i] = C * cos(i / (N - 1) * PI) for a
//   large constant C
// Template Arguments:
//   T: the type of the x and y values, as well as the coefficients
// Function Arguments:
//   P: a vector of pairs of type T representing the N points
// Return Value: a polynomial of degree N - 1 (with N coefficients) with
//   the ith coefficient associated with the term x^i
// In practice, has a small constant
// Time Complexity: O(N^2)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/angieandfunctions
template <class T>
vector<T> lagrangePolynomialInterpolation(vector<pair<T, T>> P) {
  int N = P.size(); vector<T> A(N, T()), temp(N, T()); temp[0] = T(1);
  for (int k = 0; k < N; k++) for (int i = k + 1; i < N; i++)
    P[i].second = (P[i].second - P[k].second) / (P[i].first - P[k].first);
  T last = T(); for (int k = 0; k < N; k++) for (int i = 0; i < N; i++) {
    A[i] += P[k].second * temp[i]; swap(last, temp[i]);
    temp[i] -= P[k].first * last;
  }
  return A;
}