Submit Info #67552

Problem Lang User Status Time Memory
Polynomial Interpolation cpp pikelrik AC 1914 ms 36.53 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 1914 ms 36.51 MiB
max_random_01 AC 1913 ms 36.53 MiB
random_00 AC 1744 ms 35.02 MiB
random_01 AC 913 ms 19.12 MiB
random_02 AC 658 ms 15.01 MiB

#pragma GCC optimize("Ofast,inline") #pragma GCC target("avx,avx2,f16c,fma,sse3,ssse3,sse4.1,sse4.2") #include <bits/stdc++.h> using namespace std; template<uint32_t M> struct montgomery_mint { static_assert(0 < M, "Modulus must be non-zero"); static_assert(M % 2 != 0, "Modulus must be odd"); static_assert(M < (1u << 31), "Modulus must fit in a 32-bit signed integer"); constexpr static auto get_m() { return M; } constexpr static auto get_r() { return -M % M; } constexpr static auto get_r_inverse() { constexpr auto r = get_r(); uint32_t r_inverse = 1, r_pow = r; for (uint32_t n = M - 2; n != 0; n /= 2) { if (n % 2 != 0) { r_inverse = (uint64_t) r_inverse * r_pow % M; } r_pow = (uint64_t) r_pow * r_pow % M; } return r_inverse; } constexpr static auto get_r_squared() { return (uint64_t) get_r() * get_r() % M; } constexpr static auto get_m_dash() { constexpr auto r_inverse = get_r_inverse(); return (uint32_t) ((((uint64_t) r_inverse << 32) - 1) / M); } constexpr static uint32_t reduce(uint64_t x) { constexpr auto m_dash = get_m_dash(); uint32_t y = (x + (uint64_t) ((uint32_t) x * m_dash) * M) >> 32; return y >= M ? y - M : y; } uint32_t val; constexpr montgomery_mint() : val() {} constexpr montgomery_mint(long long x) : val() { constexpr auto r_squared = get_r_squared(); x %= (int) M; val = reduce((x < 0 ? x + M : x) * r_squared); } constexpr uint32_t get_val() const { return reduce(val); } constexpr montgomery_mint pow(uint64_t n) const { montgomery_mint ans = 1, x(*this); for (; n != 0; n /= 2) { if (n & 1) ans *= x; x *= x; } return ans; } constexpr montgomery_mint inv() const { return pow(M - 2); } constexpr montgomery_mint operator+() const { montgomery_mint m; m.val = val; return m; } constexpr montgomery_mint operator-() const { montgomery_mint m; m.val = (val == 0 ? 0 : M - val); return m; } constexpr montgomery_mint &operator+=(const montgomery_mint &m) { if ((val += m.val) >= M) val -= M; return *this; } constexpr montgomery_mint &operator-=(const montgomery_mint &m) { if ((val -= m.val) >= M) val += M; return *this; } constexpr montgomery_mint &operator*=(const montgomery_mint &m) { val = reduce((uint64_t) val * m.val); return *this; } constexpr montgomery_mint &operator/=(const montgomery_mint &m) { val = reduce((uint64_t) val * m.inv().val); return *this; } constexpr friend montgomery_mint operator+(const montgomery_mint &lhs, const montgomery_mint &rhs) { return montgomery_mint(lhs) += rhs; } constexpr friend montgomery_mint operator-(const montgomery_mint &lhs, const montgomery_mint &rhs) { return montgomery_mint(lhs) -= rhs; } constexpr friend montgomery_mint operator*(const montgomery_mint &lhs, const montgomery_mint &rhs) { return montgomery_mint(lhs) *= rhs; } constexpr friend montgomery_mint operator/(const montgomery_mint &lhs, const montgomery_mint &rhs) { return montgomery_mint(lhs) /= rhs; } constexpr friend bool operator==(const montgomery_mint &lhs, const montgomery_mint &rhs) { return lhs.val == rhs.val; } constexpr friend bool operator!=(const montgomery_mint &lhs, const montgomery_mint &rhs) { return lhs.val != rhs.val; } constexpr montgomery_mint &operator++() { return *this += 1; } constexpr montgomery_mint &operator--() { return *this -= 1; } constexpr montgomery_mint operator++(int) { montgomery_mint result(*this); *this += 1; return result; } constexpr montgomery_mint operator--(int) { montgomery_mint result(*this); *this -= 1; return result; } template<typename T> constexpr explicit operator T() const { return T(val); } template <uint32_t M1> constexpr explicit operator montgomery_mint<M1>() const { return montgomery_mint<M1>(get_val()); } friend std::ostream &operator<<(std::ostream &os, const montgomery_mint &m) { return os << m.get_val(); } friend std::istream &operator>>(std::istream &is, montgomery_mint &m) { long long x; return is >> x, m = x, is; } }; template<typename> struct is_mint_helper : std::false_type { }; template<uint32_t M> struct is_mint_helper<montgomery_mint<M>> : std::true_type { }; template<typename T> struct is_mint : is_mint_helper<typename std::decay<T>::type> { }; namespace ntt { template<uint32_t P> struct prime_info { constexpr static uint32_t root = 0, root_pw = 0; }; template<> struct prime_info<7340033> { constexpr static uint32_t root = 5, root_pw = 1 << 20; }; template<> struct prime_info<998244353> { constexpr static uint32_t root = 15311432, root_pw = 1 << 23; }; template<> struct prime_info<754974721> { constexpr static uint32_t root = 739831874, root_pw = 1 << 24; }; template<> struct prime_info<167772161> { constexpr static uint32_t root = 243, root_pw = 1 << 25; }; template<> struct prime_info<469762049> { constexpr static uint32_t root = 2187, root_pw = 1 << 26; }; std::vector<int> rev = {0}; void compute_bit_reverse(int lg) { static int computed = 0; if (lg <= computed) return; rev.resize(1 << lg); for (int i = (1 << computed); i < (1 << lg); i++) { rev[i] = (rev[i >> 1] >> 1) ^ ((i & 1) << 30); } computed = lg; } template<uint32_t M> std::vector<montgomery_mint<M>> root = {0, 1}; template<uint32_t M> void compute_roots(int lg) { static int computed = 1; if (lg <= computed) return; root<M>.resize(1 << lg); for (int k = computed; k < lg; k++) { montgomery_mint<M> z(prime_info<M>::root); for (int i = (1 << (k + 1)); i < prime_info<M>::root_pw; i <<= 1) { z *= z; } for (int i = (1 << (k - 1)); i < (1 << k); i++) { root<M>[i << 1] = root<M>[i]; root<M>[i << 1 | 1] = root<M>[i] * z; } } computed = lg; } template<int M> void ntt(std::vector<montgomery_mint<M>> &a) { int n = int(a.size()), lg = 32 - __builtin_clz(n) - 1; compute_bit_reverse(lg); compute_roots<M>(lg); int shift = 31 - lg; for (int i = 0; i < n; i++) { if (i < (rev[i] >> shift)) { std::swap(a[i], a[rev[i] >> shift]); } } for (int k = 1; k < n; k <<= 1) { for (int i = 0; i < n; i += 2 * k) { for (int j = 0; j < k; j++) { auto z = root<M>[j + k] * a[i + j + k]; a[i + j + k] = a[i + j] - z; a[i + j] = a[i + j] + z; } } } } template<typename mint> std::enable_if_t<is_mint<mint>::value, std::vector<mint>> naive_convolution(const std::vector<mint> &a, const std::vector<mint> &b) { std::vector<mint> c(a.size() + b.size() - 1); for (int i = 0; i < (int) a.size(); i++) { for (int j = 0; j < (int) b.size(); j++) { c[i + j] += a[i] * b[j]; } } return c; } template<uint32_t M> std::enable_if_t<prime_info<M>::root != 0, std::vector<montgomery_mint<M>>> convolution(std::vector<montgomery_mint<M>> a, std::vector<montgomery_mint<M>> b) { if (a.empty() || b.empty()) { return {}; } if (a.size() < 8 || b.size() < 8) { return naive_convolution(a, b); } int n = 1; while (n < a.size() + b.size()) { n <<= 1; } a.resize(n), b.resize(n); ntt<M>(a); ntt<M>(b); montgomery_mint<M> n_inv = montgomery_mint<M>(n).inv(); for (int i = 0; i < n; i++) { a[i] *= b[i] * n_inv; } std::reverse(a.begin() + 1, a.end()); ntt<M>(a); return a; } template<uint32_t M> std::enable_if_t<prime_info<M>::root != 0, void> inplace_convolution(std::vector<montgomery_mint<M>> &a, std::vector<montgomery_mint<M>> b) { if (a.empty() || b.empty()) { a.clear(); return; } if (a.size() < 8 || b.size() < 8) { a = naive_convolution(a, b); return; } int n = 1; while (n < a.size() + b.size()) { n <<= 1; } a.resize(n), b.resize(n); ntt<M>(a); ntt<M>(b); montgomery_mint<M> n_inv = montgomery_mint<M>(n).inv(); for (int i = 0; i < n; i++) { a[i] *= b[i] * n_inv; } std::reverse(a.begin() + 1, a.end()); ntt<M>(a); } template<uint32_t M> montgomery_mint<M> garner(int a1, int a2, int a3) { constexpr auto M1 = 754974721, M2 = 167772161, M3 = 469762049; constexpr int R12 = montgomery_mint<M2>(M1).inv().get_val(); constexpr int R13 = montgomery_mint<M3>(M1).inv().get_val(); constexpr int R23 = montgomery_mint<M3>(M2).inv().get_val(); int x1 = a1; int x2 = (long long) (a2 - x1) * R12 % M2; if (x2 < 0) x2 += M2; int x3 = ((long long) (a3 - x1) * R13 % M3 - x2) * R23 % M3; if (x3 < 0) x3 += M3; return montgomery_mint<M>(x1) + montgomery_mint<M>(x2) * M1 + montgomery_mint<M>(x3) * M1 * M2; } template<uint32_t M> std::enable_if_t<prime_info<M>::root == 0, std::vector<montgomery_mint<M>>> convolution(std::vector<montgomery_mint<M>> a, const std::vector<montgomery_mint<M>> &b) { constexpr auto M1 = 754974721u, M2 = 167772161u, M3 = 469762049u; auto c1 = convolution(std::vector<montgomery_mint<M1>>(a.begin(), a.end()), std::vector<montgomery_mint<M1>>(b.begin(), b.end())); auto c2 = convolution(std::vector<montgomery_mint<M2>>(a.begin(), a.end()), std::vector<montgomery_mint<M2>>(b.begin(), b.end())); auto c3 = convolution(std::vector<montgomery_mint<M3>>(a.begin(), a.end()), std::vector<montgomery_mint<M3>>(b.begin(), b.end())); int n = (int) c1.size(); a.resize(n); for (int i = 0; i < n; i++) { a[i] = garner<M>(c1[i].get_val(), c2[i].get_val(), c3[i].get_val()); } return a; } template<uint32_t M = 998244353, typename T> std::enable_if_t<!is_mint<T>::value, std::vector<montgomery_mint<M>>> convolution(const std::vector<T> &a, const std::vector<T> &b) { return convolution(std::vector<montgomery_mint<M>>(a.begin(), a.end()), std::vector<montgomery_mint<M>>(b.begin(), b.end())); } int garner(int a1, int a2, int a3, uint32_t M) { constexpr auto M1 = 754974721, M2 = 167772161, M3 = 469762049; constexpr int R12 = montgomery_mint<M2>(M1).inv().get_val(); constexpr int R13 = montgomery_mint<M3>(M1).inv().get_val(); constexpr int R23 = montgomery_mint<M3>(M2).inv().get_val(); int x1 = a1; int x2 = (long long) (a2 - x1) * R12 % M2; if (x2 < 0) x2 += M2; int x3 = ((long long) (a3 - x1) * R13 % M3 - x2) * R23 % M3; if (x3 < 0) x3 += M3; return (x1 + (long long) x2 * M1 + (long long) x3 * M1 % M * M2) % M; } template<typename T> std::enable_if_t<!is_mint<T>::value, std::vector<T>> convolution(std::vector<T> a, const std::vector<T> &b, uint32_t M) { constexpr auto M1 = 754974721u, M2 = 167772161u, M3 = 469762049u; auto c1 = convolution(std::vector<montgomery_mint<M1>>(a.begin(), a.end()), std::vector<montgomery_mint<M1>>(b.begin(), b.end())); auto c2 = convolution(std::vector<montgomery_mint<M2>>(a.begin(), a.end()), std::vector<montgomery_mint<M2>>(b.begin(), b.end())); auto c3 = convolution(std::vector<montgomery_mint<M3>>(a.begin(), a.end()), std::vector<montgomery_mint<M3>>(b.begin(), b.end())); int n = (int) c1.size(); a.resize(n); for (int i = 0; i < n; i++) { a[i] = garner(c1[i].get_val(), c2[i].get_val(), c3[i].get_val(), M); } return a; } template<typename T> void normalize(const std::vector<T> &a) { for (int i = int(a.size()) - 1; i >= 0; i--) { if (a[i]) { a.resize(i + 1); return; } } a.clear(); } } template <typename T> struct polynomial : public std::vector<T> { template <typename...args> polynomial(args...A) : std::vector<T>(A...) { } polynomial(const std::initializer_list<T> &l) : std::vector<T>(l) { } int deg() const noexcept { return (int) this->size() - 1; } void normalize() { for (int i = deg(); i >= 0; i--) { if ((*this)[i]) { this->resize(i + 1); return; } } this->clear(); } polynomial &operator+=(const polynomial &q) { if (q.size() > this->size()) { this->resize(q.size()); } for (int i = 0; i < q.size(); i++) { (*this)[i] += q[i]; } normalize(); return *this; } polynomial &operator-=(const polynomial &q) { if (q.size() > this->size()) { this->resize(q.size()); } for (int i = 0; i < q.size(); i++) { (*this)[i] -= q[i]; } normalize(); return *this; } void naive_mul(polynomial &a, const polynomial &b) const { polynomial result(a.deg() + b.deg() + 1); for (int i = 0; i <= a.deg(); i++) { for (int j = 0; j <= b.deg(); j++) { result[i + j] += a[i] * b[j]; } } a.swap(result); } polynomial &operator*=(const polynomial &q) { if (this->empty() || q.empty()) { this->clear(); } else if (this->size() <= 60) { naive_mul(*this, q); } else { ntt::inplace_convolution(*this, q); } normalize(); return *this; } polynomial &operator*=(const T &val) { for (auto &x : *this) { x *= val; } return *this; } void divide(polynomial &a, polynomial b) const { assert(!b.empty()); if (a.deg() < b.deg()) { a.clear(); return; } reverse(a.begin(), a.end()); int sz = a.deg() - b.deg() + 1; a %= sz; reverse(b.begin(), b.end()); a *= b.inv(sz); a %= sz; reverse(a.begin(), a.end()); } polynomial &operator/=(const polynomial &q) { divide(*this, q); normalize(); return *this; } polynomial &operator/=(T val) { val = 1 / val; for (auto &x : *this) { x *= val; } return *this; } polynomial &operator%=(const polynomial &q) { *this -= q * (*this / q); normalize(); return *this; } polynomial &operator%=(size_t k) { if (k <= deg()) this->resize(k); return *this; } polynomial &operator<<=(size_t k) { if (this->size() <= k) { this->clear(); } else { polynomial result(this->begin() + k, this->end()); this->swap(result); } return *this; } polynomial &operator>>=(size_t k) { polynomial result(this->size() + k); std::copy(this->begin(), this->end(), result.begin() + k); this->swap(result); return *this; } T operator()(const T &val) { T ans = 0; for (int i = deg(); i >= 0; i--) { ans = (ans * val + (*this)[i]); } return ans; } friend polynomial operator+(polynomial p, const polynomial &q) { p += q; return p; } friend polynomial operator-(polynomial p, const polynomial &q) { p -= q; return p; } friend polynomial operator*(polynomial p, const polynomial &q) { p *= q; return p; } friend polynomial operator*(polynomial p, const T &val) { p *= val; return p; } friend polynomial operator*(const T &val, polynomial p) { p *= val; return p; } friend polynomial operator/(polynomial p, const polynomial &q) { p /= q; return p; } friend polynomial operator/(polynomial p, const T &val) { p /= val; return p; } friend polynomial operator%(polynomial p, const polynomial &q) { p %= q; return p; } friend polynomial operator%(polynomial p, size_t k) { p %= k; return p; } friend polynomial operator<<(polynomial p, size_t k) { p <<= k; return p; } friend polynomial operator>>(polynomial p, size_t k) { p >>= k; return p; } polynomial inv(int k = -1) const { if (k == -1) k = this->size(); assert(!this->empty() && (*this)[0] != 0); polynomial b(1, 1 / (*this)[0]); for (int i = 2; i <= (k << 1); i <<= 1) { polynomial temp = (*this) % i; temp *= b, temp %= i; temp *= T(-1), temp[0] += 2; b *= temp, b %= i; } b.resize(k); return b; } polynomial derivative() const { if (deg() < 1) { return {}; } polynomial result(this->size() - 1); for (int i = 1; i < this->size(); i++) { result[i - 1] = i * (*this)[i]; } return result; } polynomial integral() const { if (this->empty()) { return {}; } polynomial result(this->size() + 1); for (int i = 0; i < this->size(); i++) { result[i + 1] = (*this)[i] / (i + 1); } return result; } polynomial log(int k = -1) const { assert(!this->empty() && (*this)[0] == 1); if (k == -1) k = this->size(); polynomial result = ((derivative() % k) * inv(k)).integral(); result.resize(k); return result; } polynomial exp(int k = -1) const { assert(this->empty() || (*this)[0] == 0); if (k == -1) k = this->size(); polynomial q(1, 1); for (int i = 2; i <= (k << 1); i <<= 1) { polynomial temp = polynomial(1, 1) + (*this % i) - q.log(i); q *= temp, q %= i; } q.resize(k); return q; } polynomial pow(int n, int k = -1) const { if (k == -1) k = this->size(); if (this->empty()) return polynomial(k); T alpha = 0; int pw = 0; for (int i = 0; i < this->size(); i++) { if ((*this)[i]) { alpha = (*this)[i]; pw = i; break; } } if ((long long)pw * n >= k) { return polynomial(k); } polynomial<T> b = (*this) << pw; b /= alpha; b = (T(n) * b.log(k)).exp(k); b >>= pw * n; b *= alpha.pow(n); b.resize(k); return b; } }; namespace evaluation_and_interpolation { template <typename T, typename U> std::vector<polynomial<T>> product_tree(const std::vector<U> &x) { int n = (int) x.size(); std::vector<polynomial<T>> tree(2 * n); for (int i = 0; i < n; i++) { tree[i + n] = {-T(x[i]), T(1)}; } for (int i = n - 1; i > 0; i--) { tree[i] = tree[i << 1] * tree[i << 1 | 1]; } return tree; } template <typename T, typename U> std::enable_if_t<!std::is_same<polynomial<T>, U>::value, std::vector<T>> multipoint_evaluation(const polynomial<T> &p, const std::vector<U> &x) { int n = (int) x.size(); auto tree = product_tree<T>(x); tree[1] = p % tree[1]; for (int i = 1; i < n; i++) { tree[i << 1] = tree[i] % tree[i << 1]; tree[i << 1 | 1] = tree[i] % tree[i << 1 | 1]; } std::vector<T> result(n); for (int i = 0; i < n; i++) { result[i] = (tree[i + n].empty() ? T(0) : tree[i + n][0]); } return result; } template <typename T, typename U> std::enable_if_t<std::is_same<polynomial<T>, U>::value, std::vector<T>> multipoint_evaluation(const polynomial<T> &p, const std::vector<U> &tree) { int n = (int) tree.size() / 2; std::vector<T> result(n); auto recurse = [&](int i, polynomial<T> q, const auto &self) -> void { q %= tree[i]; if (i >= n) { result[i - n] = q.empty() ? T(0) : q[0]; } else { self(i << 1, q, self); self(i << 1 | 1, q, self); } }; recurse(1, p, recurse); return result; } template <typename T, typename U> std::enable_if_t<!std::is_same<polynomial<T>, U>::value, polynomial<T>> interpolation(const std::vector<U> &x, const std::vector<U> &y) { int n = (int) x.size(); auto tree = product_tree<T>(x); auto ps = multipoint_evaluation(tree[1].derivative(), tree); auto recurse = [&](int i, const auto &self) -> polynomial<T> { if (i >= n) { return {T(y[i - n]) / ps[i - n]}; } else { return self(i << 1, self) * tree[i << 1 | 1] + self(i << 1 | 1, self) * tree[i << 1]; } }; return recurse(1, recurse); } template <typename T, typename U> std::enable_if_t<std::is_same<polynomial<T>, U>::value, polynomial<T>> interpolation(const std::vector<U> &tree, const std::vector<U> &y) { int n = (int) tree.size() / 2; auto ps = multipoint_evaluation(tree[1].derivative(), tree); auto recurse = [&](int i, const auto &self) -> polynomial<T> { if (i >= n) { return {T(y[i - n]) / ps[i - n]}; } else { return self(i << 1, self) * tree[i << 1 | 1] + self(i << 1 | 1, self) * tree[i << 1]; } }; return recurse(1, recurse); } } int main() { ios::sync_with_stdio(false); cin.tie(nullptr); using mint = montgomery_mint<998244353>; using evaluation_and_interpolation::interpolation; int n; cin >> n; polynomial<mint> x(n), y(n); for (auto &z : x) cin >> z; for (auto &z : y) cin >> z; auto poly = interpolation<mint>(x, y); for (int i = 0; i < n; i++) { cout << poly[i] << ' '; } cout << '\n'; return 0; }