Submit Info #35870

Problem Lang User Status Time Memory
Polynomial Interpolation cpp (anonymous) AC 3098 ms 96.87 MiB

ケース詳細
Name Status Time Memory
example_00 AC 0 ms 0.68 MiB
example_01 AC 1 ms 0.62 MiB
max_random_00 AC 3098 ms 96.87 MiB
max_random_01 AC 3093 ms 96.83 MiB
random_00 AC 2897 ms 93.43 MiB
random_01 AC 1559 ms 50.53 MiB
random_02 AC 1126 ms 37.98 MiB

#include <bits/stdc++.h> using namespace std; #define rep(i, a, b) for(int i = a; i < (b); ++i) #define trav(i, a) for(auto& i : a) #define all(x) begin(x), end(x) #define sz(x) (int)(x).size() typedef long long ll; typedef pair<int, int> pii; typedef vector<int> vi; const ll mod = (119 << 23) + 1, root = 62; // = 998244353 ll modpow(ll b, ll e) { ll ans = 1; for (; e; b = b * b % mod, e /= 2) if (e & 1) ans = ans * b % mod; return ans; } ll inv(ll a, ll md) { return 1<a ? md - inv(md%a,a)*md/a : 1; } typedef vector<ll> vl; void ntt(vl &a) { int n = sz(a), L = 31 - __builtin_clz(n); static vl rt(2, 1); for (static int k = 2, s = 2; k < n; k *= 2, s++) { rt.resize(n); ll z[] = {1, modpow(root, mod >> s)}; rep(i,k,2*k) rt[i] = rt[i / 2] * z[i & 1] % mod; } vi rev(n); rep(i,0,n) rev[i] = (rev[i / 2] | (i & 1) << L) / 2; rep(i,0,n) if (i < rev[i]) swap(a[i], a[rev[i]]); for (int k = 1; k < n; k *= 2) for (int i = 0; i < n; i += 2 * k) rep(j,0,k) { ll z = rt[j + k] * a[i + j + k] % mod, &ai = a[i + j]; a[i + j + k] = ai - z + (z > ai ? mod : 0); ai += (ai + z >= mod ? z - mod : z); } } vl conv(const vl &a, const vl &b) { if (a.empty() || b.empty()) return {}; int s = sz(a) + sz(b) - 1, B = 32 - __builtin_clz(s), n = 1 << B; int inv = modpow(n, mod - 2); vl L(a), R(b), out(n); L.resize(n), R.resize(n); ntt(L), ntt(R); rep(i,0,n) out[-i & (n - 1)] = (ll)L[i] * R[i] % mod * inv % mod; ntt(out); return {out.begin(), out.begin() + s}; } struct Mod { ll v; static Mod as(ll v) { return Mod(v%mod); } Mod() : v(0) {} // WARNING: neg numbers explicit Mod(ll vv) : v(vv) { if(v >= mod) v -= mod; } // WARNING: neg numbers Mod operator+(Mod b) { return Mod(v + b.v); } Mod operator-(Mod b) { return Mod(v - b.v + mod); } Mod operator*(Mod b) { return Mod::as(v * b.v); } Mod operator/(Mod b) { return *this * Mod(inv(b.v, mod)); } Mod operator^(ll e) { return Mod(modpow(v,e)); } // optional explicit operator ll() const { return v; } }; // ============== BASIC ================= typedef Mod num; // or double typedef vector<num> poly; poly &operator+=(poly &a, const poly &b) { a.resize(max(sz(a), sz(b))); rep(i,0,sz(b)) a[i] = a[i] + b[i]; return a; } poly &operator-=(poly &a, const poly &b) { a.resize(max(sz(a), sz(b))); rep(i,0,sz(b)) a[i] = a[i] - b[i]; return a; } poly &operator*=(poly &a, const poly &b) { if (sz(a) + sz(b) < 1){ // optional poly res(sz(a) + sz(b) - 1); rep(i,0,sz(a)) rep(j,0,sz(b)) res[i + j] = (res[i + j] + a[i] * b[j]); return (a = res); } // use vd or convMod<mod> below depending on field auto res = conv(vl(all(a)), vl(all(b))); return (a = poly(all(res))); } poly operator*(poly a, const num b) { trav(i, a) i = i * b; return a; } #define OP(o) \ poly operator o(poly a, const poly& b) { \ poly c = a; return c o##= b; } OP(*) OP(+) OP(-); // ============== DIVISION / MOD ================= poly modK(poly a, int k) { return {a.begin(), a.begin() + min(k, sz(a))}; } poly inverse(poly A) { // assumes a[0] != 0 poly B = poly({num(1) / A[0]}); while (sz(B) < sz(A)) B = modK(B*(poly{num(2)} - modK(A,2*sz(B))*B), 2*sz(B)); return modK(B, sz(A)); } poly &operator/=(poly &a, poly b) { // assumes a[0] != 0 if (sz(a) < sz(b)) return a = {}; int s = sz(a) - sz(b) + 1; reverse(all(a)), reverse(all(b)); a.resize(s), b.resize(s); a *= inverse(b); a.resize(s), reverse(all(a)); return a; } OP(/); poly &operator%=(poly &a, const poly &b) { // assumes a[0] != 0 if (sz(a) < sz(b)) return a; poly c = (a / b) * b; a.resize(sz(b) - 1); rep(i,0,sz(a)) a[i] = a[i] - c[i]; return a; } OP(%); // ============== DERIV / INTEG ================= poly deriv(poly a) { if (a.empty()) return {}; poly b(sz(a) - 1); rep(i,1,sz(a)) b[i - 1] = a[i] * num(i); return b; } poly integr(poly a) { if (a.empty()) return {num(0)}; poly b(sz(a) + 1); b[1] = num(1); // for reals: b[i] = 1/i instead rep(i,2,sz(b)) b[i] = b[mod%i]*Mod(-mod/i+mod); rep(i,1,sz(b)) b[i] = a[i-1] * b[i]; return b; } // ============== EVAL / INTERP ================= vector<num> eval(const poly &a, const vector<num> &x) { int n = sz(x); if (!n) return {}; vector<poly> up(2 * n); rep(i, 0, n) up[i + n] = poly({num(0) - x[i], Mod(1)}); for (int i=n; --i;) up[i] = up[2*i] * up[2*i+1]; vector<poly> down(2 * n); down[1] = a % up[1]; rep(i, 2, 2*n) down[i] = down[i/2] % up[i]; vector<num> y(n); rep(i, 0, n) y[i] = down[i + n][0]; return y; } poly interp(vector<num> x, vector<num> y) { int n=sz(x); vector<poly> up(n*2); rep(i,0,n) up[i+n] = poly({num(0)-x[i], num(1)}); for(int i=n; --i;) up[i] = up[2*i]*up[2*i+1]; vector<num> a = eval(deriv(up[1]), x); vector<poly> down(2*n); rep(i,0,n) down[i+n] = poly({y[i]*(num(1)/a[i])}); for(int i=n; --i;) down[i] = down[i*2] * up[i*2+1] + down[i*2+1] * up[i*2]; return down[1]; } // ============== LOG / EXP ================= poly log(poly a) { // assumes a[0] = 1 return modK(integr(deriv(a) * inverse(a)), sz(a)); } poly exp(poly a) { // assumes a[0] = 0 poly b(1, num(1)); if (a.empty()) return b; while (sz(b) < sz(a)) { b.resize(sz(b) * 2); b *= (poly({num(1)}) + modK(a, sz(b)) - log(b)); b.resize(sz(b) / 2 + 1); } return modK(b, sz(a)); } // ============== POW ================= poly pow(poly a, ll m) { // assumes a != 0 or m != 0 int p = 0, n = sz(a); while (p < sz(a) && a[p].v == 0) ++p; if (m*p >= sz(a)) return poly(sz(a)); num j = a[p]; a = {a.begin() + p, a.end()}; a = a * (num(1) / j); a.resize(n); auto res = exp(log(a) * num(m)) * (j ^ m); res.insert(res.begin(), p*m, Mod(0)); return {res.begin(), res.begin()+n}; } int main() { cin.tie(0)->sync_with_stdio(0); cin.exceptions(cin.failbit); int n; cin >> n; vector<num> x(n), y(n); rep(i,0,n) { ll a; cin>>a; x[i] = Mod::as(a); } rep(i,0,n) { ll a; cin>>a; y[i] = Mod::as(a); } poly p = interp(x,y); rep(i,0,sz(p)) cout << ll(p[i]) << " "; cout << endl; }