pe_nt

#ifndef PE_NT_
#define PE_NT_

#include "pe_base"
#include "pe_int128"
#include "pe_type_traits"
#include "pe_int"
#include "pe_span"
#include "pe_mod"
#include "pe_nt_base"

namespace pe {
#if defined(COMPILER_MSVC) || !OS_TYPE_WIN
template <typename T>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T) SqrtI(T d) {
  const auto result = static_cast<T>(std::sqrt(static_cast<long double>(d)));
  return result * result > d ? result - 1 : result;
}
#else
template <typename T>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T) SqrtI(T v) {
  return static_cast<T>(std::sqrt(static_cast<long double>(v)));
}
#endif

namespace internal {
template <typename T>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T) RootIImpl(T v, int n) {
  // n > 2 && v > 1
  T r = static_cast<T>(
      std::pow(static_cast<long double>(v), static_cast<long double>(1.) / n));
#if defined(COMPILER_MSVC)
  // Avoid floating round error when the compiler is msvc.
  if (r > 2) {
    --r;
  }
#endif
  // If r is 0, the for loop may take a long time.
  if (r < 1) {
    r = 1;
  }
  for (T nextr = r + 1;; ++nextr) {
    T t = v;
    for (int i = 0; i < n; ++i) {
      t /= nextr;
      if (t == 0) {
        return r;
      }
    }
    r = nextr;
  }
}
}  // namespace internal

template <typename T>
SL REQUIRES((is_builtin_integer_v<T>)) RETURN(T) RootI(T v, int n) {
  PE_ASSERT(n >= 0);
  if (n == 0) {
    return 1;
  }
  if (n == 1) {
    return v;
  }
  if (v == 0 || v == 1) {
    return v;
  }
  if (v == -1) {
    PE_ASSERT(n % 2 == 1);
    return v;
  }
  if (n == 2) {
    PE_ASSERT(v > 0);
    return SqrtI(v);
  }
  return v > 0 ? internal::RootIImpl(v, n) : -internal::RootIImpl(-v, n);
}

SL int64 SquareRootSafe(int64 n) {
  if (n <= 0) return 0;
  if (n == 1) return 1;
  for (int64 test = SqrtI(n);; ++test) {
    const int64 t = test * test;
    if (t == n) return test;
    if (t > n) return 0;
  }
  return 0;
}

SL int64 SquareRootUnsafe(int64 n) {
  if (n <= 0) return 0;
  if (n == 1) return 1;
  int64 ret = SqrtI(n);
  return ret * ret == n ? ret : 0;
}

SL int64 SquareRoot(int64 n) {
  if (n <= 0) return 0;
  if (n == 1) return 1;
  // see Donald L. Graham, Donald E. Knuth, Oren Patashnik, CONCRETE
  // MATHEMATICS, pp 70-72
  // int64 ret = static_cast<int64>(std::sqrt(n + 0.5));
  int64 ret = SqrtI(n);
  return ret * ret == n ? ret : 0;
}

// Returns the maximum idx where Power(v, idx) <= n.
// Returns -1 if v == 1 and n > 1.
template <typename T1, typename T2>
SL REQUIRES((is_builtin_integer_v<T1> && is_builtin_integer_v<T2>)) RETURN(T2)
    LogI(T1 v, T2 n) {
  PE_ASSERT(v >= 1);
  PE_ASSERT(n >= 1);
  if (n == 1) {
    return 0;
  }

  if (v == 1) {
    return -1;
  }

  const T2 limit = n / v;
  int idx = 0;
  for (T2 now = 1; now <= limit; now *= v, ++idx) {
  }

  return idx;
}

template <typename T1, typename T2>
SL std::tuple<T1, int64> ExtractFactor(T1 A, T2 B) {
  int64 c = 0;
  while (!IsZero(A) && A % B == 0) {
    A = static_cast<T1>(A / B), ++c;
  }
  return std::tuple<T1, int64>{A, c};
}

// Greatest common divisor
template <typename T>
SL REQUIRES((is_builtin_or_extended_integer_v<T>)) RETURN(T) Gcd(T m, T n) {
  while (!IsZero(n)) {
    T t = m % n;
    m = n;
    n = t;
  }
  return m;
}

// #if !defined(COMPILER_MSVC) && !PE_HAS_CPP17
template <typename T0, typename T1>
SL REQUIRES((is_general_integer_v<T0> && is_general_integer_v<T1>))
    RETURN(IntPromotion(T0, T1)) Gcd(const T0& m, const T1& n) {
  using ResultType = typename BiggerIntegralN<T0, T1>::type;
  return Gcd(static_cast<ResultType>(m), static_cast<ResultType>(n));
}
// #endif

template <typename T>
SL REQUIRES((is_general_integer_v<T>)) RETURN(T) Gcd(const T& a) {
  return a;
}

template <typename TH, typename... TT>
SL REQUIRES(((sizeof...(TT) >= 2) && (is_general_integer_v<TH>) &&
             (predicate_all<is_general_integer, TT...>::value)))
    RETURN(IntPromotion(TH, TT...)) Gcd(TH a, TT... b) {
  return Gcd(a, Gcd(b...));
}

template <typename T>
SL REQUIRES((is_general_integer_v<T>)) RETURN(T)
    Gcd(std::initializer_list<T> init) {
  T ret = 0;
  for (auto& v : init) {
    if (IsZero(ret)) {
      ret = v;
    } else {
      ret = Gcd(ret, v);
    }
  }
  return ret;
}

template <typename T>
SL REQUIRES((is_general_integer_v<T>)) RETURN(T)
    Gcd(const std::vector<T>& init) {
  T ret = 0;
  for (auto& v : init) {
    if (IsZero(ret)) {
      ret = v;
    } else {
      ret = Gcd(ret, v);
    }
  }
  return ret;
}

template <typename T>
SL REQUIRES((is_general_integer_v<T>)) RETURN(T) Lcm(const T& m, const T& n) {
  return m / Gcd(m, n) * n;
}

template <typename T0, typename T1>
SL REQUIRES((is_general_integer_v<T0> && is_general_integer_v<T1>))
    RETURN(IntPromotion(T0, T1)) Lcm(const T0& m, const T1& n) {
  using ResultType = typename BiggerIntegralN<T0, T1>::type;
  return Lcm(static_cast<ResultType>(m), static_cast<ResultType>(n));
}

template <typename T>
SL REQUIRES((is_general_integer_v<T>)) RETURN(T) Lcm(const T& a) {
  return a;
}

template <typename TH, typename... TT>
SL REQUIRES(((sizeof...(TT) >= 2) && (is_general_integer_v<TH>) &&
             (predicate_all<is_general_integer, TT...>::value)))
    RETURN(IntPromotion(TH, TT...)) Lcm(TH a, TT... b) {
  return Lcm(a, Lcm(b...));
}

template <typename T>
SL REQUIRES((is_general_integer_v<T>)) RETURN(T)
    Lcm(std::initializer_list<T> init) {
  T ret = 0;
  for (auto& v : init) {
    if (IsZero(ret)) {
      ret = v;
    } else {
      ret = ret / Gcd(ret, v) * v;
    }
  }
  return ret;
}

template <typename T>
SL REQUIRES((is_general_integer_v<T>)) RETURN(T)
    Lcm(const std::vector<T>& init) {
  T ret = 0;
  for (auto& v : init) {
    if (IsZero(ret)) {
      ret = v;
    } else {
      ret = ret / Gcd(ret, v) * v;
    }
  }
  return ret;
}

// Extended greatest common divisor
namespace internal {
template <typename T>
SL REQUIRES((is_builtin_integer_v<T> && pe_is_signed_v<T>)) RETURN(T)
    ExGcdImpl(T a, T b, T& x, T& y) {
  if (b == 0) {
    x = 1;
    y = 0;
    return a;
  }

  T x1 = 0, y1 = 0;
  const T ret = ExGcdImpl(b, a % b, x1, y1);
  // b * x1 + (a%b) * y1 = ret
  // b * x1 + (a - a / b * b) * y1 = ret
  // b * x1 + a * y1 - a / b * b * y1 = ret
  // y1 * a + (x1 - a / b * y1) b = ret
  x = y1;
  y = x1 - a / b * y1;
  return ret;
}
}  // namespace internal

template <typename T>
SL REQUIRES((is_builtin_integer_v<T> && pe_is_signed_v<T>)) RETURN(T)
    ExGcd(T a, T b, T& x, T& y) {
  PE_ASSERT(a != 0 || b != 0);
  const T r = internal::ExGcdImpl(a, b, x, y);
  PE_ASSERT(a * x + b * y == r);
  return r;
}

template <typename T>
SL REQUIRES((is_builtin_integer_v<T> && pe_is_signed_v<T>))
    RETURN(std::tuple<T, T, T>) ExGcd(T a, T b) {
  PE_ASSERT(a != 0 || b != 0);
  T x = 0, y = 0;
  const T r = internal::ExGcdImpl(a, b, x, y);
  PE_ASSERT(a * x + b * y == r);
  return std::tuple<T, T, T>{r, x, y};
}

template <typename T>
SL REQUIRES((is_builtin_integer_v<T> && pe_is_signed_v<T>)) RETURN(T)
    ModInv(T a, T mod) {
  T b = mod, s = 1, t = 0, old_a = a;
  while (b > 0) {
    std::swap(s -= t * (a / b), t);
    std::swap(a %= b, b);
  }
  if (a > 1) {
    std::cerr << "mod inv error: " << old_a << " " << mod << std::endl;
    exit(1);
  }
  return T(s) < 0 ? s + mod : s;
}

// p is prime
SL int64 InvOf(int64 x, int64 p) { return PowerMod(x, p - 2, p); }

template <typename T>
struct ModValue {
  int ok = 0;
  T value;
  T mod;
};

// Solves equation a x = b (mod m)
template <typename T>
SL REQUIRES((is_builtin_integer_v<T> && pe_is_signed_v<T>)) RETURN(ModValue<T>)
    SolveLinearEquation(T a, T b, T m) {
  PE_ASSERT(m);

  if (m < 0) m = -m;

  if (sizeof(T) <= 8) {
    PE_ASSERT(m <= PE_SOI63);
  }

  a %= m;
  b %= m;

  if (a < 0) a += m;
  if (b < 0) b += m;

  T xx = 0, yy = 0, d = ExGcd(a, m, xx, yy);

  if (b % d) return ModValue<T>{0};

  xx = (xx % m + m) % m;
  xx = xx * (b / d) % m;

  T mod = m / d;
  T value = xx % mod;

  return ModValue<T>{1, value, mod};
}

// Solves equation x = a (mod m1) x = b (mod m2)
template <typename T>
SL REQUIRES((is_builtin_integer_v<T> && pe_is_signed_v<T>)) RETURN(ModValue<T>)
    Crt2(T a, T m1, T b, T m2) {
  PE_ASSERT(m1 != 0 && m2 != 0);
  if (m1 < 0) m1 = -m1;
  if (m2 < 0) m2 = -m2;
  // x = a (mod m1)
  // x = b (mod m2)
  // x = k0 m1 + a
  // k0 m1 + a = b (mod m2)
  // k0 m1 = b - a (mod m2)
  a %= m1;
  b %= m2;
  if (a < 0) a = -a;
  if (b < 0) b = -b;

  ModValue<T> mod_value = SolveLinearEquation(m1, b - a, m2);
  if (!mod_value.ok) {
    return ModValue<T>{0};
  }
  // k0 = k1 * s + t
  // x = (k1 * s + t) * m1 + a
  // x = k1 * s * m1 + (a + t * m1)

  T mod = mod_value.mod * m1;
  T value = ((a + mod_value.value * m1) % mod + mod) % mod;

  return ModValue<T>{1, value, mod};
}

template <typename T>
SL REQUIRES((is_builtin_integer_v<T> && pe_is_signed_v<T>)) RETURN(ModValue<T>)
    CrtN(const T* val, const T* mod, int n) {
  PE_ASSERT(val != NULL && mod != NULL);
  ModValue<T> ret{1, 0, 1};
  for (int i = 0; i < n; ++i) {
    ret = Crt2(ret.value, ret.mod, val[i], mod[i]);
    if (!ret.ok) {
      return ModValue<T>{0};
    }
  }
  return ret;
}

template <typename T>
SL REQUIRES((is_builtin_integer_v<T> && pe_is_signed_v<T>)) RETURN(ModValue<T>)
    CrtN(const std::vector<T>& val, const std::vector<T>& mod) {
  return CrtN<T>(&val[0], &mod[0], std::size(val));
}

// Returns the number of x in the range [1, n] such that x is a multiple of d.
template <typename T = int64>
SL T CountMultiple(const int64 n, const int64 d) {
  return n / d;
}

// Returns the sum of x in the range [1, n] such that x is a multiple of d.
template <typename T = int64>
SL REQUIRES((!IsNModNumberV<T>)) RETURN(T)
    SumMultiple(const int64 n, const int64 d) {
  const T t0 = n / d;
  const T t1 = t0 & 1 ? ((t0 + 1) >> 1) * t0 : (t0 >> 1) * (t0 + 1);
  return t1 * d;
}

template <typename T = int64>
SL REQUIRES((IsNModNumberV<T>)) RETURN(T)
    SumMultiple(const int64 n, const int64 d) {
  const auto t0 = ExtractValue(T(n / d));
  const T t1 = t0 & 1 ? T((t0 + 1) >> 1) * t0 : T(t0 >> 1) * (t0 + 1);
  return t1 * d;
}

// Returns the number of x in the range [1, n] such that x = remain mod(mod).
template <typename T = int64>
SL T CountModValue(const int64 n, int64 remain, int64 mod) {
  if (remain == 0) {
    return CountMultiple<T>(n, mod);
  } else {
    return CountMultiple<T>(n + mod - remain, mod);
  }
}

namespace internal {
// The sum of i, i + d, i + 2d, ..., i + (n-1) * d
template <typename T = int64>
SL REQUIRES((!IsNModNumberV<T>)) RETURN(T)
    ApSum(const int64 i, const int64 d, const int64 n) {
  const T t0 = n & 1 ? ((n - 1) >> 1) * n : (n >> 1) * (n - 1);
  return n * i + t0 * d;
}

template <typename T = int64>
SL REQUIRES((IsNModNumberV<T>)) RETURN(T)
    ApSum(const int64 i, const int64 d, const int64 n) {
  const T t0 = n & 1 ? T((n - 1) >> 1) * n : T(n >> 1) * (n - 1);
  return T(n) * i + t0 * d;
}
}  // namespace internal

// Returns the sum of x in the range [1, n] such that x = remain mod(mod).
template <typename T = int64>
SL T SumModValue(const int64 n, int64 remain, int64 mod) {
  if (remain == 0) {
    return SumMultiple<T>(n, mod);
  } else {
    return internal::ApSum(remain, mod, (n + mod - remain) / mod);
  }
}

// Returns the number of x in the range [1, n] such that x = remain0 mod(mod0),
// and x = remain1 mod(mod1).
template <typename T = int64>
SL T CountModValue(const int64 n, int64 remain0, int64 mod0, int64 remain1,
                   int64 mod1) {
  ModValue<int64> t = Crt2<int64>(remain0, mod0, remain1, mod1);
  if (!t.ok) return 0;
  return CountModValue<T>(n, t.value, t.mod);
}

// Returns the sum of x in the range [1, n] such that x = remain0 mod(mod0),
// and x = remain1 mod(mod1).
template <typename T = int64>
SL T SumModValue(const int64 n, int64 remain0, int64 mod0, int64 remain1,
                 int64 mod1) {
  ModValue<int64> t = Crt2<int64>(remain0, mod0, remain1, mod1);
  if (!t.ok) return 0;
  return SumModValue<T>(n, t.value, t.mod);
}

// Returns the number of x in the range [1, n] such that gcd(x, a) = 1.
// The rad factors with mu are given by rad_factors_with_mu.
template <typename T = int64>
SL T CountCoprime(
    int64 n, const std::vector<std::pair<int64, int>>& rad_factors_with_mu) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : rad_factors_with_mu) {
    ret += CountMultiple<T>(n, iter.first) * iter.second;
  }
  return ret;
}

template <typename T = int64>
SL T CountCoprime(int64 n, const IntegerFactorization& f) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : GetRadFactorsWithMu(f, n)) {
    ret += CountMultiple<T>(n, iter.first) * iter.second;
  }
  return ret;
}

template <typename T = int64>
SL T CountCoprime(int64 n, int64 a) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : GetRadFactorsWithMu(a, n)) {
    ret += CountMultiple<T>(n, iter.first) * iter.second;
  }
  return ret;
}

// Returns the number x in the range [1, n] such that gcd(x, a) = 1 and x =
// remain mod(mod). The rad factors with mu are given by rad_factors_with_mu.
template <typename T = int64>
SL T CountCoprime(int64 n,
                  const std::vector<std::pair<int64, int>>& rad_factors_with_mu,
                  int64 remain, int64 mod) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : rad_factors_with_mu) {
    ret += CountModValue<T>(n, remain, mod, 0, iter.first) * iter.second;
  }
  return ret;
}

template <typename T = int64>
SL T CountCoprime(int64 n, const IntegerFactorization& f, int64 remain,
                  int64 mod) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : GetRadFactorsWithMu(f, n)) {
    ret += CountModValue<T>(n, remain, mod, 0, iter.first) * iter.second;
  }
  return ret;
}

template <typename T = int64>
SL T CountCoprime(int64 n, int64 a, int64 remain, int64 mod) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : GetRadFactorsWithMu(a, n)) {
    ret += CountModValue<T>(n, remain, mod, 0, iter.first) * iter.second;
  }
  return ret;
}

// Returns the sum x in the range [1, n] such that gcd(x, a) = 1.
// The rad factors with mu are given by rad_factors_with_mu.
template <typename T = int64>
SL T SumCoprime(int64 n,
                const std::vector<std::pair<int64, int>>& rad_factors_with_mu) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : rad_factors_with_mu) {
    ret += SumMultiple<T>(n, iter.first) * iter.second;
  }
  return ret;
}

template <typename T = int64>
SL T SumCoprime(int64 n, const IntegerFactorization& f) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : GetRadFactorsWithMu(f, n)) {
    ret += SumMultiple<T>(n, iter.first) * iter.second;
  }
  return ret;
}

template <typename T = int64>
SL T SumCoprime(int64 n, int64 a) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : GetRadFactorsWithMu(a, n)) {
    ret += SumMultiple<T>(n, iter.first) * iter.second;
  }
  return ret;
}

// Returns the sum x in the range [1, n] such that gcd(x, a) = 1 and x =
// remain mod(mod) The rad factors with mu are given by rad_factors_with_mu.
// The rad factors with mu are given by rad_factors_with_mu.
template <typename T = int64>
SL T SumCoprime(int64 n,
                const std::vector<std::pair<int64, int>>& rad_factors_with_mu,
                int64 remain, int64 mod) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : rad_factors_with_mu) {
    ret += SumModValue<T>(n, remain, mod, 0, iter.first) * iter.second;
  }
  return ret;
}

template <typename T = int64>
SL T SumCoprime(int64 n, const IntegerFactorization& f, int64 remain,
                int64 mod) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : GetRadFactorsWithMu(f, n)) {
    ret += SumModValue<T>(n, remain, mod, 0, iter.first) * iter.second;
  }
  return ret;
}

template <typename T = int64>
SL T SumCoprime(int64 n, int64 a, int64 remain, int64 mod) {
  if (n == 0) {
    return 0;
  }
  T ret = 0;
  for (auto& iter : GetRadFactorsWithMu(a, n)) {
    ret += SumModValue<T>(n, remain, mod, 0, iter.first) * iter.second;
  }
  return ret;
}

// dest[0] = 0
// dest[1] = 1
// dest[2] = 1 / 2
// dest[maxn-1] = 1 / (maxn-1)
// dest[maxn] = 1 / maxn
//
// maxn can be greater than or equal to mod
template <typename IT, typename TM>
SL REQUIRES((!IsNModNumberV<typename std::iterator_traits<IT>::value_type>))
    RETURN(void) InitInverse(IT dest, int64 maxn, TM mod) {
  PE_ASSERT(maxn >= 0);
  auto start = dest;
  *dest++ = 0;
  if (maxn >= 1) {
    *dest++ = 1;
  }
  using T = typename std::iterator_traits<IT>::value_type;
  auto maxi = std::min(maxn, static_cast<int64>(mod) - 1);
  for (int64 i = 2; i <= maxi; ++i) {
    const auto t = mod / i;
    const auto k = mod % i;
    const auto v = MulMod(t, t, mod);
    const auto u = MulMod(start[k], start[k], mod);
    *dest++ = static_cast<T>(MulMod(i, MulMod(v, u, mod), mod));
  }
  for (auto i = mod; i <= maxn; ++i) {
    *dest++ = start[i - mod];
  }
}

template <typename T, typename TM>
SL REQUIRES((!IsNModNumberV<T>)) RETURN(void)
    InitInverse(Span<T> dest, TM mod) {
  InitInverse<T*, TM>(std::data(dest), std::size(dest) - 1, mod);
}

template <typename IT>
SL REQUIRES((IsNModNumberV<typename std::iterator_traits<IT>::value_type>))
    RETURN(void) InitInverse(IT dest, int64 maxn) {
  PE_ASSERT(maxn >= 0);
  auto start = dest;
  *dest++ = 0;
  if (maxn >= 1) {
    *dest++ = 1;
  }
  using T = typename std::iterator_traits<IT>::value_type;
  using ModType = typename T::ints;
  const ModType mod = T().Mod();
  auto maxi = std::min(maxn, static_cast<int64>(mod) - 1);
  for (int64 i = 2; i <= maxi; ++i) {
    const auto t = mod / i;
    const auto k = mod % i;
    const auto v = MulMod(t, t, mod);
    const auto u = start[k] * start[k];
    *dest++ = u * v * i;
  }
  for (auto i = mod; i <= maxn; ++i) {
    *dest++ = start[i - mod];
  }
}

template <typename T>
SL REQUIRES((IsNModNumberV<T>)) RETURN(void) InitInverse(Span<T> dest) {
  InitInverse<T*>(std::data(dest), std::size(dest) - 1);
}

// 1 <= s <= e
//
// dest[0] = 1
// dest[1] = s
// dest[2] = s(s+1)
// dest[e-s+1]=s(s+1)...e
// idest[0] = 1
// idest[1] = 1/s
// idest[2] = 1/(s(s+1))
// idest[e-s+1]=1/(s(s+1)...e)
template <typename T>
SL REQUIRES((!IsNModNumberV<T>)) RETURN(void)
    InitSeqProd2(T* dest, T* idest, int64 s, int64 e, T mod) {
  if (dest == nullptr && idest == nullptr) {
    return;
  }

  const int64 cnt = static_cast<int64>(e - s + 1);
  if (mod == 1) {
    if (dest != nullptr) {
      std::fill(dest, dest + cnt + 1, 0);
    }
    if (idest != nullptr) {
      std::fill(idest, idest + cnt + 1, 0);
    }
    return;
  }

  T lastv = 1;
  if (dest == nullptr) {
    T v = s % mod;
    for (int64 i = s; i <= e; ++i) {
      lastv = MulMod(lastv, v, mod);
      v = v == mod - 1 ? 0 : v + 1;
    }
  } else {
    dest[0] = lastv;
    T v = s % mod;
    for (int64 i = s, j = 1; i <= e; ++i, ++j) {
      dest[j] = lastv = MulMod(lastv, v, mod);
      v = v == mod - 1 ? 0 : v + 1;
    }
  }
  if (idest == nullptr) {
    return;
  }

  // Assumption 1: mod is a prime number
  // Assumption 2: lastv != 0
  lastv = PowerMod(lastv, mod - 2, mod);
  idest[0] = 1;
  idest[cnt] = lastv;

  T v = e % mod;
  for (int64 j = cnt - 1; j > 0; --j) {
    idest[j] = MulMod(idest[j + 1], v, mod);
    v = v == 0 ? mod - 1 : v - 1;
  }
}

template <typename T>
SL REQUIRES((!IsNModNumberV<T>)) RETURN(void)
    InitSeqProd2(T* dest, T* idest, int64 e, T mod) {
  InitSeqProd2<T>(dest, idest, 1, e, mod);
}

template <typename T>
SL REQUIRES((!IsNModNumberV<T>)) RETURN(void)
    InitSeqProd2(Span<T> dest, Span<T> idest, T mod) {
  const int64 size1 = std::size(dest);
  const int64 size2 = std::size(idest);
  PE_ASSERT(size1 == size2);
  InitSeqProd2<T>(std::data(dest), std::data(idest), 1, size1 - 1, mod);
}

template <typename T>
SL REQUIRES((IsNModNumberV<T>)) RETURN(void)
    InitSeqProd2(T* dest, T* idest, int64 s, int64 e) {
  if (dest == nullptr && idest == nullptr) {
    return;
  }

  using ModType = typename T::ints;
  const ModType mod = T().Mod();
  const int64 cnt = static_cast<int64>(e - s + 1);
  if (mod == 1) {
    return;
  }

  ModType lastv = 1;
  if (dest == nullptr) {
    ModType v = s % mod;
    for (int64 i = s; i <= e; ++i) {
      lastv = MulMod(lastv, v, mod);
      v = v == mod - 1 ? 0 : v + 1;
    }
  } else {
    dest[0].SetValue(lastv);
    ModType v = s % mod;
    for (int64 i = s, j = 1; i <= e; ++i, ++j) {
      lastv = MulMod(lastv, v, mod);
      dest[j].SetValue(lastv);
      v = v == mod - 1 ? 0 : v + 1;
    }
  }
  if (idest == nullptr) {
    return;
  }

  // Assumption 1: mod is a prime number
  // Assumption 2: lastv != 0
  lastv = PowerMod(lastv, mod - 2, mod);
  idest[0] = 1;
  idest[cnt].SetValue(lastv);

  ModType v = e % mod;
  for (int64 j = cnt - 1; j > 0; --j) {
    idest[j].SetValue(MulMod(idest[j + 1].value(), v, mod));
    v = v == 0 ? mod - 1 : v - 1;
  }
}

template <typename T>
SL REQUIRES((IsNModNumberV<T>)) RETURN(void)
    InitSeqProd2(T* dest, T* idest, int64 e) {
  InitSeqProd2<T>(dest, idest, 1, e);
}

template <typename T>
SL REQUIRES((IsNModNumberV<T>)) RETURN(void)
    InitSeqProd2(Span<T> dest, Span<T> idest) {
  const int64 size1 = std::size(dest);
  const int64 size2 = std::size(idest);
  PE_ASSERT(size1 == size2);
  InitSeqProd2<T>(std::data(dest), std::data(idest), 1, size1 - 1);
}

template <typename IT, typename T, typename G>
SL void GenerateSequence(IT first, int64 maxn, T init, const G& gen) {
  if (maxn < 0) {
    return;
  }

  T value = init;
  *first++ = init;

  for (int64 i = 1; i <= maxn; ++i) {
    value = gen(value);
    *first++ = value;
  }
}

template <typename T, typename G>
SL void GenerateSequence(Span<T> target, T init, const G& gen) {
  GenerateSequence(std::data(target), static_cast<int64>(std::size(target)) - 1,
                   init, gen);
}

template <typename T, int D>
SL REQUIRES((!IsNModNumberV<T>)) RETURN(void)
    InitComb(T (*comb)[D], int maxn, int64 mod) {
  const int64 ONE = 1 % mod;
  for (int i = 0; i <= maxn; ++i) {
    for (int j = 1; j < i; ++j) {
      comb[i][j] = AddMod(comb[i - 1][j - 1], comb[i - 1][j], mod);
    }
    comb[i][0] = comb[i][i] = ONE;
  }
}

template <typename T, int D>
SL REQUIRES((!IsNModNumberV<T>)) RETURN(void)
    InitComb(T (*comb)[D], int64 mod) {
  InitComb(comb, D - 1, mod);
}

template <typename T, int D>
SL REQUIRES((IsNModNumberV<T>)) RETURN(void) InitComb(T (*comb)[D], int maxn) {
  const T ONE = 1;
  for (int i = 0; i <= maxn; ++i) {
    for (int j = 1; j < i; ++j) {
      comb[i][j] = comb[i - 1][j - 1] + comb[i - 1][j];
    }
    comb[i][0] = comb[i][i] = ONE;
  }
}

template <typename T, int D>
SL REQUIRES((IsNModNumberV<T>)) RETURN(void) InitComb(T (*comb)[D]) {
  InitComb(comb, D - 1);
}

// Checks whether x is the primitive root of p while  the factorization of p -1
// is f.
SL bool IsPrimitiveRoot(int64 x, int64 p, const IntegerFactorization& f) {
  return std::none_of(std::begin(f), std::end(f),
                      [=](const std::pair<int64, int>& it) {
                        return PowerMod(x, (p - 1) / it.first, p) == 1;
                      });
}

// Checks whether x is the primitive root of p.
SL bool IsPrimitiveRoot(int64 x, int64 p) {
  return IsPrimitiveRoot(x, p, Factorize(p - 1));
}

SL int64 FindPrimitiveRoot(int64 p) {
  if (p <= 1) return -1;
  if (p == 2) return 1;
  auto fp1 = Factorize(p - 1);
  for (int64 i = 2;; ++i) {
    if (IsPrimitiveRoot(i, p, fp1)) return i;
  }
  return -1;
}

SL int64 FindPrimitiveRoot(int64 p, const IntegerFactorization& fp1) {
  if (p <= 1) return -1;
  if (p == 2) return 1;
  for (int64 i = 2;; ++i) {
    if (IsPrimitiveRoot(i, p, fp1)) return i;
  }
  return -1;
}

SL int64 Ind(int64 a, int64 pr, int64 p) {
  using map_t = std::unordered_map<int64, int>;
  PE_ASSERT(a > 0);
  int64 inv_pr = InvOf(pr, p);
  auto mx = static_cast<int64>(std::sqrt(double(p - 1) + 0.5));
  if (mx * mx != p - 1) ++mx;
  map_t mem;
  int64 t1 = 1, t2 = a;
  for (unsigned int i = 0; i <= mx - 1; ++i) {
    mem[t2] = i;
    t2 = t2 * inv_pr % p;
  }
  int64 temp = PowerMod(pr, mx, p);
  for (unsigned int i = 0; i <= mx - 1; ++i) {
    auto where = mem.find(t1);
    if (where != mem.end()) {
      return i * mx + where->second;
    }
    t1 = t1 * temp % p;
  }
  return -1;
}

class IndSolver {
  static int64 SquareRootForInd(int64 x) {
    auto ret = static_cast<int64>(std::sqrt(double(x)));
    if (ret * ret < x) ++ret;
    PE_ASSERT(ret * ret >= x);
    return ret;
  }

 public:
  IndSolver(const int64 p, int64 pr = -1)
      : p_(p),
        pr_(pr < 0 ? FindPrimitiveRoot(p) : pr),
        mx_(SquareRootForInd(p_ - 1)) {
    const int64 inv_pr = InvOf(pr_, p);
    int64 t = 1;
    for (int i = 0; i <= mx_ - 1; ++i) {
      mem_[t] = i;
      t = t * inv_pr % p_;
    }
  }

  int64 operator()(int64 a) const { return Cal(a); }

  int64 Cal(int64 a) const {
    PE_ASSERT(a > 0);
    a %= p_;

    const int64 temp = PowerMod(pr_, mx_, p_);
    const int64 inv_a = InvOf(a, p_);
    int64 t = 1;

    for (int i = 0; i <= mx_ - 1; ++i) {
      auto where = mem_.find(MulMod(t, inv_a, p_));
      if (where != mem_.end()) {
        return i * mx_ + where->second;
      }
      t = t * temp % p_;
    }

    return -1;
  }

  int64 pr() const { return pr_; }

 private:
  const int64 p_;
  const int64 pr_;
  const int64 mx_;
  std::unordered_map<int64, int> mem_;
};

namespace flint {
#if ENABLE_FLINT && GMP_LIMB_BITS == 64
std::vector<int64> SquareRootMod(int64 n, const IntegerFactorization& mod_f) {
  const int size = static_cast<int>(std::size(mod_f));
  n_factor_t ff;
  for (int i = 0; i < size; ++i) {
    ff.p[i] = mod_f[i].first;
    ff.exp[i] = mod_f[i].second;
  }
  ff.num = size;

  ulong* ans;
  int cnt = static_cast<int>(n_sqrtmodn(&ans, n, &ff));
  if (cnt == 0) return {};

  std::vector<int64> ret(ans, ans + cnt);
  flint_free(ans);

  return ret;
}

std::vector<int64> SquareRootMod(int64 n, int64 mod) {
  return SquareRootMod(n, Factorize(mod));
}

namespace pmod {
std::vector<int64> SquareRootMod(int64 n, int64 p, int e = 1) {
  ulong* ans;
  int cnt = static_cast<int>(n_sqrtmod_primepow(&ans, n, p, e));
  if (cnt == 0) return {};

  std::vector<int64> ret(ans, ans + cnt);
  flint_free(ans);

  return ret;
}
}  // namespace pmod
#define HAS_SQUARE_ROOT_MOD_FLINT 1
#else
#define HAS_SQUARE_ROOT_MOD_FLINT 0
#endif
}  // namespace flint

namespace pmod {
namespace internal {
// returns x where x^2 = n mod p
// if no solution, return -1.
// https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
SL int64 SquareRootModImpl(int64 n, int64 p) {
  if (p == 2) {
    return n & 1;
  }
  if (n == 0) {
    return 0;
  }
  if (n == 1) {
    return 1;
  }
  if (PowerMod(n, (p - 1) / 2, p) == p - 1) {
    return -1;
  }
  if (p % 4 == 3) {
    return PowerMod(n, (p + 1) >> 2, p);
  }
  int64 Q = p - 1, S = 0;
  while (Q % 2 == 0) Q >>= 1, ++S;
  int64 z = 1;
  while (PowerMod(z, (p - 1) >> 1, p) == 1) ++z;
  int64 M = S, c = PowerMod(z, Q, p);
  int64 t = PowerMod(n, Q, p), R = PowerMod(n, (Q + 1) >> 1, p);
  for (;;) {
    if (t == 0) return -1;
    if (t == 1) return R;
    int64 i = 0, tt = t;
    while (tt != 1) tt = MulMod(tt, tt, p), ++i;
    const int64 b = PowerMod(c, 1 << (M - i - 1), p);
    M = i;
    c = MulMod(b, b, p);
    t = MulMod(t, c, p);
    R = MulMod(R, b, p);
  }
}
}  // namespace internal

SL std::vector<int64> SquareRootMod(int64 n, int64 p) {
  const int64 x = internal::SquareRootModImpl(n, p);
  if (x != -1) {
    const int64 y = x == 0 ? 0 : p - x;
    return x == y ? std::vector<int64>{x} : std::vector<int64>{x, y};
  } else {
    return std::vector<int64>{};
  }
}

SL std::vector<int64> RootMod(int64 a, int64 n, int64 p,
                              const int64 returned_count = -1) {
  if (p == 2) {
    return std::vector<int64>{a};
  }
  if (a == 0) {
    return std::vector<int64>{0};
  }
  IndSolver sv(p);
  const int64 r = sv(a);
  if (r == -1) return std::vector<int64>{};
  // x^n = a mod p
  // n log_x = r mod q
  const int64 q = p - 1;
  ModValue<int64> t = SolveLinearEquation(n % q, r % q, q);
  if (!t.ok) {
    return std::vector<int64>{};
  }
  std::vector<int64> result;
  int64 coe = PowerMod(sv.pr(), t.mod, p);
  int64 x = PowerMod(sv.pr(), t.value, p);
  int64 oldx = x;
  do {
    if (returned_count < 0 ||
        static_cast<int64>(std::size(result)) < returned_count) {
      result.push_back(x);
    } else {
      break;
    }
    x = MulMod(x, coe, p);
  } while (oldx != x);
  return result;
}
}  // namespace pmod

// (x+yi)(a+bi)
SL void GaussMul(int64& x, int64& y, int64 a, int64 b) {
  const int64 xx = x * a - b * y;
  const int64 yy = x * b + a * y;
  x = xx, y = yy;
}

#if PE_HAS_INT128
namespace internal {
SL void FactorizeForTwoSquaresByPmask(int64 n, IntegerFactorization& ret) {
  while (n != 1) {
    int now = pmask[n];
    int c = 0;
    while (n % now == 0) n /= now, ++c;
    if (c) {
      if ((now & 3) == 3) {
        if (IsOdd(c)) {
          IntegerFactorization{}.swap(ret);
          return;
        }
      }
      ret.emplace_back(static_cast<int64>(now), c);
    }
  }
}

// n > 2
SL IntegerFactorization FactorizeForTwoSquares(int64 n) {
  IntegerFactorization ret;

  int p2 = 0;
  while (IsEven(n)) n >>= 1, ++p2;
  if (p2 > 0) {
    ret.emplace_back(2LL, p2);
  }

  for (int i = 1; i < pcnt; ++i) {
    if (n <= maxp) {
      FactorizeForTwoSquaresByPmask(n, ret);
      return ret;
    }
    const int64 p = plist[i];
    const int64 test = p * p;
    if (test > n) break;
    int c = 0;
    while (n % p == 0) n /= p, ++c;
    if (c) {
      if ((p & 3) == 3) {
        if (IsOdd(c)) {
          return IntegerFactorization();
        }
      }
      ret.emplace_back(p, c);
    }
  }
  if (n != 1) {
    if ((n & 3) == 3) {
      return IntegerFactorization();
    }
    ret.emplace_back(n, 1);
  }
  return ret;
}
}  // namespace internal

// 4k+1 primes p = x^2+y^2 (x >= 0, y >= 0)
// 2 = 1^2 + 1^2
// 0, 1 are handled as special case.
SL std::vector<std::pair<int64, int64>> PrimeTwoSquares(const int64 p) {
  if (p == 0) {
    return std::vector<std::pair<int64, int64>>{{0LL, 0LL}};
  } else if (p == 1) {
    return std::vector<std::pair<int64, int64>>{{0LL, 1LL}};
  } else if (p == 2) {
    return std::vector<std::pair<int64, int64>>{{1LL, 1LL}};
  }
  if ((p & 3) != 1) {
    return std::vector<std::pair<int64, int64>>{};
  }
  int128 y = 1, x = 2;
  for (x = 2; PowerMod(x, (p - 1) >> 1, p) != p - 1; ++x) {
    ;
  }
  x = PowerMod(x, (p - 1) >> 2, p);

  for (int128 s = x * x + y * y; s != p;) {
    const int128 k = s / p, k2 = k >> 1;
    int128 a = (x % k + k) % k, b = (y % k + k) % k;
    if (a > k2) a = k - a;
    if (b > k2) b = k - b;
    if ((x * a + y * b) % k) std::swap(a, b);
    const int128 u = (x * a + y * b) / k, v = (b * x - a * y) / k;
    x = u, y = v;
    s = x * x + y * y;
  }
  if (y < 0) y = -y;
  if (x < 0) x = -x;
  if (x > y) std::swap(x, y);
  return std::vector<std::pair<int64, int64>>{{x, y}};
}

namespace internal {
SL void TwoSquaresImpl(int now, const IntegerFactorization& f, int64 u, int64 v,
                       std::vector<std::pair<int64, int64>>& result) {
  if (now == sz(f)) {
    u = Abs(u), v = Abs(v);
    if (u > v) std::swap(u, v);
    result.emplace_back(u, v);
  } else {
    const auto& item = f[now];
    if (item.first == 2) {
      for (int i = 0; i < item.second; ++i) {
        GaussMul(u, v, 1, 1);
      }
      TwoSquaresImpl(now + 1, f, u, v, result);
    } else if ((item.first & 3) == 3) {
      const int64 t = item.first;
      int64 me = 1;
      for (int i = 0, j = item.second >> 1; i < j; ++i) me *= t;
      TwoSquaresImpl(now + 1, f, u * me, v * me, result);
    } else {
      auto t = PrimeTwoSquares(item.first);
      const int64 a = t[0].first, b = t[0].second;
      std::vector<std::pair<int64, int64>> U, V;
      {
        int64 x = 1, y = 0;
        for (int i = 0; i <= item.second; ++i) {
          U.emplace_back(x, y);
          GaussMul(x, y, a, b);
        }
        x = 1, y = 0;
        for (int i = 0; i <= item.second; ++i) {
          V.emplace_back(x, y);
          GaussMul(x, y, a, -b);
        }
      }
      for (int i = 0, j = item.second; i <= item.second; ++i, --j) {
        int64 uu = u, vv = v;
        GaussMul(uu, vv, U[i].first, U[i].second);
        GaussMul(uu, vv, V[j].first, V[j].second);
        TwoSquaresImpl(now + 1, f, uu, vv, result);
      }
    }
  }
}

SL std::vector<std::pair<int64, int64>> TwoSquaresImpl(
    const IntegerFactorization& f) {
  std::vector<std::pair<int64, int64>> ret;
  TwoSquaresImpl(0, f, 1, 0, ret);
  std::sort(std::begin(ret), std::end(ret));
  ret.resize(unique(ret.begin(), ret.end()) - ret.begin());
  return ret;
}
}  // namespace internal

SL std::vector<std::pair<int64, int64>> TwoSquares(
    const IntegerFactorization& f) {
  for (auto& iter : f) {
    if ((iter.first & 3) == 3) {
      if (IsOdd(iter.second)) {
        return std::vector<std::pair<int64, int64>>{};
      }
    }
  }
  return internal::TwoSquaresImpl(f);
}

// Represent n as two squares: n = x^2 + y^2 where 0 <= x <= y.
SL std::vector<std::pair<int64, int64>> TwoSquares(const int64 n) {
  if (n == 0) {
    return std::vector<std::pair<int64, int64>>{{0LL, 0LL}};
  } else if (n == 1) {
    return std::vector<std::pair<int64, int64>>{{0LL, 1LL}};
  } else if (n == 2) {
    return std::vector<std::pair<int64, int64>>{{1LL, 1LL}};
  }
  auto f = internal::FactorizeForTwoSquares(n);
  if (std::empty(f)) {
    return std::vector<std::pair<int64, int64>>{};
  }
  return internal::TwoSquaresImpl(f);
}

SL void ExpandTwoSquares(std::vector<std::pair<int64, int64>>& ret,
                         std::pair<int64, int64> p) {
  ret.emplace_back(p.first, p.second);
  if (p.first != 0) {
    ret.emplace_back(-p.first, p.second);
  }
  if (p.second != 0) {
    ret.emplace_back(p.first, -p.second);
  }
  if (p.first != 0 && p.second != 0) {
    ret.emplace_back(-p.first, -p.second);
  }
}

SL std::vector<std::pair<int64, int64>> TwoSquaresFull(const int64 n) {
  if (n == 0) {
    return std::vector<std::pair<int64, int64>>{{0LL, 0LL}};
  }
  auto tmp = TwoSquares(n);
  std::vector<std::pair<int64, int64>> ret;
  ret.reserve(std::size(tmp) * 4);
  for (auto& iter : tmp) {
    ExpandTwoSquares(ret, iter);
    if (iter.first != iter.second) {
      ExpandTwoSquares(ret, {iter.second, iter.first});
    }
  }
  return ret;
}
#endif

template <typename T>
SL std::vector<int> ToBaseK(T n, int k, int width = 0) {
  if (n == 0) {
    return width <= 1 ? std::vector<int>{0} : std::vector<int>(width, 0);
  }
  PE_ASSERT(n >= 0 || k < 0);
  PE_ASSERT(Abs(k) >= 2);
  if (width <= 0) {
    std::vector<int> ret;
    while (n != 0) {
      int m = n % k;
      n /= k;
      if (m < 0) {
        m -= k;
        ++n;
      }
      ret.push_back(m);
    }
    return ret;
  } else {
    std::vector<int> ret(width, 0);
    for (int i = 0; i < width; ++i) {
      int m = n % k;
      n /= k;
      if (m < 0) {
        m -= k;
        ++n;
      }
      ret[i] = m;
    }
    return ret;
  }
}

template <typename T>
SL T FromBaseK(const std::vector<int>& vec, int k) {
  T ret = 0;
  const int size = static_cast<int>(std::size(vec));
  for (int i = size - 1; i >= 0; --i) ret = ret * k + vec[i];
  return ret;
}

SL std::pair<std::vector<int64>,
             std::vector<std::vector<std::pair<int64, int>>>>
Sieve(int64 start, int64 end) {
  std::vector<int64> sieve(end - start + 1);
  std::vector<std::vector<std::pair<int64, int>>> factorization(
      end - start + 1, std::vector<std::pair<int64, int>>());
  for (int64 i = start; i <= end; ++i) {
    sieve[i - start] = i;
  }

  int64 sqrt_n = SqrtI(end);

  for (int i = 0; i < pcnt; ++i) {
    const int64 p = plist[i];
    if (p > sqrt_n) {
      break;
    }
    int64 pos = start / p * p;
    while (pos < start) pos += p;
    for (;; pos += p) {
      if (pos > end) {
        break;
      }
      int c = 0;
      while (sieve[pos - start] % p == 0) {
        sieve[pos - start] /= p, ++c;
      }
      if (c > 0) {
        factorization[pos - start].push_back({p, c});
      }
    }
  }

  for (int64 i = start; i <= end; ++i) {
    if (i > 1 && sieve[i - start] != 1) {
      factorization[i - start].push_back({sieve[i - start], 1});
    }
  }
  return {sieve, factorization};
}
}  // namespace pe
// 65701
#endif
/*

*/