Submit Info #27511

Problem Lang User Status Time Memory
Polynomial Interpolation cpp kcvlex AC 3503 ms 45.25 MiB

ケース詳細
Name Status Time Memory
example_00 AC 2 ms 0.67 MiB
example_01 AC 1 ms 0.68 MiB
max_random_00 AC 3503 ms 45.25 MiB
max_random_01 AC 3497 ms 45.20 MiB
random_00 AC 3394 ms 43.62 MiB
random_01 AC 1885 ms 25.09 MiB
random_02 AC 1310 ms 18.26 MiB

#include <limits> #include <initializer_list> #include <utility> #include <bitset> #include <tuple> #include <type_traits> #include <functional> #include <string> #include <array> #include <deque> #include <list> #include <queue> #include <stack> #include <vector> #include <map> #include <set> #include <unordered_map> #include <unordered_set> #include <iterator> #include <algorithm> #include <complex> #include <random> #include <numeric> #include <iostream> #include <iomanip> #include <sstream> #include <regex> #include <cassert> #include <cstddef> #include <variant> #define endl codeforces #define ALL(v) std::begin(v), std::end(v) #define ALLR(v) std::rbegin(v), std::rend(v) using ll = std::int64_t; using ull = std::uint64_t; using pii = std::pair<int, int>; using tii = std::tuple<int, int, int>; using pll = std::pair<ll, ll>; using tll = std::tuple<ll, ll, ll>; using size_type = ssize_t; template <typename T> using vec = std::vector<T>; template <typename T> using vvec = vec<vec<T>>; template <typename T> const T& var_min(const T &t) { return t; } template <typename T> const T& var_max(const T &t) { return t; } template <typename T, typename... Tail> const T& var_min(const T &t, const Tail&... tail) { return std::min(t, var_min(tail...)); } template <typename T, typename... Tail> const T& var_max(const T &t, const Tail&... tail) { return std::max(t, var_max(tail...)); } template <typename T, typename... Tail> void chmin(T &t, const Tail&... tail) { t = var_min(t, tail...); } template <typename T, typename... Tail> void chmax(T &t, const Tail&... tail) { t = var_max(t, tail...); } template <typename T, std::size_t Head, std::size_t... Tail> struct multi_dim_array { using type = std::array<typename multi_dim_array<T, Tail...>::type, Head>; }; template <typename T, std::size_t Head> struct multi_dim_array<T, Head> { using type = std::array<T, Head>; }; template <typename T, std::size_t... Args> using mdarray = typename multi_dim_array<T, Args...>::type; template <typename T, typename F, typename... Args> void fill_seq(T &t, F f, Args... args) { if constexpr (std::is_invocable<F, Args...>::value) { t = f(args...); } else { for (size_type i = 0; i < t.size(); i++) fill_seq(t[i], f, args..., i); } } template <typename T> vec<T> make_v(size_type sz) { return vec<T>(sz); } template <typename T, typename... Tail> auto make_v(size_type hs, Tail&&... ts) { auto v = std::move(make_v<T>(std::forward<Tail>(ts)...)); return vec<decltype(v)>(hs, v); } namespace init__ { struct InitIO { InitIO() { std::cin.tie(nullptr); std::ios_base::sync_with_stdio(false); std::cout << std::fixed << std::setprecision(30); } } init_io; } template <typename T> T ceil_pow2(T bound) { T ret = 1; while (ret < bound) ret *= 2; return ret; } template <typename T> T ceil_div(T a, T b) { return a / b + !!(a % b); } namespace math { template <typename T> constexpr T mul_id_ele() { if constexpr (std::is_fundamental<T>::value) { return T(1); } else { return T::mul_id_ele(); } } template <typename T> constexpr T add_id_ele() { if constexpr (std::is_fundamental<T>::value) { return T(0); } else { return T::add_id_ele(); } } template <typename T> constexpr T pow(const T &n, ll k) { T ret = mul_id_ele<T>(); T cur = n; while (k) { if (k & 1) ret *= cur; cur *= cur; k /= 2; } return ret; } template <typename T> typename std::enable_if<std::is_integral<T>::value, T>::type gcd(T a, T b) { return b ? gcd(a % b, b) : a; } } namespace math { template <ll Mod> struct Modint { constexpr Modint(ll x) noexcept : x((Mod + x % Mod) % Mod) { } constexpr Modint() noexcept : Modint(0) { } constexpr static Modint add_id_ele() { return Modint(0); } constexpr static Modint mul_id_ele() { return Modint(1); } constexpr ll& value() noexcept { return x; } constexpr ll value() const noexcept { return x; } constexpr Modint& operator+=(const Modint &oth) noexcept { x += oth.value(); if (Mod <= x) x -= Mod; return *this; } constexpr Modint& operator-=(const Modint &oth) noexcept { x += Mod - oth.value(); if (Mod <= x) x -= Mod; return *this; } constexpr Modint& operator*=(const Modint &oth) noexcept { x *= oth.value(); x %= Mod; return *this; } constexpr Modint& operator/=(const Modint &oth) noexcept { x *= oth.inv().value(); x %= Mod; return *this; } constexpr Modint operator+(const Modint &oth) const noexcept { return Modint(x) += oth; } constexpr Modint operator-(const Modint &oth) const noexcept { return Modint(x) -= oth; } constexpr Modint operator*(const Modint &oth) const noexcept { return Modint(x) *= oth; } constexpr Modint operator/(const Modint &oth) const noexcept { return Modint(x) /= oth; } constexpr Modint operator-() const noexcept { return Modint((x != 0) * (Mod - x)); } constexpr bool operator==(const Modint &oth) const noexcept { return value() == oth.value(); } template <typename T> constexpr typename std::enable_if<std::is_integral<T>::value, const Modint&>::type operator=(T t) noexcept { (*this) = Modint(std::forward<T>(t)); return *this; } constexpr Modint inv() const noexcept { return ::math::pow(*this, Mod - 2); } constexpr ll mod() const noexcept { return Mod; } private: ll x; }; template <ll Mod> Modint<Mod> inv(Modint<Mod> m) { return m.inv(); } template <ll Mod> std::istream& operator>>(std::istream &is, Modint<Mod> &m) { ll v; is >> v; m = v; return is; } template <ll Mod> std::ostream& operator<<(std::ostream &os, Modint<Mod> m) { os << m.value(); return os; } } namespace poly { template <typename T, typename Impl> struct convolution_interface { void multiply(vec<T> &a, vec<T> b) { static_cast<Impl*>(this)->multiply(a, b); } // TODO : value_type template <typename Container> void multiply_sparse(vec<T> &a, Container terms) { std::sort(ALLR(terms), [&](const auto &p, const auto &q) { return p.second < q.second; }); a.resize(a.size() + terms[0].second); for (size_type i = size_type(a.size() - 1); 0 <= i; i--) { auto cur = std::move(a[i]); a[i] = T(); for (const auto &[ v, deg ] : terms) if (i + deg < a.size()) a[i + deg] += cur * v; } } template <typename Container> void divide_sparse(vec<T> &a, Container terms) { std::sort(ALL(terms), [&](const auto &p, const auto &q) { return p.second < q.second; }); assert(terms[0].second == 0); for (size_type i = 0; i < size_type(a.size()); i++) { a[i] *= terms[0].first; for (const auto &[ v, deg ] : terms) if (deg && 0 <= i - deg) a[i] -= v * a[i - deg]; } } }; template <typename T> struct SimpleConvolution : convolution_interface<T, SimpleConvolution<T>> { void multiply(vec<T> &a, vec<T> b) { a.resize(a.size() + b.size() - 1); for (size_type i = size_type(a.size()) - 1; 0 <= i; i--) { auto cur = std::move(a[i]); a[i] = T(); for (size_type j = size_type(b.size()); 0 <= j; j--) { if (i + j < a.size()) a[i + j] += cur * b[j]; } } } }; } namespace poly { template <typename T, typename ConvolutionImpl> struct FPS : vec<T> { using value_type = T; using vec<T>::vec; using conv_type = convolution_interface<T, ConvolutionImpl>; static conv_type** get_conv() { static conv_type *conv = nullptr; return &conv; } static void set_conv(conv_type *conv) { *(get_conv()) = conv; } size_type degree() const { return size_type(this->size()) - 1; } FPS& reverse() { std::reverse(ALL(*this)); return *this; } FPS& extend(size_type deg) { this->resize(deg + 1); return *this; } FPS& prefix(size_type deg) { this->resize(deg + 1); return *this; } FPS& shift(size_type s) { size_type deg = degree(); for (size_type i = 0; i + s <= deg; i++) { (*this)[i] = (*this)[i + s]; } return prefix(deg - s); } FPS& middle(size_type l, size_type r) { return this->shift(l).prefix(r - l); } FPS& add(const FPS &rhs) { if (degree() < rhs.degree()) extend(rhs.degree()); for (size_type i = 0; i <= degree(); i++) (*this)[i] += rhs[i]; return *this; } FPS& sub(const FPS &rhs) { if (degree() < rhs.degree()) extend(rhs.degree()); for (size_type i = 0; i <= degree(); i++) (*this)[i] -= rhs[i]; return *this; } FPS& mul(FPS rhs) { (*get_conv())->multiply(*this, std::move(rhs)); return *this; } FPS& mul(T t) { for (auto &e : (*this)) e *= t; return *this; } FPS& shrink() { while (1 <= degree() && this->back() == T(0)) this->pop_back(); return *this; } // f(0) != 0 FPS& inv(size_type deg) { FPS orig(*this); prefix(0); (*this)[0] = (*this)[0].inv(); for (size_type i = 1; i <= deg; i *= 2) { FPS d(i * 2); for (size_type j = 0; j <= std::min(d.degree(), orig.degree()); j++) d[j] = orig[j]; d.mul(*this).mul(*this).prefix(i * 2 - 1); this->mul(T(2)) .sub(std::move(d)) .prefix(i * 2 - 1); } return prefix(deg); } // den[0] != 0 FPS& quo(FPS den) { if (degree() < den.degree()) { this->clear(); this->push_back(T(0)); return *this; } size_type sz = degree() - den.degree(); den.reverse().inv(sz); return this->reverse() .prefix(sz) .mul(std::move(den)) .prefix(sz) .reverse(); } // den[0] != 0 FPS& mod(FPS den) { auto q = *this; q.quo(den).mul(std::move(den)); return this->sub(std::move(q)) .shrink(); } FPS& diff() { for (ll i = 0; i + 1 <= degree() ; i++) (*this)[i] = (*this)[i + 1] * T(i + 1); if (1 <= degree()) prefix(degree() - 1); return *this; } FPS& integrate() { this->resize(this->size() + 1); for (ll i = degree(); 1 <= i; i--) (*this)[i] = (*this)[i - 1] * inv(i); (*this)[0] = 0; return *this; } template <typename Container> FPS& multiply_sparse(Container terms) { (*get_conv())->multiply_sparse(*this, terms); return *this; } template <typename Container> FPS& divide_sparse(Container terms) { (*get_conv())->divide_sparse(*this, terms); return *this; } }; template <typename Poly> Poly prefix(Poly f, size_type sz) { return f.prefix(sz); } template <typename Poly> Poly inv(Poly f, size_type sz) { return f.inv(sz); } template <typename Poly> Poly add(Poly lhs, const Poly &rhs) { return lhs.add(rhs); } template <typename Poly> Poly sub(Poly lhs, const Poly &rhs) { return lhs.sub(rhs); } template <typename Poly> Poly mul(Poly lhs, const Poly &rhs) { return lhs.mul(rhs); } template <typename Poly> Poly quo(Poly lhs, const Poly &rhs) { return lhs.quo(rhs); } template <typename Poly> Poly mod(Poly lhs, const Poly &rhs) { return lhs.mod(rhs); } template <typename Poly> Poly mul(Poly lhs, Poly &&rhs) { return lhs.mul(std::move(rhs)); } template <typename Poly> Poly quo(Poly lhs, Poly &&rhs) { return lhs.quo(std::move(rhs)); } template <typename Poly> Poly mod(Poly lhs, Poly &&rhs) { return lhs.mod(std::move(rhs)); } template <typename Poly1, typename Poly2> std::pair<Poly1, Poly1> div(Poly1 lhs, Poly2 &&rhs) { auto q = std::move(quo(lhs, rhs)); auto m = std::move(mul(q, std::forward<Poly2>(rhs))); lhs.sub(std::move(m)).shrink(); return std::make_pair(std::move(q), std::move(lhs)); } // http://q.c.titech.ac.jp/docs/progs/polynomial_division.html // [x^N](P/Q) template <typename Poly> typename Poly::value_type coef_of_div(Poly P, Poly Q, size_type N) { assert(Q[0] == 1); assert(P.size() + 1 == Q.size()); for (; N; N /= 2) { FPS Qn = Q; for (size_type i = 1; i < Q.size(); i += 2) Qn[i] *= -1; auto PQ = std::move(mul(P, Qn)); auto QQ = std::move(mul(Q, Qn)); for (size_type i = 0, offset = N % 2; i < P.size(); i++) P[i] = PQ[2 * i + offset]; for (size_type i = 0; i < Q.size(); i++) Q[i] = QQ[2 * i]; } return P[0]; } } namespace poly { namespace ntt_helper { template <ll Mod> constexpr bool is_primitive_root(ll r) { math::Modint<Mod> mr(r); for (ll d = 2; d * d <= Mod; d++) { if ((Mod - 1) % d == 0) { if (pow(mr, d).value() == 1) return false; if (pow(mr, (Mod - 1) / d).value() == 1) return false; } } return true; } template <ll Mod> constexpr ll find_primitive_root(ll r) { return (is_primitive_root<Mod>(r) ? r : find_primitive_root<Mod>(r + 1)); } template <ll Mod> constexpr ll find_primitive_root() { return find_primitive_root<Mod>(2); } constexpr auto calc_max_base(ll m) { ll ret = 0; for (; m % 2 == 0; ret++, m /= 2); return ret; } } template <ll Mod, ll PrimitiveRoot> struct ntt__ : convolution_interface<math::Modint<Mod>, ntt__<Mod, PrimitiveRoot>> { using mint = math::Modint<Mod>; using value_type = mint; constexpr ntt__() : root_lis(max_size_log), iroot_lis(max_size_log) { for (size_type i = 0; i < static_cast<size_type>(root_lis.size()); i++) { root_lis[i] = pow(mint(PrimitiveRoot), (Mod - 1) >> (i + 2)) * -1; iroot_lis[i] = root_lis[i].inv(); } } void multiply(vec<value_type> &a, vec<value_type> b) { size_type m = a.size(), n = b.size(); size_type sz = 1, res_sz = n + m - 1; while (sz < res_sz) sz *= 2; a.resize(sz); b.resize(sz); ntt(a, sz, false); ntt(b, sz, false); auto isz = mint(sz).inv(); for (size_type i = 0; i < sz; i++) a[i] *= b[i] * isz; ntt(a, sz, true); a.resize(res_sz); } private: static constexpr size_type max_size_log = ntt_helper::calc_max_base(Mod - 1); static constexpr size_type max_size = 1ll << max_size_log; static constexpr size_type max_conv_size = max_size * 2; vec<mint> root_lis, iroot_lis; using buf_iterator = typename vec<mint>::iterator; void ntt(vec<value_type> &v, size_type sz, bool inv) { if (!inv) { for (int m = sz / 2; m; m /= 2) { mint mul = 1; for(int s = 0, k = 0; s < sz; s += 2 * m) { for(int i = s, j = s + m; i < s + m; ++i, ++j) { auto x = v[i], y = v[j] * mul; v[i] = x + y; v[j] = x - y; } mul *= root_lis[__builtin_ctz(++k)]; } } } else { for (int m = 1; m < sz; m *= 2) { mint mul = 1; for (int s = 0, k = 0; s < sz; s += 2 * m) { for (int i = s, j = s + m; i < s + m; i++, j++) { auto l = v[i], r = v[j]; v[i] = l + r; v[j] = (l - r) * mul; } mul *= iroot_lis[__builtin_ctz(++k)]; } } } } }; template <ll Mod> using NTT = ntt__<Mod, ntt_helper::find_primitive_root<Mod>()>; } namespace poly { template <typename Poly> struct SubproductTree { using value_type = typename Poly::value_type; template <typename Container> SubproductTree(const Container &cs) { rebuild(cs); } size_type size() const noexcept { return sz; } size_type height() const noexcept { return rank[0]; } size_type get_len(size_type r) const { return len_sum[r + 1] - len_sum[r]; } // rank, index size_type get_poly_idx(size_type r, size_type idx) const { return idx + (len_sum.back() - len_sum[r + 1]); } const Poly& get_poly(size_type r, size_type idx) const { return nodes[get_poly_idx(r, idx)]; } private: vec<Poly> nodes; vec<ll> rank; vec<size_type> len_sum; size_type sz, nsz; template <typename Container> void rebuild(const Container &cs) { resize(cs.size()); size_type offset = nsz - sz; for (size_type i = 0; i < sz; i++) { nodes[i + offset] = Poly { value_type(-cs[i]), value_type(1) }; rank[i + offset] = 0; } size_type len = sz; size_type end = offset; while (1 < len) { len = (len + 1) / 2; size_type begin = end - len; for (size_type i = 0; i < len; i++) { size_type c1 = end + 2 * i; size_type c2 = end + 2 * i + 1; if (nsz <= c2 || rank[c1] != rank[c2]) { nodes[begin + i] = nodes[c1]; } else { nodes[begin + i] = std::move(mul(nodes[c1], nodes[c2])); } rank[begin + i] = rank[c1] + 1; } end = begin; } len_sum.resize(rank[0] + 2); for (size_type i = 0; i < nsz; i++) len_sum[rank[i] + 1]++; for (size_type i = 0; i + 1 < rank[0] + 2; i++) len_sum[i + 1] += len_sum[i]; } void resize(size_type sz_arg) { sz = sz_arg; nsz = sz_arg; while (1 < sz_arg) { auto nxt = (sz_arg + 1) / 2; sz_arg = nxt; nsz += nxt; } nodes.resize(nsz); rank.resize(nsz); } }; } namespace poly { namespace internal { template <typename Poly> struct MultipointEvaluationSolver { using value_type = typename Poly::value_type; MultipointEvaluationSolver(SubproductTree<Poly> *ptree) : ptree(ptree), ret(ptree->size()) { } const vec<value_type>& solve(Poly f) { dfs(std::move(f)); return ret; } private: const SubproductTree<Poly> *ptree; vec<value_type> ret; void dfs(const Poly &f) { auto r = ceil_pow2(ret.size()); dfs(f, 0, r, ptree->height(), 0); } void dfs(Poly f, size_type l, size_type r, size_type rank, size_type idx) { if (ret.size() <= l) return; if (r - l == 1) { ret[l] = f[0]; return; } auto m = (l + r) / 2; auto rec = [&](ll nl, ll nr, ll nidx) { if (ptree->get_len(rank - 1) <= nidx) return; auto g = std::move(mod(f, ptree->get_poly(rank - 1, nidx))); dfs(std::move(g), nl, nr, rank - 1, nidx); }; rec(l, m, 2 * idx + 0); rec(m, r, 2 * idx + 1); } }; } template <typename Poly> struct MultipointEvaluation { using value_type = typename Poly::value_type; template <typename Container> MultipointEvaluation(const Container &cs) : ptree(cs), solver(&ptree) { } const vec<value_type>& eval(Poly f) { return solver.solve(std::move(f)); } private: SubproductTree<Poly> ptree; internal::MultipointEvaluationSolver<Poly> solver; }; } namespace poly { template <typename Poly> Poly interpolation(const vec<typename Poly::value_type> &xv, const vec<typename Poly::value_type> &yv) { SubproductTree<Poly> ptree(xv); auto f = ptree.get_poly(ptree.height(), 0); f.diff(); auto z = internal::MultipointEvaluationSolver(&ptree).solve(std::move(f)); vec<Poly> b(ptree.size()); for (size_type i = 0; i < size_type(z.size()); i++) b[i] = Poly { yv[i] / z[i] }; for (size_type i = 0; i < ptree.height(); i++) { auto len = ptree.get_len(i); for (size_type j = 0; j < len / 2; j++) { auto m1 = std::move(mul(ptree.get_poly(i, 2 * j), std::move(b[2 * j + 1]))); auto m2 = std::move(mul(ptree.get_poly(i, 2 * j + 1), std::move(b[2 * j]))); b[j] = std::move(add(std::move(m1), std::move(m2))); } if (len & 1) b[len / 2] = std::move(b[len - 1]); } return b[0]; } } constexpr ll mod = 998244353; using mint = math::Modint<mod>; using ntt_type = poly::NTT<mod>; using poly_type = poly::FPS<mint, ntt_type>; ntt_type ntt; int main() { poly_type::set_conv(&ntt); ll n; std::cin >> n; vec<mint> xv(n), yv(n); for (auto &&e : xv) std::cin >> e; for (auto &&e : yv) std::cin >> e; auto ans = poly::interpolation<poly_type>(xv, yv); for (ll i = 0; i < n; i++) std::cout << ans[i] << " \n"[i + 1 == n]; return 0; }