Submit Info #26541

Problem Lang User Status Time Memory
Polynomial Interpolation cpp hiiragi4000 AC 7438 ms 39.20 MiB

ケース詳細
Name Status Time Memory
example_00 AC 2 ms 0.67 MiB
example_01 AC 1 ms 0.67 MiB
max_random_00 AC 7438 ms 39.20 MiB
max_random_01 AC 7430 ms 39.18 MiB
random_00 AC 6026 ms 33.69 MiB
random_01 AC 5339 ms 21.41 MiB
random_02 AC 2352 ms 14.95 MiB

#ifndef POLYNOMIAL_HH #define POLYNOMIAL_HH #ifndef CONVOLUTION_HH #define CONVOLUTION_HH #ifndef NUMBER_THEORY_HH #define NUMBER_THEORY_HH #include<array> #include<numeric> #include<utility> #include<vector> #define U32 unsigned #define I64 long long #define U64 unsigned long long inline std::pair<std::vector<bool>, std::vector<int>> prime_table_leq(int n){ // assert(n >= 1); std::vector<bool> isp(n+1u, true); std::vector<int> table; isp[0] = isp[1] = false; for(I64 i=2; i<=n; ++i){ if(isp[i]) table.push_back(i); for(int p: table){ if(p*i > n) break; isp[p*i] = false; if(i%p == 0) break; } } return {isp, table}; } inline std::vector<int> phi_table_leq(int n){ // assert(n >= -1); std::vector<int> phi(n+1u); std::vector<int> primes; iota(phi.begin(), phi.end(), 0u); for(I64 i=2; i<=n; ++i){ if(phi[i] == i){ --phi[i]; primes.push_back(i); } for(int p: primes){ if(p*i > n) break; if(i%p){ phi[p*i] = (p-1)*phi[i]; }else{ phi[p*i] = p*phi[i]; break; } } } return phi; } inline std::vector<signed char> mu_table_leq(int n){ // assert(n >= 1); std::vector<signed char> mu(n+1u, -2); std::vector<int> primes; mu[1] = 1; for(I64 i=2; i<=n; ++i){ if(mu[i] == -2){ mu[i] = -1; primes.push_back(i); } for(int p: primes){ if(p*i > n) break; if(i%p){ mu[p*i] = -mu[i]; }else{ mu[p*i] = 0; break; } } } return mu; } constexpr std::array<I64, 2> extgcd(I64 a, I64 b) noexcept{ // assert(a>=0 && b>=0); if(b == 0) return {1, 0}; auto [x, y] = extgcd(b, a%b); return {y, x-a/b*y}; } constexpr I64 mod_inv(I64 a, I64 n) noexcept{ // assert(a>0 && n>0); auto [x, y] = extgcd(a, n); // assert(a*x+n*y == 1); return (x%=n)<0? x+n: x; } constexpr int pow_mod(int a, I64 n, int m) noexcept{ // assert(m > 0); if(n == 0) return 1%m; a = ((I64)a%m+m)%m; if(n < 0) a = mod_inv(a, m); I64 res = 1, t = a; while(n){ if(n%2) res = res*t%m; if(n/=2) t = t*t%m; } return res; } struct RingZn{ RingZn() = default; constexpr RingZn(int t) noexcept: t(t){} template<typename T> constexpr RingZn(U32 n, T a) noexcept: n(n), a(a<0? (a%(I64)n+n)%n: a%n){} constexpr U32 order() const noexcept{ return n; } template<typename T> constexpr explicit operator T() const noexcept{ return n? (I64)a: (I64)t; } constexpr RingZn operator+() const noexcept{ return *this; } constexpr RingZn operator-() const noexcept{ return n? RingZn(n, a? n-a: 0): RingZn(-t); } constexpr RingZn &operator+=(RingZn rhs) noexcept{ if(!n){ if(!rhs.n){ t += rhs.t; return *this; } *this = RingZn(rhs.n, t); }else if(!rhs.n){ rhs = RingZn(n, rhs.t); } // else assert(n == rhs.n); a = ((U64)a+rhs.a)%n; return *this; } constexpr RingZn &operator-=(RingZn rhs) noexcept{ if(!n){ if(!rhs.n){ t -= rhs.t; return *this; } *this = RingZn(rhs.n, t); }else if(!rhs.n){ rhs = RingZn(n, rhs.t); } // else assert(n == rhs.n); a = ((U64)a+n-rhs.a)%n; return *this; } constexpr RingZn &operator*=(RingZn rhs) noexcept{ if(!n){ if(!rhs.n){ t *= rhs.t; return *this; } *this = RingZn(rhs.n, t); }else if(!rhs.n){ rhs = RingZn(n, rhs.t); } // else assert(n == rhs.n); a = (U64)a*rhs.a%n; return *this; } constexpr RingZn &operator/=(RingZn rhs) noexcept{ if(!n){ if(!rhs.n){ t *= rhs.t; return *this; } *this = RingZn(rhs.n, t); }else if(!rhs.n){ rhs = RingZn(n, rhs.t); } // else assert(n == rhs.n); a = (U64)a*mod_inv(rhs.a, n)%n; return *this; } constexpr RingZn &operator++() noexcept{ if(n) a = a==n-1? 0: a+1; else ++t; return *this; } constexpr RingZn operator++(int) noexcept{ RingZn res = *this; ++*this; return res; } constexpr RingZn &operator--() noexcept{ if(n) a = a==0? n-1: a-1; else --t; return *this; } constexpr RingZn operator--(int) noexcept{ RingZn res = *this; --*this; return res; } private: U32 n = 0; union{ int t = 0; U32 a; }; }; #define DEF_BIOP(op)\ constexpr RingZn operator op(RingZn lhs, RingZn rhs) noexcept{\ return lhs op##= rhs;\ } DEF_BIOP(+) DEF_BIOP(-) DEF_BIOP(*) DEF_BIOP(/) #undef DEF_BIOP constexpr bool operator==(RingZn lhs, RingZn rhs) noexcept{ if(!lhs.order()){ if(!rhs.order()){ return (int)lhs == (int)rhs; } lhs = RingZn(rhs.order(), (int)lhs); }else if(!rhs.order()){ rhs = RingZn(lhs.order(), (int)rhs); }else if(lhs.order() != rhs.order()){ return false; } return (U32)lhs == (U32)rhs; } constexpr bool operator!=(RingZn lhs, RingZn rhs) noexcept{ return !(lhs == rhs); } #undef U64 #undef I64 #undef U32 #endif #include<algorithm> #include<complex> #include<cmath> #define I64 long long template<int A, int B, int R> struct BasicNtt{ static constexpr int max_size() noexcept{ return 1<<B; } static constexpr int prime() noexcept{ return P; } static constexpr int primitive_root() noexcept{ return R; } template<typename It> static void transform_in_place(bool inverse, It x, int n){ // assert(n>=0 && (n&-n)==n && n<=max_size()); // assert(*std::min_element(x, x+n) >= 0); for(int i=1, j=0; i<n; ++i){ for(int k=n>>1; !((j^=k)&k); k>>=1); if(i < j) std::swap(x[i], x[j]); } int z2[B+1] = {pow_mod(R, inverse? (P-1)/n: -(P-1)/n, P)}, lg = 1; for(; z2[lg-1]!=1; ++lg){ z2[lg] = (I64)z2[lg-1]*z2[lg-1]%P; } std::reverse(z2, z2+lg); for(int i=1; 1<<i<=n; ++i){ for(int j=0; j<n; j+=1<<i){ for(int k=0, t=1; k<1<<(i-1); ++k, t=(I64)t*z2[i]%P){ I64 u = x[j|k], v = (I64)t*x[j|1<<(i-1)|k]%P; x[j|k] = (u+v)%P; x[j|1<<(i-1)|k] = (u+P-v)%P; } } } if(inverse){ for(int i=0; i<n; ++i){ x[i] = x[i]*(P-(P-1ll)/n)%P; } } } template<typename It> static std::vector<int> transform(bool inverse, It b, It e, int n){ std::vector<int> x(n); int d = e-b<n? e-b: n; for(int i=0; i<d; ++i){ x[i] = (b[i]%(I64)P+P)%P; } transform_in_place(inverse, x.begin(), n); return x; } private: static constexpr int P = (A<<B)+1; }; using Ntt = BasicNtt<15, 27, 31>; template<typename T> struct BasicFft{ template<typename It> static void transform_in_place(bool inverse, It x, int n){ // assert(n>=0 && (n&-n)==n); for(int i=1, j=0; i<n; ++i){ for(int k=n>>1; !((j^=k)&k); k>>=1); if(i < j) std::swap(x[i], x[j]); } auto exp_table = calc_exp(n); for(int i=1; 1<<i<=n; ++i){ int d = inverse? n>>i: n-(n>>i); for(int j=0; j<n; j+=1<<i){ for(int k=0, a=0; k<1<<(i-1); ++k, a=a+d&n-1){ std::complex<T> u = x[j|k], v = exp_table[a]*x[j|1<<(i-1)|k]; x[j|k] = u+v; x[j|1<<(i-1)|k] = u-v; } } } if(inverse){ for(int i=0; i<n; ++i){ x[i] /= n; } } } template<typename It> static std::vector<std::complex<T>> transform(bool inverse, It b, It e, int n){ std::vector<std::complex<T>> x(n); copy_n(b, e-b<n? e-b: n, x.begin()); transform_in_place(inverse, x.begin(), n); return x; } private: static constexpr T PI = 3.14159265358979323846l; static std::vector<std::complex<T>> calc_exp(int n){ std::vector<std::complex<T>> table; if(n == 1){ table = {1}; }else if(n == 2){ table = {1, -1}; }else if(n == 4){ table = {1, {0, 1}, -1, {0, -1}}; }else if(n >= 8){ table.resize(n); table[0] = {1, 0}; for(int i=1; i<n/8; ++i){ table[i] = {cos(2*PI*i/n), sin(2*PI*i/n)}; } T cos45 = sqrt((T).5); table[n/8] = {cos45, cos45}; for(int i=n/8+1; i<=n/4; ++i){ table[i] = {table[n/4-i].imag(), table[n/4-i].real()}; } for(int i=n/4+1; i<=n/2; ++i){ table[i] = {-table[n/2-i].real(), table[n/2-i].imag()}; } for(int i=n/2+1; i<n; ++i){ table[i] = -table[i-n/2]; } } return table; } }; using Fft = BasicFft<double>; using Fftf = BasicFft<float>; using Fftl = BasicFft<long double>; #undef I64 #endif #include<initializer_list> #include<iterator> #define I64 long long template<typename NttT> struct PolynomialNtt; template<typename NttT> bool operator==(PolynomialNtt<NttT> const&, PolynomialNtt<NttT> const&); template<typename NttT> struct PolynomialNtt: private std::vector<int>{ static constexpr int max_size() noexcept{ return NttT::max_size(); } static constexpr int mod() noexcept{ return NttT::prime(); } using Base = std::vector<int>; using Base::capacity; using Base::clear; using Base::get_allocator; using Base::reserve; using Base::resize; using Base::shrink_to_fit; using Base::size; PolynomialNtt() = default; PolynomialNtt(int c): PolynomialNtt(Base{c}){} PolynomialNtt(std::initializer_list<int> coef): PolynomialNtt(Base(coef)){} template<typename It> PolynomialNtt(It b, It e): PolynomialNtt(Base(b, e)){} explicit PolynomialNtt(Base coef): Base(move(coef)){ for(size_t i=0; i<size(); ++i){ data()[i] = ((I64)data()[i]%mod() + mod()) % mod(); } } int deg() const noexcept{ auto it = find_if(crbegin(), crend(), [](int ci){return ci;}); return it.base()-cbegin()-1; } explicit operator bool() const noexcept{ return deg() >= 0; } int coef(int i) const noexcept{ return i>=(int)size()? 0: data()[i]; } template<typename T> void set_coef(int i, T ci){ if(i >= (int)size()) resize(i+1); data()[i] = (int)RingZn(mod(), ci); } int operator()(int x) const noexcept{ RingZn res(mod(), 0); for(size_t i=size(); i-->0; ){ res = x*res + data()[i]; } return (int)res; } void truncate() noexcept{ resize(deg()+1); } PolynomialNtt rev() const{ int d = deg(); vector<int> res(d+1); reverse_copy(cbegin(), cbegin()+(d+1), res.begin()); return PolynomialNtt(move(res)); } PolynomialNtt operator+() const{ return *this; } PolynomialNtt operator-() const{ auto res = *this; for(size_t i=0; i<size(); ++i) if(res[i]){ res[i] = mod() - res[i]; } return res; } PolynomialNtt &operator+=(PolynomialNtt const &rhs){ int d2 = rhs.deg(); if((int)size() < d2+1) resize(d2+1); for(int i=0; i<=d2; ++i){ data()[i] = ((I64)data()[i] + rhs[i]) % mod(); } return *this; } PolynomialNtt &operator-=(PolynomialNtt const &rhs){ int d2 = rhs.deg(); if((int)size() < d2+1) resize(d2+1); for(int i=0; i<=d2; ++i){ data()[i] = ((I64)data()[i] + mod() - rhs[i]) % mod(); } return *this; } PolynomialNtt &operator*=(PolynomialNtt const &rhs){ int d1 = deg(), d2 = rhs.deg(), n = d1+d2+1; if(d1==-1 || d2==-1){ clear(); return *this; } if(d1<32 || d2<32){ return *this = PolynomialNtt(naive_multiply(data(), d1, rhs.data(), d2, 0, d1+d2)); } while((n&-n) != n) n += n&-n; resize(n); NttT::transform_in_place(false, begin(), n); if(this == &rhs){ for(int i=0; i<n; ++i){ data()[i] = (I64)data()[i] * data()[i] % mod(); } }else{ auto b = NttT::transform(false, rhs.cbegin(), rhs.cend(), n); for(int i=0; i<n; ++i){ data()[i] = (I64)data()[i] * b[i] % mod(); } } NttT::transform_in_place(true, begin(), n); return *this; } PolynomialNtt &operator/=(PolynomialNtt const &rhs){ auto rf = rev(), rg = rhs.rev(); int n = (int)rf.size()-1, d = (int)rg.size()-1; if(n < d){ clear(); return *this; } rf.resize(n-d+1); (*this = rf*rg.reciprocal(n-d+1)).resize(n-d+1); reverse(begin(), end()); return *this; } PolynomialNtt &operator%=(PolynomialNtt const &rhs){ int n = deg(), d = rhs.deg(); if(n < d) return *this; auto q = *this; q /= rhs; q *= rhs; resize(d); q.resize(d); return *this -= q; } PolynomialNtt reciprocal(int n_terms) const{ if(n_terms <= 0) return {}; int d = deg(), n = n_terms; while((n&-n) != n) n += n&-n; std::vector<int> res(n); res[0] = mod_inv(coef(0), mod()); for(int l=1; l<n; l<<=1){ std::vector<int> g2; if(l <= 32){ PolynomialNtt h(naive_multiply(data(), std::min(d, l-1), res.data(), l-1, l, 2*l-1)); if(d >= l){ h += PolynomialNtt(naive_multiply(data()+l, std::min(d-l, l-1), res.data(), l-1, 0, l-1)); } g2 = naive_multiply(res.data(), l-1, h.data(), l-1, 0, l-1); for(int i=0; i<l; ++i) if(g2[i]){ g2[i] = mod()-g2[i]; } }else{ auto f1 = NttT::transform(false, cbegin(), cbegin()+std::min(d+1, l), 2*l); auto g1 = NttT::transform(false, res.cbegin(), res.cbegin()+l, 2*l); std::vector<int> h(2*l); for(int i=0; i<2*l; ++i){ h[i] = ((I64)f1[i]*g1[i] + mod() - 1) % mod(); if(i&1 && h[i]) h[i] = mod()-h[i]; } if(d >= l){ auto f2 = NttT::transform(false, cbegin()+l, cbegin()+std::min(d+1, 2*l), 2*l); for(int i=0; i<2*l; ++i){ f2[i] = (I64)f2[i] * g1[i] % mod(); } NttT::transform_in_place(true, f2.begin(), 2*l); fill_n(f2.begin()+l, l, 0); NttT::transform_in_place(false, f2.begin(), 2*l); for(int i=0; i<2*l; ++i){ h[i] = ((I64)h[i] + f2[i]) % mod(); } } g2.resize(2*l); for(int i=0; i<2*l; ++i){ g2[i] = (I64)(mod()-g1[i]) * h[i] % mod(); } NttT::transform_in_place(true, g2.begin(), 2*l); } copy_n(g2.cbegin(), l, res.begin()+l); } res.resize(n_terms); return PolynomialNtt(move(res)); } friend bool operator==<>(PolynomialNtt const &lhs, PolynomialNtt const &rhs); private: static std::vector<int> naive_multiply(int const *f, int deg_f, int const *g, int deg_g, int lo, int hi){ std::vector<int> res(hi-lo+1); for(int i=std::max(lo-deg_g, 0); i<=std::min(deg_f, hi); ++i){ for(int j=std::max(lo-i, 0); j<=std::min(deg_g, hi-i); ++j){ res[i+j-lo] = (res[i+j-lo] + (I64)f[i]*g[j]) % mod(); } } return res; } }; #define DEF_BIOP(op) template<typename NttT>\ PolynomialNtt<NttT> operator op(PolynomialNtt<NttT> lhs, PolynomialNtt<NttT> const &rhs){\ return lhs op##= rhs;\ } DEF_BIOP(+) DEF_BIOP(-) DEF_BIOP(*) DEF_BIOP(/) DEF_BIOP(%) #undef DEF_BIOP template<typename NttT> bool operator==(PolynomialNtt<NttT> const &lhs, PolynomialNtt<NttT> const &rhs){ int d = lhs.deg(); return rhs.deg()==d && std::equal(lhs.cbegin(), lhs.cend(), rhs.cbegin()); } template<typename NttT> bool operator!=(PolynomialNtt<NttT> const &lhs, PolynomialNtt<NttT> const &rhs){ return !(lhs == rhs); } #undef I64 #endif #include<cstdio> using namespace std; using fx_t = PolynomialNtt<BasicNtt<119, 23, 3>>; void st_build(int const *x, int i, int l, int r, fx_t *st){ if(l == r){ st[i].resize(2); st[i].set_coef(1, 1); st[i].set_coef(0, -x[l]); return; } int m = l+(r-l)/2; st_build(x, 2*i, l, m, st); st_build(x, 2*i+1, m+1, r, st); if(i > 1){ st[i] = st[2*i] * st[2*i+1]; } } void st_eval(int const *x, fx_t const *st, fx_t f, int i, int l, int r, int *res){ if(r-l < 32){ for(int j=l; j<=r; ++j){ res[j] = f(x[j]); } return; } int m = l+(r-l)/2; st_eval(x, st, f%st[2*i], 2*i, l, m, res); st_eval(x, st, move(f%=st[2*i+1]), 2*i+1, m+1, r, res); } fx_t st_interpolate(int const *x, int *y, int *z, fx_t const *st, int i, int l, int r){ if(l == r){ return y[l]; } int m = l+(r-l)/2; st_eval(x, st, st[2*i+1], 2*i, l, m, z); st_eval(x, st, st[2*i], 2*i+1, m+1, r, z); for(int j=l; j<=r; ++j){ y[j] = y[j] * mod_inv(z[j], fx_t::mod()) % fx_t::mod(); } auto fl = st_interpolate(x, y, z, st, 2*i, l, m); auto fr = st_interpolate(x, y, z, st, 2*i+1, m+1, r); return fl*st[2*i+1]+fr*st[2*i]; } int main(){ int n; scanf("%d", &n); vector<int> x(n), y(n), z(n); for(int i=0; i<n; ++i){ scanf("%d", x.data()+i); } for(int i=0; i<n; ++i){ scanf("%d", y.data()+i); } vector<fx_t> st(4*n); st_build(x.data(), 1, 0, n-1, st.data()); auto f = st_interpolate(x.data(), y.data(), z.data(), st.data(), 1, 0, n-1); for(int i=0; i<n; ++i){ printf("%d%c", f.coef(i), " \n"[i==n-1]); } }