Submit Info #62412

Problem Lang User Status Time Memory
Multipoint Evaluation cpp brunovsky AC 1758 ms 39.60 MiB

ケース詳細
Name Status Time Memory
example_00 AC 1 ms 0.45 MiB
example_01 AC 1 ms 0.45 MiB
max_random_00 AC 1755 ms 39.60 MiB
max_random_01 AC 1758 ms 39.60 MiB
random_00 AC 275 ms 15.11 MiB
random_01 AC 176 ms 9.83 MiB
random_02 AC 833 ms 27.21 MiB
zero_00 AC 1 ms 0.45 MiB

#pragma GCC optimize("O3,unroll-loops") #pragma GCC target("avx,avx2,sse,sse2,ssse3") #include <bits/stdc++.h> #ifdef LOCAL #include "code/formatting.hpp" #else #define debug(...) (void)0 #endif using namespace std; static_assert(sizeof(int) == 4 && sizeof(long) == 8); template <typename Fun> class y_combinator_result { Fun fun_; public: template <typename T> explicit y_combinator_result(T&& fun) : fun_(std::forward<T>(fun)) {} template <typename... Args> decltype(auto) operator()(Args&&... args) { return fun_(std::ref(*this), std::forward<Args>(args)...); } }; template <typename Fun> auto y_combinator(Fun&& fun) { return y_combinator_result<std::decay_t<Fun>>(std::forward<Fun>(fun)); } template <uint32_t mod> struct modnum { // change these if you need another size of integers static constexpr inline uint32_t MOD = mod; using u32 = uint32_t; using u64 = uint64_t; using i32 = int32_t; using i64 = int64_t; static_assert(mod > 0 && mod < UINT_MAX / 2); u32 n; constexpr modnum() : n(0) {} constexpr modnum(u64 v) : n(v >= mod ? v % mod : v) {} constexpr modnum(u32 v) : n(v >= mod ? v % mod : v) {} constexpr modnum(i64 v) : modnum(v >= 0 ? u64(v) : u64(mod + v % int(mod))) {} constexpr modnum(i32 v) : modnum(v >= 0 ? u32(v) : u32(mod + v % int(mod))) {} explicit constexpr operator int() const { return n; } explicit constexpr operator bool() const { return n != 0; } static constexpr u32 fit(u32 x) { return x >= mod ? x - mod : x; } static constexpr int modinv(u32 x) { int nx = 1, ny = 0; u32 y = mod; while (x) { auto k = y / x; y = y % x; ny = ny - k * nx; swap(x, y), swap(nx, ny); } return ny < 0 ? mod + ny : ny; } friend constexpr modnum modpow(modnum b, int64_t e) { modnum p = 1; while (e > 0) { if (e & 1) p = p * b; if (e >>= 1) b = b * b; } return p; } constexpr modnum inv() const { return modinv(n); } constexpr modnum operator-() const { return n == 0 ? n : mod - n; } constexpr modnum operator+() const { return *this; } constexpr modnum operator++(int) { return n = fit(n + 1), *this - 1; } constexpr modnum operator--(int) { return n = fit(mod + n - 1), *this + 1; } constexpr modnum& operator++() { return n = fit(n + 1), *this; } constexpr modnum& operator--() { return n = fit(mod + n - 1), *this; } constexpr modnum& operator+=(modnum v) { return n = fit(n + v.n), *this; } constexpr modnum& operator-=(modnum v) { return n = fit(mod + n - v.n), *this; } constexpr modnum& operator*=(modnum v) { return n = (u64(n) * v.n) % mod, *this; } constexpr modnum& operator/=(modnum v) { return *this *= v.inv(); } friend constexpr modnum operator+(modnum lhs, modnum rhs) { return lhs += rhs; } friend constexpr modnum operator-(modnum lhs, modnum rhs) { return lhs -= rhs; } friend constexpr modnum operator*(modnum lhs, modnum rhs) { return lhs *= rhs; } friend constexpr modnum operator/(modnum lhs, modnum rhs) { return lhs /= rhs; } friend string to_string(modnum v) { return to_string(v.n); } friend constexpr bool operator==(modnum lhs, modnum rhs) { return lhs.n == rhs.n; } friend constexpr bool operator!=(modnum lhs, modnum rhs) { return lhs.n != rhs.n; } friend ostream& operator<<(ostream& out, modnum v) { return out << v.n; } friend istream& operator>>(istream& in, modnum& v) { i64 n; return in >> n, v = modnum(n), in; } }; // Base include namespace fft { template <typename T> struct my_complex { using self = my_complex<T>; T x, y; constexpr my_complex(T x = T(0), T y = T(0)) : x(x), y(y) {} constexpr T& real() { return x; } constexpr T& imag() { return y; } constexpr const T& real() const { return x; } constexpr const T& imag() const { return y; } constexpr void real(T v) { x = v; } constexpr void imag(T v) { y = v; } constexpr friend auto real(self a) { return a.x; } constexpr friend auto imag(self a) { return a.y; } constexpr self rot_ccw(self a) { return self(-a.y, a.x); } constexpr self rot_cw(self a) { return self(a.y, -a.x); } constexpr friend auto abs(self a) { return sqrt(norm(a)); } constexpr friend auto arg(self a) { return atan2(a.y, a.x); } constexpr friend auto norm(self a) { return a.x * a.x + a.y * a.y; } constexpr friend auto conj(self a) { return self(a.x, -a.y); } constexpr friend auto inv(self a) { return self(a.x / norm(a), -a.y / norm(a)); } constexpr friend auto polar(T r, T theta = T()) { return self(r * cos(theta), r * sin(theta)); } constexpr T& operator[](int i) { assert(i == 0 || i == 1), *(&x + i); } constexpr const T& operator[](int i) const { assert(i == 0 || i == 1), *(&x + i); } constexpr self& operator+=(self b) { return *this = *this + b; } constexpr self& operator-=(self b) { return *this = *this - b; } constexpr self& operator*=(self b) { return *this = *this * b; } constexpr self& operator/=(self b) { return *this = *this / b; } constexpr friend self operator*(self a, T b) { return self(a.x * b, a.y * b); } constexpr friend self operator*(T a, self b) { return self(a * b.x, a * b.y); } constexpr friend self operator/(self a, T b) { return self(a.x / b, a.y / b); } constexpr friend self operator/(T a, self b) { return a * inv(b); } constexpr friend self operator+(self a) { return self(a.x, a.y); } constexpr friend self operator-(self a) { return self(-a.x, -a.y); } constexpr friend self operator+(self a, self b) { return self(a.x + b.x, a.y + b.y); } constexpr friend self operator-(self a, self b) { return self(a.x - b.x, a.y - b.y); } constexpr friend self operator/(self a, self b) { return a * inv(b); } constexpr friend self operator*(self a, self b) { return self(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); } }; using default_complex = my_complex<double>; constexpr double TAU = 6.283185307179586476925286766559; int next_two(int32_t N) { return N > 1 ? 8 * sizeof(N) - __builtin_clz(N - 1) : 0; } template <typename T, typename D> constexpr T fft_round(D coefficient) { return is_integral<T>::value ? llround(coefficient) : coefficient; } template <typename T> void trim_vector(vector<T>& v) { if constexpr (is_floating_point<T>::value) while (!v.empty() && abs(v.back()) <= 30 * numeric_limits<T>::epsilon()) v.pop_back(); else while (!v.empty() && v.back() == T(0)) v.pop_back(); } template <typename T> auto naive_multiply(const vector<T>& a, const vector<T>& b) { int A = a.size(), B = b.size(), C = A && B ? A + B - 1 : 0; vector<T> c(C); for (int i = 0; i < A && B; i++) for (int j = 0; j < B; j++) c[i + j] += a[i] * b[j]; trim_vector(c); return c; } struct fft_reverse_cache { static inline vector<vector<int>> rev; static const int* get(int N) { int n = next_two(N), r = rev.size(); rev.resize(max(r, n + 1)); if (rev[n].empty()) { int R = 1 << n; rev[n].assign(R, 0); for (int i = 0; i < N; i++) { rev[n][i] = (rev[n][i >> 1] | ((i & 1) << n)) >> 1; } } return rev[n].data(); } }; template <typename T> struct root_of_unity {}; template <typename D> struct root_of_unity<my_complex<D>> { static auto get(int n) { assert(n > 0); return my_complex<D>(cos(TAU / n), sin(TAU / n)); } }; template <typename C> struct fft_roots_cache { static inline vector<C> root = vector<C>(2, C(1)); static inline vector<C> invroot = vector<C>(2, C(1)); static inline vector<C> scratch_a, scratch_b; static auto get(int N) { for (int k = root.size(); k < N; k *= 2) { root.resize(2 * k); invroot.resize(2 * k); auto z = root_of_unity<C>::get(2 * k); auto iz = C(1) / z; for (int i = k / 2; i < k; i++) { root[2 * i] = root[i]; root[2 * i + 1] = root[i] * z; invroot[2 * i] = invroot[i]; invroot[2 * i + 1] = invroot[i] * iz; } } return make_pair(cref(root), cref(invroot)); } static auto get_scratch(int N) { if (int(scratch_a.size()) < N) { scratch_a.resize(N); scratch_b.resize(N); } return make_pair(ref(scratch_a), ref(scratch_b)); } }; template <typename T> void fft_bit_reverse(vector<T>& a, int N) { auto rev = fft_reverse_cache::get(N); for (int i = 0; i < N; i++) { if (i < rev[i]) { swap(a[i], a[rev[i]]); } } } template <bool inverse, typename T> void fft_transform(vector<T>& a, int N, bool reverse = true) { if (reverse) { fft_bit_reverse(a, N); } auto [root, invroot] = fft_roots_cache<T>::get(N); for (int k = 1; k < N; k *= 2) { for (int i = 0; i < N; i += 2 * k) { for (int l = i, r = i + k, j = 0; j < k; j++, l++, r++) { auto z = inverse ? invroot[j + k] : root[j + k]; auto t = a[r] * z; a[r] = a[l] - t; a[l] = a[l] + t; } } } if constexpr (inverse) { auto inv = T(1) / T(N); for (int i = 0; i < N; i++) { a[i] *= inv; } } } } // namespace fft // NTT with modnums namespace fft { template <uint32_t MOD> struct root_of_unity<modnum<MOD>> { using type = modnum<MOD>; static constexpr int ntt_primitive_root(int p) { if (p == 2 || p == 3) return p - 1; if (p == 167772161 || p == 469762049 || p == 998244353) return 3; if (p == 1000000007) return 5; assert(p % 2 == 1); auto modpow = [](int base, int64_t e, int mod) { int x = 1; base = base % mod; while (e > 0) { if (e & 1) x = int64_t(x) * int64_t(base) % mod; if (e >>= 1) base = int64_t(base) * int64_t(base) % mod; } return x; }; int divs[20] = {2}; int cnt = 1; int x = (p - 1) / 2; // phi(p) while (x % 2 == 0) { x /= 2; } for (int i = 3; i * i <= x; i += 2) { if (x % i == 0) { divs[cnt++] = i; while (x % i == 0) { x /= i; } } } if (x > 1) { divs[cnt++] = x; } for (int g = 2;; g++) { bool ok = true; for (int i = 0; i < cnt; i++) { if (modpow(g, (p - 1) / divs[i], p) == 1) { ok = false; break; } } if (ok) { return g; } } __builtin_unreachable(); } static type get(int n) { modnum<MOD> g = ntt_primitive_root(MOD); assert(n > 0 && (MOD - 1) % n == 0 && "Modulus cannot handle NTT this large"); return modpow(g, int64_t(MOD - 1) / n); } }; template <uint32_t MOD> auto ntt_multiply(const vector<modnum<MOD>>& a, const vector<modnum<MOD>>& b) { using T = modnum<MOD>; if (a.empty() || b.empty()) { return vector<T>(); } int A = a.size(), B = b.size(); int S = A + B - 1, s = next_two(S), N = 1 << s; if (1.0 * A * B <= 2.5 * N * s) { return naive_multiply(a, b); } auto fb = fft_roots_cache<T>::get_scratch(N).second; vector<T> c(N); copy_n(begin(a), A, begin(c)); copy_n(begin(b), B, begin(fb)); fill_n(begin(fb) + B, N - B, T(0)); fft_transform<0>(c, N); fft_transform<0>(fb, N); for (int i = 0; i < N; i++) { c[i] = c[i] * fb[i]; } fft_transform<1>(c, N); trim_vector(c); return c; } } // namespace fft namespace polymath { static mt19937 rng(random_device{}()); using T = modnum<998244353>; using poly = vector<T>; // TODO: Fill in core FFT functions void bitr(poly& u) { fft::fft_bit_reverse(u, size(u)); } void fft(poly& u) { fft::fft_transform<0>(u, size(u)); } void ifft(poly& u) { fft::fft_transform<1>(u, size(u)); } void sfft(poly& u) { fft::fft_transform<0>(u, size(u)), bitr(u); } void sifft(poly& u) { fft::fft_transform<1>(u, size(u), false); } poly convolve(const poly& a, const poly& b) { return fft::ntt_multiply(a, b); } poly rand_poly(int n) { static uniform_int_distribution<int> dist(1, T::MOD - 1); poly v(n); for (int i = 0; i < n; i++) { v[i] = dist(rng); } return v; } // Utility stuff, then elementary operations, then FFT and power series operations int size(const poly& u) { return u.size(); } bool empty(const poly& u) { return u.empty(); } auto& assign(poly& u, int n, T val) { return u.assign(n, val), u; } auto& clear(poly& u) { return u.clear(), u; } auto& pop(poly& u) { return empty(u) ? u : (u.pop_back(), u); } auto& push(poly& u, T add) { return u.push_back(add), u; } auto& reverse(poly& u) { return reverse(begin(u), end(u)), u; } auto& resize(poly& u, int n) { return u.resize(n), u; } auto& shrink(poly& u, int n) { return u.resize(min(n, size(u))), u; } auto& grow(poly& u, int n) { return u.resize(max(n, size(u))), u; } auto& reserve(poly& u, int n) { return u.reserve(n), u; } auto& trim(poly& u, int upto = 0) { while (size(u) > upto && u.back() == T(0)) u.pop_back(); return u; } auto shrunk(poly u, int n) { return u.resize(min(n, size(u))), u; } auto resized(poly u, int n) { return u.resize(n), u; } auto trimmed(poly u) { return trim(u), u; } auto& operator+=(poly& u, T value) { return grow(u, 1), u[0] += value, u; } auto& operator-=(poly& u, T value) { return grow(u, 1), u[0] -= value, u; } auto& operator*=(poly& u, T value) { for (int i = 0; i < size(u); i++) u[i] *= value; return u; } auto& operator/=(poly& u, T value) { return u *= T(1) / value; } auto& operator+=(poly& u, const poly& v) { grow(u, size(v)); for (int i = 0; i < size(v); i++) u[i] += v[i]; return u; } auto& operator-=(poly& u, const poly& v) { grow(u, size(v)); for (int i = 0; i < size(v); i++) u[i] -= v[i]; return u; } auto& operator*=(poly& u, const poly& v) { return u = convolve(u, v), u; } auto operator-(poly u) { for (int i = 0; i < size(u); i++) u[i] = -u[i]; return u; } auto operator+(poly u, T value) { return u += value, u; } auto operator-(poly u, T value) { return u -= value, u; } auto operator*(poly u, T value) { return u *= value, u; } auto operator+(T value, poly u) { return u += value, u; } auto operator-(T value, poly u) { return u -= value, u; } auto operator*(T value, poly u) { return u *= value, u; } auto operator/(poly u, T value) { return u /= value, u; } auto operator+(poly u, const poly& v) { return u += v, u; } auto operator-(poly u, const poly& v) { return u -= v, u; } auto operator*(const poly& u, const poly& v) { return convolve(u, v); } auto& pointwisein(poly& u, const poly& o) { grow(u, size(o)); for (int i = 0, n = o.size(); i < n; i++) u[i] *= o[i]; return u; } auto pointwise(poly u, const poly& o) { return pointwisein(u, o), u; } auto& rem_shift(poly& u, int s) { return u.erase(begin(u), begin(u) + min(size(u), s)), u; } auto& add_shift(poly& u, int s, T add = T(0)) { return u.insert(begin(u), s, add), u; } auto& operator>>=(poly& u, int s) { return rem_shift(u, s); } auto& operator<<=(poly& u, int s) { return add_shift(u, s); } auto operator>>(poly u, int s) { return u >>= s, u; } auto operator<<(poly u, int s) { return u <<= s, u; } auto& derivativein(poly& u) { for (int i = 1; i < size(u); i++) u[i - 1] = T(i) * u[i]; return pop(u); } auto derivative(poly u) { return derivativein(u), u; } auto& integralin(poly& u, T zero = 0) { push(u, T(0)); for (int i = size(u) - 1; i > 0; i--) u[i] = u[i - 1] / T(i); return u[0] = zero, u; } auto integral(poly u) { return integralin(u), u; } auto eval(const poly& u, T x) { T sum = 0; for (int i = size(u) - 1; i >= 0; i--) sum = u[i] + sum * x; return sum; } } // namespace polymath // Series: inverse, exp, log, pow, division, composition namespace polymath { auto inverse(const poly& u, int n = -1) { assert(!empty(u) && u[0] != T(0)); n = n < 0 ? size(u) : n; poly ans(n, T(1) / u[0]); for (int d = 1; d < n; d *= 2) { poly f(2 * d), g(2 * d); for (int j = 0; j < min(2 * d, size(u)); j++) f[j] = u[j]; for (int j = 0; j < d; j++) g[j] = ans[j]; fft(f), fft(g); for (int j = 0; j < 2 * d; j++) f[j] *= g[j]; ifft(f); for (int j = 0; j < d; j++) f[j] = 0; fft(f); for (int j = 0; j < 2 * d; j++) f[j] *= g[j]; ifft(f); for (int j = d; j < min(2 * d, n); j++) ans[j] = -f[j]; } return ans; } auto exp(const poly& u, int n = -1) { assert(empty(u) || u[0] == T(0)); n = n < 0 ? size(u) : n; poly b{1, size(u) > 1 ? u[1] : 0}, c{1}, z1, z2{1, 1}; b.reserve(n); for (int m = 2; m < n; m *= 2) { auto y = b; resize(y, 2 * m); sfft(y); z1 = move(z2); poly z(m); for (int i = 0; i < m; ++i) z[i] = y[i] * z1[i]; sifft(z); fill(begin(z), begin(z) + m / 2, T(0)); sfft(z); for (int i = 0; i < m; ++i) z[i] *= -z1[i]; sifft(z); c.insert(end(c), begin(z) + m / 2, end(z)); z2 = c; resize(z2, 2 * m); sfft(z2); poly x(begin(u), begin(u) + min(size(u), m)); resize(x, m), derivativein(x), push(x, 0); sfft(x); for (int i = 0; i < m; ++i) x[i] *= y[i]; sifft(x); x -= derivative(b); resize(x, 2 * m); for (int i = 0; i < m - 1; ++i) x[m + i] = x[i], x[i] = T(0); sfft(x); for (int i = 0; i < 2 * m; ++i) x[i] *= z2[i]; sifft(x); integralin(pop(x)); for (int i = m; i < min(size(u), 2 * m); ++i) x[i] += u[i]; fill(begin(x), begin(x) + m, T(0)); sfft(x); for (int i = 0; i < 2 * m; ++i) x[i] *= y[i]; sifft(x); b.insert(end(b), begin(x) + m, end(x)); } resize(b, n); return b; } auto log(const poly& u, int n = -1) { assert(u[0] == 1); n = n < 0 ? size(u) : n; auto v = derivative(u) * inverse(u, n); integralin(resize(v, n - 1)); return v; } auto pow(const poly& u, int64_t k, int n = -1) { n = n < 0 ? size(u) : n; for (int i = 0; i < size(u) && i * k < n; i++) { if (u[i] != 0) { poly ans = exp(k * log((u >> i) / u[i])) * modpow(u[i], k); resize(ans <<= i * k, n); return ans; } } return poly(n, 0); } auto quotient(poly u, poly v, int n = -1) { n = n < 0 ? min(size(u), size(v)) : n; trim(u), trim(v); assert(!empty(v) && v[0] != T(0)); if (size(u) < size(v) || size(u) == 1) { u /= v[0]; } else { shrink(u, n); u *= inverse(v, n); } resize(u, n); return u; } auto operator/(poly u, poly v) { trim(u), trim(v); assert(!empty(v)); if (size(u) < size(v)) { return poly(); } else { int n = size(u) - size(v) + 1; reverse(u), reverse(v); auto d = u * inverse(v, n); shrink(d, n), reverse(d), trim(d); return d; } } auto operator%(poly u, poly v) { return u -= v * (u / v), trim(u), u; } auto& operator/=(poly& u, poly v) { return u = u / move(v), u; } auto& operator%=(poly& u, poly v) { return u = u % move(v), u; } auto divrem(poly u, poly v) { auto d = u / v; u -= d * v, trim(u); return make_pair(move(d), move(u)); } // Composition p(q(x)), naively quadratic, complexity O(n² log n) auto naive_composition(const poly& p, poly q, int n = -1) { n = n < 0 ? (size(p) - 1) * (size(q) - 1) + 1 : n; if (n == 0) { return poly(); } trim(q); poly ans(n), u{1}; u.reserve(n); for (int i = 0; i < min(n, size(p)); i++) { ans += p[i] * u, u = shrunk(u * q, n); } return ans; } // Composition p(q(x)), square-root decomposition, complexity O((n log n)^3/2) auto sqrt_composition(const poly& p, poly q, int n = -1) { // Source: https://judge.yosupo.jp/submission/43112 n = n < 0 ? (size(p) - 1) * (size(q) - 1) + 1 : n; if (n == 0) { return poly(); } trim(q); int k = sqrt(n); int d = n / k + (k * k < n); vector<poly> small(d + 1); small[0] = {1}, small[1] = shrunk(q, n); for (int i = 2; i <= d; i++) { small[i] = small[i / 2] * small[(i + 1) / 2]; shrink(small[i], n); } vector<poly> fi(k); for (int i = 0; i < k; i++) { for (int j = 0; j < d && i * d + j < min(n, size(p)); j++) { int x = i * d + j; fi[i] += p[x] * small[j]; } } poly ans(n), big = {1}; for (int i = 0; i < k; i++) { fi[i] = shrunk(fi[i] * big, n); ans += fi[i]; big = shrunk(big * small[d], n); } return ans; } auto composition(const poly& u, const poly& v, int n = -1) { return sqrt_composition(u, v, n); } } // namespace polymath // Series: gcd, resultant namespace polymath { struct polymat { bool identity = false; poly a0, a1, b0, b1; polymat() : identity(true) {} explicit polymat(poly v) : a1({1}), b0({1}), b1(move(v)) {} polymat(poly a0, poly a1, poly b0, poly b1) : identity(false), a0(move(a0)), a1(move(a1)), b0(move(b0)), b1(move(b1)) {} friend polymat operator*(const polymat& lhs, const polymat& rhs) { if (lhs.identity || rhs.identity) return lhs.identity ? rhs : lhs; return polymat{ trimmed(lhs.a0 * rhs.a0 + lhs.a1 * rhs.b0), // trimmed(lhs.a0 * rhs.a1 + lhs.a1 * rhs.b1), // trimmed(lhs.b0 * rhs.a0 + lhs.b1 * rhs.b0), // trimmed(lhs.b0 * rhs.a1 + lhs.b1 * rhs.b1), // }; } auto multiply(const poly& a, const poly& b) { if (identity) return make_pair(a, b); return make_pair(trimmed(a0 * a + a1 * b), trimmed(b0 * a + b1 * b)); } }; auto halfgcd(const poly& p0, const poly& p1) { int m = size(p0) >> 1; if (size(p1) <= m) { return polymat(); } polymat R = halfgcd(p0 >> m, p1 >> m); auto [a, b] = R.multiply(p0, p1); if (size(b) <= m) { return R; } auto [q, r] = divrem(a, b); polymat Q(-q); int k = (m << 1) - size(b) + 1; if (size(r) <= k) { return Q * R; } return halfgcd(b >> k, r >> k) * Q * R; } auto cogcd(const poly& p0, const poly& p1) { polymat M = halfgcd(p0, p1); auto [a, b] = M.multiply(p0, p1); if (empty(b)) { return M; } auto [q, r] = divrem(a, b); polymat Q(-q); if (empty(r)) { return Q * M; } return cogcd(b, r) * Q * M; } auto exgcd(poly a, poly b, poly& x, poly& y) { // ax + by = monic gcd(a,b) trim(a), trim(b); if (empty(a) || empty(b)) { x = y = {1}; return empty(b) ? a : b; } if (size(a) > size(b)) { auto c = cogcd(a, b); x = move(c.a0), y = move(c.a1); } else if (size(a) < size(b)) { auto c = cogcd(b, a); x = move(c.a1), y = move(c.a0); } else { auto [q, r] = divrem(a, b); auto c = cogcd(b, r) * polymat(-q); x = move(c.a0), y = move(c.a1); } auto g = trimmed(a * x + b * y); if (auto v = g.back(); v != T(1)) { x /= v, y /= v, g /= v; } return g; } auto gcd(poly a, poly b) { while (size(a) > 0) { b %= a; swap(a, b); } return empty(b) ? b : b / b.back(); } auto naive_gcd(poly a, poly b, poly& x, poly& y) { if (empty(a) && empty(b)) { x = y = {}; return poly(); } poly xn = {1}, yn = {}; x = {}, y = {1}; while (size(a) > 0) { auto [q, r] = divrem(b, a); b = move(r); x = x - q * xn, trim(x); y = y - q * yn, trim(y); swap(a, b); swap(x, xn); swap(y, yn); } if (auto v = b.back(); v != T(1)) { x /= v, y /= v, b /= v; } return b; } auto resultant(const poly& a, const poly& b) { int A = size(a), B = size(b); if (B == 0) { return T(); } else if (B == 1) { return modpow(b[0], A - 1); } else { auto c = a % b; A -= size(c); int sign = ((A - 1) & (B - 1) & 1) ? -1 : 1; auto mul = modpow(b[0], A - 1) * sign; return mul * resultant(b, c); } } } // namespace polymath // Series: multipoint evaluation, interpolation namespace polymath { constexpr int MULTIEVAL_THRESHOLD = 50; auto multieval(const poly& p, const poly& x) { // Source: https://github.com/Aeren1564/Algorithms int S = size(x); if (S == 0) { return poly(); } if (size(p) == 0) { return poly(S); } vector<poly> st(2 * S); y_combinator([&](auto self, int u, int l, int r) -> void { if (r - l == 1) { st[u] = {-x[l], 1}; } else { int m = l + ((r - l) >> 1), v = u + ((m - l) << 1); self(u + 1, l, m), self(v, m, r); st[u] = st[u + 1] * st[v]; } })(0, 0, S); poly ans(S); y_combinator([&](auto self, int u, int l, int r, poly f) -> void { f %= st[u]; if (size(f) <= MULTIEVAL_THRESHOLD) { for (int i = l; i < r; i++) { ans[i] = eval(f, x[i]); } } else if (r - l == 1) { ans[l] = f[0]; } else { int m = l + ((r - l) >> 1), v = u + ((m - l) << 1); self(u + 1, l, m, f), self(v, m, r, move(f)); } })(0, 0, S, p); return ans; } auto interpolate(const poly& x, const poly& y) { assert(size(x) == size(y)); int S = size(x); if (S == 0) { return poly(); } vector<poly> st(2 * S); y_combinator([&](auto self, int u, int l, int r) -> void { if (r - l == 1) { st[u] = {-x[l], 1}; } else { int m = l + ((r - l) >> 1), v = u + ((m - l) << 1); self(u + 1, l, m), self(v, m, r); st[u] = st[u + 1] * st[v]; } })(0, 0, S); poly val(S); y_combinator([&](auto self, int u, int l, int r, poly f) -> void { f %= st[u]; if (size(f) <= MULTIEVAL_THRESHOLD) { for (int i = l; i < r; i++) { val[i] = eval(f, x[i]); } } else if (r - l == 1) { val[l] = f[0]; } else { int m = l + ((r - l) >> 1), v = u + ((m - l) << 1); self(u + 1, l, m, f), self(v, m, r, move(f)); } })(0, 0, S, derivative(st[0])); for (int i = 0; i < S; i++) { val[i] = y[i] / val[i]; } return y_combinator([&](auto self, int u, int l, int r) -> poly { if (r - l == 1) { return poly{val[l]}; } else { int m = l + ((r - l) >> 1), v = u + ((m - l) << 1); return self(u + 1, l, m) * st[v] + self(v, m, r) * st[u + 1]; } })(0, 0, S); } } // namespace polymath int main() { ios::sync_with_stdio(false), cin.tie(nullptr); using namespace polymath; int N, M; cin >> N >> M; vector<T> p(N), x(M); for (int i = 0; i < N; i++) { cin >> p[i]; } for (int i = 0; i < M; i++) { cin >> x[i]; } auto ans = multieval(p, x); for (int i = 0; i < M; i++) { cout << ans[i] << " \n"[i + 1 == M]; } return 0; }