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

// Computes the non negative remainder of a by mod
// Template Arguments:
//   T: the type of a and mod
// Function Arguments:
//   a: the value to find the non negative remainder
//   mod: the modulo
// Return Value: the non negative remainder of a by mod
// In practice, has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/bts18p5
template <class T> T posMod(T a, T mod) {
  T ret = -mod < a && a < mod ? a : a % mod; return 0 <= ret ? ret : ret + mod;
}

// Adds two integers a and b modulo mod
// Template Arguments:
//   T: the type of a, b, and mod
// Function Arguments:
//   a: the first value, must be in the range [0, mod)
//   b: the second value, must be in the range [0, mod)
//   mod: the modulo
// Return Value: a + b modulo mod, a value in the range [0, mod)
// In practice, has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/modulararithmetic
//   https://dmoj.ca/problem/bts18p5
template <class T> T addMod(T a, T b, T mod) {
  T ret = a + b; return ret < mod ? ret : ret - mod;
}

// Subtracts two integers a and b modulo mod
// Template Arguments:
//   T: the type of a, b, and mod
// Function Arguments:
//   a: the first value, must be in the range [0, mod)
//   b: the second value, must be in the range [0, mod)
//   mod: the modulo
// Return Value: a - b modulo mod, a value in the range [0, mod)
// In practice, has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/modulararithmetic
template <class T> T subMod(T a, T b, T mod) {
  return a >= b ? a - b : a + mod - b;
}

// Multiplies two integers a and b modulo mod where a * b does not overflow
// Template Arguments:
//   T: the type of a, b, and mod
// Function Arguments:
//   a: the first value, must be in the range [0, mod)
//   b: the second value, must be in the range [0, mod)
//   mod: the modulo
// Return Value: a * b modulo mod, a value in the range [0, mod)
// In practice, has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/modulararithmetic
//   https://atcoder.jp/contests/m-solutions2019/tasks/m_solutions2019_e
//   https://dmoj.ca/problem/bts18p5
template <class T> T mulMod(T a, T b, T mod) { return a * b % mod; }

// Computes base to the power pow modulo mod where mod * mod does not overflow
// Assumes 0^0 = 1
// Template Arguments:
//   T: the type of base and mod
//   U: the type of pow
// Function Arguments:
//   base: the base, must be in the range [0, mod)
//   pow: the power, must be non negative and integral
//   mod: the modulo
// Return Value: base to the power pow modulo mod,
//   a value in the range [0, mod)
// In practice, has a small constant
// Time Complexity: O(log pow)
// Memory Complexity: O(1)
// Tested:
//   https://atcoder.jp/contests/m-solutions2019/tasks/m_solutions2019_e
template <class T, class U> T powMod(T base, U pow, T mod) {
  T x = 1; while (true) {
    if (pow % 2 == 1) x = mulMod(x, base, mod);
    if ((pow /= 2) == 0) break;
    base = mulMod(base, base, mod);
  }
  return x;
}

// Computes the multiplicative inverse of a for a prime mod p where p * p
//   does not overflow
// Guaranteed to exist
// Template Arguments:
//   T: the type of a and p
// Function Arguments:
//   a: the value to find the inverse of, must be in the range [0, p)
//   p: the prime modulo
// Return Value: the multiplicative inverse of a for a prime mod p
// In practice, has a small constant
// Time Complexity: O(log p)
// Memory Complexity: O(1)
// Tested:
//   https://atcoder.jp/contests/m-solutions2019/tasks/m_solutions2019_e
template <class T> T mulInvModPrime(T a, T p) { return powMod(a, p - 2, p); }

// Divides a by b modulo a prime p where p * p does not overflow
// Guaranteed to be valid
// Template Arguments:
//   T: the type of a, b and p
// Function Arguments:
//   a: the dividend, must be in the range [0, p)
//   b: the divisor, must be in the range [1, p)
//   p: the prime modulo
// Return Value: a / b modulo p
// In practice, has a small constant
// Time Complexity: O(log p)
// Memory Complexity: O(1)
// Tested:
//   https://atcoder.jp/contests/m-solutions2019/tasks/m_solutions2019_e
template <class T> T divModPrime(T a, T b, T p) {
  return mulMod(a, mulInvModPrime(b, p), p);
}

// Divides a by b modulo mod
// Only valid if gcd(b, mod) == 1
// Template Arguments:
//   T: the type of a, b and mod
// Function Arguments:
//   a: the dividend, must be in the range [0, mod)
//   b: the divisor, must be in the range [1, mod)
//   mod: the modulo
// Return Value: a / b modulo mod, or -1 if invalid
// In practice, has a small constant
// Time Complexity: O(log p)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/modulararithmetic
template <class T> T divMod(T a, T b, T mod) {
  T inv = mulInv(b, mod); return inv == -1 ? -1 : mulMod(a, inv, mod);
}

// The following functions are useful when mod * mod overflow

// Multiplies two integers a and b modulo mod where a * b does overflow, but
//   mod + mod does not overflow
// Template Arguments:
//   T: the type of a, b, and mod
// Function Arguments:
//   a: the first value, must be in the range [0, mod)
//   b: the second value, must be in the range [0, mod)
//   mod: the modulo
// Return Value: a * b modulo mod, a value in the range [0, mod)
// In practice, has a small constant
// Time Complexity: O(log b)
// Memory Complexity: O(1)
// Tested:
//   https://atcoder.jp/contests/m-solutions2019/tasks/m_solutions2019_e
template <class T> T mulModOvf(T a, T b, T mod) {
  T x = 0; while (true) {
    if (b % 2 == 1) x = addMod(x, a, mod);
    if ((b /= 2) == 0) break;
    a = addMod(a, a, mod);
  }
  return x;
}

// Computes base to the power pow modulo mod where mod * mod overflows, but
//   mod + mod does not overflow
// Template Arguments:
//   T: the type of base and mod
//   U: the type of pow
// Function Arguments:
//   base: the base, must be in the range [0, mod)
//   pow: the power, must be non negative and integral
//   mod: the modulo
// Return Value: base to the power pow modulo mod,
//   a value in the range [0, mod)
// In practice, has a small constant
// Time Complexity: O(log pow log mod)
// Memory Complexity: O(1)
// Tested:
//   https://atcoder.jp/contests/m-solutions2019/tasks/m_solutions2019_e
template <class T, class U> T powModOvf(T base, U pow, T mod) {
  T x = 1; while (true) {
    if (pow % 2 == 1) x = mulModOvf(x, base, mod);
    if ((pow /= 2) == 0) break;
    base = mulModOvf(base, base, mod);
  }
  return x;
}

// Computes the multiplicative inverse of a for a prime mod p where p * p
//   overflows, but mod + mod does not overflow
// Guaranteed to exist
// Template Arguments:
//   T: the type of a and p
// Function Arguments:
//   a: the value to find the inverse of, must be in the range [0, p)
//   p: the prime modulo
// Return Value: the multiplicative inverse of a for a prime mod p
// In practice, has a small constant
// Time Complexity: O(log p log mod)
// Memory Complexity: O(1)
// Tested:
//   https://atcoder.jp/contests/m-solutions2019/tasks/m_solutions2019_e
template <class T> T mulInvModPrimeOvf(T a, T p) {
  return powModOvf(a, p - 2, p);
}

// Divides a by b modulo a prime p where p * p overflows, but mod + mod does
//   not overflow
// Guaranteed to be valid
// Template Arguments:
//   T: the type of a, b and p
// Function Arguments:
//   a: the dividend, must be in the range [0, p)
//   b: the divisor, must be in the range [1, p)
//   p: the prime modulo
// Return Value: a / b modulo p
// In practice, has a small constant
// Time Complexity: O(log p log mod)
// Memory Complexity: O(1)
// Tested:
//   https://atcoder.jp/contests/m-solutions2019/tasks/m_solutions2019_e
template <class T> T divModPrimeOvf(T a, T b, T p) {
  return mulModOvf(a, mulInvModPrimeOvf(b, p), p);
}

// Struct supporting operations in Montgomery form
// Constructor Arguments:
//   mod: the modulo of the space
// Functions:
//   init(x): transforms a number into Montgomery form
//   reduce(x): transforms a number from Montgomery form
//   mul(a, b): multiplies the numbers a and b in Montgomery form and returns
//     their product modulo mod in Montgomery form
// In practice, has a small constant
// Time Complexity:
//   constructor, init, reduce, mul: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://loj.ac/p/6466
//   https://www.spoj.com/problems/FACT2/
struct Montgomery {
  using u64 = uint64_t; using u128 = __uint128_t; using s128 = __int128_t;
  struct u256 {
    static u128 HI(u128 x) { return x >> 64; }
    static u128 LO(u128 x) { return u64(x); }
    u128 hi, lo; u256(u128 lo = 0) : hi(0), lo(lo) {}
    u256(u128 hi, u128 lo) : hi(hi), lo(lo) {}
    static u256 mul(u128 x, u128 y) {
      u128 t1 = LO(x) * LO(y), t2 = HI(x) * LO(y) + HI(y) * LO(x) + HI(t1);
      return u256(HI(x) * HI(y) + HI(t2), (t2 << 64) + LO(t1));
    }
  };
  u128 mod, inv, r2;
  Montgomery(u128 mod = 1) : mod(mod), inv(1), r2(-mod % mod) {
    for (int i = 0; i < 7; i++) inv *= 2 - mod * inv;
    for (int i = 0; i < 4; i++) if ((r2 <<= 1) >= mod) r2 -= mod;
    for (int i = 0; i < 5; i++) r2 = mul(r2, r2);
  }
  u128 init(u128 x) { return mul(x, r2); }
  u128 reduce(u256 x) {
    u128 q = x.lo * inv; s128 a = x.hi - u256::mul(q, mod).hi;
    return a < 0 ? a + mod : a;
  }
  u128 mul(u128 a, u128 b) { return reduce(u256::mul(a, b)); }
};

// Multiplies two unsigned 128-bit integers a and b modulo mod even if a * b
//   overflows, using Montgomery reduction
// Function Arguments:
//   a: the first value, must be in the range [0, mod)
//   b: the second value, must be in the range [0, mod)
//   mod: the modulo
// Return Value: a * b modulo mod, a value in the range [0, mod)
// In practice, has a small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://loj.ac/p/6466
//   https://www.spoj.com/problems/FACT2/
__uint128_t mulMod(__uint128_t a, __uint128_t b, __uint128_t mod) {
  static Montgomery mont; if (mont.mod != mod) mont = Montgomery(mod);
  return mont.reduce(mont.mul(mont.init(a), mont.init(b)));
}

// Barrett Reduction for fast modulo
// Constructor Arguments:
//   mod: the modulo of the space
// Functions:
//   reduce(a): returns a value congruent to a modulo mod in the range
//     [0, 2 mod)
// In practice, has a very small constant
// Time Complexity:
//   constructor, reduce: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/modulararithmetic
struct Barrett {
  uint64_t mod, inv;
  Barrett(uint64_t mod = 1) : mod(mod), inv(uint64_t(-1) / mod) {}
  uint64_t reduce(uint64_t a) {
    return a - uint64_t((__uint128_t(inv) * a) >> 64) * mod;
  }
};

// Fast modulo using Barrett Reduction
// Function Arguments:
//   a: the value to modulo
//   mod: the modulo
// Return Value: a modulo mod
// In practice, has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/modulararithmetic
uint64_t bmod(uint64_t a, uint64_t mod) {
  static Barrett b; if (b.mod != mod) b = Barrett(mod);
  uint64_t ret = b.reduce(a); if (ret >= mod) ret -= mod;
  return ret;
}


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

// Big Integer
// Constructor Arguments:
//   s: a string of digits
//   v: a long long
// Functions:
//   read(s): reads a BigInt from a string
//   write(s): write a BigInt to a string
//   >>, <<: input and output operators
//   <, <=, >, >=, ==, !=: comparison operators
//   isZero(): returns true if the big integer is 0
//   abs(): returns the absolute value
//   ++, --, +, +=, -, -=, *, *=, /, /=, %, %=: standard arithmetic operators
//   divMod(a1, b1): returns a pair with the quotient and remainder of
//     a1 divided by b1
// In practice, multiplication has a moderate constant, all other functions
//   have a small constant
// Time Complexity:
//   constructor: O(N)
//   isZero: O(1)
//   read, write, >>, <<, <, <=, >, >=, ==, !=, abs, ++, --, +, +=, -, -=: O(N)
//   *, *= (long long variants): O(N)
//   *, *= (BigInt variants): O(N log N)
//   /, /=, % (long long variants): O(N)
//   /, /=, %, %= (BigInt variants): O(N^2)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/aplusb2
//   https://dmoj.ca/problem/atimesb
//   https://dmoj.ca/problem/ccc96s2
//   https://dmoj.ca/problem/ccc97s5
struct BigInt {
  using T = long long; static constexpr const int DIG = 4;
  static constexpr const T BASE = pow(10, DIG);
  vector<T> a; int sign;
  void trim() {
    while (!a.empty() && a.back() == 0) a.pop_back();
    if (a.empty()) sign = 1;
  }
  void read(const string &s) {
    sign = 1; a.clear(); int pos = 0;
    for (; pos < int(s.size()) && (s[pos] == '-' || s[pos] == '+'); ++pos)
      if (s[pos] == '-') sign = -sign;
    for (int i = int(s.size()) - 1; i >= pos; i -= DIG) {
      T x = 0;
      for (int j = max(pos, i - DIG + 1); j <= i; j++) x = x * 10 + s[j] - '0';
      a.push_back(x);
    }
    trim();
  }
  friend istream& operator >> (istream &stream, BigInt &v) {
    string s; stream >> s; v.read(s); return stream;
  }
  string write() const {
    string ret = ""; if (sign == -1) ret.push_back('-');
    for (char c : to_string(a.empty() ? 0 : a.back())) ret.push_back(c);
    for (int i = int(a.size()) - 2; i >= 0; i--) {
      string s = to_string(a[i]);
      for (int j = int(s.size()); j < DIG; j++) ret.push_back('0');
      for (char c : s) ret.push_back(c);
    }
    return ret;
  }
  friend ostream& operator << (ostream &stream, const BigInt &v) {
    stream << v.write(); return stream;
  }
  bool operator < (const BigInt &v) const {
    if (sign != v.sign) return sign < v.sign;
    if (int(a.size()) != int(v.a.size()))
      return int(a.size()) * sign < int(v.a.size()) * v.sign;
    for (int i = int(a.size()) - 1; i >= 0; i--) if (a[i] != v.a[i])
      return a[i] * sign < v.a[i] * sign;
    return false;
  }
  bool operator <= (const BigInt &v) const { return !(v < *this); }
  bool operator > (const BigInt &v) const { return v < *this; }
  bool operator >= (const BigInt &v) const { return !(*this < v); }
  bool operator == (const BigInt &v) const {
    return !(*this < v) && !(v < *this);
  }
  bool operator != (const BigInt &v) const { return *this < v || v < *this; }
  BigInt() : sign(1) {}
  BigInt(const string &s) { read(s); }
  BigInt(T v) : BigInt(to_string(v)) {}
  bool isZero() const { return a.empty() || (int(a.size()) == 1 && !a[0]); }
  BigInt operator + () const { return *this; }
  BigInt operator - () const {
    BigInt res = *this; res.sign = -sign; return res;
  }
  BigInt abs() const { BigInt res = *this; res.sign *= res.sign; return res; }
  T value() const {
    T res = 0;
    for (int i = int(a.size()) - 1; i >= 0; i--) res = res * BASE + a[i];
    return res * sign;
  }
  BigInt operator ++ () { return *this += BigInt(1); }
  BigInt operator ++ (int) {
    BigInt ret = *this; *this += BigInt(1); return ret;
  }
  BigInt operator -- () { return *this -= BigInt(1); }
  BigInt operator -- (int) {
    BigInt ret = *this; *this -= BigInt(1); return ret;
  }
  BigInt operator + (const BigInt &v) const {
    if (sign == v.sign) {
      BigInt res = v; T carry = 0;
      for (int i = 0; i < max(int(a.size()), int(v.a.size())) || carry; i++) {
        if (i == int(res.a.size())) res.a.push_back(0);
        res.a[i] += carry + (i < int(a.size()) ? a[i] : 0);
        carry = res.a[i] >= BASE; if (carry) res.a[i] -= BASE;
      }
      return res;
    }
    return *this - (-v);
  }
  BigInt &operator += (const BigInt &v) { return *this = *this + v; }
  BigInt operator - (const BigInt &v) const {
    if (sign == v.sign) {
      if (abs() >= v.abs()) {
        BigInt res = *this; T carry = 0;
        for (int i = 0; i < int(v.a.size()) || carry; i++) {
          res.a[i] -= carry + (i < int(v.a.size()) ? v.a[i] : 0);
          carry = res.a[i] < 0; if (carry) res.a[i] += BASE;
        }
        res.trim(); return res;
      }
      return -(v - *this);
    }
    return *this + (-v);
  }
  BigInt &operator -= (const BigInt &v) { return *this = *this - v; }
  BigInt operator * (T v) const { BigInt res = *this; res *= v; return res; }
  BigInt &operator *= (T v) {
    if (v < 0) { sign = -sign; v = -v; }
    T carry = 0; for (int i = 0; i < int(a.size()) || carry; i++) {
      if (i == int(a.size())) a.push_back(0);
      T cur = a[i] * v + carry; carry = cur / BASE; a[i] = cur % BASE;
    }
    trim(); return *this;
  }
  BigInt operator * (const BigInt &v) const {
    BigInt res; res.a = mulIntFFT(a, v.a, BASE); res.sign = sign * v.sign;
    res.trim(); return res;
  }
  BigInt &operator *= (const BigInt &v) { return *this = *this * v; }
  friend pair<BigInt, BigInt> divmod(const BigInt &a1, const BigInt &b1) {
    T norm = BASE / (b1.a.back() + 1);
    BigInt a = a1.abs() * norm, b = b1.abs() * norm, q, r;
    q.a.resize(int(a.a.size()));
    for (int i = int(a.a.size()) - 1; i >= 0; i--) {
      r *= BASE; r += a.a[i];
      T s1 = int(r.a.size()) <= int(b.a.size()) ? 0 : r.a[int(b.a.size())];
      T s2 = int(r.a.size()) <= int(b.a.size()) - 1 ? 0
                                                    : r.a[int(b.a.size()) - 1];
      T d = (BASE * s1 + s2) / b.a.back(); r -= b * d;
      for (; r < 0; r += b) --d;
      q.a[i] = d;
    }
    q.sign = a1.sign * b1.sign; r.sign = a1.sign; q.trim(); r.trim();
    return make_pair(q, r / norm);
  }
  BigInt operator / (T v) const { BigInt res = *this; res /= v; return res; }
  BigInt &operator /= (T v) {
    if (v < 0) { sign = -sign; v = -v; }
    for (int i = int(a.size()) - 1, rem = 0; i >= 0; i--) {
      T cur = a[i] + rem * BASE; a[i] = cur / v; rem = cur % v;
    }
    trim(); return *this;
  }
  T operator % (T v) const {
    if (v < 0) v = -v;
    T m = 0;
    for (int i = int(a.size()) - 1; i >= 0; i--) m = (a[i] + m * BASE) % v;
    return m * sign;
  }
  BigInt operator / (const BigInt &v) const { return divmod(*this, v).first; }
  BigInt operator % (const BigInt &v) const { return divmod(*this, v).second; }
  BigInt &operator /= (const BigInt &v) { return *this = *this / v; }
  BigInt &operator %= (const BigInt &v) { return *this = *this % v; }
};

repo_name: wesley-a-leung/Resources

mime_type: text/x-c++


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

// Type templated greatest common divisor function
// Always returns a non negative integer
// Template Arguments:
//   T: the type of a and b
// Function Arguments:
//   a: the first value
//   b: the second value
// Return Value: the greatest common divisor of a and b
// In practice, has a small constant
// Time Complexity: O(log(min(a, b)))
// Memory Complexity: O(1)
// Tested:
//   https://hackerrank.com/contests/projecteuler/challenges/euler005/problem
template <class T> T gcd(T a, T b) {
  return b == 0 ? (a >= 0 ? a : -a) : gcd(b, a % b);
}

// Type templated lowest common multiple function
// Always returns a non negative integer
// Template Arguments:
//   T: the type of a and b
// Function Arguments:
//   a: the first value
//   b: the second value
// Return Value: the lowest common multiple of a and b
// In practice, has a small constant
// Time Complexity: O(log(min(a, b)))
// Memory Complexity: O(1)
// Tested:
//   https://hackerrank.com/contests/projecteuler/challenges/euler005/problem
template <class T> T lcm(T a, T b) {
  if (a == 0 || b == 0) return 0;
  if (a < 0) a = -a;
  if (b < 0) b = -b;
  return a / gcd(a, b) * b;
}

// Extended Euclidean Algorithm to compute x and y, where ax + by = gcd(a, b)
// Template Arguments:
//   T: the type of a and b
// Function Arguments:
//   a: the first value
//   b: the second value
//   x: a reference to a variable to store x
//   y: a reference to a variable to store y
// Return Value: the greatest common divisor of a and b
// In practice, has a small constant
// Time Complexity: O(log(min(a, b)))
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/modinv
//   https://www.spoj.com/problems/CEQU/
//   https://codeforces.com/problemsets/acmsguru/problem/99999/106
template <class T> T EEA(T a, T b, T &x, T &y) {
  T xx = y = 0, yy = x = 1; while (b != 0) {
    T q = a / b; a %= b; swap(a, b);
    x -= q * xx; swap(x, xx); y -= q * yy; swap(y, yy);
  }
  if (a < 0) { a = -a; x = -x; y = -y; }
  return a;
}

// Computes the multiplicative inverse of a in Zm
// Inverse only exists if gcd(a, m) == 1
// Template Arguments:
//   T: the type of a and m
// Function Arguments:
//   a: the value to find the inverse of, must be in the range [0, m)
//   m: the mod, must be positive
// Return Value: the multiplicative inverse of a in Zm, -1 if no inverse
// In practice, has a small constant
// Time Complexity: O(log m)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/modinv
template <class T> T mulInv(T a, T m) {
  if (a == 0) return -1;
  T x, y; if (EEA(a, m, x, y) != 1) return -1;
  x %= m; return x < 0 ? x + m : x;
}

// Solves the linear congruence ax = c mod m
// Template Arguments:
//   T: the type of a, c, and m
// Function Arguments:
//   a: must be in the range [0, m)
//   c: must be in the range [0, m)
//   m: the mod, must be positive
// Return Value: a pair containing x, and the modulus of the answer (equal to
//   m / gcd(a, m)), x is -1 if there is no solution
// In practice, has a small constant
// Time Complexity: O(log m)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> pair<T, T> solveCongruence(T a, T c, T m) {
  T x, y, g = EEA(a, m, x, y);
  if (c % g != 0) return make_pair(-1, m / g);
  x = (x % m + m) % m; x = (x * c / g) % (m / g); return make_pair(x, m / g);
}

// Solves the Linear Diophantine Equation ax + by = c
// All pairs of integers (s, t) where s = x.first + k * x.second
//   and t = y.first + k * y.second for all integers k are solutions
// Edge cases:
//   1. a == 0 && b == 0 is satisfied by all integer pairs (s, t) if c == 0,
//     no solutions otherwise
//   2. a == 0 is satisfied by all integer pairs (s, t) with t = c / b if c is
//     divisible by b, no solutions otherwise
//   3. b == 0 is satisfied by all integer pairs (s, t) with s = c / a if c is
//     divisible by a, no solutions otherwise
// Template Arguments:
//   T: the type of a, b, c
// Function Arguments:
//   a: the first value
//   b: the second value
//   c: the value of ax + by
//   x: a reference to a pair storing x and its mod (equal to b / gcd(a, b))
//   y: a reference to a pair storing y and its mod (equal to -a / gcd(a, b))
// Return Value: true if there is a solution, false otherwise
// In practice, has a small constant
// Time Complexity: O(log(min(a, b)))
// Memory Complexity: O(1)
// Tested:
//   https://www.spoj.com/problems/CEQU/
//   https://codeforces.com/problemsets/acmsguru/problem/99999/106
template <class T> bool LDE(T a, T b, T c, pair<T, T> &x, pair<T, T> &y) {
  assert(a != 0 && b != 0); T xg, yg, g = EEA(a, b, xg, yg);
  if (c % g != 0) return false;
  x = make_pair(xg * (c / g), b / g); y = make_pair(yg * (c / g), -a / g);
  return true;
}

// Generalized Chinese Remainder Theorem to find the solution to
//   x mod lcm(a.second, b.second) given x = a.first mod a.second
//   and x = b.first mod b.second
// Template Arguments:
//   T: the type of x
// Function Arguments:
//   a: the first value and its associated mod (0 <= a.first < a.second)
//   b: the second value and its associated mod (0 <= b.first < b.second)
// Return Value: the pair x and lcm(a.second, b.second) where
//   0 <= x < lcm(a.second, b.second), x is -1 if there is no solution
// In practice, has a small constant
// Time Complexity: O(log(min(a.second, b.second)))
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/generalchineseremainder
template <class T> pair<T, T> CRT(pair<T, T> a, pair<T, T> b) {
  T g = gcd(a.second, b.second), l = a.second / g * b.second;
  if ((b.first - a.first) % g != 0) return make_pair(-1, l);
  T A = a.second / g, B = b.second / g;
  T mul = (b.first - a.first) / g * mulInv(A % B, B) % B;
  return make_pair(((mul * a.second + a.first) % l + l) % l, l);
}


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

// Struct representing integers modulo MOD where MOD can change during runtime
// Template Arguments:
//   T: the type of the integer, must be integral
// Constructor Arguments:
//   x: the value to initialize the struct with
// Fields:
//   static MOD: the value to mod by, must be non negative
//   static PRIME_MOD: a boolean indicating whether MOD is prime
//   static MUL_OVERFLOW: a boolean indicating whether MOD * MOD overflows
//   v: the value
// Functions:
//   statis setMod(mod, primeMod, mulOverflow): sets the mod to MOD, with
//     primeMod indicating whether the mod is prime and mulOverflow indicating
//     whether multiplication overflows
//   <, <=, >, >=, ==, !=: comparison operators
//   ++, --, +, +=, -, -=, *, *=: standard arithmetic operators modulo MOD
//   pow(p): returns this value raises this to the power of p
//   hasMulInv(): returns whether a multiplicative inverse of this value exists
//   mulInv(): returns the multiplicative inverse of this value, only exists
//     if gcd(v, MOD) is 1
//   /, /=: division in modular arithmetic, only valid if a multiplicative
//     inverse of the divisor exists
//   >>, <<: input and output operators
// Time Complexity:
//   constructor: O(1)
//   setMod: O(1)
//   <, <=, >, >=, ==, !=, ++, --, +, +=, -, -=, >>, <<: O(1)
//   *, *=: O(1) if MUL_OVERFLOW is false, O(log MOD) otherwise
//   pow: O(log MOD) if MUL_OVERFLOW is false, O((log MOD)^2) otherwise
//   hasMulInv, mulInv, /, /=: O(log MOD)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/modulararithmetic
template <class T> struct DynamicIntMod {
  static_assert(is_integral<T>::value, "T must be an integral type");
  static_assert(is_signed<T>::value, "T must be a signed type");
  using IM = DynamicIntMod<T>;
  static T MOD; static bool PRIME_MOD, MUL_OVERFLOW;
  static void setMod(T mod, bool primeMod = false, bool mulOverflow = false) {
    MOD = mod; PRIME_MOD = primeMod; MUL_OVERFLOW = mulOverflow;
  }
  T v; DynamicIntMod() : v(0) {}
  DynamicIntMod(const T &x) {
    v = -MOD < x && x < MOD ? x : x % MOD; if (v < 0) v += MOD;
  }
  bool operator < (const IM &i) const { return v < i.v; }
  bool operator <= (const IM &i) const { return v <= i.v; }
  bool operator > (const IM &i) const { return v > i.v; }
  bool operator >= (const IM &i) const { return v >= i.v; }
  bool operator == (const IM &i) const { return v == i.v; }
  bool operator != (const IM &i) const { return v != i.v; }
  IM operator ++ () {
    if (++v == MOD) v = 0;
    return *this;
  }
  IM operator ++ (int) {
    IM ret = *this; if (++v == MOD) v = 0;
    return ret;
  }
  IM operator -- () {
    if (v-- == 0) v = MOD - 1;
    return *this;
  }
  IM operator -- (int) {
    IM ret = *this; if (v-- == 0) v = MOD - 1; return ret;
  }
  IM operator + (const IM &i) const { return IM(*this) += i; }
  IM &operator += (const IM &i) {
    if ((v += i.v) >= MOD) v -= MOD;
    return *this;
  }
  IM operator - (const IM &i) const { return IM(*this) -= i; }
  IM &operator -= (const IM &i) {
    if ((v -= i.v) < 0) v += MOD;
    return *this;
  }
  IM operator + () const { return *this; }
  IM operator - () const { return IM(-v); }
  IM operator * (const IM &i) const { return IM(*this) *= i; }
  IM &operator *= (const IM &i) {
    if (MUL_OVERFLOW) {
      IM x = 0, y = *this; T z = i.v;
      for (; z > 0; z /= 2, y += y) if (z % 2 == 1) x += y;
      *this = x;
    } else v = v * i.v % MOD;
    return *this;
  }
  bool hasMulInv() const {
    if (v == 0) return false;
    if (!PRIME_MOD || MUL_OVERFLOW) {
      T g = MOD, r = v; while (r != 0) { g %= r; swap(g, r); }
      return g == 1;
    } else return true;
  }
  IM mulInv() const {
    assert(v != 0);
    if (!PRIME_MOD || MUL_OVERFLOW) {
      T g = MOD, r = v, x = 0, y = 1; while (r != 0) {
        T q = g / r; g %= r; swap(g, r); x -= q * y; swap(x, y);
      }
      assert(g == 1); return IM(x);
    } else return pow2(*this, MOD - 2);
  }
  IM operator / (const IM &i) const { return IM(*this) /= i; }
  IM &operator /= (const IM &i) { return *this *= i.mulInv(); }
  friend istream &operator >> (istream &stream, IM &i) {
    T v; stream >> v; i = IM(v); return stream;
  }
  friend ostream &operator << (ostream &stream, const IM &i) {
    return stream << i.v;
  }
};

template <class T> T DynamicIntMod<T>::MOD = T(1);
template <class T> bool DynamicIntMod<T>::PRIME_MOD = false;
template <class T> bool DynamicIntMod<T>::MUL_OVERFLOW = false;


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

// Computes the Kth term of a linear recurrence of length N
//   A[i] = sum (C[j] * A[i - 1 - j]) for 0 <= j < N
// Template Arguments:
//   T: the type of the elements in the sequence
//   U: the type of K
// Function Arguments:
//   A: a vector of the first N terms of the sequence
//   C: a vector of the the N recurrence coefficients of the sequence
//   K: the term of the sequence to find
// Return Value: the Kth term (0-indexed) of the linear recurrence
// Time Complexity: O(N^2 log K)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/ddrp6
template <class T, class U>
T linearRecurrence(const vector<T> &A, const vector<T> &C, U K) {
  assert(A.size() == C.size()); int N = A.size(); if (K < N) return A[K];
  auto comb = [&] (const vector<T> &a, const vector<T> &b) {
    vector<T> res(N * 2 + 1, T());
    for (int i = 0; i <= N; i++) for (int j = 0; j <= N; j++)
      res[i + j] += a[i] * b[j];
    for (int i = N * 2; i > N; i--) for (int j = 0; j < N; j++)
      res[i - 1 - j] += res[i] * C[j];
    return res;
  };
  vector<T> P(N + 1, T()), E = P; P[0] = E[1] = T(1); K++; while (true) {
    if (K % 2 == 1) P = comb(P, E);
    if ((K /= 2) == 0) break;
    E = comb(E, E);
  }
  T res = T(); for (int i = 0; i < N; i++) res += P[i + 1] * A[i];
  return res;
}


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

// Adds two nimbers a and b
// Template Arguments:
//   T: the type of a, b, must be an unsigned integral type
// Function Arguments:
//   a: the first value
//   b: the second value
// Return Value: a + b with nimber addition
// In practice, has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   https://judge.yosupo.jp/problem/nim_product_64
template <class T> T addNim(T a, T b) {
  static_assert(is_unsigned<T>::value, "T must be an unsigned integral type");
  return a ^ b;
}

// Helper function to initialize the product cache array
template <class T, const int BITS> T initNim(T prod[BITS][BITS]) {
  for (int i = 0; i < BITS; i++) fill(prod[i], prod[i] + BITS, T());
  return T();
}

// Helper function to multiply 2^i and 2^j
template <class T, const int BITS>
T mulNimPow2(T prod[BITS][BITS], int i, int j) {
  static_assert(is_unsigned<T>::value, "T must be an unsigned integral type");
  T &res = prod[i][j]; if (res) return res;
  if (!(i & j)) return res = T(1) << (i | j);
  int a = (i & j) & -(i & j); return res = mulNimPow2(prod, i ^ a, j)
      ^ mulNimPow2(prod, (i ^ a) | (a - 1), (j ^ a) | (i & (a - 1)));
}

// Multiplies two nimbers a and b
// Template Arguments:
//   T: the type of a, b, must be an unsigned integral type
// Function Arguments:
//   a: the first value
//   b: the second value
// Return Value: a * b with nimber multiplication
// In practice, has a very small constant
// Time Complexity: O(BITS^2)
// Memory Complexity: O(1)
// Tested:
//   https://judge.yosupo.jp/problem/nim_product_64
template <class T> T mulNim(T a, T b) {
  static_assert(is_unsigned<T>::value, "T must be an unsigned integral type");
  static constexpr const int BITS = 8 * sizeof(T);
  static T prod[BITS][BITS]; static T ZERO = initNim(prod);
  T res = ZERO; for (int i = 0; i < BITS; i++) if ((a >> i) & 1)
    for (int j = 0; j < BITS; j++) if ((b >> j) & 1)
      res = addNim(res, mulNimPow2(prod, i, j));
  return res;
}

// Computes base to the power pow in nimber arithmetic
// Template Arguments:
//   T: the type of base, must be an unsigned integral type
//   U: the type of pow
// Function Arguments:
//   base: the base
//   pow: the power, must be non negative and integral
// Return Value: base to the power pow in nimber arithmetic
// In practice, has a very small constant
// Time Complexity: O(BITS^2 log pow)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T, class U> T powNim(T base, U pow) {
  static_assert(is_unsigned<T>::value, "T must be an unsigned integral type");
  T x = 1; while (true) {
    if (pow % 2 == 1) x = mulNim(x, base);
    if ((pow /= 2) == 0) break;
    base = mulNim(base, base);
  }
  return x;
}

// Computes the multiplicative inverse of a in nimber arithmetic
// Guaranteed to exist
// Template Arguments:
//   T: the type of a, must be an unsigned integral type
// Function Arguments:
//   a: the value to find the inverse of
// Return Value: the multiplicative inverse of a in nimber arithmetic
// In practice, has a small constant
// Time Complexity: O(BITS^3)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> T mulInvNim(T a) {
  static_assert(is_unsigned<T>::value, "T must be an unsigned integral type");
  return powNim(a, T(-2));
}


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

// Helper functions
template <class F, class T> T S(F &f, T a, T b) {
  return (f(a) + 4 * f((a + b) / 2) + f(b)) * (b - a) / 6;
}

template <class F, class T> T rec(F &f, T a, T b, T eps, T s) {
  T c = (a + b) / 2, S1 = S(f, a, c), S2 = S(f, c, b), SS = S1 + S2;
  if (abs(SS - s) <= 15 * eps || b - a < 1e-10) return SS + (SS - s) / 15;
  return rec(f, a, c, eps / 2, S1) + rec(f, c, b, eps / 2, S2);
}

// Integrates the function f over the range [a, b]
// Template Arguments:
//   F: the type of f
//   T: the type of the bounds and the return value of the function f
// Function Arguments:
//   a: the lower bound
//   b: the upper bound
//   f(x): the function returning the y value at x
//   eps: a value for epsilon
// Return Value: the integral of f over the range [a, b]
// Time Complexity: O(b - a) * (cost to compute f(x))
// Memory Complexity: O(1)
// Tested:
//   https://codeforces.com/problemsets/acmsguru/problem/99999/217
template <class F, class T>
T integrate(F f, T a, T b, T eps) { return rec(f, a, b, eps, S(f, a, b)); }


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

// Given a sequence A of length N, compute a sequence C of minimum length M,
//   such that A[i] = sum (C[j] * A[i - 1 - j]) for 0 <= j < M
// Template Arguments:
//   T: the type of the elements in the sequence
// Function Arguments:
//   A: the sequence to find the recurrence
// Return Value: the shortest vector C of length M, representing the
//   recurrence coefficients
// Time Complexity: O(N^2)
// Memory Complexity: O(N)
// Tested:
//   https://judge.yosupo.jp/problem/find_linear_recurrence
//   https://dmoj.ca/problem/ddrp6
template <class T> vector<T> berlekampMassey(const vector<T> &A) {
  if (A.empty()) return vector<T>();
  int N = A.size(), L = 0; vector<T> C(N, T()), B = C, tmp;
  T b = C[0] = B[0] = T(1); for (int i = 0, m = 1; i < N; i++, m++) {
    T d = A[i]; for (int j = 1; j <= L; j++) d += C[j] * A[i - j];
    if (d == T()) continue;
    tmp = C; T c = d / b; for (int j = m; j < N; j++) C[j] -= c * B[j - m];
    if (2 * L > i) continue;
    L = i + 1 - L; B = tmp; b = d; m = 0;
  }
  C.resize(L + 1); C.erase(C.begin()); for (auto &&c : C) c = -c;
  return C;
}


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

// Graphic Matroid for spanning forests
// Constructor Arguments:
//   V: the number of vertices in the simple graph
//   edges: a vector of pairs in the form (v, w) representing
//     an undirected edge in the graph between vertices v and w
// Functions:
//   clear(): clears all edges in the independent set
//   add(i): adds the ith edge to the independent set
//   independent(i): returns whether adding edge i to the current
//     existing independent set would still result in an independent set
// In practice, has a small constant
// Time Complexity:
//   constructor: O(V + E)
//   clear: O(K) for K edges added
//   add, independent: O(log V)
// Memory Complexity: O(V + E)
// Tested:
//   https://www.spoj.com/problems/COIN/
struct GraphicMatroid {
  int V; UnionFindUndo uf; vector<pair<int, int>> edges;
  GraphicMatroid(int V, const vector<pair<int, int>> &edges)
      : V(V), uf(V), edges(edges) {}
  void clear() { while (!uf.history.empty()) uf.undo(); }
  void add(int i) { uf.join(edges[i].first, edges[i].second); }
  bool independent(int i) {
    return !uf.connected(edges[i].first, edges[i].second);
  }
};

// Colorful Matroid
// Constructor Arguments:
//   color: the colors of each element in the ground set
// Functions:
//   clear(): clears all elements in the independent set
//   add(i): adds the ith element to the independent set
//   independent(i): returns whether adding element i to the current
//     existing independent set would still result in an independent set
// In practice, has a small constant
// Time Complexity:
//   constructor: O(N) for N elements in the ground set
//   clear: O(C) for C colors
//   add, independent: O(1)
// Memory Complexity: O(N + C) for N elements in the ground set and C colors
// Tested:
//   https://codeforces.com/gym/102156/problem/D
//   https://www.spoj.com/problems/COIN/
struct ColorfulMatroid {
  vector<bool> in; vector<int> color;
  ColorfulMatroid(const vector<int> &color)
      : in(*max_element(color.begin(), color.end()) + 1, false),
        color(color) {}
  void clear() { fill(in.begin(), in.end(), false); }
  void add(int i) { in[color[i]] = true; }
  bool independent(int i) { return !in[color[i]]; }
};

// Matroid for vectors in Z2
// Template Arguments:
//   BITS: the number of bits
// Constructor Arguments:
//   vec: the vectors in the ground set
// Functions:
//   clear(): clears all vector in the independent set
//   add(i): adds the ith vector to the independent set
//   independent(i): returns whether adding vector i to the current
//     existing independent set would still result in an independent set
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(N BITS) for N vectors in the ground set
//   clear: O(BITS)
//   add, independent: O(BITS^2 / 64)
// Memory Complexity: O(N BITS) for N vectors in the ground set
// Tested:
//   https://codeforces.com/gym/102156/problem/D
template <const int BITS> struct Z2Matroid {
  vector<bitset<BITS>> basis, vec;
  Z2Matroid(const vector<bitset<BITS>> &vec) : basis(BITS), vec(vec) {}
  void clear() { fill(basis.begin(), basis.end(), bitset<BITS>()); }
  void add(int i) {
    bitset<BITS> v = vec[i]; for (int j = 0; j < BITS; j++) if (v[j]) {
      if (basis[j].none()) { basis[j] = v; return; }
      v ^= basis[j];
    }
  }
  bool independent(int i) {
    bitset<BITS> v = vec[i]; for (int j = 0; j < BITS; j++) if (v[j]) {
      if (basis[j].none()) return true;
      v ^= basis[j];
    }
    return false;
  }
};

// Linear Matroid for vectors
// Template Arguments:
//   T: the type of each element
// Constructor Arguments:
//   vec: the vectors in the ground set
//   EPS: a value for epsilon
// Functions:
//   clear(): clears all vector in the independent set
//   add(i): adds the ith vector to the independent set
//   independent(i): returns whether adding vector i to the current
//     existing independent set would still result in an independent set
// In practice, has a small constant
// Time Complexity:
//   constructor: O(NV) for N vectors of size V in the ground set
//   clear: O(V) for vectors of size V
//   add, independent: O(V^2) for vectors of size V
// Memory Complexity: O(NV) for N vectors of size V in the ground set
template <class T> struct LinearMatroid {
  static T abs(T a) { return a >= 0 ? a : -a; }
  int vars; T EPS; vector<vector<T>> basis, vec;
  LinearMatroid(const vector<vector<T>> &vec, T EPS = T(1e-9))
      : vars(vec.empty() ? 0 : vec[0].size()), EPS(EPS),
        basis(vars), vec(vec) {}
  void clear() { fill(basis.begin(), basis.end(), vector<T>()); }
  void add(int i) {
    vector<T> v = vec[i]; for (int j = 0; j < vars; j++) if (abs(v[j]) > EPS) {
      if (basis[j].empty()) { basis[j] = v; return; }
      T alpha = v[j] / basis[j][j];
      for (int k = j; k < vars; k++) v[j] -= alpha * basis[j][k];
    }
  }
  bool independent(int i) {
    vector<T> v = vec[i]; for (int j = 0; j < vars; j++) if (abs(v[j]) > EPS) {
      if (basis[j].empty()) return true;
      T alpha = v[j] / basis[j][j];
      for (int k = j; k < vars; k++) v[j] -= alpha * basis[j][k];
    }
    return false;
  }
};

// Find the largest subset of a ground set of N elements that is
//   independent in both matroids
// A matroid must have the basic properties of:
//   - the empty set is independent
//   - any subset of an independent set is empty
//   - if independent set A is smaller than independent set B, then there is
//       at least 1 element of B can be added to A without loss of independency
// Elements are 0-indexed
// Template Arguments:
//   Matroid1: the type of the first matroid
//     Required Functions:
//       clear(): clears all elements in the independent set of this matroid
//       add(i): adds element i to the independent set of this matroid
//       independent(i): returns whether adding element i to the current
//         existing independent set would still result in an independent set
//   Matroid2: the type of the second matroid
//     Required Functions: same as Matroid1
// Constructor Arguments:
//   N: the number of elements in the ground set
//   m1: an instance of the first matroid
//   m2: an instance of the second matroid
// Fields:
//   N: the number of elements in the ground set
//   inSet: a vector of booleans indicating whether each element is in the
//     largest independent set or not
//   independentSet: a vector of the indices of the elements in the largest
//     independent set
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(NI sqrt I) * (time complexity of m1.add, m2.add,
//     m1.independent, and m2.independent) where I is the size of the largest
//     independent set
// Memory Complexity: O(N + M1 + M2) where M1 and M2 are the memory
//   complexities of m1 and m2
// Tested:
//   https://codeforces.com/gym/102156/problem/D
//   https://www.spoj.com/problems/COIN/
template <class Matroid1, class Matroid2> struct MatroidIntersection {
  int N; vector<bool> inSet; vector<int> independentSet;
  bool augment(Matroid1 &m1, Matroid2 &m2) {
    vector<int> par(N + 1, -1); queue<int> q; q.push(N); while (!q.empty()) {
      int v = q.front(); q.pop(); if (inSet[v]) {
        m1.clear(); for (int i = 0; i < N; i++) if (inSet[i] && i != v)
          m1.add(i);
        for (int i = 0; i < N; i++)
          if (!inSet[i] && par[i] == -1 && m1.independent(i)) {
            par[i] = v; q.push(i);
          }
      } else {
        auto backE = [&] {
          m2.clear(); for (int c = 0; c < 2; c++) for (int i = 0; i < N; i++)
            if ((v == i || inSet[i]) && (par[i] == -1) == c) {
              if (!m2.independent(i)) {
                if (c) { par[i] = v; q.push(i); return i; }
                else return -1;
              }
              m2.add(i);
            }
          return N;
        };
        for (int w = backE(); w != -1; w = backE()) if (w == N) {
          for (; v != N; v = par[v]) inSet[v] = !inSet[v];
          return true;
        }
      }
    }
    return false;
  }
  MatroidIntersection(int N, Matroid1 m1, Matroid2 m2)
      : N(N), inSet(N + 1, false) {
    m1.clear(); m2.clear(); inSet[N] = true; for (int i = N - 1; i >= 0; i--)
      if (m1.independent(i) && m2.independent(i)) {
        inSet[i] = true; m1.add(i); m2.add(i);
      }
    while (augment(m1, m2));
    inSet.pop_back();
    for (int i = 0; i < N; i++) if (inSet[i]) independentSet.push_back(i);
  }
};


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

// Struct representing integers modulo MOD
// Template Arguments:
//   T: the type of the integer, must be integral
//   MOD: the value to mod by, must be non negative
//   PRIME_MOD: a boolean indicating whether MOD is prime, can be computed
//     automatically if C++ 14 or above, and if MOD is small
//   MUL_OVERFLOW: a boolean indicating whether MOD * MOD overflows, can be
//     computed automatically if numeric_limits<T>::max() exists
// Constructor Arguments:
//   x: the value to initialize the struct with
// Fields:
//   v: the value
// Functions:
//   <, <=, >, >=, ==, !=: comparison operators
//   ++, --, +, +=, -, -=, *, *=: standard arithmetic operators modulo MOD
//   pow(p): returns this value raises this to the power of p
//   hasMulInv(): returns whether a multiplicative inverse of this value exists
//   mulInv(): returns the multiplicative inverse of this value, only exists
//     if gcd(v, MOD) is 1
//   /, /=: division in modular arithmetic, only valid if a multiplicative
//     inverse of the divisor exists
//   >>, <<: input and output operators
// Time Complexity:
//   constructor: O(1)
//   <, <=, >, >=, ==, !=, ++, --, +, +=, -, -=, >>, <<: O(1)
//   *, *=: O(1) if MUL_OVERFLOW is false, O(log MOD) otherwise
//   pow: O(log MOD) if MUL_OVERFLOW is false, O((log MOD)^2) otherwise
//   hasMulInv, mulInv, /, /=: O(log MOD)
// Memory Complexity: O(1)
// Tested:
//   https://dmoj.ca/problem/angieandfunctions
template <class T, const T MOD,
#if __cplusplus < 201402L
          const bool PRIME_MOD,
#else
          const bool PRIME_MOD = isPrime(MOD),
#endif
          const bool MUL_OVERFLOW = (numeric_limits<T>::max() / MOD < MOD)>
    struct IntMod {
  static_assert(is_integral<T>::value, "T must be an integral type");
  static_assert(is_signed<T>::value, "T must be a signed type");
  static_assert(0 < MOD, "MOD must be a positive integer");
  using IM = IntMod<T, MOD, PRIME_MOD, MUL_OVERFLOW>;
  T v; IntMod() : v(0) {}
  IntMod(const T &x) {
    v = -MOD < x && x < MOD ? x : x % MOD; if (v < 0) v += MOD;
  }
  bool operator < (const IM &i) const { return v < i.v; }
  bool operator <= (const IM &i) const { return v <= i.v; }
  bool operator > (const IM &i) const { return v > i.v; }
  bool operator >= (const IM &i) const { return v >= i.v; }
  bool operator == (const IM &i) const { return v == i.v; }
  bool operator != (const IM &i) const { return v != i.v; }
  IM operator ++ () {
    if (++v == MOD) v = 0;
    return *this;
  }
  IM operator ++ (int) {
    IM ret = *this; if (++v == MOD) v = 0;
    return ret;
  }
  IM operator -- () {
    if (v-- == 0) v = MOD - 1;
    return *this;
  }
  IM operator -- (int) {
    IM ret = *this; if (v-- == 0) v = MOD - 1; return ret;
  }
  IM operator + (const IM &i) const { return IM(*this) += i; }
  IM &operator += (const IM &i) {
    if ((v += i.v) >= MOD) v -= MOD;
    return *this;
  }
  IM operator - (const IM &i) const { return IM(*this) -= i; }
  IM &operator -= (const IM &i) {
    if ((v -= i.v) < 0) v += MOD;
    return *this;
  }
  IM operator + () const { return *this; }
  IM operator - () const { return IM(-v); }
  IM operator * (const IM &i) const { return IM(*this) *= i; }
  IM &operator *= (const IM &i) {
    if (MUL_OVERFLOW) {
      IM x = 0, y = *this; T z = i.v;
      for (; z > 0; z /= 2, y += y) if (z % 2 == 1) x += y;
      *this = x;
    } else v = v * i.v % MOD;
    return *this;
  }
  bool hasMulInv() const {
    if (v == 0) return false;
    if (!PRIME_MOD || MUL_OVERFLOW) {
      T g = MOD, r = v; while (r != 0) { g %= r; swap(g, r); }
      return g == 1;
    } else return true;
  }
  IM mulInv() const {
    if (!PRIME_MOD || MUL_OVERFLOW) {
      T g = MOD, r = v, x = 0, y = 1; while (r != 0) {
        T q = g / r; g %= r; swap(g, r); x -= q * y; swap(x, y);
      }
      assert(g == 1); return IM(x);
    } else return pow2(*this, MOD - 2);
  }
  IM operator / (const IM &i) const { return IM(*this) /= i; }
  IM &operator /= (const IM &i) { return *this *= i.mulInv(); }
  friend istream &operator >> (istream &stream, IM &i) {
    T v; stream >> v; i = IM(v); return stream;
  }
  friend ostream &operator << (ostream &stream, const IM &i) {
    return stream << i.v;
  }
};


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

// Finds a candidate for the majority element with constant space
// A second pass is required to ensure the candidate is actually a
//   majority element
// Template Arguments:
//   F: the type of the function generating the elements
// Function Arguments:
//   N: the number of elements
//   f: a function that returns the ith element on the ith call
// Return Value: returns an element that is a candidate for the
//   majority element
// In practice, has a very small constant
// Time Complexity: O(N)
// Memory Complexity: O(1)
// Tested:
//   https://open.kattis.com/problems/farmingmars
template <class F>
auto boyerMooreMajority(int N, F f) -> typename decay<decltype(f())>::type {
  typename decay<decltype(f())>::type ret = f();
  for (int cnt = 1, i = 1; i < N; i++) {
    auto v = f(); if (cnt == 0) { ret = v; cnt++; }
    else if (v == ret) cnt++;
    else cnt--;
  }
  return ret;
}


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

// Helper function to find a primitive root for the mod p
template <class T> T primitiveRoot(T p) {
  T phi = p - 1, n = phi; vector<T> pf; for (T i = 2; i * i <= n; i++) {
    if (n % i == 0) pf.push_back(i);
    while (n % i == 0) n /= i;
  }
  for (T res = 2; res <= p; res++) {
    bool ok = true; for (int i = 0; i < int(pf.size()) && ok; i++)
      ok &= powMod(res, phi / pf[i], p) != 1;
    if (ok) return res;
  }
  return -1;
}

// Number Theoretic Transform
// Template Arguments:
//   T: the type of the integer, must be integral
//   MOD: the value to mod by, must be in the form C * 2^K + 1 for some
//       integers C, K where size(a) <= 2^K
//     Sample Mods:
//       998244353 = 119 * 2^23 + 1, primitiveRoot = 3
//       754974721 = 45 * 2^24 + 1, primitiveRoot = 11
//       167772161 = 5 * 2^25 + 1, primitiveRoot = 3
//       469762049 = 7 * 2^26 + 1, primitiveRoot = 3
//       1004535809 = 479 * 2^21 + 1, primitiveRoot = 3
//       1012924417 = 483 * 2^21 + 1, primitiveRoot = 5
// Function Arguments:
//   a: a reference to the vector of type T in the range [0, MOD)
// In practice, has a moderate constant, faster than FFT and Karatsuba
// Time Complexity: O(N log N) where N = size(a)
// Memory Complexity: O(N)
// Tested:
//   https://judge.yosupo.jp/problem/convolution_mod
//   https://judge.yosupo.jp/problem/convolution_mod_1000000007
template <class T, const T MOD> void ntt(vector<T> &a) {
  static_assert(is_integral<T>::value, "T must be an integral type");
  static_assert(is_signed<T>::value, "T must be a signed type");
#if __cplusplus >= 201402L
  static_assert(isPrime(MOD), "MOD must be prime");
#endif
  static T pk = 1, PK = 0;
  for (; pk < MOD; pk *= 2) if ((MOD - 1) % pk == 0) PK = pk;
  assert(PK != 0 && "MOD must be of the form C * 2^K + 1");
  static T ROOT = powMod(primitiveRoot(MOD), (MOD - 1) / PK, MOD);
  static vector<T> rt(2, 1); static vector<int> ord; static int k = 2, len = 1;
  int N = a.size(); assert(N <= PK); assert(!(N & (N - 1)));
  for (; k < N; k <<= 1, len++) {
    rt.resize(N, T()); T x = powMod(ROOT, PK >> (len + 1), MOD);
    for (int i = k; i < (k << 1); i++)
      rt[i] = i & 1 ? mulMod(rt[i >> 1], x, MOD) : rt[i >> 1];
  }
  if (int(ord.size()) != N) {
    ord.assign(N, 0); int len = __builtin_ctz(N); for (int i = 0; i < N; i++)
      ord[i] = (ord[i >> 1] >> 1) + ((i & 1) << (len - 1));
  }
  for (int i = 0; i < N; i++) if (i < ord[i]) swap(a[i], a[ord[i]]);
  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 = mulMod(a[len + i + j], rt[len + j], MOD);
      a[i + j] = addMod(u, v, MOD); a[len + i + j] = subMod(u, v, MOD);
    }
}

// Polynomial Multiplication
// Template Arguments:
//   T: the type of the integer, must be integral
//   MOD: the value to mod by, must be in the form C * 2^K + 1 for some
//     integers C, K where size(a) + size(b) - 1 <= 2^K
// Functions Arguments:
//   a: a vector of type T representing the first polynomial, a[i] stores
//     the coefficient of x^i, must be in the range [0, MOD)
//   b: a vector of type T representing the second polynomial, b[i] stores
//     the coefficient of x^i, must be in the range [0, MOD)
//   eq: a boolean indicating whether a and b are equal (to reduce the number
//     of ntt calls)
// Return Value: a vector of type T representing the polynomial a times b with
//   no trailing zeros, all value in the range [0, MOD)
// In practice, has a moderate constant
// Time Complexity: O(N log N) where N = size(a) + size(b)
// Memory Complexity: O(N)
// Tested:
//   https://judge.yosupo.jp/problem/convolution_mod
//   https://judge.yosupo.jp/problem/convolution_mod_1000000007
const int NTT_CUTOFF = 30000;
template <class T, const T MOD>
vector<T> mulPolyNTT(const vector<T> &a, const vector<T> &b, bool eq = false) {
  int N = max(0, int(a.size()) + int(b.size()) - 1);
  if ((long long)(a.size()) * (long long)(b.size()) <= NTT_CUTOFF) {
    vector<T> res(N, T()); for (int i = 0; i < int(a.size()); i++)
      for (int j = 0; j < int(b.size()); j++)
        res[i + j] = addMod(res[i + j], mulMod(a[i], b[j], MOD), MOD);
    while (int(res.size()) > 1 && res.back() == T()) res.pop_back();
    return res;
  }
  while (N & (N - 1)) N++;
  vector<T> fa(N, T()), fb, res(N, T());
  copy(a.begin(), a.end(), fa.begin()); ntt<T, MOD>(fa); if (eq) fb = fa;
  else {
    fb.assign(N, T()); copy(b.begin(), b.end(), fb.begin()); ntt<T, MOD>(fb);
  }
  T invN = mulInvModPrime(T(N), MOD);
  res[0] = mulMod(mulMod(fa[0], fb[0], MOD), invN, MOD);
  for (int i = 1; i < N; i++)
    res[N - i] = mulMod(mulMod(fa[i], fb[i], MOD), invN, MOD);
  ntt<T, MOD>(res);
  while (int(res.size()) > 1 && res.back() == T()) res.pop_back();
  return res;
}


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

// constexpr function (in C++14 and above) to determine if x is prime
// Template Arguments:
//   T: the type of x
// Function Arguments:
//   x: the value to check if prime
// Return Value: true if x is prime, false otherwise
// In practice, has a small constant
// Time Complexity: O(sqrt x)
// Memory Complexity: O(1)
// Tested:
//   https://www.spoj.com/SHORTEN/problems/ISPRIME/
//   https://dmoj.ca/problem/angieandfunctions
template <class T> constexpr bool isPrime(T x) {
  if (x < 2) return false;
  for (T i = 2; i * i <= x; i++) if (x % i == 0) return false;
  return true;
}

// Finds the prime factors of x
// Template Arguments:
//   T: the type of x
// Function Arguments:
//   x: the value to prime factor
// Return Value: a sorted vector of T with the prime factorization of x
// In practice, has a small constant
// Time Complexity: O(sqrt x)
// Memory Complexity: O(log x)
// Tested:
//   https://www.spoj.com/problems/FACT0/
template <class T> vector<T> primeFactor(T x) {
  vector<T> ret; for (T i = 2; i * i <= x; i++)
    while (x % i == 0) { ret.push_back(i); x /= i; }
  if (x > 1) ret.push_back(x);
  return ret;
}

// Finds the prime factors of x and its exponent
// Template Arguments:
//   T: the type of x
// Function Arguments:
//   x: the value to prime factor
// Return Value: a sorted vector of pairs of T and int with each
//   prime factor of x and its exponent
// In practice, has a small constant
// Time Complexity: O(sqrt x)
// Memory Complexity: O(log x)
// Tested:
//   https://www.spoj.com/problems/FACT0/
template <class T> vector<pair<T, int>> primeFactorWithCount(T x) {
  vector<pair<T, int>> ret; for (T i = 2; i * i <= x; i++) if (x % i == 0)
    for (ret.emplace_back(i, 0); x % i == 0; x /= i) ret.back().second++;
  if (x > 1) ret.emplace_back(x, 1);
  return ret;
}

// Computes phi(x) which is the number of positive integers less than or equal
//   to x that are relatively prime to x
// Template Arguments:
//   T: the type of x
// Function Arguments:
//   x: the value to find phi(x)
// Return Value: phi(x)
// In practice, has a small constant
// Time Complexity: O(sqrt x)
// Memory Complexity: O(1)
// Tested:
//   https://www.spoj.com/problems/ETF/
template <class T> T phi(T x) {
  T ret = x; for (T i = 2; i * i <= x; i++) if (x % i == 0)
    for (ret -= ret / i; x % i == 0; x /= i);
  if (x > 1) ret -= ret / x;
  return ret;
}

// Computes mobius(x), which is 0 if x contains a squared prime factor,
//   -1 if x has an odd number of prime factors with no squares,
//   1 if x has an even number of prime factors with no squares
// Template Arguments:
//   T: the type of x
// Function Arguments:
//   x: the value to find mobius(x)
// Return Value: the value of mobius(x)
// In practice, has a small constant
// Time Complexity: O(sqrt x)
// Memory Complexity: O(1)
template <class T> int mobius(T x) {
  int ret = 1; for (T i = 2; i * i <= x; i++)
    if (x % i == 0) { x /= i; ret *= -1; if (x % i == 0) return 0; }
  if (x > 1) ret *= -1;
  return ret;
}

// Finds the factors of x
// Template Arguments:
//   T: the type of x
// Function Arguments:
//   x: the value to find the factors
// Return Value: a vector of T with the factors of x,
//   not guaranteed to be sorted
// In practice, has a small constant
// Time Complexity: O(sqrt x)
// Memory Complexity: O(sqrt x)
// Tested:
//   https://www.spoj.com/problems/CZ_PROB2/
template <class T> vector<T> factor(T x) {
  vector<T> ret; for (T i = 1; i * i <= x; i++)
    if (x % i == 0) { ret.push_back(i); if (x / i != i) ret.push_back(x / i); }
  return ret;
}

// Determine if x is prime using the Miller Rabin Primality Test
// Template Arguments:
//   T: the type of x
// Function Arguments:
//   x: the value to check if prime, mulMod(x, x, x) should not overflow
//   iterations: number of iterations to run, should be at least 7
// Return Value: true if x is prime, false otherwise, can return
//   false positives for x > 7e18, otherwise is guaranteed to be correct
// In practice, has a small constant
// Time Complexity: O(time complexity of powMod)
// Memory Complexity: O(iterations)
// Tested:
//   https://dmoj.ca/problem/bf3hard
template <class T> bool millerRabin(T x, int iterations = 7) {
  if (x < 2 || x % 6 % 4 != 1) return (x | 1) == 3;
  vector<T> A = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
  while (int(A.size()) < iterations)
    A.push_back(uniform_int_distribution<long long>(
        1795265023, numeric_limits<long long>::max())(rng64));
  int s = 0; while (!(((x - 1) >> s) & 1)) s++;
  T d = x >> s; for (T a : A) {
    T p = powMod(a % x, d, x); int i = s;
    while (p != 1 && p != x - 1 && a % x != 0 && i--) p = mulMod(p, p, x);
    if (p != x - 1 && i != s) return false;
  }
  return true;
}

// Sieve of Erathosthenes to identify primes less than or equal to N
// Constructor Arguments:
//   N: the maximum value
//   f(i): the function to run a callback on for each prime i
// Fields:
//   p: a vector of booleans with the ith indicating whether the integer
//     i * 2 + 1 is prime or not
// Functions:
//   isPrime(x): returns whether x is prime or not
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(N log log N)
//   isPrime: O(1)
// Memory Complexity: O(N / 128)
// Tested:
//   https://judge.yosupo.jp/problem/enumerate_primes
//   https://codeforces.com/contest/1149/problem/A
struct Sieve {
  vector<bool> p;
  template <class F = function<void(int)>>
  Sieve(int N, F f = [] (int) {}) : p(N / 2 + 1, true) {
    p[0] = false; if (N >= 2) f(2);
    for (int i = 3; i <= N; i += 2) if (p[i / 2]) {
      f(i); if (i <= (N + i - 1) / i) for (int j = i * i; j <= N; j += i * 2)
        p[j / 2] = false;
    }
  }
  bool isPrime(int x) { return x == 2 || (x % 2 == 1 && p[x / 2]); }
};

// Sieve of Erathosthenes to identify primes, the smallest prime factor of
//   each number from 1 to N, the phi function, and the mobius function
// Constructor Arguments:
//   N: the maximum value
// Fields:
//   isPrime: a vector of N + 1 booleans indicating whether each integer is
//     prime or not
//   primes: a vector of all primes less than or equal to N
//   SPF: a vector of N + 1 integers representing the smallest prime factor
//     less than or equal to each integer
//   phi: a vector of integers with phi[i] equal to the number of positive
//     integers less than or equal to i that are relatively prime to i
//   mobius: a vector of chars with mobius[i] equal to 0 if i contains a
//     squared prime factor, -1 if it has an odd number of prime factors with
//     no squares, 1 if it has an even number of prime factors with no squares
// Functions:
//   primeFactor(x): returns a sorted vector of integers with the
//     prime factorization of x
//   primeFactorWithCount(x): returns a sorted vector of pairs of integers
//     with each prime factor of x and its exponent
// In practice, has a small constant
// Time Complexity:
//   constructor: O(N)
//   primeFactor, primeFactorWithCount: O(log x)
// Memory Complexity: O(N)
// Tested:
//   https://tlx.toki.id/problems/troc-16/E
//   https://dmoj.ca/problem/nccc7j5
struct SievePrimeFactoring {
  vector<bool> isPrime; vector<int> primes, SPF, phi; vector<char> mobius;
  SievePrimeFactoring(int N)
      : isPrime(N + 1, true), SPF(N + 1, 0), phi(N + 1, 0), mobius(N + 1, 1) {
    isPrime[0] = isPrime[1] = false; phi[1] = 1; for (int i = 2; i <= N; i++) {
      if (isPrime[i]) {
        primes.push_back(SPF[i] = i); phi[i] = i - 1; mobius[i] = -1;
      }
      for (int p : primes) {
        if (i * p > N) break;
        isPrime[i * p] = false; SPF[i * p] = p;
        if (i % p == 0) { phi[i * p] = phi[i] * p; mobius[i * p] = 0; break; }
        phi[i * p] = phi[i] * phi[p]; mobius[i * p] = mobius[i] * -1;
      }
    }
  }
  vector<int> primeFactor(int x) {
    vector<int> ret; while (x > 1) { ret.push_back(SPF[x]); x /= SPF[x]; }
    return ret;
  }
  vector<pair<int, int>> primeFactorWithCount(int x) {
    vector<pair<int, int>> ret; while (x > 1) {
      ret.emplace_back(SPF[x], 0);
      for (int spf = SPF[x]; x % spf == 0; x /= spf) ret.back().second++;
    }
    return ret;
  }
};

// Segmented Sieve of Eratosthenes to identify primes between lo and hi
// Template Arguments:
//   T: the type of lo and hi
//   F: the type of f
// Functions Arguments:
//   lo: the inclusive lower bound
//   hi: the inclusive upper bound
//   f(i): the function to run a callback on for each prime i
// In practice, has a small constant
// Time Complexity: O(sqrt(hi) log(log(hi)) + (hi - lo))
// Memory Complexity: O(sqrt(hi) + (hi - lo))
// Tested:
//   https://dmoj.ca/problem/phantom3
template <class T, class F> void segmentedSieve(T lo, T hi, F f) {
  lo = max(lo, T(2)); T sqrtHi = sqrtl(hi) + 2;
  while (sqrtHi * sqrtHi > hi) sqrtHi--;
  vector<bool> p1(sqrtHi + 1, false), p2(hi - lo + 1, false);
  p1[0] = p1[1] = true; for (T i = 2; i <= sqrtHi; i++) if (!p1[i]) {
    for (T j = i * i; j <= sqrtHi; j += i) p1[j] = true;
    for (T j = (lo + i - 1) / i * i; j <= hi; j += i)
      if (j != i && !p2[j - lo]) p2[j - lo] = true;
  }
  for (T i = 0; i < hi - lo + 1; i++) if (!p2[i]) f(lo + i);
}

// Returns an arbitrary divisor of non prime N
// Template Arguments:
//   T: the type of N
// Function Arguments:
//   N: the non prime value to find a divisor, mulMod(N, N, N) should not
//     overflow
//   iterations: number of iterations to run, should be at least 40
// Return Value: an arbitrary divisor of N
// In practice, has a moderate constant
// Time Complexity: O(N^(1/4))
// Memory Complexity: O(1)
// Tested:
//   https://loj.ac/p/6466
//   https://www.spoj.com/problems/FACT2/
//   https://judge.yosupo.jp/problem/factorize
template <class T> T pollardsRho(T N, int iterations = 40) {
  if (N == 1) return 1;
  auto f = [&] (T x) { return mulMod(x, x, N) + 1; };
  T x = 0, y = 0, p = 2, q; int t = 0, i = 1;
  while (t++ % iterations != 0 || gcd(p, N) == 1) {
    if (x == y) y = f(x = ++i);
    if ((q = mulMod(p, max(x, y) - min(x, y), N)) != 0) p = q;
    x = f(x); y = f(f(y));
  }
  return gcd(p, N);
}

// Prime factors a large integer x
// Template Arguments:
//   T: the type of N
// Function Arguments:
//   x: the value to prime factor, mulMod(x, x, x) should not overflow
//   pollardsRhoIters: number of iterations to run pollard rho, should be
//     at least 40
//   millerRabinIters: number of iterations to run miller rabin, should be
//     at least 7
// Return Value: a sorted vector of T with the prime factorization of x
// In practice, has a moderate constant
// Time Complexity: O(x^(1/4) log(x))
// Memory Complexity: O(log x)
// Tested:
//   https://www.spoj.com/problems/FACT2/
//   https://judge.yosupo.jp/problem/factorize
template <class T>
vector<T> pollardsRhoPrimeFactor(T x, int pollardsRhoIters = 40,
                                 int millerRabinIters = 7) {
  if (x == 1) return vector<T>();
  vector<T> ret; queue<T> q; q.push(x); while (!q.empty()) {
    T y = q.front(); q.pop();
    if (millerRabin(y, millerRabinIters)) ret.push_back(y);
    else { q.push(pollardsRho(y, pollardsRhoIters)); q.push(y / q.back()); }
  }
  return ret;
}

// Prime factors a large integer x and its exponent
// Template Arguments:
//   T: the type of N
// Function Arguments:
//   x: the value to prime factor, mulMod(x, x, x) should not overflow
//   pollardsRhoIters: number of iterations to run pollard rho, should be
//     at least 40
//   millerRabinIters: number of iterations to run miller rabin, should be
//     at least 7
// Return Value: a sorted vector of pairs of T and int with each
//   prime factor of x and its exponent
// In practice, has a moderate constant
// Time Complexity: O(x^(1/4) log(x))
// Memory Complexity: O(log x)
// Tested:
//   https://www.spoj.com/problems/FACT2/
template <class T>
vector<pair<T, int>> pollardsRhoPrimeFactorWithCount(
    T x, int pollardsRhoIters = 40, int millerRabinIters = 7) {
  vector<T> pf = pollardsRhoPrimeFactor(x, pollardsRhoIters, millerRabinIters);
  sort(pf.begin(), pf.end()); vector<pair<T, int>> ret;
  for (int i = 0, cnt = 0; i < int(pf.size()); i++) {
    cnt++; if (i + 1 == int(pf.size()) || pf[i] != pf[i + 1]) {
      ret.emplace_back(pf[i], cnt); cnt = 0;
    }
  }
  return ret;
}

// Determines the factors of all numbers from 1 to N
// Constructor Arguments:
//   N: the maximum value
// Fields:
//   factors: a vector of vectors of integers with factors[i] being the sorted
//     factors of i
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(N log N)
// Memory Complexity: O(N log N)
// Tested:
//   https://codeforces.com/contest/111/problem/B
struct Factors {
  vector<vector<int>> factors;
  Factors(int N) : factors(N + 1) {
    for (int i = 1; i <= N; i++) for (int j = i; j <= N; j += i)
      factors[j].push_back(i);
  }
};


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

// Computes base to the power of pow
// Template Arguments:
//   T: the type of base
//   U: the type of pow
// Function Arguments:
//   base: the base
//   pow: the power, must be non negative and integral
// Return Value: returns base to the power of pow
// In practice, has a small constant
// Time Complexity: O(log pow)
// Memory Complexity: O(1)
// Tested:
//   https://cses.fi/problemset/task/1095
template <class T, class U> T pow2(T base, U pow) {
  T x = 1; while (true) {
    if (pow % 2 == 1) x = x * base;
    if ((pow /= 2) == 0) break;
    base = base * base;
  }
  return x;
}


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

// Solve the equation Ax = b for a matrix A of size M x N, and vectors b of
//   size M, and vector x of size N
// Template Arguments:
//   T: the type of each element
// Constructor Arguments:
//   A: a M x N matrix
//   b: a vector, with dimension M
//   EPS: a value for epsilon
// Fields:
//   NO_SOLUTION: static const int representing no solution
//   ONE_SOLUTION: static const int representing one solution
//   INF_SOLUTIONS: static const int representing infinite solution
//   M: the number of equations
//   N: the number of variables
//   solutions: the number of solution to the equation, equal to one of
//     NO_SOLUTION, ONE_SOLUTION, or INF_SOLUTIONS
//   EPS: the value for epsilon
//   x: the solution vector with dimension N, can be any solution if there are
//     infinite solutions, undefined if no solutions
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(min(M, N) MN)
// Memory Complexity: O(MN)
// Tested:
//   https://open.kattis.com/problems/equationsolver
template <class T> struct GaussianElimination {
  static const int NO_SOLUTION = 0, ONE_SOLUTION = 1, INF_SOLUTIONS = 2;
  T abs(T a) { return a >= 0 ? a : -a; }
  int M, N, solutions; T EPS; vector<T> x;
  GaussianElimination(vector<vector<T>> A, vector<T> b, T EPS = T(1e-9))
      : M(A.size()), N(M == 0 ? 0 : A[0].size()), solutions(ONE_SOLUTION),
        EPS(EPS), x(N, T()) {
    vector<int> ind(N, -1); for (int r = 0, c = 0; c < N && r < M; c++) {
      int mx = r; for (int i = r + 1; i < M; i++)
        if (abs(A[i][c]) > abs(A[mx][c])) mx = i;
      if (abs(A[mx][c]) <= EPS) continue;
      if (r != mx) { A[r].swap(A[mx]); swap(b[r], b[mx]); }
      ind[c] = r; T inv = T(1) / A[r][c];
      for (int i = 0; i < M; i++) if (i != r) {
        T alpha = A[i][c] * inv; b[i] -= b[r] * alpha;
        for (int j = c; j < N; j++) A[i][j] -= A[r][j] * alpha;
      }
      r++;
    }
    for (int i = 0; i < N; i++) if (ind[i] != -1)
      x[i] = b[ind[i]] / A[ind[i]][i];
    for (int i = 0; i < M; i++) {
      T sm = T(); for (int j = 0; j < N; j++) sm += x[j] * A[i][j];
      if (abs(sm - b[i]) > EPS) {
        solutions = NO_SOLUTION; x.clear(); return;
      }
    }
    for (int i = 0; i < N; i++) if (ind[i] == -1) solutions = INF_SOLUTIONS;
  }
};


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

// Solve the xor boolean satisfiability problem for a boolean matrix A of
//   size M x N, and vectors b of size M, and vector x of size N, so that
//   b_i = xor(A[i][j] & x[j]) for 0 <= j < N; equivalent to solving the
//   equation Ax = b mod 2 for a matrix A, and vectors b, x where all values
//   are 0 or 1
// Template Arguments:
//   MAXN: the maximmum value of N
// Constructor Arguments:
//   A: a M x N boolean matrix
//   b: a boolean vector, with dimension M
// Fields:
//   NO_SOLUTION: static const int representing no solution
//   ONE_SOLUTION: static const int representing one solution
//   INF_SOLUTIONS: static const int representing infinite solution
//   M: the number of equations
//   N: the number of variables
//   solutions: the number of solution to the equation, equal to one of
//     NO_SOLUTION, ONE_SOLUTION, or INF_SOLUTIONS
//   x: the solution vector with dimension N, can be any solution if there are
//     infinite solutions, undefined if no solutions
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(min(M, N) M MAXN / 64)
// Memory Complexity: O(M MAXN / 64)
// Tested:
//   https://dmoj.ca/problem/tle16c1p6
//   https://dmoj.ca/problem/dmopc20c4p6
template <const int MAXN> struct XorSat {
  static const int NO_SOLUTION = 0, ONE_SOLUTION = 1, INF_SOLUTIONS = 2;
  int M, N, solutions; bitset<MAXN> x;
  XorSat(vector<bitset<MAXN>> A, vector<bool> b)
      : M(A.size()), N(M == 0 ? 0 : A[0].size()), solutions(ONE_SOLUTION) {
    vector<int> ind(N, -1); for (int r = 0, c = 0; c < N && r < M; c++) {
      int mx = r; for (int i = r + 1; i < M; i++)
        if (A[i][c]) { mx = i; break; }
      if (!A[mx][c]) continue;
      if (r != mx) { swap(A[r], A[mx]); swap(b[r], b[mx]); }
      ind[c] = r; for (int i = 0; i < M; i++)
        if (i != r && A[i][c]) { b[i] = b[i] ^ b[r]; A[i] ^= A[r]; }
      r++;
    }
    for (int i = 0; i < N; i++) if (ind[i] != -1) x[i] = b[ind[i]];
    for (int i = 0; i < M; i++) {
      bool val = false; for (int j = 0; j < N; j++) val ^= x[j] & A[i][j];
      if (val != b[i]) { solutions = NO_SOLUTION; x.reset(); return; }
    }
    for (int i = 0; i < N; i++) if (ind[i] == -1) solutions = INF_SOLUTIONS;
  }
};


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

// Computes the basis of vectors in base 2
// Template Arguments:
//   BITS: the number of bits
// Fields:
//   basis: a vector of bitsets of length BITS, with the ith element being the
//     reduced vector with the most significant bit of i, or empty if there is
//     no such vector
//   ids: a vector of integers representing the id of the ith vector in basis,
//     or -1 if basis[i] is empty
// Functions:
//   addVector(v, id): adds the vector v to the set with and id of id
//   inSpan(v): returns whether a vector v is in the span of the basis
// In practice, has a small constant
// Time Complexity:
//   constructor, addVector, inSpan: O(BITS^2 / 64)
// Memory Complexity: O(BITS^2 / 64)
// Tested:
//   https://dmoj.ca/problem/dmopc19c5p5
template <const int BITS> struct XorBasis {
  vector<bitset<BITS>> basis; vector<int> ids;
  XorBasis() : basis(BITS), ids(BITS, -1) {}
  void addVector(bitset<BITS> v, int id) {
    for (int i = 0; i < BITS; i++) if (v[i]) {
      if (basis[i].none()) { basis[i] = v; ids[i] = id; return; }
      v ^= basis[i];
    }
  }
  bool inSpan(bitset<BITS> v) {
    for (int i = 0; i < BITS; i++) if (v[i]) v ^= basis[i];
    return v.none();
  }
};


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

// Fast Fourier Transform
// Template Arguments:
//   F: the floating point type of the complex number components
// Function Arguments:
//   a: a reference to the vector of complex numbers (of type complex<F>)
//     to convolute
// In practice, has a moderate constant, faster than using Karatsuba, slower
//   than NTT
// Time Complexity: O(N log N) where N = size(a)
// Memory Complexity: O(N)
// Tested:
//   https://open.kattis.com/problems/polymul2
//   https://dmoj.ca/problem/atimesb
template <class F> void fft(vector<complex<F>> &a) {
  static_assert(is_floating_point<F>::value,
                "F must be a floating point type");
  static vector<complex<F>> rt(2, 1); static vector<int> ord; static int k = 2;
  static F PI = acos(F(-1)); int N = a.size(); assert(!(N & (N - 1)));
  for (; k < N; k <<= 1) {
    rt.resize(N); complex<F> x = polar(F(1), PI / k);
    for (int i = k; i < (k << 1); i++)
      rt[i] = i & 1 ? rt[i >> 1] * x : rt[i >> 1];
  }
  if (int(ord.size()) != N) {
    ord.assign(N, 0); int len = __builtin_ctz(N); for (int i = 0; i < N; i++)
      ord[i] = (ord[i >> 1] >> 1) + ((i & 1) << (len - 1));
  }
  for (int i = 0; i < N; i++) if (i < ord[i]) swap(a[i], a[ord[i]]);
  for (int len = 1; len < N; len <<= 1) for (int i = 0; i < N; i += len << 1)
    for (int j = 0; j < len; j++) {
      complex<F> u = a[i + j], x = a[len + i + j], y = rt[len + j];
      complex<F> v(real(x) * real(y) - imag(x) * imag(y),
                   real(x) * imag(y) + imag(x) * real(y));
      a[i + j] = u + v; a[len + i + j] = u - v;
    }
}

// Polynomial Multiplication
// Template Arguments:
//   T: the type of each element
//   F: the floating point type to use for FFT
// 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 moderate constant
// Time Complexity: O(N log N) where N = size(a) + size(b)
// Memory Complexity: O(N)
// Tested:
//   https://open.kattis.com/problems/polymul2
const int FFT_CUTOFF = 30000;
template <class T, class F = long double>
vector<T> mulPolyFFT(const vector<T> &a, const vector<T> &b) {
  int N = max(0, int(a.size()) + int(b.size()) - 1);
  if ((long long)(a.size()) * (long long)(b.size()) <= FFT_CUTOFF) {
    vector<T> res(N, T()); for (int i = 0; i < int(a.size()); i++)
      for (int j = 0; j < int(b.size()); j++) res[i + j] += a[i] * b[j];
    while (int(res.size()) > 1 && res.back() == T()) res.pop_back();
    return res;
  }
  while (N & (N - 1)) N++;
  vector<complex<F>> f(N, complex<F>(0, 0));
  for (int i = 0; i < int(a.size()); i++) f[i].real(a[i]);
  for (int i = 0; i < int(b.size()); i++) f[i].imag(b[i]);
  fft(f); complex<F> r(0, F(-0.25) / N); for (int i = 0; i <= N / 2; i++) {
    int j = (N - i) & (N - 1);
    complex<F> prod = (f[j] * f[j] - conj(f[i] * f[i])) * r;
    f[i] = prod; f[j] = conj(prod);
  }
  fft(f); vector<T> res(N, T()); bool isIntegral = is_integral<T>::value;
  for (int i = 0; i < N; i++)
    res[i] = isIntegral ? round(real(f[i])) : real(f[i]);
  while (int(res.size()) > 1 && res.back() == T()) res.pop_back();
  return res;
}

// Integer Multiplication
// Template Arguments:
//   T: the type of each digit
//   F: the floating point type to use for FFT
// 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 moderate constant
// Time Complexity: O(n log n) where n = size(a) + size(b)
// Memory Complexity: O(n)
// Tested:
//   https://dmoj.ca/problem/atimesb
template <class T, class F = long double>
vector<T> mulIntFFT(const vector<T> &a, const vector<T> &b, T BASE = T(10)) {
  vector<T> res = mulPolyFFT<T, F>(a, 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;

// Two Phase Simplex to solve a linear programming problem with N variables
//   and M equations in canonical form:
//     max c^T x
//     subject to Ax <= b and x >= 0
//     where A is an M x N matrix; b is a vector with dimension M;
//     c, x are vectors with dimension N
// Template Arguments:
//   F: a floating point type
// Constructor Arguments:
//   A: a M x N coefficient matrix, memory is saved if this is moved and has
//     a capacity of M + 2 x N + 2
//   b: a constraint vector, with dimension M
//   c: a cost vector, with dimension N
//   INF: a value for infinity
//   EPS: a value for epsilon
// Fields:
//   M: the number of equations
//   N: the number of variables
//   INF: the value for infinity
//   EPS: the value for epsilon
//   val: the value of the maximum cost of the objective function c^T x,
//     INF if unbounded, -INF if infeasible
//   x: a solution vector of dimension N producing the maximum cost, can be
//     any solution if unbounded, empty if infeasible
// In practice, has a very small constant
// Time Complexity:
//   constructor: O(2^N), worst case, much faster in practice
// Memory Complexity: O(MN)
// Tested:
//   https://open.kattis.com/problems/cheeseifyouplease
//   https://open.kattis.com/problems/goatropes
//   https://codeforces.com/contest/375/problem/E
//   https://codeforces.com/contest/605/problem/C
//   https://codeforces.com/contest/1061/problem/E
//   https://codeforces.com/contest/1288/problem/F
//   https://www.spoj.com/problems/BABY/
//   https://dmoj.ca/problem/noi08p3
template <class F> struct Simplex {
  static_assert(is_floating_point<F>::value,
                "F must be a floating point type");
  int M, N; F INF, EPS, val; vector<int> IN, OUT; vector<F> x;
  vector<vector<F>> T;
  bool cmp(F a, int b, F c, int d) { return abs(a - c) > EPS ? a < c : b < d; }
  void pivot(int r, int s) {
    auto &a1 = T[r]; F inv1 = F(1) / a1[s];
    for (int i = 0; i <= M + 1; i++) if (i != r && abs(T[i][s]) > EPS) {
      auto &a2 = T[i]; F inv2 = a2[s] * inv1;
      for (int j = 0; j <= N + 1; j++) a2[j] -= a1[j] * inv2;
      a2[s] = a1[s] * inv2;
    }
    for (int j = 0; j <= N + 1; j++) if (j != s) a1[j] *= inv1;
    for (int i = 0; i <= M + 1; i++) if (i != r) T[i][s] *= -inv1;
    a1[s] = inv1; swap(IN[r], OUT[s]);
  }
  bool simplex(int phase) {
    int x = M + phase - 1; while (true) {
      int s = -1; for (int j = 0; j <= N; j++)
        if (OUT[j] != -phase
            && (s == -1 || cmp(T[x][j], OUT[j], T[x][s], OUT[s]))) s = j;
      if (T[x][s] >= -EPS) return true;
      int r = -1; for (int i = 0; i < M; i++)
        if (T[i][s] > EPS
            && (r == -1 || cmp(T[i][N + 1] * T[r][s], IN[i],
                               T[r][N + 1] * T[i][s], IN[r]))) r = i;
      if (r == -1) return false;
      pivot(r, s);
    }
  }
  Simplex(vector<vector<F>> A, const vector<F> &b, const vector<F> &c,
          F INF = numeric_limits<F>::infinity(), F EPS = F(1e-9))
      : M(b.size()), N(c.size()), INF(INF), EPS(EPS), IN(M), OUT(N + 1),
        T(move(A)) {
    T.reserve(M + 2); for (int i = 0; i < M; i++) {
      T[i].resize(N + 2, F());
      IN[i] = N + i; T[i][N] = F(-1); T[i][N + 1] = b[i];
    }
    T.emplace_back(N + 2, F()); T.emplace_back(N + 2, F());
    for (int j = 0; j < N; j++) { OUT[j] = j; T[M][j] = -c[j]; }
    OUT[N] = -1; T[M + 1][N] = F(1); int r = 0;
    for (int i = 1; i < M; i++) if (T[i][N + 1] < T[r][N + 1]) r = i;
    if (T[r][N + 1] < -EPS) {
      pivot(r, N);
      if (!simplex(2) || T[M + 1][N + 1] < -EPS) { val = -INF; return; }
      for (int i = 0; i < M; i++) if (IN[i] == -1) {
        int s = 0; for (int j = 1; j <= N; j++)
          if (s == -1 || cmp(T[i][j], OUT[j], T[i][s], OUT[s])) s = j;
        pivot(i, s);
      }
    }
    bool unbounded = !simplex(1); x.assign(N, F());
    for (int i = 0; i < M; i++) if (IN[i] < N) x[IN[i]] = T[i][N + 1];
    val = unbounded ? INF : T[M][N + 1];
  }
};


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

template <class T> using Matrix = vector<vector<T>>;

// Returns the size of the first dimension of a matrix
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   A: the matrix
// Return Value: the size of the first dimension of the matrix A
// In practice, has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T> int getN(const Matrix<T> &A) { return A.size(); }

// Returns the size of the second dimension of a matrix
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   A: the matrix
// Return Value: the size of the second dimension of the matrix A
// In practice, has a very small constant
// Time Complexity: O(1)
// Memory Complexity: O(1)
// Tested:
//   Fuzz Tested
template <class T>
int getM(const Matrix<T> &A) { return A.empty() ? 0 : A[0].size(); }

// Creates a matrix of size N x M filled with the default value of T
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   N: the size of the first dimension
//   M: the size of the second dimension
// Return Value: a matrix of size N x M filled with the default value of T
// In practice, has a very small constant
// Time Complexity: O(NM)
// Memory Complexity: O(NM)
// Tested:
//   Fuzz Tested
template <class T> Matrix<T> makeMatrix(int N, int M) {
  return vector<vector<T>>(N, vector<T>(M, T()));
}

// Creates an identity matrix of size N x N
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   N: the size of both dimensions
// Return Value: an identity matrix of size N x N
// In practice, has a very small constant
// Time Complexity: O(N^2)
// Memory Complexity: O(N^2)
// Tested:
//   https://www.spoj.com/problems/MPOW/
template <class T> Matrix<T> identityMatrix(int N) {
  Matrix<T> A = makeMatrix<T>(N, N);
  for (int i = 0; i < N; i++) A[i][i] = T(1);
  return A;
}

// Returns the transpose of a matrix
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   A: the matrix
// Return Value: the transpose of the matrix A
// In practice, has a very small constant
// Time Complexity: O(NM)
// Memory Complexity: O(NM)
// Tested:
//   Fuzz Tested
template <class T> Matrix<T> transpose(const Matrix<T> &A) {
  Matrix<T> C = makeMatrix<T>(getM(A), getN(A));
  for (int i = 0; i < getN(C); i++) for (int j = 0; j < getM(C); j++)
    C[i][j] = A[j][i];
  return C;
}

// Adds B to the matrix A
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   A: a reference to the first matrix
//   B: a constant reference to the matrix to add
// Return Value: a reference to the matrix A after B is added to it
// In practice, has a very small constant
// Time Complexity: O(NM)
// Memory Complexity: O(NM)
// Tested:
//   Fuzz Tested
template <class T> Matrix<T> &operator += (Matrix<T> &A, const Matrix<T> &B) {
  assert(getN(A) == getN(B) && getM(A) == getM(B));
  for (int i = 0; i < getN(A); i++) for (int j = 0; j < getM(A); j++)
    A[i][j] += B[i][j];
  return A;
}

// Adds the matrices A and B
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   A: the first matrix
//   B: a constant reference to the matrix to add
// Return Value: the matrix A plus B
// In practice, has a very small constant
// Time Complexity: O(NM)
// Memory Complexity: O(NM)
// Tested:
//   Fuzz Tested
template <class T>
Matrix<T> operator + (Matrix<T> A, const Matrix<T> &B) { return A += B; }

// Subtracts B from the matrix A
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   A: a reference to the first matrix
//   B: a constant reference to the matrix to subtract
// Return Value: a reference to the matrix A after B is subtracted from it
// Time Complexity: O(NM)
// Memory Complexity: O(NM)
// Tested:
//   Fuzz Tested
template <class T> Matrix<T> &operator -= (Matrix<T> &A, const Matrix<T> &B) {
  assert(getN(A) == getN(B) && getM(A) == getM(B));
  for (int i = 0; i < getN(A); i++) for (int j = 0; j < getM(A); j++)
    A[i][j] -= B[i][j];
  return A;
}

// Subtracts the matrices A and B
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   A: the first matrix
//   B: a constant reference to the matrix to add
// Return Value: the matrix A minus B
// In practice, has a very small constant
// Time Complexity: O(NM)
// Memory Complexity: O(NM)
// Tested:
//   Fuzz Tested
template <class T>
Matrix<T> operator - (Matrix<T> A, const Matrix<T> &B) { return A -= B; }

// Multiplies the matrices A and B, getM(A) must equal getN(B)
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   A: a constant reference to the first matrix
//   B: the matrix to multiply
// Return Value: the matrix A times B
// In practice, has a very small constant
// Time Complexity: O(N(A) M(B) M(A))
// Memory Complexity: O((N(A) + N(B)) M(B))
// Tested:
//   https://www.spoj.com/problems/MPOW/
template <class T> Matrix<T> operator * (const Matrix<T> &A, Matrix<T> B) {
  assert(getM(A) == getN(B));
  B = transpose(B); Matrix<T> C = makeMatrix<T>(getN(A), getN(B));
  for (int i = 0; i < getN(C); i++) for (int j = 0; j < getM(C); j++) {
    T tmp = T(); for (int k = 0; k < getM(A); k++) tmp += A[i][k] * B[j][k];
    C[i][j] = tmp;
  }
  return C;
}

// Multiplies the matrix B into the matrix A, getM(A) must equal getN(B)
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   A: a reference to the first matrix
//   B: a constant reference to the matrix to multiply
// Return Value: a reference to the matrix A after it is multipled by B
// In practice, has a very small constant
// Time Complexity: O(N(A) M(B) M(A))
// Memory Complexity: O((N(A) + N(B)) M(B))
// Tested:
//   https://www.spoj.com/problems/MPOW/
template <class T>
Matrix<T> &operator *= (Matrix<T> &A, const Matrix<T> &B) { return A = A * B; }

// Raises the square matrix A to the power of pow, getN(A) must equal getM(A)
// Template Arguments:
//   T: the type of each element in the matrix
//   U: the type of pow
// Function Arguments:
//   A: the base square matrix
// Return Value: the matrix A raised to the power pow
// In practice, has a very small constant
// Time Complexity: O(N^3 log pow)
// Memory Complexity: O(N^2)
// Tested:
//   https://www.spoj.com/problems/MPOW/
template <class T, class U> Matrix<T> powMat(Matrix<T> A, U pow) {
  assert(getN(A) == getM(A));
  Matrix<T> x = identityMatrix<T>(getN(A)); while (true) {
    if (pow % 2 == 1) x *= A;
    if ((pow /= 2) == 0) break;
    A *= A;
  }
  return x;
}

// Computes the determinant of a square matrix A, getN(A) must equal getM(A)
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   A: the matrix
// Return Value: the determinant of the matrix A
// In practice, has a very small constant
// Time Complexity: O(N^3)
// Memory Complexity: O(N^2)
// Tested:
//   https://judge.yosupo.jp/problem/matrix_det
//   https://dmoj.ca/problem/det
template <class T> T det(Matrix<T> A) {
  auto abs = [&] (T a) { return a >= 0 ? a : -a; };
  int N = getN(A); assert(N == getM(A)); T ret = T(1);
  for (int i = 0; i < N; i++) {
    int mx = i; for (int j = i + 1; j < N; j++)
      if (abs(A[j][i]) > abs(A[mx][i])) mx = j;
    if (i != mx) { ret = -ret; A[i].swap(A[mx]); }
    T inv = T(1) / A[i][i]; for (int j = i + 1; j < N; j++) {
      T alpha = A[j][i] * inv;
      for (int k = i + 1; k < N; k++) A[j][k] -= alpha * A[i][k];
    }
  }
  for (int i = 0; i < N; i++) ret *= A[i][i];
  return ret;
}

// Computes the inverse of a square matrix A, getN(A) must equal getM(A)
// Template Arguments:
//   T: the type of each element in the matrix
// Function Arguments:
//   A: the matrix
//   EPS: a value for epsilon
// Return Value: the inverse of the matrix A,
//   or a 0 x 0 matrix if it does not exist
// In practice, has a very small constant
// Time Complexity: O(N^3)
// Memory Complexity: O(N^2)
// Tested:
//   https://www.spoj.com/problems/MIFF/
//   https://judge.yosupo.jp/problem/inverse_matrix
template <class T> Matrix<T> invMat(Matrix<T> A, T EPS = T(1e-9)) {
  auto abs = [&] (T a) { return a >= 0 ? a : -a; };
  int N = getN(A); assert(N == getM(A)); Matrix<T> I = identityMatrix<T>(N);
  for (int i = 0; i < N; i++) {
    int mx = i; for (int j = i; j < N; j++)
      if (A[j][i] != 0) { mx = j; break; }
    if (abs(A[mx][i]) <= EPS) return makeMatrix<T>(0, 0);
    if (i != mx) { A[i].swap(A[mx]); I[i].swap(I[mx]); }
    T inv = T(1) / A[i][i]; for (int k = i; k < N; k++) A[i][k] *= inv;
    for (int k = 0; k < N; k++) I[i][k] *= inv;
    for (int j = 0; j < N; j++) if (j != i) {
      T alpha = A[j][i]; if (abs(alpha) <= EPS) continue;
      for (int k = i; k < N; k++) A[j][k] -= alpha * A[i][k];
      for (int k = 0; k < N; k++) I[j][k] -= alpha * I[i][k];
    }
  }
  return I;
}


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

// Computes the maximum axis aligned rectangular area of a histogram where
//   all bars have width 1
// Template Arguments:
//   T: the type of the heights
// Function Arguments:
//   A: a vector of the heights of each bar
// Return Value: the maximum axis aligned rectangular area of the histogram
// In practice, has a moderate constant
// Time Complexity: O(N)
// Memory Complexity: O(N)
// Tested:
//   https://www.acmicpc.net/problem/6549
template <class T> T maxRectAreaHistogram(vector<T> A) {
  T ret = T(); int N = A.size(), top = 0; vector<int> stk(N);
  for (int i = 0; i < int(A.size()); i++) {
    int j = i; while (top > 0 && A[stk[top - 1]] >= A[i]) {
      ret = max(ret, (i - stk[top - 1]) * A[stk[top - 1]]);
      A[j = stk[--top]] = A[i];
    }
    stk[top++] = j;
  }
  for (; top > 0; top--) ret = max(ret, (N - stk[top - 1]) * A[stk[top - 1]]);
  return ret;
}


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

// Computes the maximum submatrix of a boolean matrix where
//   all elements are true
// Function Arguments:
//   A: a vector of vectors of booleans
// Return Value: the maximum submatrix of the boolean matrix where
//   all elements are true
// In practice, has a moderate constant
// Time Complexity: O(NM)
// Memory Complexity: O(M)
// Tested:
//   https://dmoj.ca/problem/ccoprep16q1
int maxOneSubmatrix(const vector<vector<bool>> &A) {
  int N = A.size(), M = N == 0 ? 0 : A[0].size(), ret = 0; vector<int> H(M, 0);
  for (int i = N - 1; i >= 0; i--) {
    for (int j = 0; j < M; j++) H[j] = A[i][j] ? i == N - 1 ? 1 : H[j] + 1 : 0;
    ret = max(ret, maxRectAreaHistogram(H));
  }
  return ret;
}


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

// Computes an array z, where z[i] is the length of the longest subarray
//   starting at S[i], which is also a prefix of S (maximum k such
//   that S[i + j] == S[j] for all 0 <= j < k)
// Indices are 0-indexed and ranges are inclusive
//   T: the type of each element in the array
// Constructor Arguments:
//   S: a vector of type T
// Fields:
//   N: the length of the array
//   S: a vector of type T representing the array
//   z: a vector of integers where the ith element is the length of the
//     longest subarray starting at index i which is also a prefix of S
// In practice, has a moderate constant
// Time Complexity:
//   constructor: O(N)
// Memory Complexity: O(N)
// Tested:
//   https://judge.yosupo.jp/problem/zalgorithm
//   https://open.kattis.com/problems/stringmatching
template <class T> struct ZAlgorithm {
  int N; vector<int> z;
  ZAlgorithm(const vector<T> &S) : N(S.size()), z(N, 0) {
    z[0] = N; for (int i = 1, l = 0, r = 0; i < N; i++) {
      if (i <= r) z[i] = min(r - i + 1, z[i - l]);
      while (i + z[i] < N && S[z[i]] == S[i + z[i]]) z[i]++;
      if (i + z[i] - 1 > r) r = (l = i) + z[i] - 1;
    }
  }
};

// Finds all starting indices of a pattern array in a text array
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element in the array
// Constructor Arguments:
//   pat: the pattern array
//   txt: the text array
// Fields:
//   Z: the associated ZAlgorithm object constructed from the concatenation
//     of the pattern and the text array
//   matches: a vector of all starting indices of the pattern array in
//     the text array
// In practice, has a moderate constant, slightly slower than KMP
// Time Complexity:
//   constructor: O(N + M)
// Memory Complexity: O(N + M)
// Tested:
//   https://open.kattis.com/problems/stringmatching
template <class T> struct StringMatching {
  ZAlgorithm<T> Z; vector<int> matches;
  vector<T> init(vector<T> pat, const vector<T> &txt) {
    pat.reserve(pat.size() + txt.size());
    pat.insert(pat.end(), txt.begin(), txt.end()); return pat;
  }
  StringMatching(const vector<T> &pat, const vector<T> &txt)
      : Z(init(pat, txt)) {
    assert(int(pat.size()) >= 1); for (int i = 0; i < int(txt.size()); i++)
      if (Z.z[pat.size() + i] >= int(pat.size())) matches.push_back(i);
  }
};


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

// KMP array searching
// Indices are 0-indexed and ranges are inclusive
// Template Arguments:
//   T: the type of each element in the array
// Constructor Arguments:
//   pat: the pattern array
// Fields:
//   N: the length of the pattern
//   fail: a vector of integers where the ith value is the fallback index when
//     a match failure occures
// Functions
//   search(txt): returns the starting index of the first occurrence of pat in
//     txt, or -1 if there are none
//   multisearch(txt): returns a vector of all starting indices of the
//     pattern array in the text txt
// In practice, has a moderate constant, slightly faster than ZAlgorithm
// Time Complexity:
//   constructor: O(N)
//   search, multisearch: O(M)
// Memory Complexity: O(N)
// Tested:
//   https://dmoj.ca/problem/bf4
//   https://open.kattis.com/problems/stringmatching
//   https://codeforces.com/contest/1466/problem/G
template <class T> struct KMP {
  int N; vector<T> pat; vector<int> fail;
  KMP(const vector<T> &pat) : N(pat.size()), pat(pat), fail(N + 1, -1) {
    assert(N >= 1); for (int i = 0, j = -1; i < N;) {
      while (j >= 0 && pat[i] != pat[j]) j = fail[j];
      i++; j++; fail[i] = (i != N && pat[i] == pat[j] ? fail[j] : j);
    }
  }
  int search(const vector<T> &txt) {
    for (int i = 0, j = 0; i < int(txt.size()); i++, j++) {
      while (j >= 0 && txt[i] != pat[j]) j = fail[j];
      if (j == N - 1) return i - j;
    }
    return -1;
  }
  vector<int> multiSearch(const vector<T> &txt) {
    vector<int> ret; for (int i = 0, j = 0; i < int(txt.size()); i++, j++) {
      while (j >= 0 && (j == N || txt[i] != pat[j])) j = fail[j];
      if (j == N - 1) ret.push_back(i - j);
    }
    return ret;
  }
};