pe_nt_base

#ifndef PE_NT_BASE_
#define PE_NT_BASE_

#include "pe_base"
#include "pe_int128"
#include "pe_type_traits"
#include "pe_int"
#include "pe_mod"
#include "pe_range"

namespace pe {
// About primes
static int64 maxp;
static int64 maxp2;
static int pcnt;
static int* pmask = nullptr;
static int* plist = nullptr;
static int* phi = nullptr;
static int* mu = nullptr;
static int* rad = nullptr;
static int* sigma0 = nullptr;
static int64* sigma1 = nullptr;

SL void InitMaxp(int64 v) {
  maxp = v;
  maxp2 = maxp * maxp;
  assert(maxp >= 100000);
  assert(maxp <= INT_MAX);
}

// Prime pi
// http://oeis.org/A006880
// https://en.wikipedia.org/wiki/Prime-counting_function
// More results can be found here
// http://sweet.ua.pt/tos/primes.html
// http://numbers.computation.free.fr/Constants/Primes/Pix/resultstable.php
// http://www.primefan.ru/stuff/primes/table.html
// http://www.primefan.ru/stuff/primes/17(12)18.txt
// http://www.primefan.ru/stuff/primes/17(13)18.txt
// http://www.primefan.ru/stuff/primes/17(14)18.txt
// A fast implementation:
// https://github.com/kimwalisch/primecount
constexpr int64 kPrimePi[] = {
    0LL,
    4LL,                    // 10
    25LL,                   // 100
    168LL,                  // 1000
    1229LL,                 // 1e4
    9592LL,                 // 1e5
    78498LL,                // 1e6
    664579LL,               // 1e7
    5761455LL,              // 1e8
    50847534LL,             // 1e9
    455052511LL,            // 1e10
    4118054813LL,           // 1e11
    37607912018LL,          // 1e12
    346065536839LL,         // 1e13
    3204941750802LL,        // 1e14
    29844570422669LL,       // 1e15
    279238341033925LL,      // 1e16
    2623557157654233LL,     // 1e17
    24739954287740860LL,    // 1e18
    234057667276344607LL,   // 1e19
    2220819602560918840LL,  // 1e20
};

SL int64 EstimatePrimePi(int64 n) {
  if (n <= 1) return 0;
  auto ret = static_cast<int64>(1. * n / std::log(1. * n));
  return std::max(ret, static_cast<int64>(1LL));
}

SL void InitPmaskPlist(int*& pmask, int*& plist) {
  pmask = new int[maxp + 1];
  int64 size =
      std::max(static_cast<int64>((EstimatePrimePi(maxp + 1) + 1) * 1.1),
               static_cast<int64>(100000LL));
  plist = new int[size];
}

SL void DeinitPrimes() {
  pcnt = 0;
  maxp2 = maxp = 0;
  if (pmask) {
    delete[] pmask;
    pmask = nullptr;
  }
  if (plist) {
    delete[] plist;
    plist = nullptr;
  }
  if (phi) {
    delete[] phi;
    phi = nullptr;
  }
  if (mu) {
    delete[] mu;
    mu = nullptr;
  }
  if (rad) {
    delete[] rad;
    rad = nullptr;
  }
  if (sigma0) {
    delete[] sigma0;
    sigma0 = nullptr;
  }
  if (sigma1) {
    delete[] sigma1;
    sigma1 = nullptr;
  }
}

SL void InitPrimes() {
  if (maxp == 0) {
    InitMaxp(1000000);
  }

  InitPmaskPlist(pmask, plist);

  pcnt = 0;
  for (int i = 1; i <= maxp; ++i) pmask[i] = i;
  for (int i = 2; i <= maxp; ++i) {
    if (pmask[i] == i) {
      plist[pcnt++] = i;
    }
    for (int j = 0; j < pcnt; ++j) {
      const int64 t = static_cast<int64>(plist[j]) * i;
      if (t > maxp) break;
      pmask[t] = plist[j];
      if (static_cast<unsigned int>(i) % static_cast<unsigned int>(plist[j]) ==
          0) {
        break;
      }
    }
  }
}

SL void InitSigma0(int* sigma0) {
  sigma0[0] = 0;
  sigma0[1] = 1;
  for (int i = 2; i <= maxp; ++i) {
    if (pmask[i] == i) {
      sigma0[i] = 2;
    } else {
      const int p = pmask[i];
      int e = 2;
      int j = i / p;
      for (; j % p == 0; j /= p) {
        ++e;
      }
      sigma0[i] = e * sigma0[j];
    }
  }
}

SL void InitSigma1(int64* sigma1) {
  sigma1[0] = 0;
  sigma1[1] = 1;
  for (int i = 2; i <= maxp; ++i) {
    if (pmask[i] == i) {
      sigma1[i] = i + 1;
    } else {
      const int p = pmask[i];
      int d = p;
      int64 s = 1 + p;
      int j = i / p;
      for (; j % p == 0; j /= p) {
        d *= p;
        s += d;
      }
      sigma1[i] = s * sigma1[j];
    }
  }
}

SL void InitPrimes(int cal_phi, int cal_mu, int cal_rad, int cal_sigma0,
                   int cal_sigma1) {
  if (maxp == 0) {
    InitMaxp(1000000);
  }

  InitPmaskPlist(pmask, plist);

  pcnt = 0;
  if (cal_phi) phi = new int[maxp + 1];
  if (cal_mu) mu = new int[maxp + 1];
  if (cal_rad) rad = new int[maxp + 1];

  for (int i = 1; i <= maxp; ++i) pmask[i] = i;
  if (phi) phi[0] = 0, phi[1] = 1;
  if (mu) mu[0] = 0, mu[1] = 1;
  if (rad) rad[0] = 0, rad[1] = 1;

  for (int i = 2; i <= maxp; ++i) {
    if (pmask[i] == i) {
      if (mu) mu[i] = -1;
      if (phi) phi[i] = i - 1;
      if (rad) rad[i] = i;
      plist[pcnt++] = i;
    }
    for (int j = 0; j < pcnt; ++j) {
      const int64 t = static_cast<int64>(plist[j]) * i;
      if (t > maxp) break;
      pmask[t] = plist[j];
      if (static_cast<unsigned int>(i) % static_cast<unsigned int>(plist[j]) ==
          0) {
        if (mu) mu[t] = 0;
        if (phi) phi[t] = phi[i] * plist[j];
        if (rad) rad[t] = 0;
        break;
      } else {
        if (mu) mu[t] = -mu[i];
        if (phi) phi[t] = phi[i] * (plist[j] - 1);
        if (rad) rad[t] = rad[i] * plist[j];
      }
    }
  }

  if (cal_sigma0) {
    sigma0 = new int[maxp + 1];
    InitSigma0(sigma0);
  }
  if (cal_sigma1) {
    sigma1 = new int64[maxp + 1];
    InitSigma1(sigma1);
  }
}

struct IntegerFactorization;
SL IntegerFactorization Factorize(int64 n);

struct IntegerFactorization : public std::vector<std::pair<int64, int>> {
  using base = std::vector<std::pair<int64, int>>;
  // using base::base;

  IntegerFactorization() = default;

  IntegerFactorization(int64 t) : IntegerFactorization(Factorize(t)) {}

  IntegerFactorization(const IntegerFactorization& o) = default;

  IntegerFactorization(IntegerFactorization&& o) noexcept
      : base(std::move(o)) {}

  IntegerFactorization(const base& o) : base(o) {}

  IntegerFactorization(base&& o) : base(std::move(o)) {}

  IntegerFactorization(const std::map<int64, int>& o)
      : base(o.begin(), o.end()) {}

  IntegerFactorization& operator=(const IntegerFactorization& o) {
    this->base::operator=(o);
    return *this;
  }

  IntegerFactorization& operator=(IntegerFactorization&& o) noexcept {
    this->base::operator=(std::move(o));
    return *this;
  }

  IntegerFactorization& operator=(const base& o) {
    this->base::operator=(o);
    return *this;
  }

  IntegerFactorization& operator=(base&& o) {
    this->base::operator=(std::move(o));
    return *this;
  }

  std::map<int64, int> ToMap() const {
    return std::map<int64, int>(base::begin(), base::end());
  }
};

static const IntegerFactorization IF_ONE;

SL IntegerFactorization operator*(const IntegerFactorization& a,
                                  const IntegerFactorization& b) {
  if (std::empty(a)) {
    return b;
  }
  if (std::empty(b)) {
    return a;
  }
  auto ret = a.ToMap();
  for (auto& iter : b) {
    ret[iter.first] += iter.second;
  }
  return ret;
}

SL std::ostream& operator<<(std::ostream& o, const IntegerFactorization& v) {
  const int64 n = static_cast<int64>(std::size(v));
  for (int64 i = 0; i < n; ++i) {
    o << v[i].first;
    if (v[i].second > 1) o << "^" << v[i].second;
    if (i < n - 1) o << " * ";
  }
  return o;
}

SL void FactorizeByPmask(uint64 n, IntegerFactorization& ret) {
  while (n != 1) {
    unsigned int now = pmask[n];
    int c = 0;
    while (n >= now && n % now == 0) n /= now, ++c;
    if (c) ret.emplace_back(static_cast<int64>(now), c);
  }
}

namespace internal {
SL void FactorizeWarningOnce(int64 n) {
  static std::atomic_int WARNING_ONCE_FLAG{0};
  if (WARNING_ONCE_FLAG.exchange(1) == 0) {
    std::cerr << "Factorize is slow." << std::endl;
    std::cerr << std::setw(40) << "maxp = " << maxp << std::endl;
    std::cerr << std::setw(40) << "maxp2 = " << maxp2 << std::endl;
    if (pcnt > 0) {
      std::cerr << std::setw(40) << "plist[pcnt-1] = " << plist[pcnt - 1]
                << std::endl;
      std::cerr << std::setw(40) << "plist[pcnt-1]^2 = "
                << (int64)plist[pcnt - 1] * plist[pcnt - 1] << std::endl;
    }
    std::cerr << std::setw(40) << "n = " << n << std::endl;
  }
}
}  // namespace internal

SL IntegerFactorization Factorize(int64 n) {
  IntegerFactorization ret;
  if (n <= 1) {
    return ret;
  }

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

  if (i == pcnt) {
    internal::FactorizeWarningOnce(n);
    {
      int p2 = 0;
      while ((n & 1) == 0) n >>= 1, ++p2;
      if (p2) ret.emplace_back(2, p2);
    }
    const int64 start = pcnt <= 1 ? 3 : plist[pcnt - 1] + 2;
    for (int64 p = start;; p += 2) {
      const int64 test = p * p;
      if (test > n) break;
      int c = 0;
      while (n >= p && n % p == 0) n /= p, ++c;
      if (c) ret.emplace_back(p, c);
    }
  }

  if (n != 1) ret.emplace_back(n, 1);

  return ret;
}

SL IntegerFactorization Factorize(int64 n, const std::vector<int64>& hint) {
  IntegerFactorization ret;
  if (n <= 1) {
    return ret;
  }

  for (auto& h : hint) {
    for (auto iter : Factorize(h)) {
      const int64 p = iter.first;
      int c = 0;
      while (n % p == 0) n /= p, ++c;
      if (c) ret.emplace_back(p, c);
    }
  }

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

  if (i == pcnt) {
    internal::FactorizeWarningOnce(n);
    {
      int p2 = 0;
      while ((n & 1) == 0) n >>= 1, ++p2;
      if (p2) ret.emplace_back(2, p2);
    }
    const int64 start = pcnt <= 1 ? 3 : plist[pcnt - 1] + 2;
    for (int64 p = start;; p += 2) {
      const int64 test = p * p;
      if (test > n) break;
      int c = 0;
      while (n >= p && n % p == 0) n /= p, ++c;
      if (c) ret.emplace_back(p, c);
    }
  }

  if (n != 1) ret.emplace_back(n, 1);

  return ret;
}

SL void FactorizePowerByPmask(int64 n, int r, IntegerFactorization& ret) {
  while (n != 1) {
    int now = pmask[n];
    int c = 0;
    while (n % now == 0) n /= now, ++c;
    if (c) ret.emplace_back(static_cast<int64>(now), c * r);
  }
}

SL IntegerFactorization FactorizePower(int64 n, int r) {
  IntegerFactorization ret;
  if (n <= 1) {
    return ret;
  }

  for (int i = 0; i < pcnt; ++i) {
    if (n <= maxp) {
      FactorizePowerByPmask(n, r, 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) ret.emplace_back(p, c * r);
  }
  if (n != 1) ret.emplace_back(n, r);
  return ret;
}

SL int IsPrime(int64 n) {
  if (n <= 1) return 0;
  if (n == 2) return 1;
  if ((n & 1) == 0) return 0;
  if (n <= maxp) return pmask[n] == n;
  PE_ASSERT(n <= maxp2);
  for (int i = 0; i < pcnt; ++i) {
    const int64 p = plist[i];
    const int64 test = p * p;
    if (test > n) break;
    if (n % p == 0) return 0;
  }
  return 1;
}

// prime test for n <= 2^63 - 1
// http://oeis.org/A014233
constexpr int64 sopp[] = {
    2047LL,                 // 2
    1373653LL,              // 3
    25326001LL,             // 5
    3215031751LL,           // 7
    2152302898747LL,        // 11
    3474749660383LL,        // 13
    341550071728321LL,      // 17
    341550071728321LL,      // 19
    3825123056546413051LL,  // 23
    3825123056546413051LL,  // 29
    3825123056546413051LL,  // 31
                            // bigger than 2^64             // 37
                            // 2^63 = 9223372036854775808   //
};
constexpr int64 kSoppSize = std::size(sopp);

namespace internal {
SL int MrTestImpl(int s, int64 t, int64 n, int x) {
  // PE_ASSERT((n&1) == 0)

  t = PowerMod<int64, int64, int64>(x, t, n);
  if (t == 1 || t == n - 1) return 1;

  for (int i = 1; i < s; ++i) {
    t = MulMod(t, t, n);
    if (t == n - 1) return 1;
  }

  return 0;
}
}  // namespace internal

SL int MrTest(int64 n, int x) {
  if (n <= 1) return 0;
  if ((n & 1) == 0) return 0;

  int s = 0;
  int64 t = n - 1;
  while ((t & 1) == 0) ++s, t >>= 1;

  return internal::MrTestImpl(s, t, n, x);
}

SL int IsPrimeEx(int64 n) {
  if (n <= 1) return 0;
  if (n == 2) return 1;
  if ((n & 1) == 0) return 0;
  if (n <= maxp) return pmask[n] == n;

  for (int i = 1; i < 20; ++i) {
    if (n % plist[i] == 0) return 0;
  }

  int s = 0;
  int64 t = n - 1;
  while ((t & 1) == 0) ++s, t >>= 1;

  for (int i = 0; i < kSoppSize; ++i) {
    const int64 val = sopp[i];
    if (n == val) return 0;
    if (!internal::MrTestImpl(s, t, n, plist[i])) return 0;
    if (n < val) break;
  }
  return 1;
}

template <typename T>
struct FactorizedInteger {
  FactorizedInteger() : n_(1) {}

  FactorizedInteger(T n) : n_(n), factorization_(n) {}

  FactorizedInteger(T n, const IntegerFactorization& f)
      : n_(n), factorization_(f) {}
  FactorizedInteger(T n, IntegerFactorization&& f)
      : n_(n), factorization_(std::move(f)) {}

  FactorizedInteger(const IntegerFactorization& f)
      : n_(-1), factorization_(f) {}
  FactorizedInteger(IntegerFactorization&& f)
      : n_(-1), factorization_(std::move(f)) {}

  FactorizedInteger& operator=(const FactorizedInteger& o) {
    n_ = o.n_;
    factorization_ = o.factorization_;
    return *this;
  }

  FactorizedInteger& operator=(FactorizedInteger&& o) noexcept {
    n_ = o.n_;
    factorization_ = std::move(o.factorization_);
    return *this;
  }

  T GetValue() const {
    if (n_ < 0) {
      T tmp = 1;
      for (auto iter : factorization_) {
        for (int i = 0; i < iter.second; ++i) {
          tmp *= iter.first;
        }
      }
      n_ = tmp;
    }
    return n_;
  }

  const IntegerFactorization& GetFactorization() const {
    return factorization_;
  }

  std::vector<T> GetFactors(const T maxv = -1) const {
    if (maxv == 0) {
      return std::vector<T>();
    }

    int max_size = 1;
    for (auto& iter : factorization_) {
      max_size *= iter.second + 1;
    }

    std::vector<T> result(max_size);
    result[0] = 1;

    int top = 1;
    for (auto& iter : factorization_) {
      const int last = top;
      const T p = iter.first;
      T v = 1;
      for (int i = 0; i < iter.second; ++i) {
        v *= p;
        if (maxv > 0 && v > maxv) break;
        for (int j = 0; j < last; ++j) {
          const T tmp = result[j] * v;
          if (maxv < 0 || tmp <= maxv) {
            result[top++] = tmp;
          }
        }
      }
    }

    result.resize(top);
    return result;
  }

  std::vector<T> GetPrimeFactors() const {
    if (std::empty(factorization_)) {
      return {};
    }

    std::vector<T> ans;
    ans.reserve(std::size(factorization_));

    for (auto& iter : factorization_) {
      ans.push_back(iter.first);
    }
    return ans;
  }

  std::vector<T> GetRadFactors(const T maxv = -1) const {
    if (maxv == 0) {
      return std::vector<T>();
    }

    const int max_size = 1 << std::size(factorization_);
    std::vector<T> result(max_size);
    result[0] = 1;

    int top = 1;
    for (auto& iter : factorization_) {
      const int last = top;
      const int64 p = iter.first;
      T v = 1;
      for (int i = 0; i < 1; ++i) {
        v *= p;
        if (maxv > 0 && v > maxv) break;
        for (int j = 0; j < last; ++j) {
          const T tmp = result[j] * v;
          if (maxv < 0 || tmp <= maxv) {
            result[top++] = tmp;
          }
        }
      }
    }

    result.resize(top);
    return result;
  }

  std::vector<std::pair<T, int>> GetRadFactorsWithMu(const T maxv = -1) const {
    if (maxv == 0) {
      return std::vector<std::pair<T, int>>();
    }

    const int max_size = 1 << std::size(factorization_);
    std::vector<std::pair<T, int>> result(max_size);
    result[0] = {1LL, 1};

    int top = 1;
    for (auto& iter : factorization_) {
      const int last = top;
      const int64 p = iter.first;
      T v = 1;
      for (int i = 0; i < 1; ++i) {
        v *= p;
        if (maxv > 0 && v > maxv) break;
        for (int j = 0; j < last; ++j) {
          const T tmp = result[j].first * v;
          const int mu = -result[j].second;
          if (maxv < 0 || tmp <= maxv) {
            result[top++] = {tmp, mu};
          }
        }
      }
    }

    result.resize(top);
    return result;
  }

  T CalPhi() const {
    if (std::empty(factorization_)) {
      return 1;
    }
    if (phi && n_ >= 0 && n_ <= maxp) {
      return phi[n_];
    }
    T phi_n = GetValue();
    for (auto& it : factorization_) {
      phi_n -= phi_n / it.first;
    }
    return phi_n;
  }

  int CalMu() const {
    if (std::empty(factorization_)) {
      return 1;
    }
    if (mu && n_ >= 0 && n_ <= maxp) {
      return mu[n_];
    }
    for (auto& iter : factorization_) {
      if (iter.second > 0) {
        return 0;
      }
    }
    return std::size(factorization_) & 1 ? -1 : 1;
  }

  T CalRad() const {
    if (std::empty(factorization_)) {
      return 1;
    }
    if (rad && n_ >= 0 && n_ <= maxp) {
      return rad[n_];
    }

    T ret = 1;
    for (auto& iter : factorization_) {
      ret *= iter.first;
    }
    return ret;
  }

  int CalSmallOmega() const {
    return static_cast<int>(std::size(factorization_));
  }
  int CalBigOmega() const {
    int ret = 0;
    for (auto& iter : factorization_) {
      ret += iter.second;
    }
    return ret;
  }

  T CalSigma0() const {
    if (std::empty(factorization_)) {
      return 1;
    }
    if (sigma0 && n_ >= 0 && n_ <= maxp) {
      return sigma0[n_];
    }

    T ret = 1;
    for (auto& iter : factorization_) {
      ret *= iter.second + 1;
    }
    return ret;
  }

  T CalSigma1() const {
    if (std::empty(factorization_)) {
      return 1;
    }
    if (sigma1 && n_ >= 0 && n_ <= maxp) {
      return sigma1[n_];
    }

    T ret = 1;
    for (auto& iter : factorization_) {
      T s = 1 + iter.first;
      T d = iter.first;
      for (int i = 1; i < iter.second; ++i) {
        d *= iter.first;
        s += d;
      }
      ret *= s;
    }
    return ret;
  }

  int IsSquareFree() const {
    for (auto& iter : factorization_) {
      if (iter.second > 1) {
        return 0;
      }
    }
    return 1;
  }

 private:
  mutable T n_;  // It can be -1
  IntegerFactorization factorization_;
};

using FI64 = FactorizedInteger<int64>;

SL std::vector<int64> GetFactors(const FactorizedInteger<int64>& fi,
                                 int64 maxv = -1) {
  return fi.GetFactors(maxv);
}

template <typename T>
SL std::vector<T> GetFactors(const FactorizedInteger<T>& fi, T maxv = -1) {
  return fi.GetFactors(maxv);
}

SL std::vector<int64> GetPrimeFactors(const FactorizedInteger<int64>& fi) {
  return fi.GetPrimeFactors();
}

template <typename T>
SL std::vector<T> GetPrimeFactors(const FactorizedInteger<T>& fi) {
  return fi.GetPrimeFactors();
}

SL std::vector<int64> GetRadFactors(const FactorizedInteger<int64>& fi,
                                    int64 maxv = -1) {
  return fi.GetRadFactors(maxv);
}

template <typename T>
SL std::vector<T> GetRadFactors(const FactorizedInteger<T>& fi, T maxv = -1) {
  return fi.GetRadFactors(maxv);
}

SL std::vector<std::pair<int64, int>> GetRadFactorsWithMu(
    const FactorizedInteger<int64>& fi, int64 maxv = -1) {
  return fi.GetRadFactorsWithMu(maxv);
}

template <typename T>
SL std::vector<std::pair<T, int>> GetRadFactorsWithMu(
    const FactorizedInteger<T>& fi, T maxv = -1) {
  return fi.GetRadFactorsWithMu(maxv);
}

SL int64 CalPhi(int64 n) {
  if (n <= 0) return 0;
  if (phi && n <= maxp) return phi[n];
  auto fn = Factorize(n);
  int64 phi_n = n;
  for (auto& it : fn) {
    phi_n -= phi_n / it.first;
  }
  return phi_n;
}

SL int64 CalPhi(const FactorizedInteger<int64>& fi) { return fi.CalPhi(); }

template <typename T>
SL T CalPhi(const FactorizedInteger<T>& fi) {
  return fi.CalPhi();
}

namespace internal {
SL int64 CalMuImplByPmask(int64 n, int64 v = 1) {
  while (n != 1) {
    const int now = pmask[n];
    int c = 0;
    while (n % now == 0) n /= now, ++c;
    if (c > 1) {
      return 0;
    } else if (c == 1) {
      v = -v;
    }
  }
  return v;
}

SL int64 CalMuImpl(int64 n) {
  int64 v = 1;

  for (int i = 0; i < pcnt; ++i) {
    if (n <= maxp) {
      return CalMuImplByPmask(n, v);
    }
    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 > 1) {
      return 0;
    } else if (c == 1) {
      v = -v;
    }
  }
  if (n != 1) v = -v;
  return v;
}
}  // namespace internal

SL int64 CalMu(int64 n) {
  if (n <= 0) return 0;
  if (n == 1) return 1;
  if (mu && n <= maxp) return mu[n];
  return internal::CalMuImpl(n);
}

SL int64 CalMu(const FactorizedInteger<int64>& fi) { return fi.CalMu(); }

template <typename T>
SL T CalMu(const FactorizedInteger<T>& fi) {
  return fi.CalMu();
}

namespace internal {
SL int64 CalRadImplByPmask(int64 n, int64 v = 1) {
  while (n != 1) {
    const int now = pmask[n];
    int c = 0;
    while (n % now == 0) n /= now, ++c;
    if (c >= 1) v *= now;
  }
  return v;
}

SL int64 CalRadImpl(int64 n) {
  int64 v = 1;

  for (int i = 0; i < pcnt; ++i) {
    if (n <= maxp) {
      return CalRadImplByPmask(n, v);
    }
    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 >= 1) v *= p;
  }
  if (n != 1) v *= n;
  return v;
}
}  // namespace internal

SL int64 CalRad(int64 n) {
  if (n <= 0) return 0;
  if (n == 1) return 1;
  if (rad && n <= maxp) return rad[n];
  return internal::CalRadImpl(n);
}

SL int64 CalRad(const FactorizedInteger<int64>& fi) { return fi.CalRad(); }

template <typename T>
SL T CalRad(const FactorizedInteger<T>& fi) {
  return fi.CalRad();
}

SL int CalSmallOmega(const FactorizedInteger<int64>& fi) {
  return fi.CalSmallOmega();
}

template <typename T>
SL int CalSmallOmega(const FactorizedInteger<T>& fi) {
  return fi.CalSmallOmega();
}

SL int CalBigOmega(const FactorizedInteger<int64>& fi) {
  return fi.CalBigOmega();
}

template <typename T>
SL int CalBigOmega(const FactorizedInteger<T>& fi) {
  return fi.CalBigOmega();
}

namespace internal {
SL int64 CalSigma0ByPmask(int64 n, int64 v = 1) {
  while (n != 1) {
    const int p = pmask[n];
    int c = 0;
    while (n % p == 0) n /= p, ++c;
    if (c >= 1) v *= c + 1;
  }
  return v;
}

SL int64 CalSigma0Impl(int64 n) {
  int64 v = 1;

  for (int i = 0; i < pcnt; ++i) {
    if (n <= maxp) {
      return CalSigma0ByPmask(n, v);
    }
    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 >= 1) v *= c + 1;
  }
  if (n != 1) v *= 2;
  return v;
}
}  // namespace internal

SL int64 CalSigma0(int64 n) {
  if (n <= 0) return 0;
  if (n == 1) return 1;
  if (sigma0 && n <= maxp) return sigma0[n];

  return internal::CalSigma0Impl(n);
}

SL int64 CalSigma0(const FactorizedInteger<int64>& fi) {
  return fi.CalSigma0();
}

template <typename T>
SL T CalSigma0(const FactorizedInteger<T>& fi) {
  return fi.CalSigma0();
}

namespace internal {
SL int64 CalSigma1ByPmask(int64 n, int64 v = 1) {
  while (n != 1) {
    const int p = pmask[n];
    int64 d = 1;
    int64 s = 1;
    while (n % p == 0) n /= p, d *= p, s += d;
    if (d > 1) v *= s;
  }
  return v;
}

SL int64 CalSigma1Impl(int64 n) {
  int64 v = 1;

  for (int i = 0; i < pcnt; ++i) {
    if (n <= maxp) {
      return CalSigma1ByPmask(n, v);
    }
    const int64 p = plist[i];
    const int64 test = p * p;
    if (test > n) break;
    int64 d = 1;
    int64 s = 1;
    while (n % p == 0) n /= p, d *= p, s += d;
    if (d > 1) v *= s;
  }
  if (n != 1) v *= n + 1;
  return v;
}
}  // namespace internal

SL int64 CalSigma1(int64 n) {
  if (n <= 0) return 0;
  if (n == 1) return 1;
  if (sigma1 && n <= maxp) return sigma1[n];

  return internal::CalSigma1Impl(n);
}

SL int64 CalSigma1(const FactorizedInteger<int64>& fi) {
  return fi.CalSigma1();
}

template <typename T>
SL T CalSigma1(const FactorizedInteger<T>& fi) {
  return fi.CalSigma1();
}

namespace internal {
SL int IsSquareFreeByPmask(int64 n) {
  while (n != 1) {
    const int now = pmask[n];
    int c = 0;
    while (n % now == 0) n /= now, ++c;
    if (c > 1) return 0;
  }
  return 1;
}
}  // namespace internal

SL int IsSquareFree(int64 n) {
  if (n <= 0) {
    return 0;
  }

  for (int i = 0; i < pcnt; ++i) {
    if (n <= maxp) {
      return internal::IsSquareFreeByPmask(n);
    }
    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 > 1) return 0;
  }
  return 1;
}

SL int IsSquareFree(const FactorizedInteger<int64>& fi) {
  return fi.IsSquareFree();
}

template <typename T>
SL int IsSquareFree(const FactorizedInteger<T>& fi) {
  return fi.IsSquareFree();
}

template <typename T>
SL T CalMultiplicativeOrder(T n, T mod, T phi_mod = -1) {
  if (phi_mod == -1) {
    phi_mod = CalPhi(mod);
  }
  FactorizedInteger<T> fi(phi_mod);
  T ret = phi_mod;
  for (auto p : fi.GetPrimeFactors()) {
    if (PowerMod(n, phi_mod / p, mod) == 1) {
      ret /= p;
    }
  }
  return ret;
}

template <typename T>
SL T CalMultiplicativeOrder(T n, T mod,
                            const FactorizedInteger<T>& phi_mod_pi) {
  T phi_mod = phi_mod_pi.GetValue();
  T ret = phi_mod;
  for (auto p : phi_mod_pi.GetPrimeFactors()) {
    if (PowerMod(n, phi_mod / p, mod) == 1) {
      ret /= p;
    }
  }
  return ret;
}

SL int FindPrimeIdx(const std::function<int(int64)>& f) {
  for (int i = 0; i < pcnt; ++i) {
    if (f(plist[i])) return i;
  }
  return pcnt;
}

// sg: square greater
SL int FindPrimeIdxSg(int64 n) {
  int l = 0, r = pcnt - 1;
  while (l <= r) {
    const int mid = (l + r) >> 1;
    if (static_cast<int64>(plist[mid]) * plist[mid] > n) {
      r = mid - 1;
    } else {
      l = mid + 1;
    }
  }
  return l;
}

template <typename T = int64>
struct PrimeEnumeratorPe {
  PrimeEnumeratorPe(T start, T end = -1)
      : start_(start <= 1 ? 2 : start), end_(end) {}

  struct Iterator {
    using reference = T;
    using value_type = T;

    Iterator(T now, T end) : now_(now), end_(end) {}

    int operator==(const Iterator& o) const { return now_ == o.now_; }

    int operator!=(const Iterator& o) const { return !(this->operator==(o)); }

    T operator*() const { return now_; }

    Iterator& operator++() {
      if (now_ > 0) {
        ++now_;
        SeekNext();
      }
      return *this;
    }

    Iterator operator++(int) {
      Iterator r{now_, end_};
      ++*this;
      return r;
    }

    void SeekNext() {
      if (now_ > 0) {
        if (end_ > 0) {
          while (now_ <= end_ && !IsPrimeEx(now_)) {
            now_ += now_ & 1 ? 2 : 1;
          }
          if (now_ > end_) {
            now_ = -1;
          }
        } else {
          while (!IsPrimeEx(now_)) {
            now_ += now_ & 1 ? 2 : 1;
          }
        }
      }
    }

   private:
    T now_;
    T end_;
  };

  using iterator = Iterator;
  using const_iterator = Iterator;

  Iterator begin() const {
    Iterator ret{start_, end_};
    ret.SeekNext();
    return ret;
  }

  Iterator end() const { return {-1, end_}; }

  ContainerRange<Iterator> ToRange() const {
    return ContainerRange<Iterator>(begin(), end());
  }

 private:
  T start_;
  T end_;
};

template <typename T = int64>
SL std::vector<T> GetPrimesInRangePe(int64 start, int64 end) {
  if (start <= 1) {
    start = 2;
  }
  if (start > end) {
    return {};
  }
  std::vector<T> result;
  if (start == 2) {
    result.emplace_back(2);
    ++start;
  }
  start |= 1;
  for (int64 p = start; p <= end; p += 2) {
    if (IsPrimeEx(p)) {
      result.emplace_back(p);
    }
  }
  return result;
}

#if ENABLE_PRIME_SIEVE

template <typename T = int64>
struct PrimeEnumeratorPs {
  PrimeEnumeratorPs(T start, T end = -1)
      : start_(start <= 1 ? 2 : start),
        end_(end),
        it_(std::max(start_ - 1, 1LL),
            end < 0 ? primesieve::get_max_stop() : end) {}

  struct Iterator {
    using reference = T;
    using value_type = T;

    Iterator(T now, T end, primesieve::iterator* it)
        : now_(now), end_(end), it_(it) {}

    int operator==(const Iterator& o) const { return now_ == o.now_; }

    int operator!=(const Iterator& o) const { return !(this->operator==(o)); }

    T operator*() const { return now_; }

    Iterator& operator++() {
      if (now_ > 0) {
        ++now_;
        SeekNext();
      }
      return *this;
    }

    Iterator operator++(int) {
      Iterator r{now_, end_, it_};
      ++*this;
      return r;
    }

    void SeekNext() {
      if (now_ > 0) {
        now_ = it_->next_prime();
        if (end_ > 0 && now_ > end_) {
          now_ = -1;
        }
      }
    }

   private:
    T now_;
    T end_;
    primesieve::iterator* it_;
  };

  using iterator = Iterator;
  using const_iterator = Iterator;

  Iterator begin() const {
    Iterator ret{start_, end_, const_cast<primesieve::iterator*>(&it_)};
    ret.SeekNext();
    return ret;
  }

  Iterator end() const {
    return {-1, end_, const_cast<primesieve::iterator*>(&it_)};
  }

  ContainerRange<Iterator> ToRange() const {
    return ContainerRange<Iterator>(begin(), end());
  }

 private:
  T start_;
  T end_;
  primesieve::iterator it_;
};

template <typename T = int64>
SL std::vector<T> GetPrimesInRangePs(int64 start, int64 end) {
  if (start <= 1) {
    start = 2;
  }
  if (start > end) {
    return {};
  }

  std::vector<T> result;
  primesieve::iterator it(std::max(start - 1, 1LL), end);
  for (auto p = it.next_prime(); p <= end; p = it.next_prime()) {
    result.emplace_back(p);
  }
  return result;
}
#endif

#if ENABLE_PRIME_SIEVE
template <typename T = int64>
using PrimeEnumerator = PrimeEnumeratorPs<T>;

template <typename T = int64>
SL std::vector<T> GetPrimesInRange(int64 start, int64 end) {
  return GetPrimesInRangePs(start, end);
}
#else
template <typename T = int64>
using PrimeEnumerator = PrimeEnumeratorPe<T>;

template <typename T = int64>
SL std::vector<T> GetPrimesInRange(int64 start, int64 end) {
  return GetPrimesInRangePe(start, end);
}
#endif
}  // namespace pe
// 65701
#endif
/*

*/