目录
dsu
struct DSU {
std::vector<int> f, siz;
DSU(int n) : f(n), siz(n, 1) { std::iota(f.begin(), f.end(), 0); }
int leader(int x) {
while (x != f[x]) x = f[x] = f[f[x]];
return x;
}
bool same(int x, int y) { return leader(x) == leader(y); }
bool merge(int x, int y) {
x = leader(x);
y = leader(y);
if (x == y) return false;
siz[x] += siz[y];
f[y] = x;
return true;
}
int size(int x) { return siz[leader(x)]; }
};
组合数
const int N = 2e5 + 7;
const int mod = 998244353;
int fac[N], inv[N];
void init(int n) {
fac[0] = 1;
for (int i = 1; i <= n; ++i) fac[i] = fac[i - 1] * i % mod;
inv[n] = qpow(fac[n], mod - 2);
for (int i = n - 1; ~i; --i) inv[i] = inv[i + 1] * (i + 1) % mod;
}
int C(int n, int m) {
return n < m ? 0 : fac[n] * inv[m] % mod * inv[n - m] % mod;
}
取模
constexpr int P = 1000000007;
// assume -P <= x < 2P
int norm(int x) {
if (x < 0) {
x += P;
}
if (x >= P) {
x -= P;
}
return x;
}
template<class T>
T power(T a, int b) {
T res = 1;
for (; b; b /= 2, a *= a) {
if (b % 2) {
res *= a;
}
}
return res;
}
struct Z {
int x;
Z(int x = 0) : x(norm(x% P)) {}
int val() const {
return x;
}
Z operator-() const {
return Z(norm(P - x));
}
Z inv() const {
assert(x != 0);
return power(*this, P - 2);
}
Z& operator*=(const Z& rhs) {
x = x * rhs.x % P;
return *this;
}
Z& operator+=(const Z& rhs) {
x = norm(x + rhs.x);
return *this;
}
Z& operator-=(const Z& rhs) {
x = norm(x - rhs.x);
return *this;
}
Z& operator/=(const Z& rhs) {
return *this *= rhs.inv();
}
friend Z operator*(const Z& lhs, const Z& rhs) {
Z res = lhs;
res *= rhs;
return res;
}
friend Z operator+(const Z& lhs, const Z& rhs) {
Z res = lhs;
res += rhs;
return res;
}
friend Z operator-(const Z& lhs, const Z& rhs) {
Z res = lhs;
res -= rhs;
return res;
}
friend Z operator/(const Z& lhs, const Z& rhs) {
Z res = lhs;
res /= rhs;
return res;
}
friend std::istream& operator>>(std::istream& is, Z& a) {
int v;
is >> v;
a = Z(v);
return is;
}
friend std::ostream& operator<<(std::ostream& os, const Z& a) {
return os << a.val();
}
};
std::vector<Z> fac(N + 1), invfac(N + 1);
fac[0] = 1;
for (int i = 1; i <= N; i++) {
fac[i] = fac[i - 1] * i;
}
invfac[N] = fac[N].inv();
for (int i = N; i; i--) {
invfac[i - 1] = invfac[i] * i;
}
线段树
struct Max {
int x;
Max(int x = 0) : x(x) {}
};
Max operator+(const Max& a, const Max& b) {
return std::max(a.x, b.x);
}
template<class Info,
class Merge = std::plus<Info>>
struct SegmentTree {
const int n;
const Merge merge;
std::vector<Info> info;
SegmentTree(int n) : n(n), merge(Merge()), info(4 << std::__lg(n)) {}
SegmentTree(std::vector<Info> init) : SegmentTree(init.size()) {
std::function<void(int, int, int)> build = [&](int p, int l, int r) {
if (r - l == 1) {
info[p] = init[l];
return;
}
int m = (l + r) / 2;
build(2 * p, l, m);
build(2 * p + 1, m, r);
pull(p);
};
build(1, 0, n);
}
void pull(int p) {
info[p] = merge(info[2 * p], info[2 * p + 1]);
}
void modify(int p, int l, int r, int x, const Info& v) {
if (r - l == 1) {
info[p] = v;
return;
}
int m = (l + r) / 2;
if (x < m) {
modify(2 * p, l, m, x, v);
} else {
modify(2 * p + 1, m, r, x, v);
}
pull(p);
}
void modify(int p, const Info& v) {
modify(1, 0, n, p, v);
}
Info rangeQuery(int p, int l, int r, int x, int y) {
if (l >= y || r <= x) {
return Info();
}
if (l >= x && r <= y) {
return info[p];
}
int m = (l + r) / 2;
return merge(rangeQuery(2 * p, l, m, x, y), rangeQuery(2 * p + 1, m, r, x, y));
}
Info rangeQuery(int l, int r) {
return rangeQuery(1, 0, n, l, r);
}
};
lazy线段树
struct Max {
int x;
Max(int x = -inf) : x(x) {}
};
Max operator+(const Max& a, const Max& b) {
return std::max(a.x, b.x);
}
void apply(Max& a, const Max& b) {
a.x = std::max(a.x, b.x);
}
template<class Info, class Tag,
class Merge = std::plus<Info>>
struct LazySegmentTree {
const int n;
const Merge merge;
std::vector<Info> info;
std::vector<Tag> tag;
LazySegmentTree(int n) : n(n), merge(Merge()), info(4 << std::__lg(n)), tag(4 << std::__lg(n)) {}
LazySegmentTree(std::vector<Info> init) : LazySegmentTree(init.size()) {
std::function<void(int, int, int)> build = [&](int p, int l, int r) {
if (r - l == 1) {
info[p] = init[l];
return;
}
int m = (l + r) / 2;
build(2 * p, l, m);
build(2 * p + 1, m, r);
pull(p);
};
build(1, 0, n);
}
void pull(int p) {
info[p] = merge(info[2 * p], info[2 * p + 1]);
}
void apply(int p, const Tag& v) {
::apply(info[p], v);
::apply(tag[p], v);
}
void push(int p) {
apply(2 * p, tag[p]);
apply(2 * p + 1, tag[p]);
tag[p] = Tag();
}
void modify(int p, int l, int r, int x, const Info& v) {
if (r - l == 1) {
info[p] = v;
return;
}
int m = (l + r) / 2;
push(p);
if (x < m) {
modify(2 * p, l, m, x, v);
} else {
modify(2 * p + 1, m, r, x, v);
}
pull(p);
}
void modify(int p, const Info& v) {
modify(1, 0, n, p, v);
}
Info rangeQuery(int p, int l, int r, int x, int y) {
if (l >= y || r <= x) {
return Info();
}
if (l >= x && r <= y) {
return info[p];
}
int m = (l + r) / 2;
push(p);
return merge(rangeQuery(2 * p, l, m, x, y), rangeQuery(2 * p + 1, m, r, x, y));
}
Info rangeQuery(int l, int r) {
return rangeQuery(1, 0, n, l, r);
}
void rangeApply(int p, int l, int r, int x, int y, const Tag& v) {
if (l >= y || r <= x) {
return;
}
if (l >= x && r <= y) {
apply(p, v);
return;
}
int m = (l + r) / 2;
push(p);
rangeApply(2 * p, l, m, x, y, v);
rangeApply(2 * p + 1, m, r, x, y, v);
pull(p);
}
void rangeApply(int l, int r, const Tag& v) {
return rangeApply(1, 0, n, l, r, v);
}
};
2-sat
struct TwoSat {
int n;
std::vector<std::vector<int>> e;
std::vector<bool> ans;
TwoSat(int n) : n(n), e(2 * n), ans(n) {}
void addClause(int u, bool f, int v, bool g) {
e[2 * u + !f].push_back(2 * v + g);
e[2 * v + !g].push_back(2 * u + f);
}
bool satisfiable() {
std::vector<int> id(2 * n, -1), dfn(2 * n, -1), low(2 * n, -1);
std::vector<int> stk;
int now = 0, cnt = 0;
std::function<void(int)> tarjan = [&](int u) {
stk.push_back(u);
dfn[u] = low[u] = now++;
for (auto v : e[u]) {
if (dfn[v] == -1) {
tarjan(v);
low[u] = std::min(low[u], low[v]);
} else if (id[v] == -1) {
low[u] = std::min(low[u], dfn[v]);
}
}
if (dfn[u] == low[u]) {
int v;
do {
v = stk.back();
stk.pop_back();
id[v] = cnt;
} while (v != u);
++cnt;
}
};
for (int i = 0; i < 2 * n; ++i) if (dfn[i] == -1) tarjan(i);
for (int i = 0; i < n; ++i) {
if (id[2 * i] == id[2 * i + 1]) return false;
ans[i] = id[2 * i] > id[2 * i + 1];
}
return true;
}
std::vector<bool> answer() { return ans; }
};
树状数组
template<typename T>
struct Fenwick {
int n;
std::vector<T> a;
Fenwick(int n) : n(n), a(n + 1) {}
void add(int x, T v) {
for (int i = x; i <= n; i += i & -i) {
a[i] += v;
}
}
T sum(int x) {
T ans = 0;
for (int i = x; i > 0; i -= i & -i) {
ans += a[i];
}
return ans;
}
T rangeSum(int l, int r) {
return sum(r) - sum(l);
}
T find_kth(int k) {
T ans = 0, cnt = 0;
for (int i = 20; i >= 0; i--) {
ans += (1 << i);
if (ans > n || cnt + a[ans] >= k) {
ans -= (1 << i);
} else {
cnt += a[ans];
}
}
return ++ans;
}
};
FFT
constexpr int P = 998244353;
// assume -P <= x < 2P
int norm(int x) {
if (x < 0) {
x += P;
}
if (x >= P) {
x -= P;
}
return x;
}
template<class T>
T power(T a, int b) {
T res = 1;
for (; b; b /= 2, a *= a) {
if (b % 2) {
res *= a;
}
}
return res;
}
struct Z {
int x;
Z(int x = 0) : x(norm(x)) {}
int val() const {
return x;
}
Z operator-() const {
return Z(norm(P - x));
}
Z inv() const {
assert(x != 0);
return power(*this, P - 2);
}
Z& operator*=(const Z& rhs) {
x = x * rhs.x % P;
return *this;
}
Z& operator+=(const Z& rhs) {
x = norm(x + rhs.x);
return *this;
}
Z& operator-=(const Z& rhs) {
x = norm(x - rhs.x);
return *this;
}
Z& operator/=(const Z& rhs) {
return *this *= rhs.inv();
}
friend Z operator*(const Z& lhs, const Z& rhs) {
Z res = lhs;
res *= rhs;
return res;
}
friend Z operator+(const Z& lhs, const Z& rhs) {
Z res = lhs;
res += rhs;
return res;
}
friend Z operator-(const Z& lhs, const Z& rhs) {
Z res = lhs;
res -= rhs;
return res;
}
friend Z operator/(const Z& lhs, const Z& rhs) {
Z res = lhs;
res /= rhs;
return res;
}
friend std::istream& operator>>(std::istream& is, Z& a) {
int v;
is >> v;
a = Z(v);
return is;
}
friend std::ostream& operator<<(std::ostream& os, const Z& a) {
return os << a.val();
}
};
std::vector<int> rev;
std::vector<Z> roots{ 0, 1 };
void dft(std::vector<Z>& a) {
int n = a.size();
if ((int) rev.size() != n) {
int k = __builtin_ctz(n) - 1;
rev.resize(n);
for (int i = 0; i < n; i++) {
rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
}
}
for (int i = 0; i < n; i++) {
if (rev[i] < i) {
std::swap(a[i], a[rev[i]]);
}
}
if ((int) roots.size() < n) {
int k = __builtin_ctz(roots.size());
roots.resize(n);
while ((1 << k) < n) {
Z e = power(Z(3), (P - 1) >> (k + 1));
for (int i = 1 << (k - 1); i < (1 << k); i++) {
roots[2 * i] = roots[i];
roots[2 * i + 1] = roots[i] * e;
}
k++;
}
}
for (int k = 1; k < n; k *= 2) {
for (int i = 0; i < n; i += 2 * k) {
for (int j = 0; j < k; j++) {
Z u = a[i + j];
Z v = a[i + j + k] * roots[k + j];
a[i + j] = u + v;
a[i + j + k] = u - v;
}
}
}
}
void idft(std::vector<Z>& a) {
int n = a.size();
std::reverse(a.begin() + 1, a.end());
dft(a);
Z inv = (1 - P) / n;
for (int i = 0; i < n; i++) {
a[i] *= inv;
}
}
struct Poly {
std::vector<Z> a;
Poly() {}
Poly(const std::vector<Z>& a) : a(a) {}
Poly(const std::initializer_list<Z>& a) : a(a) {}
int size() const {
return a.size();
}
void resize(int n) {
a.resize(n);
}
Z operator[](int idx) const {
if (idx < size()) {
return a[idx];
} else {
return 0;
}
}
Z& operator[](int idx) {
return a[idx];
}
Poly mulxk(int k) const {
auto b = a;
b.insert(b.begin(), k, 0);
return Poly(b);
}
Poly modxk(int k) const {
k = std::min(k, size());
return Poly(std::vector<Z>(a.begin(), a.begin() + k));
}
Poly divxk(int k) const {
if (size() <= k) {
return Poly();
}
return Poly(std::vector<Z>(a.begin() + k, a.end()));
}
friend Poly operator+(const Poly& a, const Poly& b) {
std::vector<Z> res(std::max(a.size(), b.size()));
for (int i = 0; i < (int) res.size(); i++) {
res[i] = a[i] + b[i];
}
return Poly(res);
}
friend Poly operator-(const Poly& a, const Poly& b) {
std::vector<Z> res(std::max(a.size(), b.size()));
for (int i = 0; i < (int) res.size(); i++) {
res[i] = a[i] - b[i];
}
return Poly(res);
}
friend Poly operator*(Poly a, Poly b) {
if (a.size() == 0 || b.size() == 0) {
return Poly();
}
int sz = 1, tot = a.size() + b.size() - 1;
while (sz < tot) {
sz *= 2;
}
a.a.resize(sz);
b.a.resize(sz);
dft(a.a);
dft(b.a);
for (int i = 0; i < sz; ++i) {
a.a[i] = a[i] * b[i];
}
idft(a.a);
a.resize(tot);
return a;
}
friend Poly operator*(Z a, Poly b) {
for (int i = 0; i < (int) b.size(); i++) {
b[i] *= a;
}
return b;
}
friend Poly operator*(Poly a, Z b) {
for (int i = 0; i < (int) a.size(); i++) {
a[i] *= b;
}
return a;
}
Poly& operator+=(Poly b) {
return (*this) = (*this) + b;
}
Poly& operator-=(Poly b) {
return (*this) = (*this) - b;
}
Poly& operator*=(Poly b) {
return (*this) = (*this) * b;
}
Poly deriv() const {
if (a.empty()) {
return Poly();
}
std::vector<Z> res(size() - 1);
for (int i = 0; i < size() - 1; ++i) {
res[i] = (i + 1) * a[i + 1];
}
return Poly(res);
}
Poly integr() const {
std::vector<Z> res(size() + 1);
for (int i = 0; i < size(); ++i) {
res[i + 1] = a[i] / (i + 1);
}
return Poly(res);
}
Poly inv(int m) const {
Poly x{ a[0].inv() };
int k = 1;
while (k < m) {
k *= 2;
x = (x * (Poly{ 2 } - modxk(k) * x)).modxk(k);
}
return x.modxk(m);
}
Poly log(int m) const {
return (deriv() * inv(m)).integr().modxk(m);
}
Poly exp(int m) const {
Poly x{ 1 };
int k = 1;
while (k < m) {
k *= 2;
x = (x * (Poly{ 1 } - x.log(k) + modxk(k))).modxk(k);
}
return x.modxk(m);
}
Poly pow(int k, int m) const {
int i = 0;
while (i < size() && a[i].val() == 0) {
i++;
}
if (i == size() || 1LL * i * k >= m) {
return Poly(std::vector<Z>(m));
}
Z v = a[i];
auto f = divxk(i) * v.inv();
return (f.log(m - i * k) * k).exp(m - i * k).mulxk(i * k) * power(v, k);
}
Poly sqrt(int m) const {
Poly x{ 1 };
int k = 1;
while (k < m) {
k *= 2;
x = (x + (modxk(k) * x.inv(k)).modxk(k)) * ((P + 1) / 2);
}
return x.modxk(m);
}
Poly mulT(Poly b) const {
if (b.size() == 0) {
return Poly();
}
int n = b.size();
std::reverse(b.a.begin(), b.a.end());
return ((*this) * b).divxk(n - 1);
}
std::vector<Z> eval(std::vector<Z> x) const {
if (size() == 0) {
return std::vector<Z>(x.size(), 0);
}
const int n = std::max((int) x.size(), size());
std::vector<Poly> q(4 * n);
std::vector<Z> ans(x.size());
x.resize(n);
std::function<void(int, int, int)> build = [&](int p, int l, int r) {
if (r - l == 1) {
q[p] = Poly{ 1, -x[l] };
} else {
int m = (l + r) / 2;
build(2 * p, l, m);
build(2 * p + 1, m, r);
q[p] = q[2 * p] * q[2 * p + 1];
}
};
build(1, 0, n);
std::function<void(int, int, int, const Poly&)> work = [&](int p, int l, int r, const Poly& num) {
if (r - l == 1) {
if (l < (int) ans.size()) {
ans[l] = num[0];
}
} else {
int m = (l + r) / 2;
work(2 * p, l, m, num.mulT(q[2 * p + 1]).modxk(m - l));
work(2 * p + 1, m, r, num.mulT(q[2 * p]).modxk(r - m));
}
};
work(1, 0, n, mulT(q[1].inv(n)));
return ans;
}
};
//cdq + fft
//vector<Z> f, c, p, invp, w, sumw, invsumw;
// void cdq(int l, int r) {
// if (l == r) {
// f[l] += (f[l - 1] + c[l - 1]) * invp[l - 1];
// return;
// }
// int mid = l + r >> 1;
// cdq(l, mid);
// vector<Z> a(r - l + 2), b(r - l + 2);
// for (int i = l; i <= r; i++) {
// a[i - l] = w[i - l];
// b[i - l + 1] = (i <= mid) ? f[i] : 0;
// }
// auto res = Poly(a) * Poly(b);
// for (int i = mid + 1; i <= r; i++) {
// Z tmp = invp[i - 1] * invsumw[i - 1];
// tmp *= (1 - p[i - 1] + P);
// tmp *= res[i - l];
// f[i] -= tmp;
// }
// cdq(mid + 1, r);
// }
std::vector<Z> fac(N + 1), invfac(N + 1);
fac[0] = 1;
for (int i = 1; i <= N; i++) {
fac[i] = fac[i - 1] * i;
}
invfac[N] = fac[N].inv();
for (int i = N; i; i--) {
invfac[i - 1] = invfac[i] * i;
}
SAM
struct SuffixAutomaton {
static constexpr int ALPHABET_SIZE = 26, N = 1e5;
struct Node {
int len;
int link;
int next[ALPHABET_SIZE];
Node() : len(0), link(0), next{} {}
} t[2 * N];
int cntNodes;
SuffixAutomaton() {
cntNodes = 1;
std::fill(t[0].next, t[0].next + ALPHABET_SIZE, 1);
t[0].len = -1;
}
int extend(int p, int c) {
if (t[p].next[c]) {
int q = t[p].next[c];
if (t[q].len == t[p].len + 1)
return q;
int r = ++cntNodes;
t[r].len = t[p].len + 1;
t[r].link = t[q].link;
std::copy(t[q].next, t[q].next + ALPHABET_SIZE, t[r].next);
t[q].link = r;
while (t[p].next[c] == q) {
t[p].next[c] = r;
p = t[p].link;
}
return r;
}
int cur = ++cntNodes;
t[cur].len = t[p].len + 1;
while (!t[p].next[c]) {
t[p].next[c] = cur;
p = t[p].link;
}
t[cur].link = extend(p, c);
return cur;
}
};
Hash
const int L = 1e6 + 5;
const int HASH_CNT = 2;
int hashBase[HASH_CNT] = { 29, 31 };
int hashMod[HASH_CNT] = { (int) 1e9 + 9, 998244353 };
struct StringWithHash {
char s[L];
int ls;
int hsh[HASH_CNT][L];
int pwMod[HASH_CNT][L];
void init() {
ls = 0;
for (int i = 0; i < HASH_CNT; ++i) {
hsh[i][0] = 0;
pwMod[i][0] = 1;
}
}
StringWithHash() { init(); }
void extend(char c) {
s[++ls] = c;
for (int i = 0; i < HASH_CNT; ++i) {
pwMod[i][ls] =
1ll * pwMod[i][ls - 1] * hashBase[i] % hashMod[i];
hsh[i][ls] = (1ll * hsh[i][ls - 1] * hashBase[i] + c) % hashMod[i];
}
}
vector<int> getHash(int l, int r) {
vector<int> res(HASH_CNT, 0);
for (int i = 0; i < HASH_CNT; ++i) {
int t =
(hsh[i][r] - 1ll * hsh[i][l - 1] * pwMod[i][r - l + 1]) % hashMod[i];
t = (t + hashMod[i]) % hashMod[i];
res[i] = t;
}
return res;
}
};
bool equal(const vector<int>& h1, const vector<int>& h2) {
assert(h1.size() == h2.size());
for (unsigned i = 0; i < h1.size(); i++)
if (h1[i] != h2[i]) return false;
return true;
}
int128
using i128 = __int128;
std::ostream &operator<<(std::ostream &o, i128 n) {
if (!n) {
return o << 0;
} else {
std::string s;
while (n) {
s += n % 10 + '0';
n /= 10;
}
std::reverse(s.begin(), s.end());
return o << s;
}
}
i128 gcd(i128 a, i128 b) {
while (b) {
a %= b;
std::swap(a, b);
}
return a;
}
NTT
const int MOD = 998244353;
struct mod_int {
int val;
mod_int(long long v = 0) {
if (v < 0)
v = v % MOD + MOD;
if (v >= MOD)
v %= MOD;
val = v;
}
static int mod_inv(int a, int m = MOD) {
// https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Example
int g = m, r = a, x = 0, y = 1;
while (r != 0) {
int q = g / r;
g %= r; swap(g, r);
x -= q * y; swap(x, y);
}
return x < 0 ? x + m : x;
}
explicit operator int() const {
return val;
}
mod_int& operator+=(const mod_int& other) {
val += other.val;
if (val >= MOD) val -= MOD;
return *this;
}
mod_int& operator-=(const mod_int& other) {
val -= other.val;
if (val < 0) val += MOD;
return *this;
}
static unsigned fast_mod(uint64_t x, unsigned m = MOD) {
/*#if !defined(_WIN32) || defined(_WIN64)
return x % m;
#endif
// Optimized mod for Codeforces 32-bit machines.
// x must be less than 2^32 * m for this to work, so that x / m fits in a 32-bit integer.
unsigned x_high = x >> 32, x_low = (unsigned) x;
unsigned quot, rem;
asm("divl %4\n"
: "=a" (quot), "=d" (rem)
: "d" (x_high), "a" (x_low), "r" (m));
return rem;*/
return x % m;
}
mod_int& operator*=(const mod_int& other) {
val = fast_mod((uint64_t) val * other.val);
return *this;
}
mod_int& operator/=(const mod_int& other) {
return *this *= other.inv();
}
friend mod_int operator+(const mod_int& a, const mod_int& b) { return mod_int(a) += b; }
friend mod_int operator-(const mod_int& a, const mod_int& b) { return mod_int(a) -= b; }
friend mod_int operator*(const mod_int& a, const mod_int& b) { return mod_int(a) *= b; }
friend mod_int operator/(const mod_int& a, const mod_int& b) { return mod_int(a) /= b; }
mod_int& operator++() {
val = val == MOD - 1 ? 0 : val + 1;
return *this;
}
mod_int& operator--() {
val = val == 0 ? MOD - 1 : val - 1;
return *this;
}
mod_int operator++(int) { mod_int before = *this; ++* this; return before; }
mod_int operator--(int) { mod_int before = *this; --* this; return before; }
mod_int operator-() const {
return val == 0 ? 0 : MOD - val;
}
bool operator==(const mod_int& other) const { return val == other.val; }
bool operator!=(const mod_int& other) const { return val != other.val; }
mod_int inv() const {
return mod_inv(val);
}
mod_int pow(long long p) const {
assert(p >= 0);
mod_int a = *this, result = 1;
while (p > 0) {
if (p & 1)
result *= a;
a *= a;
p >>= 1;
}
return result;
}
friend ostream& operator<<(ostream& stream, const mod_int& m) {
return stream << m.val;
}
};
void debug(vector<mod_int> a) {
for (auto i : a)cerr << i << ' ';
cerr << endl;
}
void r_debug(vector<mod_int> a) {
for (auto i : a)cerr << 1 / i << ' ';
cerr << endl;
}
vector<mod_int> inv, factorial, inv_factorial;
void prepare_factorials(int maximum) {
inv.assign(maximum + 1, 1);
// Make sure MOD is prime, which is necessary for the inverse algorithm below.
for (int p = 2; p * p <= MOD; p++)
assert(MOD % p != 0);
for (int i = 2; i <= maximum; i++)
inv[i] = inv[MOD % i] * (MOD - MOD / i);
factorial.resize(maximum + 1);
inv_factorial.resize(maximum + 1);
factorial[0] = inv_factorial[0] = 1;
for (int i = 1; i <= maximum; i++) {
factorial[i] = i * factorial[i - 1];
inv_factorial[i] = inv[i] * inv_factorial[i - 1];
}
}
inline mod_int C(int n, int m) {
assert(n >= m); assert(n < factorial.size());
return factorial[n] * inv_factorial[m] * inv_factorial[n - m];
}
namespace NTT {
vector<mod_int> roots = { 0, 1 };
vector<int> bit_reverse;
int max_size = -1;
mod_int root;
bool is_power_of_two(int n) {
return (n & (n - 1)) == 0;
}
int round_up_power_two(int n) {
assert(n > 0);
while (n & (n - 1))
n = (n | (n - 1)) + 1;
return n;
}
// Given n (a power of two), finds k such that n == 1 << k.
int get_length(int n) {
assert(is_power_of_two(n));
return __builtin_ctz(n);
}
// Rearranges the indices to be sorted by lowest bit first, then second lowest, etc., rather than highest bit first.
// This makes even-odd div-conquer much easier.
void bit_reorder(int n, vector<mod_int>& values) {
if ((int) bit_reverse.size() != n) {
bit_reverse.assign(n, 0);
int length = get_length(n);
for (int i = 0; i < n; i++)
bit_reverse[i] = (bit_reverse[i >> 1] >> 1) + ((i & 1) << (length - 1));
}
for (int i = 0; i < n; i++)
if (i < bit_reverse[i])
swap(values[i], values[bit_reverse[i]]);
}
void find_root() {
int order = MOD - 1;
max_size = 1;
while (order % 2 == 0) {
order /= 2;
max_size *= 2;
}
root = 2;
// Find a max_size-th primitive root of MOD.
while (!(root.pow(max_size) == 1 && root.pow(max_size / 2) != 1))
root++;
}
void prepare_roots(int n) {
if (max_size < 0)
find_root();
assert(n <= max_size);
if ((int) roots.size() >= n)
return;
int length = get_length(roots.size());
roots.resize(n);
// The roots array is set up such that for a given power of two n >= 2, roots[n / 2] through roots[n - 1] are
// the first half of the n-th primitive roots of MOD.
while (1 << length < n) {
// z is a 2^(length + 1)-th primitive root of MOD.
mod_int z = root.pow(max_size >> (length + 1));
for (int i = 1 << (length - 1); i < 1 << length; i++) {
roots[2 * i] = roots[i];
roots[2 * i + 1] = roots[i] * z;
}
length++;
}
}
void fft_iterative(int N, vector<mod_int>& values) {
assert(is_power_of_two(N));
prepare_roots(N);
bit_reorder(N, values);
for (int n = 1; n < N; n *= 2)
for (int start = 0; start < N; start += 2 * n)
for (int i = 0; i < n; i++) {
mod_int even = values[start + i];
mod_int odd = values[start + n + i] * roots[n + i];
values[start + n + i] = even - odd;
values[start + i] = even + odd;
}
}
const int FFT_CUTOFF = 150;
vector<mod_int> mod_multiply(vector<mod_int> left, vector<mod_int> right) {
int n = left.size();
int m = right.size();
// Brute force when either n or m is small enough.
if (min(n, m) < FFT_CUTOFF) {
const uint64_t ULL_BOUND = numeric_limits<uint64_t>::max() - (uint64_t) MOD * MOD;
vector<uint64_t> result(n + m - 1);
for (int i = 0; i < n; i++)
for (int j = 0; j < m; j++) {
result[i + j] += (uint64_t) ((int) left[i]) * ((int) right[j]);
if (result[i + j] > ULL_BOUND)
result[i + j] %= MOD;
}
for (uint64_t& x : result)
if (x >= MOD)
x %= MOD;
return vector<mod_int>(result.begin(), result.end());
}
int N = round_up_power_two(n + m - 1);
left.resize(N);
right.resize(N);
bool equal = left == right;
fft_iterative(N, left);
if (equal)
right = left;
else
fft_iterative(N, right);
mod_int inv_N = mod_int(N).inv();
for (int i = 0; i < N; i++)
left[i] *= right[i] * inv_N;
reverse(left.begin() + 1, left.end());
fft_iterative(N, left);
left.resize(n + m - 1);
return left;
}
vector<mod_int> mod_multiply_all(const vector<vector<mod_int>>& polynomials) {
if (polynomials.empty())
return { 1 };
struct compare_size {
bool operator()(const vector<mod_int>& x, const vector<mod_int>& y) {
return x.size() > y.size();
}
};
priority_queue<vector<mod_int>, vector<vector<mod_int>>, compare_size> pq(polynomials.begin(), polynomials.end());
while (pq.size() > 1) {
vector<mod_int> a = pq.top(); pq.pop();
vector<mod_int> b = pq.top(); pq.pop();
pq.push(mod_multiply(a, b));
}
return pq.top();
}
}