#include #include "../numeric/fft.h" #include "modint.h" using namespace std; // https://cp-algorithms.com/algebra/polynomial.html template vector &operator+=(vector &a, const vector &b) { if (a.size() < b.size()) { a.resize(b.size()); } for (size_t i = 0; i < b.size(); i++) { a[i] += b[i]; } return a; } template vector operator+(vector a, const vector &b) { a += b; return a; } template vector &operator-=(vector &a, const vector &b) { if (a.size() < b.size()) { a.resize(b.size()); } for (int i = 0; i < b.size(); i++) { a[i] -= b[i]; } return a; } template vector operator-(vector a, const vector &b) { a -= b; return a; } template vector operator-(vector a) { for (int i = 0; i < a.size(); i++) { a[i] = -a[i]; } return a; } // fast multiplication for modint in O(n*log(n)) template vector> operator*(const vector> &a, const vector> &b) { if (a.empty() || b.empty()) { return {}; } if (min(a.size(), b.size()) < 150) { vector> c(a.size() + b.size() - 1, 0); for (size_t i = 0; i < a.size(); i++) { for (size_t j = 0; j < b.size(); j++) { c[i + j] += a[i] * b[j]; } } return c; } vector a_int(a.size()); for (size_t i = 0; i < a.size(); i++) { a_int[i] = static_cast(a[i]); } vector b_int(b.size()); for (size_t i = 0; i < b.size(); i++) { b_int[i] = static_cast(b[i]); } vector c_int = multiply_mod(a_int, b_int, mod); vector> c(c_int.size()); for (size_t i = 0; i < c.size(); i++) { c[i] = c_int[i]; } return c; } // fallback multiplication in O(n*m) template vector operator*(const vector &a, const vector &b) { if (a.empty() || b.empty()) { return {}; } vector c(a.size() + b.size() - 1, 0); for (int i = 0; i < a.size(); i++) { for (int j = 0; j < b.size(); j++) { c[i + j] += a[i] * b[j]; } } return c; } template vector &operator*=(vector &a, const vector &b) { a = a * b; return a; } template vector inverse(const vector &a) { assert(!a.empty()); int n = a.size(); vector b = {1 / a[0]}; while (b.size() < n) { vector a_cut(a.begin(), a.begin() + min(a.size(), b.size() << 1)); vector x = b * b * a_cut; b.resize(b.size() << 1); for (size_t i = b.size() >> 1; i < min(x.size(), b.size()); i++) { b[i] = -x[i]; } } b.resize(n); return b; } template vector &operator/=(vector &a, vector b) { int n = a.size(); int m = b.size(); if (n < m) { a.clear(); } else { reverse(a.begin(), a.end()); reverse(b.begin(), b.end()); b.resize(n - m + 1); a *= inverse(b); a.erase(a.begin() + n - m + 1, a.end()); reverse(a.begin(), a.end()); } return a; } template vector operator/(vector a, const vector &b) { a /= b; return a; } template vector &operator%=(vector &a, const vector &b) { int n = a.size(); int m = b.size(); if (n >= m) { vector c = (a / b) * b; a.resize(m - 1); for (int i = 0; i < m - 1; i++) { a[i] -= c[i]; } } return a; } template vector operator%(vector a, const vector &b) { a %= b; return a; } template vector power(const vector &a, long long b, const vector &mod) { assert(b >= 0); vector res = vector{1} % mod; int highest_one_bit = -1; while (1LL << (highest_one_bit + 1) <= b) ++highest_one_bit; for (int i = highest_one_bit; i >= 0; i--) { res = res * res % mod; if (b >> i & 1) { res = res * a % mod; } } return res; } template vector derivative(vector a) { for (size_t i = 0; i < a.size(); i++) { a[i] *= i; } if (!a.empty()) { a.erase(a.begin()); } return a; } template vector integrate(vector a) { a.insert(a.begin(), 0); for (int i = 1; i < a.size(); i++) { a[i] /= i; } return a; } template vector logarithm(const vector &a) { assert(!a.empty() && a[0] == 1); vector res = integrate(derivative(a) * inverse(a)); res.resize(a.size()); return res; } template vector exponent(const vector &a) { assert(!a.empty() && a[0] == 0); int n = a.size(); vector b = {1}; while (b.size() < n) { vector x(a.begin(), a.begin() + min(a.size(), b.size() << 1)); x[0] += 1; vector old_b = b; b.resize(b.size() << 1); x -= logarithm(b); x *= old_b; for (int i = b.size() >> 1; i < min(x.size(), b.size()); i++) { b[i] = x[i]; } } b.resize(n); return b; } template vector sqrt(const vector &a) { assert(!a.empty() && a[0] == 1); int n = a.size(); vector b = {1}; while (b.size() < n) { vector x(a.begin(), a.begin() + min(a.size(), b.size() << 1)); b.resize(b.size() << 1); x *= inverse(b); T inv2 = 1 / static_cast(2); for (int i = b.size() >> 1; i < min(x.size(), b.size()); i++) { b[i] = x[i] * inv2; } } b.resize(n); return b; } template vector multiply(const vector> &a) { if (a.empty()) { return {0}; } function(int, int)> mult = [&](int l, int r) { if (l == r) { return a[l]; } int m = (l + r) >> 1; return mult(l, m) * mult(m + 1, r); }; return mult(0, (int)a.size() - 1); } template T evaluate(const vector &a, const T &x) { T res = 0; for (int i = (int)a.size() - 1; i >= 0; i--) { res = res * x + a[i]; } return res; } template vector evaluate(const vector &a, const vector &x) { if (x.empty()) { return {}; } if (a.empty()) { return vector(x.size()); } int n = x.size(); vector> t(2 * n - 1); function build = [&](int v, int l, int r) { if (l == r) { t[v] = vector{-x[l], 1}; } else { int m = (l + r) >> 1; int y = v + ((m - l + 1) << 1); build(v + 1, l, m); build(y, m + 1, r); t[v] = t[v + 1] * t[y]; } }; build(0, 0, n - 1); vector res(n); function)> eval = [&](int v, int l, int r, vector p) { p %= t[v]; if (p.size() < 150) { for (int i = l; i <= r; i++) { res[i] = evaluate(p, x[i]); } return; } if (l == r) { res[l] = p[0]; } else { int m = (l + r) >> 1; int z = v + ((m - l + 1) << 1); eval(v + 1, l, m, p); eval(z, m + 1, r, p); } }; eval(0, 0, n - 1, a); return res; } template vector interpolate(const vector &x, const vector &y) { if (x.empty()) { return {}; } assert(x.size() == y.size()); int n = x.size(); vector> t(2 * n - 1); function build = [&](int v, int l, int r) { if (l == r) { t[v] = vector{-x[l], 1}; } else { int m = (l + r) >> 1; int z = v + ((m - l + 1) << 1); build(v + 1, l, m); build(z, m + 1, r); t[v] = t[v + 1] * t[z]; } }; build(0, 0, n - 1); vector val(n); function)> eval = [&](int v, int l, int r, vector p) { p %= t[v]; if (p.size() < 150) { for (int i = l; i <= r; i++) { val[i] = evaluate(p, x[i]); } return; } if (l == r) { val[l] = p[0]; } else { int m = (l + r) >> 1; int z = v + ((m - l + 1) << 1); eval(v + 1, l, m, p); eval(z, m + 1, r, p); } }; vector d = derivative(t[0]); eval(0, 0, n - 1, d); for (int i = 0; i < n; i++) { val[i] = y[i] / val[i]; } function(int, int, int)> calc = [&](int v, int l, int r) { if (l == r) { return vector{val[l]}; } int m = (l + r) >> 1; int z = v + ((m - l + 1) << 1); return calc(v + 1, l, m) * t[z] + calc(z, m + 1, r) * t[v + 1]; }; return calc(0, 0, n - 1); } // f[i] = 1^n + 2^n + ... + up^n // O(n*log(n)) complexity template vector faulhaber(const T &up, int n) { vector ex(n + 1); T e = 1; for (int i = 0; i <= n; i++) { ex[i] = e; e /= i + 1; } vector den = ex; den.erase(den.begin()); for (auto &d : den) { d = -d; } vector num(n); T p = 1; for (int i = 0; i < n; i++) { p *= up + 1; num[i] = ex[i + 1] * (1 - p); } vector res = num * inverse(den); res.resize(n); T f = 1; for (int i = 0; i < n; i++) { res[i] *= f; f *= i + 1; } return res; } // (x + 1) * (x + 2) * ... * (x + n) // (can be optimized with precomputed inverses) template vector sequence(int n) { if (n == 0) { return {1}; } if (n % 2 == 1) { return sequence(n - 1) * vector{n, 1}; } vector c = sequence(n / 2); vector a = c; reverse(a.begin(), a.end()); T f = 1; for (int i = n / 2 - 1; i >= 0; i--) { f *= n / 2 - i; a[i] *= f; } vector b(n / 2 + 1); b[0] = 1; for (int i = 1; i <= n / 2; i++) { b[i] = b[i - 1] * (n / 2) / i; } vector h = a * b; h.resize(n / 2 + 1); reverse(h.begin(), h.end()); f = 1; for (int i = 1; i <= n / 2; i++) { f /= i; h[i] *= f; } vector res = c * h; return res; } template T factorial(long long n) { if (n == 0) return 1; int m = min((long long)(sqrt(n) * 2), n); vector a = sequence(m); vector x(n / m); for (size_t i = 0; i < x.size(); ++i) { x[i] = i; x[i] *= m; } vector b = evaluate(a, x); T res = 1; for (auto v : b) res *= v; for (long long i = n / m * m + 1; i <= n; ++i) { res *= i; } return res; } template T binomial(int n, int m) { return factorial(n) / factorial(m) / factorial(n - m); } template class OnlineProduct { public: const vector a; vector b; vector c; OnlineProduct(const vector &a_) : a(a_) {} T add(const T &val) { int i = (int)b.size(); b.push_back(val); if ((int)c.size() <= i) { c.resize(i + 1); } c[i] += a[0] * b[i]; int z = 1; while ((i & (z - 1)) == z - 1 && (int)a.size() > z) { vector a_mul(a.begin() + z, a.begin() + min(z << 1, (int)a.size())); vector b_mul(b.end() - z, b.end()); vector c_mul = a_mul * b_mul; if ((int)c.size() <= i + (int)c_mul.size()) { c.resize(i + c_mul.size() + 1); } for (int j = 0; j < (int)c_mul.size(); j++) { c[i + 1 + j] += c_mul[j]; } z <<= 1; } return c[i]; } };