namespace mathtool {
// modint referenced from jiangly
template<typename T>
constexpr T power(T a, u64 b) {
T res {1};
for (; b != 0; b /= 2, a *= a) {
if (b % 2 == 1) {
res *= a;
}
}
return res;
}
template<u32 P>
constexpr u32 mulMod(u32 a, u32 b) {
return 1ULL * a * b % P;
}
template<u64 P>
constexpr u64 mulMod(u64 a, u64 b) {
u64 res = a * b - u64(1.L * a * b / P - 0.5L) * P;
res %= P;
return res;
}
template<typename U, U P>
requires std::unsigned_integral<U>
struct ModIntBase {
public:
constexpr ModIntBase() : x {0} {}
template<typename T>
requires std::integral<T>
constexpr ModIntBase(T x_) : x {norm(x_ % T {P})} {}
constexpr static U norm(U x) {
if ((x >> (8 * sizeof(U) - 1) & 1) == 1) {
x += P;
}
if (x >= P) {
x -= P;
}
return x;
}
constexpr U val() const {
return x;
}
constexpr ModIntBase operator-() const {
ModIntBase res;
res.x = norm(P - x);
return res;
}
constexpr ModIntBase inv() const {
return power(*this, P - 2);
}
constexpr ModIntBase &operator*=(const ModIntBase &rhs) & {
x = mulMod<P>(x, rhs.val());
return *this;
}
constexpr ModIntBase &operator+=(const ModIntBase &rhs) & {
x = norm(x + rhs.x);
return *this;
}
constexpr ModIntBase &operator-=(const ModIntBase &rhs) & {
x = norm(x - rhs.x);
return *this;
}
constexpr ModIntBase &operator/=(const ModIntBase &rhs) & {
return *this *= rhs.inv();
}
friend constexpr ModIntBase operator*(ModIntBase lhs, const ModIntBase
&rhs) {
lhs *= rhs;
return lhs;
}
friend constexpr ModIntBase operator+(ModIntBase lhs, const ModIntBase
&rhs) {
lhs += rhs;
return lhs;
}
friend constexpr ModIntBase operator-(ModIntBase lhs, const ModIntBase
&rhs) {
lhs -= rhs;
return lhs;
}
friend constexpr ModIntBase operator/(ModIntBase lhs, const ModIntBase
&rhs) {
lhs /= rhs;
return lhs;
}
friend constexpr std::ostream &operator<<(std::ostream &os, const
ModIntBase &a) {
return os << a.val();
}
friend constexpr bool operator==(ModIntBase lhs, ModIntBase rhs) {
return lhs.val() == rhs.val();
}
friend constexpr bool operator!=(ModIntBase lhs, ModIntBase rhs) {
return lhs.val() != rhs.val();
}
friend constexpr bool operator<(ModIntBase lhs, ModIntBase rhs) {
return lhs.val() < rhs.val();
}
private:
U x;
};
template<u32 P>
using ModInt = ModIntBase<u32, P>;
template<u64 P>
using ModInt64 = ModIntBase<u64, P>;
constexpr u32 P = 1000000007;
using Z = ModInt<P>;
template<class T>
T fpow(T a, u64 b, const LL MOD) {
T res{1};
while (b) {
if (b & 1) res = (res * a) % MOD;
a = (a * a) % MOD;
b >>= 1;
}
return res;
}
void printi28(i28 x) {
string s;
while (x > 0) {
s += (x % 10 + '0');
x /= 10;
}
reverse(all(s));
W(s);
}
struct chinese_remainder_theorem {
i28 solve(vector<i28>& a, vector<i28>& M) {
i28 m = 1;
for (const i28& p : M) {
m *= p;
}
//printi28(m);
i28 res = 0;
for (int i = 0; i < si(a); ++i) {
i28 mi = m / M[i];
res = (res + a[i] * mi % m * fpow(mi, M[i] - 2, M[i])) % m;
}
return res;
}
} CRT;
// comb referenced from jiangly
struct Comb {
int n;
std::vector<Z> _fac;
std::vector<Z> _invfac;
std::vector<Z> _inv;
Comb() : n{0}, _fac{1}, _invfac{1}, _inv{0} {}
Comb(int n) : Comb() {
init(n);
}
void init(int m) {
// m = std::min<i64>(m, Z::getMod() - 1);
if (m <= n) return;
_fac.resize(m + 1);
_invfac.resize(m + 1);
_inv.resize(m + 1);
for (int i = n + 1; i <= m; i++) {
_fac[i] = _fac[i - 1] * i;
}
_invfac[m] = _fac[m].inv();
for (int i = m; i > n; i--) {
_invfac[i - 1] = _invfac[i] * i;
_inv[i] = _invfac[i] * _fac[i - 1];
}
n = m;
}
Z fac(int m) {
if (m > n) init(2 * m);
return _fac[m];
}
Z invfac(int m) {
if (m > n) init(2 * m);
return _invfac[m];
}
Z inv(int m) {
if (m > n) init(2 * m);
return _inv[m];
}
Z binom(int n, int m) {
if (n < m || m < 0) return 0;
return fac(n) * invfac(m) * invfac(n - m);
}
} comb;
// need to be initialized
struct PrimTable {
int n;
vector<int> pr, lp;
inline PrimTable(int n_ = 0) {
n = n_;
lp.resize(n + 1);
init();
}
inline void init() {
if (n < 1) return;
lp[1] = 1;
for (int i = 2; i <= n; ++i) {
if (lp[i] == 0) {
lp[i] = i;
pr.push_back(i);
}
for (int j = 0; i * pr[j] <= n; ++j) {
lp[i * pr[j]] = pr[j];
if (pr[j] == lp[i]) {
break;
}
}
}
}
inline bool isp(int x) {
return x < 2 ? false : (lp[x] == x);
}
inline vector<int> div(int k) {
assert(k <= n);
if (k < 2) return vector<int>();
vector<int> res;
while (k > 1) {
const int f = lp[k];
do {
k /= f;
} while (k % f == 0);
res.push_back(f);
}
return res;
}
};
}
using namespace mathtool;
