# Submit Info #29851

Problem Lang User Status Time Memory
Polynomial Interpolation cpp hiiragi4000 AC 1958 ms 33.50 MiB

ケース詳細
Name Status Time Memory
example_00 AC 2 ms 0.67 MiB
example_01 AC 1 ms 0.57 MiB
max_random_00 AC 1957 ms 33.50 MiB
max_random_01 AC 1958 ms 33.39 MiB
random_00 AC 1812 ms 32.84 MiB
random_01 AC 1714 ms 19.40 MiB
random_02 AC 779 ms 14.32 MiB

#ifndef POLYNOMIAL_HH #define POLYNOMIAL_HH #ifndef CONVOLUTION_HH #define CONVOLUTION_HH #ifndef NUMBER_THEORY_HH #define NUMBER_THEORY_HH #include<array> #include<functional> #include<numeric> #include<ostream> #include<unordered_map> #include<utility> #include<vector> #include<cmath> #define U32 unsigned #define I64 long long #define U64 unsigned long long struct PrimeEnumerator{ static std::pair<std::vector<int>, std::function<bool(int)>> leq(int n){ // assert(n >= 0); int b = n/30+1; std::vector<bool> isp(b<<3, true); std::vector<int> primes{2, 3, 5}; primes.reserve(UB*30*b/std::log(30*b)); isp[0] = false; for(int i=0; i<b; ++i){ for(int j=0; j<8; ++j){ I64 d = 30ll*i + r[j]; if(d == 1) continue; if(isp[i<<3|j]) primes.push_back(d); for(int k=3; ; ++k){ I64 pd = primes[k]*d; if(pd >= 30*b) break; isp[pd/30<<3|rinv[pd%30]] = false; if(d%primes[k] == 0) break; } } } while(!primes.empty() && primes.back()>n){ primes.pop_back(); } auto is_prime = [isp=move(isp)](int i){ if(i <= 1) return false; if(i==2 || i==3 || i==5) return true; if(i%2==0 || i%3==0 || i%5==0) return false; return (bool)isp[i/30<<3|rinv[i%30]]; }; return {primes, is_prime}; } private: static constexpr int r[] = {1, 7, 11, 13, 17, 19, 23, 29}; static constexpr signed char rinv[] = { -1, 0, -1, -1, -1, -1, -1, 1, -1, -1, -1, 2, -1, 3, -1, -1, -1, 4, -1, 5, -1, -1, -1, 6, -1, -1, -1, -1, -1, 7 }; static constexpr double UB = 1.25506; }; 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){ // 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 inv_mod(I64 a, I64 mod){ // assert(a>0 && mod>0 && std::gcd(a, mod)==1); auto [x, y] = extgcd(a, mod); return (x%=mod)<0? x+mod: x; } constexpr U32 pow_mod(I64 a, I64 n, U32 mod){ // assert(mod > 0); if(mod == 1) return 0; if((a%=mod) < 0) a += mod; if(n < 0) a = inv_mod(a, mod); U64 res = 1, t = a; while(n){ if(n%2) res = res*t%mod; if(n/=2) t = t*t%mod; } return res; } constexpr bool miller_rabin(U32 n, U32 witness){ // assert(n>=3 && n%2==0); if((witness%=n) == 0) return true; U32 d = n-1; int r = 0; while(d%2 == 0){ d>>=1, ++r; } U32 x = pow_mod(witness, d, n); if(x==1 || x==n-1) return true; for(int cnt=r-1; cnt-->0; ){ x = (U64)x*x%n; if(x == n-1) return true; } return false; } constexpr bool primality_test(U32 n) noexcept{ if(n <= 1) return false; if(n == 2) return true; if(n%2 == 0) return false; return miller_rabin(n, 2) && (n < 2047 || miller_rabin(n, 7) && miller_rabin(n, 61)); } constexpr bool primality_test(int n) noexcept{ return n<0? false: primality_test((U32)n); } constexpr int jacobi_symbol(I64 a, I64 b){ // assert(b>=1 && b%2==1); if((a%=b) < 0) a += b; bool ok = true; while(a){ bool n2 = false; while(a%2 == 0){ a>>=1, n2=!n2; } if(n2 && (b%8==3 || b%8==5)){ ok = !ok; } if(a%4==3 && b%4==3){ ok = !ok; } int t = b%a; b = a; a = t; } return b>1? 0: ok? 1: -1; } constexpr int sqrt_modp(int a, int p){ // assert(primality_test(p)); a = ((I64)a%p+p)%p; if(p == 2) return a; if(a == 0) return 0; if(jacobi_symbol(a, p) == -1) return -1; if(p%4 == 3) return pow_mod(a, (p+1ll)/4, p); I64 k = 1; while(jacobi_symbol(k*k-a, p) == 1) ++k; if(k*k%p == a) return k; I64 w = ((k*k-a)%p+p)%p; I64 rc = 1, rx = 0, tc = k, tx = 1; for(I64 i=(p+1ll)/2; i>0; ){ if(i%2){ I64 rc2 = (rc*tc + rx*tx%p*w)%p; rx = (rc*tx + rx*tc)%p; rc = rc2; } if(i/=2){ I64 tc2 = (tc*tc + tx*tx%p*w)%p; tx = 2*tc*tx%p; tc = tc2; } } if(rc > p-rc) rc = p-rc; return rc; } constexpr bool primitive_root_modp_check(int r, int p){ // assert(primality_test(p)); if(r%p == 0) return false; int d = p-1; for(I64 q=2; q*q<=d; ++q) if(d%q == 0){ do d/=q; while(d%q==0); if(pow_mod(r, (p-1)/q, p) == 1) return false; } if(d>1 && pow_mod(r, (p-1)/d, p) == 1){ return false; } return true; } constexpr int primitive_root_modp(int p){ // assert(primality_test(p)); for(int i=1; ; ++i){ if(primitive_root_modp_check(i, p)){ return i; } } } inline int discrete_log(int base, int power, int mod){ // assert(mod > 0); if((base%=mod) < 0) base += mod; if((power%=mod) < 0) power += mod; if(mod==1 || power==1) return 0; int s = 1, d = 1, res = 0; while(1){ s = (I64)base*s%mod; ++res; if(s == power) return res; int d2 = std::gcd(s, mod); if(d2 == d){ if(power%d) return -1; s /= d; power /= d; mod /= d; break; } d = d2; } int m = std::ceil(std::sqrt(mod)); std::unordered_map<int, int> bs(m); for(int i=0; i<m; ++i){ bs.emplace(s, i); s = (I64)base*s%mod; } int gs = pow_mod(base, -m, mod); for(int i=0; i<m; ++i){ if(auto it=bs.find(power); it!=bs.end()){ return res + i*m + it->second; } power = (I64)gs*power%mod; } return -1; } template<U32 N> struct RingZn{ RingZn() = default; template<typename INT> constexpr RingZn(INT a) noexcept: a((a%=(I64)N)<0? (I64)a+N: a){} constexpr U32 mod() const noexcept{ return N; } template<typename INT> constexpr explicit operator INT() const noexcept{ return a; } constexpr RingZn operator+() const noexcept{ return *this; } constexpr RingZn operator-() const noexcept{ return a? N-a: 0; } constexpr RingZn &operator+=(RingZn rhs) noexcept{ a = a<N-rhs.a? a+rhs.a: a-(N-rhs.a); return *this; } constexpr RingZn &operator-=(RingZn rhs) noexcept{ a = a<rhs.a? a+(N-rhs.a): a-rhs.a; return *this; } constexpr RingZn &operator*=(RingZn rhs) noexcept{ a = (U64)a * rhs.a % N; return *this; } constexpr RingZn &operator/=(RingZn rhs) noexcept{ a = (U64)a * inv_mod(rhs.a, N) % N; return *this; } constexpr RingZn &operator++() noexcept{ a = a==N-1? 0: a+1; return *this; } constexpr RingZn operator++(int) noexcept{ RingZn res = *this; ++*this; return res; } constexpr RingZn &operator--() noexcept{ a = a? a-1: N-1; return *this; } constexpr RingZn operator--(int) noexcept{ RingZn res = *this; --*this; return res; } #define DEF_BIOP(op)\ friend constexpr RingZn operator op(RingZn lhs, RingZn rhs) noexcept{\ return lhs op##= rhs;\ } DEF_BIOP(+) DEF_BIOP(-) DEF_BIOP(*) DEF_BIOP(/) #undef DEF_BIOP friend constexpr bool operator==(RingZn lhs, RingZn rhs) noexcept{ return (U32)lhs == (U32)rhs; } friend constexpr bool operator!=(RingZn lhs, RingZn rhs) noexcept{ return !(lhs == rhs); } friend std::ostream &operator<<(std::ostream &os, RingZn a){ return os << (U32)a; } private: U32 a = 0; }; struct RingZnDyn{ RingZnDyn() = default; constexpr RingZnDyn(int t) noexcept: t(t){} template<typename INT> constexpr RingZnDyn(U32 n, INT a) noexcept: n(n), a((a%=(I64)n)<0? a+n: a){} constexpr U32 mod() const noexcept{ return n; } template<typename INT> constexpr explicit operator INT() const noexcept{ return n? (I64)a: (I64)t; } constexpr RingZnDyn operator+() const noexcept{ return *this; } constexpr RingZnDyn operator-() const noexcept{ RingZnDyn res = *this; if(n){ res.a = a? n-a: 0; }else res.t = -t; return res; } constexpr RingZnDyn &operator+=(RingZnDyn rhs){ if(!n){ if(!rhs.n){ t += rhs.t; return *this; } *this = RingZnDyn(rhs.n, t); }else if(!rhs.n){ rhs = RingZnDyn(n, rhs.t); } // else assert(n == rhs.n); a += a<n-rhs.a? rhs.a: n-rhs.a; return *this; } constexpr RingZnDyn &operator-=(RingZnDyn rhs){ if(!n){ if(!rhs.n){ t -= rhs.t; return *this; } *this = RingZnDyn(rhs.n, t); }else if(!rhs.n){ rhs = RingZnDyn(n, rhs.t); } // else assert(n == rhs.n); a = a<rhs.a? a+(n-rhs.a): a-rhs.a; return *this; } constexpr RingZnDyn &operator*=(RingZnDyn rhs){ if(!n){ if(!rhs.n){ t *= rhs.t; return *this; } *this = RingZnDyn(rhs.n, t); }else if(!rhs.n){ rhs = RingZnDyn(n, rhs.t); } // else assert(n == rhs.n); a = (U64)a * rhs.a % n; return *this; } constexpr RingZnDyn &operator/=(RingZnDyn rhs){ if(!n){ if(!rhs.n){ t /= rhs.t; return *this; } *this = RingZnDyn(rhs.n, t); }else if(!rhs.n){ rhs = RingZnDyn(n, rhs.t); } // else assert(n == rhs.n); a = (U64)a * inv_mod(rhs.a, n) % n; return *this; } constexpr RingZnDyn &operator++() noexcept{ if(n) a = a==n-1? 0: a+1; else ++t; return *this; } constexpr RingZnDyn operator++(int) noexcept{ RingZnDyn res = *this; ++*this; return res; } constexpr RingZnDyn &operator--() noexcept{ if(n) a = a==0? n-1: a-1; else --t; return *this; } constexpr RingZnDyn operator--(int) noexcept{ RingZnDyn res = *this; --*this; return res; } private: U32 n = 0; union{ int t = 0; U32 a; }; }; #define DEF_BIOP(op)\ constexpr RingZnDyn operator op(RingZnDyn lhs, RingZnDyn rhs){\ return lhs op##= rhs;\ } DEF_BIOP(+) DEF_BIOP(-) DEF_BIOP(*) DEF_BIOP(/) #undef DEF_BIOP constexpr bool operator==(RingZnDyn lhs, RingZnDyn rhs) noexcept{ if(!lhs.mod()){ if(!rhs.mod()){ return (int)lhs == (int)rhs; } lhs = RingZnDyn(rhs.mod(), (int)lhs); }else if(!rhs.mod()){ rhs = RingZnDyn(lhs.mod(), (int)rhs); }else if(lhs.mod() != rhs.mod()){ return false; } return (U32)lhs == (U32)rhs; } constexpr bool operator!=(RingZnDyn lhs, RingZnDyn rhs) noexcept{ return !(lhs == rhs); } #undef U64 #undef I64 #undef U32 #endif #include<algorithm> #include<complex> #include<iterator> #define I64 long long template<int P, int R=primitive_root_modp(P)> struct BasicNtt{ static_assert(primality_test(P)); static_assert(primitive_root_modp_check(R, P)); static constexpr int prime() noexcept{ return P; } static constexpr int primitive_root() noexcept{ return R; } static constexpr int lg_max_size() noexcept{ int d = P-1, res = 0; if((d&65535) == 0) d>>=16, res|=16; if((d&255) == 0) d>>=8, res|=8; if((d&15) == 0) d>>=4, res|=4; if((d&3) == 0) d>>=2, res|=2; if((d&1) == 0) res|=1; return res; } static constexpr int max_size() noexcept{ return 1<<lg_max_size(); } template<typename It> static void transform_in_place(bool inverse, It x, int n){ // assert(n>=0 && (n&-n)==n && n<=max_size()); using ResT = typename std::iterator_traits<It>::value_type; 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]); } RingZn<P> z2[lg_max_size()+1] = {pow_mod(R, inverse? (P-1)/n: -(P-1)/n, P)}; int lg = 1; for(; z2[lg-1]!=1; ++lg){ z2[lg] = (I64)z2[lg-1]*z2[lg-1]; } std::reverse(z2, z2+lg); for(int i=1; 1<<i<=n; ++i){ for(int j=0; j<n; j+=1<<i){ RingZn<P> t = 1; for(int k=0; k<1<<(i-1); ++k){ RingZn<P> u = x[j|k], v = t*x[j|1<<(i-1)|k]; x[j|k] = static_cast<ResT>(u+v); x[j|1<<(i-1)|k] = static_cast<ResT>(u-v); t *= z2[i]; } } } if(inverse){ RingZn<P> factor = P-(P-1)/n; for(int i=0; i<n; ++i){ x[i] = static_cast<ResT>(factor*x[i]); } } } template<typename It> static std::vector<typename std::iterator_traits<It>::value_type> transform(bool inverse, It b, It e, int n){ std::vector<typename std::iterator_traits<It>::value_type> x(n); int d = e-b<n? e-b: n; copy(b, b+d, x.begin()); transform_in_place(inverse, x.begin(), n); return x; } }; using Ntt = BasicNtt<(15<<27|1)>; 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; for(int i=1; i<n/8; ++i){ table[i] = {std::cos(2*PI*i/n), std::sin(2*PI*i/n)}; } T cos45 = std::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>; template<typename InputIt1, typename InputIt2, typename OutputIt> void convolution_int(int n, InputIt1 b1, InputIt1 e1, InputIt2 b2, InputIt2 e2, OutputIt res){ if(n <= 0) return; using ResT = typename std::iterator_traits<OutputIt>::value_type; int n1 = e1-b1, n2 = e2-b2; n1 = std::min(n1, n); n2 = std::min(n2, n); if(n1<=32 || n2<=32){ std::vector<ResT> t(n); for(int i=0; i<n1; ++i) for(int j=0; j<n2; ++j){ t[(i+j)%n] += static_cast<ResT>(b1[i])*static_cast<ResT>(b2[j]); } std::copy(t.cbegin(), t.cend(), res); return; } int n_fft = n; while((n_fft & -n_fft) != n_fft){ n_fft += n_fft & -n_fft; } if(n_fft>n && n<n1+n2-1){ n_fft <<= 1; } std::vector<std::complex<double>> x(n_fft); for(int i=0; i<n1; ++i){ x[i] = b1[i]; } for(int i=0; i<n2; ++i){ x[i] += std::complex<double>(0, b2[i]); } if(n_fft > 2*n){ for(int i=1; i<n1; ++i){ x[n_fft-n+i] = b1[i]; } } Fft::transform_in_place(false, x.begin(), n_fft); for(int i=0; i<n_fft; ++i){ x[i] *= x[i]; } Fft::transform_in_place(true, x.begin(), n_fft); for(int i=0; i<n; ++i){ res[i] = static_cast<ResT>(round(x[i].imag()/2)); } } #undef I64 #endif #include<algorithm> #include<functional> #include<initializer_list> #include<memory> #include<optional> #define U32 unsigned #define U64 unsigned long long namespace impl{ template<typename RingT> int fps_deg(std::vector<RingT> const &f) noexcept{ auto it = find_if(f.crbegin(), f.crend(), [](RingT const &ci){return ci != 0;}); return it.base()-f.cbegin()-1; } template<typename RingT> std::vector<RingT> &fps_add_assign(std::vector<RingT> &f, std::vector<RingT> const &g){ int d2 = fps_deg(g); if((int)f.size() < d2+1) f.resize(d2+1); for(int i=0; i<=d2; ++i){ f[i] += g[i]; } return f; } template<typename RingT> std::vector<RingT> &fps_sub_assign(std::vector<RingT> &f, std::vector<RingT> const &g){ int d2 = fps_deg(g); if((int)f.size() < d2+1) f.resize(d2+1); for(int i=0; i<=d2; ++i){ f[i] -= g[i]; } return f; } template<typename RingT> std::vector<RingT> fps_naive_mul(RingT const *f, int deg_f, RingT const *g, int deg_g, int lo, int hi){ std::vector<RingT> res(hi-lo+1, 0); 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] += f[i]*g[j]; } } return res; } template<typename RingT> struct FpsMulAssign{ std::vector<RingT> &operator()(std::vector<RingT> &f, std::vector<RingT> const &g){ int d1 = fps_deg(f), d2 = fps_deg(g); if(d1==-1 || d2==-1){ f.clear(); return f; } if(d1<32 || d2<32){ return f = fps_naive_mul(f.data(), d1, g.data(), d2, 0, d1+d2); } int n = std::max(d1, d2)/2+1; auto mid = f.cbegin() + std::min(n, (int)f.size()); std::vector<RingT> a(f.cbegin(), mid), b(mid, f.cend()); mid = g.cbegin() + std::min(n, (int)g.size()); std::vector<RingT> c(g.cbegin(), mid), d(mid, g.cend()); auto ac = a; (*this)(ac, c); auto bd = b; (*this)(bd, d); auto ab = a; fps_add_assign(ab, b); auto cd = c; fps_add_assign(cd, d); auto &cross = fps_sub_assign(fps_sub_assign((*this)(ab, cd), ac), bd); cross.insert(cross.cbegin(), n, 0); bd.insert(bd.cbegin(), 2*n, 0); return f = fps_add_assign(fps_add_assign(bd, cross), ac); } }; template<unsigned N> struct FpsMulAssign<RingZn<N>>{ std::vector<RingZn<N>> &operator()(std::vector<RingZn<N>> &f, std::vector<RingZn<N>> const &g){ int d1 = fps_deg(f), d2 = fps_deg(g); if(d1==-1 || d2==-1){ f.clear(); return f; } if(d1<32 || d2<32){ return f = fps_naive_mul(f.data(), d1, g.data(), d2, 0, d1+d2); } int n = d1+d2+1; while(n != (n&-n)) n += n&-n; if constexpr(primality_test(N)){ if(n <= BasicNtt<N>::max_size()){ f.resize(n); BasicNtt<N>::transform_in_place(false, f.begin(), n); if(f.data() == g.data()){ for(int i=0; i<n; ++i){ f[i] *= f[i]; } }else{ auto b = BasicNtt<N>::transform(false, g.cbegin(), g.cend(), n); for(int i=0; i<n; ++i){ f[i] *= b[i]; } } BasicNtt<N>::transform_in_place(true, f.begin(), n); f.resize(d1+d2+1); return f; } } std::vector<std::complex<double>> x(n); for(int i=0; i<=d1; ++i) x[i] = (double)f[i]; for(int i=0; i<=d2; ++i) x[i] += std::complex<double>(0, (double)g[i]); Fft::transform_in_place(false, x.begin(), n); for(int i=0; i<n; ++i){ x[i] *= x[i]; } Fft::transform_in_place(true, x.begin(), n); f.resize(d1+d2+1); if((double)n*N*N < 2e15){ for(int i=0; i<=d1+d2; ++i){ f[i] = (U64)round(x[i].imag()/2); } return f; } std::vector<RingZn<Ntt::prime()>> y(n), z(n); for(int i=0; i<=d1; ++i) y[i] = (U32)f[i]; for(int i=0; i<=d2; ++i) z[i] = (U32)g[i]; Ntt::transform_in_place(false, y.begin(), n); Ntt::transform_in_place(false, z.begin(), n); for(int i=0; i<n; ++i){ y[i] *= z[i]; } Ntt::transform_in_place(true, y.begin(), n); for(int i=0; i<=d1+d2; ++i){ U64 q = round((x[i].imag()/2-(U64)y[i])/Ntt::prime()); f[i] = (q%N)*Ntt::prime()+(U64)y[i]; } return f; } }; template<typename RingT> struct FpsReciprocal{ std::vector<RingT> operator()(std::vector<RingT> const &f, int n_terms){ if(n_terms <= 0) return {}; std::vector<RingT> res{1/f[0]}; for(int i=1; i<n_terms; i*=2){ std::vector<RingT> two{2}; std::vector<RingT> first_2i(f.cbegin(), f.cbegin()+std::min(2*i, (int)f.size())); FpsMulAssign<RingT> mul_assign; mul_assign(res, fps_sub_assign(two, mul_assign(first_2i, res))); if((int)res.size() > 2*i) res.resize(2*i); } if((int)res.size() > n_terms) res.resize(n_terms); return res; } }; } // namespace impl template<typename RingT> struct BasicFps: private std::vector<RingT>{ using Base = std::vector<RingT>; using Base::capacity; using Base::clear; using Base::get_allocator; using Base::reserve; using Base::resize; using Base::shrink_to_fit; using Base::size; BasicFps() = default; BasicFps(RingT const &c): BasicFps(c==0? Base{}: Base{c}){} BasicFps(std::initializer_list<RingT> coef): BasicFps(Base(coef)){} template<typename It> BasicFps(It b, It e): BasicFps(Base(b, e)){} explicit BasicFps(Base coef): Base(move(coef)){} std::vector<RingT> into_vec(){ std::vector<RingT> res; res.swap(*this); return res; } int deg() const noexcept{ return impl::fps_deg(*this); } explicit operator bool() const noexcept{ return deg() >= 0; } RingT coef(int i) const noexcept{ return i>=(int)size()? 0: (*this)[i]; } template<typename T> BasicFps &set_coef(int i, T &&ci){ if(i >= (int)size()) resize(i+1); (*this)[i] = std::forward<T>(ci); return *this; } RingT operator()(RingT const &x) const noexcept{ RingT res = 0; for(size_t i=size(); i-->0; ){ res = x*res + (*this)[i]; } return res; } BasicFps &operator<<=(int n_terms){ this->insert(this->cbegin(), n_terms, 0); return *this; } BasicFps operator<<(int n_terms) const{ auto res = *this; return res <<= n_terms; } BasicFps &operator>>=(int n_terms){ erase(this->cbegin(), this->cbegin()+std::min(n_terms, (int)size())); return *this; } BasicFps operator>>(int n_terms) const{ if(n_terms < (int)size()){ return BasicFps(this->cbegin()+n_terms, this->cend()); } return {}; } BasicFps &truncate(int n_terms) noexcept{ if(n_terms < (int)size()){ resize(n_terms); } return *this; } BasicFps first_n_terms(int n) const{ return BasicFps(this->cbegin(), this->cbegin()+std::min(n, (int)size())); } BasicFps rev() const{ int d = deg(); std::vector<RingT> res; res.reserve(d+1); reverse_copy(this->cbegin(), this->cbegin()+(d+1), back_inserter(res)); return BasicFps(move(res)); } BasicFps operator+() const{ return *this; } BasicFps operator-() const{ auto res = *this; transform(res.begin(), res.end(), res.begin(), std::negate<RingT>()); return res; } BasicFps &operator+=(BasicFps const &rhs){ impl::fps_add_assign(*this, rhs); return *this; } BasicFps &operator-=(BasicFps const &rhs){ impl::fps_sub_assign(*this, rhs); return *this; } BasicFps &operator*=(BasicFps const &rhs){ impl::FpsMulAssign<RingT>()(*this, rhs); return *this; } BasicFps reciprocal(int n_terms) const{ return BasicFps(impl::FpsReciprocal<RingT>()(*this, n_terms)); } BasicFps &operator/=(BasicFps const &rhs){ int n = deg(), d = rhs.deg(); if(d<32 || n-d<32){ resize(n+1); return *this = BasicFps(naive_division(std::move(*this), rhs.data(), d).first); } auto rf = rev(), rg = rhs.rev(); *this = rf.truncate(n-d+1)*rg.reciprocal(n-d+1); resize(n-d+1); reverse(this->begin(), this->end()); return *this; } BasicFps &operator%=(BasicFps const &rhs){ int n = deg(), d = rhs.deg(); if(d<32 || n-d<32){ resize(n+1); return *this = BasicFps(naive_division(std::move(*this), rhs.data(), d).second); } auto q = *this; q /= rhs; q *= rhs; truncate(d); q.truncate(d); return *this -= q; } BasicFps derivative() const{ BasicFps res; if(int d=deg(); d>0){ res.resize(d); for(int i=1; i<=d; ++i){ res[i-1] = i*(*this)[i]; } } return res; } BasicFps integral() const{ BasicFps res; if(int d=deg(); d>=0){ res.resize(d+2); for(int i=0; i<=d; ++i){ res[i+1] = (*this)[i]/(i+1); } } return res; } friend BasicFps operator*(BasicFps const &lhs, BasicFps const &rhs){ BasicFps res = lhs; return res *= std::addressof(lhs)==std::addressof(rhs)? res: rhs; } #define DEF_BIOP(op)\ friend BasicFps operator op(BasicFps lhs, BasicFps const &rhs){\ return lhs op##= rhs;\ } DEF_BIOP(+) DEF_BIOP(-) DEF_BIOP(/) DEF_BIOP(%) #undef DEF_BIOP friend bool operator==(BasicFps const &lhs, BasicFps const &rhs){ int d = lhs.deg(); return rhs.deg()==d && equal(lhs.cbegin(), lhs.cbegin()+(d+1), rhs.cbegin()); } friend bool operator!=(BasicFps const &lhs, BasicFps const &rhs){ return !(lhs == rhs); } private: static std::pair<std::vector<RingT>, std::vector<RingT>> naive_division(std::vector<RingT> f, RingT const *g, int deg_g){ int deg_f = (int)f.size()-1; if(deg_f < deg_g){ return {{}, move(f)}; } std::vector<RingT> q(deg_f-deg_g+1); for(int i=deg_f-deg_g; i>=0; --i){ RingT qi = f[i+deg_g]/g[deg_g]; f[i+deg_g] = 0; for(int j=0; j<deg_g; ++j){ f[i+j] -= qi*g[j]; } q[i] = qi; } f.resize(deg_g); return {move(q), move(f)}; } }; template<unsigned N> using FpsMod = BasicFps<RingZn<N>>; #undef U64 #undef U32 #endif #include<cstdio> using namespace std; constexpr int P = 998'244'353; void st_build(RingZn<P> const *x, int i, int l, int r, FpsMod<P> *st){ if(l == r){ st[i] = {-x[l], 1}; 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); st[i] = st[2*i] * st[2*i+1]; } void st_eval(RingZn<P> const *x, FpsMod<P> const *st, FpsMod<P> f, int i, int l, int r, RingZn<P> *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); } FpsMod<P> st_interpolate(RingZn<P> const *y, FpsMod<P> const *st, int i, int l, int r){ if(l == r){ return y[l]; } int m = l+(r-l)/2; auto g = st_interpolate(y, st, 2*i, l, m); auto h = st_interpolate(y, st, 2*i+1, m+1, r); return g*st[2*i+1]+h*st[2*i]; } int main(){ int n; scanf("%d", &n); vector<RingZn<P>> x(n), y(n), z(n); for(int i=0; i<n; ++i){ int xi; scanf("%d", &xi); x[i] = xi; } for(int i=0; i<n; ++i){ int yi; scanf("%d", &yi); y[i] = yi; } vector<FpsMod<P>> st(4*n); st_build(x.data(), 1, 0, n-1, st.data()); st_eval(x.data(), st.data(), st[1].derivative(), 1, 0, n-1, z.data()); for(int i=0; i<n; ++i){ y[i] /= z[i]; } auto f = st_interpolate(y.data(), st.data(), 1, 0, n-1); for(int i=0; i<n; ++i){ printf("%d%c", (int)f.coef(i), " \n"[i==n-1]); } }