Submit Info #21866

Problem Lang User Status Time Memory
Polynomial Interpolation cpp FlowerOfSorrow AC 1952 ms 27.32 MiB

ケース詳細
Name Status Time Memory
example_00 AC 1 ms 0.67 MiB
example_01 AC 0 ms 0.62 MiB
max_random_00 AC 1951 ms 27.32 MiB
max_random_01 AC 1952 ms 27.32 MiB
random_00 AC 1820 ms 26.88 MiB
random_01 AC 1525 ms 15.38 MiB
random_02 AC 755 ms 11.40 MiB

#include <bits/stdc++.h> using namespace std; template<typename T> struct Z_p{ using Type = typename decay<decltype(T::value)>::type; static vector<Type> mod_inv; constexpr Z_p(): value(){ } template<typename U> Z_p(const U &x){ value = normalize(x); } template<typename U> static Type normalize(const U &x){ Type v; if(-mod() <= x && x < mod()) v = static_cast<Type>(x); else v = static_cast<Type>(x % mod()); if(v < 0) v += mod(); return v; } const Type& operator()() const{ return value; } template<typename U> explicit operator U() const{ return static_cast<U>(value); } constexpr static Type mod(){ return T::value; } Z_p &operator+=(const Z_p &otr){ if((value += otr.value) >= mod()) value -= mod(); return *this; } Z_p &operator-=(const Z_p &otr){ if((value -= otr.value) < 0) value += mod(); return *this; } template<typename U> Z_p &operator+=(const U &otr){ return *this += Z_p(otr); } template<typename U> Z_p &operator-=(const U &otr){ return *this -= Z_p(otr); } Z_p &operator++(){ return *this += 1; } Z_p &operator--(){ return *this -= 1; } Z_p operator++(int){ Z_p result(*this); *this += 1; return result; } Z_p operator--(int){ Z_p result(*this); *this -= 1; return result; } Z_p operator-() const{ return Z_p(-value); } template<typename U = T> typename enable_if<is_same<typename Z_p<U>::Type, int>::value, Z_p>::type &operator*=(const Z_p& rhs){ #ifdef _WIN32 uint64_t x = static_cast<int64_t>(value) * static_cast<int64_t>(rhs.value); uint32_t xh = static_cast<uint32_t>(x >> 32), xl = static_cast<uint32_t>(x), d, m; asm( "divl %4; \n\t" : "=a" (d), "=d" (m) : "d" (xh), "a" (xl), "r" (mod()) ); value = m; #else value = normalize(static_cast<int64_t>(value) * static_cast<int64_t>(rhs.value)); #endif return *this; } template<typename U = T> typename enable_if<is_same<typename Z_p<U>::Type, int64_t>::value, Z_p>::type &operator*=(const Z_p &rhs){ int64_t q = static_cast<int64_t>(static_cast<long double>(value) * rhs.value / mod()); value = normalize(value * rhs.value - q * mod()); return *this; } template<typename U = T> typename enable_if<!is_integral<typename Z_p<U>::Type>::value, Z_p>::type &operator*=(const Z_p &rhs){ value = normalize(value * rhs.value); return *this; } Z_p operator^(long long e) const{ Z_p b = *this, res = 1; if(e < 0) b = 1 / b, e = -e; for(; e; b *= b, e >>= 1) if(e & 1) res *= b; return res; } Z_p &operator^=(const long long &e){ return *this = *this ^ e; } Z_p &operator/=(const Z_p &otr){ Type a = otr.value, m = mod(), u = 0, v = 1; if(a < int(mod_inv.size())) return *this *= mod_inv[a]; while(a){ Type t = m / a; m -= t * a; swap(a, m); u -= t * v; swap(u, v); } assert(m == 1); return *this *= u; } template<typename U> friend const Z_p<U> &abs(const Z_p<U> &v){ return v; } template<typename U> friend bool operator==(const Z_p<U> &lhs, const Z_p<U> &rhs); template<typename U> friend bool operator<(const Z_p<U> &lhs, const Z_p<U> &rhs); template<typename U> friend istream &operator>>(istream &in, Z_p<U> &number); Type value; }; template<typename T> bool operator==(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value == rhs.value; } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator==(const Z_p<T>& lhs, U rhs){ return lhs == Z_p<T>(rhs); } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator==(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) == rhs; } template<typename T> bool operator!=(const Z_p<T> &lhs, const Z_p<T> &rhs){ return !(lhs == rhs); } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator!=(const Z_p<T> &lhs, U rhs){ return !(lhs == rhs); } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator!=(U lhs, const Z_p<T> &rhs){ return !(lhs == rhs); } template<typename T> bool operator<(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value < rhs.value; } template<typename T> Z_p<T> operator+(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) += rhs; } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator+(const Z_p<T> &lhs, U rhs){ return Z_p<T>(lhs) += rhs; } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator+(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) += rhs; } template<typename T> Z_p<T> operator-(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) -= rhs; } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator-(const Z_p<T>& lhs, U rhs){ return Z_p<T>(lhs) -= rhs; } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator-(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) -= rhs; } template<typename T> Z_p<T> operator*(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) *= rhs; } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator*(const Z_p<T>& lhs, U rhs){ return Z_p<T>(lhs) *= rhs; } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator*(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) *= rhs; } template<typename T> Z_p<T> operator/(const Z_p<T> &lhs, const Z_p<T> &rhs) { return Z_p<T>(lhs) /= rhs; } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator/(const Z_p<T>& lhs, U rhs) { return Z_p<T>(lhs) /= rhs; } template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator/(U lhs, const Z_p<T> &rhs) { return Z_p<T>(lhs) /= rhs; } template<typename T> istream &operator>>(istream &in, Z_p<T> &number){ typename common_type<typename Z_p<T>::Type, int64_t>::type x; in >> x; number.value = Z_p<T>::normalize(x); return in; } template<typename T> ostream &operator<<(ostream &out, const Z_p<T> &number){ return out << number(); } /* using ModType = int; struct VarMod{ static ModType value; }; ModType VarMod::value; ModType &mod = VarMod::value; using Zp = Z_p<VarMod>; */ constexpr int mod = 998244353; using Zp = Z_p<integral_constant<decay<decltype(mod)>::type, mod>>; template<typename T> vector<typename Z_p<T>::Type> Z_p<T>::mod_inv; template<typename T = integral_constant<decay<decltype(mod)>::type, mod>> void precalc_inverse(size_t SZ){ auto &inv = Z_p<T>::mod_inv; if(inv.empty()) inv.assign(2, 1); for(; inv.size() <= SZ; ) inv.push_back((mod - 1LL * mod / int(inv.size()) * inv[mod % int(inv.size())] % mod) % mod); } template<int root = 15311432, int root_pw = 1 << 23, int inv_root = 469870224, typename IT = vector<Zp>::iterator> void number_theoric_transform(IT begin, IT end, const bool invert = false){ int n = end - begin; for(int i = 1, j = 0; i < n; ++ i){ int bit = n >> 1; for(; j & bit; bit >>= 1) j ^= bit; j ^= bit; if(i < j) swap(*(begin + i), *(begin + j)); } for(int len = 1; len < n; len <<= 1){ typename iterator_traits<IT>::value_type wlen = invert ? inv_root : root; for(int i = len << 1; i < root_pw; i <<= 1) wlen *= wlen; for(int i = 0; i < n; i += len << 1){ typename iterator_traits<IT>::value_type w = 1; for(int j = 0; j < len; ++ j){ auto u = *(begin + i + j), v = *(begin + i + j + len) * w; *(begin + i + j) = u + v; *(begin + i + j + len) = u - v; w *= wlen; } } } if(invert){ auto inv_n = typename iterator_traits<IT>::value_type(1) / n; for(auto it = begin; it != end; ++ it) *it *= inv_n; } } const size_t magic_constant = 250; template<typename Poly> void multiply(Poly &a, Poly b){ if(min(a.size(), b.size()) < magic_constant){ Poly res((int)a.size() + (int)b.size() - 1); for(size_t i = 0; i < a.size(); ++ i) for(size_t j = 0; j < b.size(); ++ j) res[i + j] += a[i] * b[j]; a = move(res); return; } int n = max((int)a.size() + (int)b.size() - 1, 1); if(__builtin_popcount(n) != 1) n = 1 << __lg(n) + 1; a.resize(n), b.resize(n); number_theoric_transform(a.begin(), a.end()), number_theoric_transform(b.begin(), b.end()); for(int i = 0; i < n; ++ i) a[i] *= b[i]; number_theoric_transform(a.begin(), a.end(), 1); while(!a.empty() && !a.back()) a.pop_back(); } template<typename T> T bpow(T x, size_t n){ T res = 1; for(; n; x *= x, n >>= 1) if(n & 1) res *= x; return res; } template<typename T> struct polynomial{ vector<T> a; void normalize(){ // get rid of leading zeroes while(!a.empty() && a.back() == T(0)) a.pop_back(); } polynomial(){} polynomial(T a0): a{a0}{ normalize(); } polynomial(const vector<T> &t): a(t){ normalize(); } polynomial &operator=(const polynomial &t){ a = t.a; return *this; } polynomial &operator+=(const polynomial &t){ a.resize(max(a.size(), t.a.size())); for(size_t i = 0; i < t.a.size(); ++ i) a[i] += t.a[i]; normalize(); return *this; } polynomial &operator-=(const polynomial &t){ a.resize(max(a.size(), t.a.size())); for(size_t i = 0; i < t.a.size(); ++ i) a[i] -= t.a[i]; normalize(); return *this; } polynomial operator+(const polynomial &t) const{ return polynomial(*this) += t; } polynomial operator-(const polynomial &t) const{ return polynomial(*this) -= t; } polynomial mod_xk(size_t k) const{ // get same polynomialnomial mod x^k return vector<T>(begin(a), begin(a) + min(k, a.size())); } polynomial mul_xk(size_t k) const{ // multiply by x^k polynomial res(*this); res.a.insert(begin(res.a), k, 0); return res; } polynomial div_xk(size_t k) const{ // divide by x^k, dropping coefficients return vector<T>(begin(a) + min(k, a.size()), end(a)); } polynomial substr(size_t l, size_t r) const{ // return mod_xk(r).div_xk(l) return vector<T>(begin(a) + min(l, a.size()), begin(a) + min(r, a.size())); } polynomial inv(size_t n) const{ // get inverse series mod x^n assert(!is_zero()); polynomial ans = 1 / a[0]; for(size_t i = 1; i < n; i <<= 1) ans = (ans * 2 - ans * ans * mod_xk(i << 1)).mod_xk(i << 1); return ans.mod_xk(n); } polynomial &operator*=(const polynomial &t){ multiply(a, t.a); normalize(); return *this; } polynomial operator*(const polynomial &t) const{ return polynomial(*this) *= t; } polynomial reverse(size_t n, bool rev = 0) const{ // reverses and leaves only n terms polynomial res(*this); if(rev) res.a.resize(max(n, res.a.size())); // If rev = 1 then tail goes to head std::reverse(res.a.begin(), res.a.end()); return res.mod_xk(n); } pair<polynomial, polynomial> divmod_slow(const polynomial &b) const{ // when divisor or quotient is small vector<T> A(a), res; for(T invx = 1 / b.a.back(); A.size() >= b.a.size(); ){ res.push_back(A.back() * invx); if(res.back() != T(0)) for(size_t i = 0; i < b.a.size(); ++ i) A[(int)A.size() - i - 1] -= res.back() * b.a[(int)b.a.size() - i - 1]; A.pop_back(); } std::reverse(begin(res), end(res)); return {res, A}; } pair<polynomial, polynomial> divmod(const polynomial &b) const{ // returns quotiend and remainder of a mod b if(deg() < b.deg()) return {polynomial{0}, *this}; int d = deg() - b.deg(); if(min(d, b.deg()) < magic_constant) return divmod_slow(b); polynomial D = (reverse(d + 1) * b.reverse(d + 1).inv(d + 1)).mod_xk(d + 1).reverse(d + 1, 1); return {D, *this - D * b}; } polynomial operator/(const polynomial &t) const{ return divmod(t).first; } polynomial operator%(const polynomial &t) const{ return divmod(t).second; } polynomial &operator/=(const polynomial &t){ return *this = divmod(t).first; } polynomial &operator%=(const polynomial &t){ return *this = divmod(t).second; } polynomial &operator*=(const T &x){ for(auto &it: a) it *= x; normalize(); return *this; } polynomial &operator/=(const T &x){ T invx = 1 / x; for(auto &it: a) it *= invx; normalize(); return *this; } polynomial operator*(const T &x) const{ return polynomial(*this) *= x; } polynomial operator/(const T &x) const{ return polynomial(*this) /= x; } void print() const{ for(auto it: a) cout << it << ' '; cout << "\n"; } T eval(T x) const{ // evaluates in single point x T res(0); for(int i = (int)a.size() - 1; i >= 0; -- i) res = res * x + a[i]; return res; } T &lead(){ // leading coefficient return a.back(); } int deg() const{ // degree return a.empty() ? numeric_limits<int>::min() / 2 : (int)a.size() - 1; } bool is_zero() const{ // is polynomialnomial zero return a.empty(); } T coef(int idx) const{ return idx >= (int)a.size() || idx < 0 ? T(0) : a[idx]; } T &operator[](int idx){ // mutable reference at idx return a[idx]; } bool operator==(const polynomial &t) const{ return a == t.a; } bool operator!=(const polynomial &t) const{ return a != t.a; } polynomial derivative() const{ // calculate derivative static vector<T> res; res.clear(); for(int i = 1; i <= deg(); ++ i) res.push_back(i * a[i]); return res; } polynomial antiderivative() const{ // calculate integral with C = 0 static vector<T> res; res.assign(1, 0); for(int i = 0; i <= deg(); ++ i) res.push_back(a[i] / (i + 1)); return res; } size_t leading_xk() const{ // Let p(x) = x^k * t(x), return k if(is_zero()) return numeric_limits<int>::max() / 2; int res = 0; while(a[res] == T(0)) ++ res; return res; } polynomial log(size_t n){ // calculate log p(x) mod x^n assert(a[0] == T(1)); return (derivative().mod_xk(n) * inv(n)).antiderivative().mod_xk(n); } polynomial exp(size_t n){ // calculate exp p(x) mod x^n if(is_zero()) return T(1); assert(a[0] == T(0)); polynomial ans = 1; for(size_t a = 1; a < n; a <<= 1){ polynomial C = ans.log(a << 1).div_xk(a) - substr(a, a << 1); ans -= (ans * C).mod_xk(a).mul_xk(a); } return ans.mod_xk(n); } polynomial pow_slow(size_t k, size_t n){ // if k is small return k ? k % 2 ? (*this * pow_slow(k - 1, n)).mod_xk(n) : (*this * *this).mod_xk(n).pow_slow(k / 2, n) : T(1); } polynomial pow(size_t k, size_t n){ // calculate p^k(n) mod x^n if(is_zero())return *this; if(k < magic_constant) return pow_slow(k, n); int i = leading_xk(); T j = a[i]; polynomial t = div_xk(i) / j; return bpow(j, k) * (t.log(n) * T(k)).exp(n).mul_xk(i * k).mod_xk(n); } polynomial mulx(T x){ // component-wise multiplication with x^k T cur = 1; polynomial res(*this); for(int i = 0; i <= deg(); ++ i) res.coef(i) *= cur, cur *= x; return res; } polynomial mulx_sq(T x){ // component-wise multiplication with x^{k^2} T cur = x, total = 1, xx = x * x; polynomial res(*this); for(int i = 0; i <= deg(); ++ i) res.coef(i) *= total, total *= cur, cur *= xx; return res; } vector<T> chirpz_even(T z, int n){ // P(1), P(z^2), P(z^4), ..., P(z^2(n-1)) int m = deg(); if(is_zero()) return vector<T>(n, 0); vector<T> vv(m + n); T zi = 1 / z, zz = zi * zi, cur = zi, total = 1; for(int i = 0; i <= max(n - 1, m); ++ i){ if(i <= m){ vv[m - i] = total; } if(i < n){ vv[m + i] = total; } total *= cur, cur *= zz; } polynomial w = (mulx_sq(z) * vv).substr(m, m + n).mulx_sq(z); vector<T> res(n); for(int i = 0; i < n; ++ i) res[i] = w[i]; return res; } vector<T> chirpz(T z, int n){ // P(1), P(z), P(z^2), ..., P(z^(n-1)) auto even = chirpz_even(z, n + 1 >> 1); auto odd = mulx(z).chirpz_even(z, n >> 1); vector<T> ans(n); for(int i = 0; i < n >> 1; ++ i) ans[i << 1] = even[i], ans[i << 1 | 1] = odd[i]; if(n & 1) ans[n - 1] = even.back(); return ans; } template<typename iter> vector<T> eval(vector<polynomial> &tree, int u, iter l, iter r){ // auxiliary evaluation function if(r - l == 1) return {eval(*l)}; else{ auto m = l + (r - l >> 1); int v = u + 1, w = u + (m - l << 1); auto A = (*this % tree[v]).eval(tree, v, l, m); auto B = (*this % tree[w]).eval(tree, w, m, r); A.insert(end(A), begin(B), end(B)); return A; } } template<typename iter> vector<T> eval(iter begin, iter end){ // evaluate polynomialnomial in (x1, ..., xn) int n = end - begin; if(is_zero()) return vector<T>(n, T(0)); vector<polynomial> tree(n << 1); build(tree, 0, begin, end); return eval(tree, 0, begin, end); } template<typename iter> polynomial interpolate(const vector<polynomial> &tree, int u, iter l, iter r, iter ly, iter ry){ // auxiliary interpolation function if(r - l == 1) return {*ly / a[0]}; else{ auto m = l + (r - l >> 1); auto my = ly + (ry - ly >> 1); int v = u + 1, w = u + (m - l << 1); auto A = (*this % tree[v]).interpolate(tree, v, l, m, ly, my); auto B = (*this % tree[w]).interpolate(tree, w, m, r, my, ry); return A * tree[w] + B * tree[v]; } } }; template<typename T> ostream &operator<<(ostream &out, const polynomial<T> &p){ for(auto x: p.a) out << x << " "; return out; } template<typename T> polynomial<T> operator*(const T& a, const polynomial<T> &b){ auto res(b); for(auto &it: res.a) it = a * it; res.normalize(); return res; } template<typename T> polynomial<T> xk(int k){ // return x^k return polynomial<T>{1}.mul_xk(k); } template<typename T> T resultant(polynomial<T> a, polynomial<T> b){ // computes resultant of a and b if(b.is_zero()) return 0; else if(b.deg() == 0) return bpow(b.lead(), a.deg()); else{ int pw = a.deg(); a %= b; pw -= a.deg(); T mul = bpow(b.lead(), pw) * T((b.deg() & a.deg() & 1) ? -1 : 1); T ans = resultant(b, a); return ans * mul; } } template<typename iter> polynomial<typename iter::value_type> generate(iter L, iter R){ // computes (x-a1)(x-a2)...(x-an) without building tree if(R - L == 1) return vector<typename iter::value_type>{-*L, 1}; else{ iter M = L + (R - L >> 1); return generate(L, M) * generate(M, R); } } template<typename T, typename iter> polynomial<T> &build(vector<polynomial<T>> &res, int u, iter L, iter R){ // builds evaluation tree for (x-a1)(x-a2)...(x-an) if(R - L == 1) return res[u] = vector<T>{-*L, 1}; else{ iter M = L + (R - L >> 1); return res[u] = build(res, u + 1, L, M) * build(res, u + (M - L << 1), M, R); } } template<typename T> polynomial<T> interpolate(const vector<T> &x, const vector<T> &y){ // interpolates minimum polynomialnomial from (xi, yi) pairs int n = x.size(); vector<polynomial<T>> tree(n << 1); return build(tree, 0, begin(x), end(x)).derivative().interpolate(tree, 0, begin(x), end(x), begin(y), end(y)); } using poly = polynomial<Zp>; // using poly = polynomial<double>; int main(){ cin.tie(0)->sync_with_stdio(0); cin.exceptions(ios::badbit | ios::failbit); int n; cin >> n; vector<Zp> x(n), y(n); for(auto i = 0; i < n; ++ i){ cin >> x[i]; } for(auto i = 0; i < n; ++ i){ cin >> y[i]; } auto res = interpolate(x, y); res.a.resize(n); cout << res << "\n"; return 0; } /* */ //////////////////////////////////////////////////////////////////////////////////////// // // // Coded by Aeren // // // ////////////////////////////////////////////////////////////////////////////////////////