Submit Info #44909

Problem Lang User Status Time Memory
Stirling Number of the First Kind cpp-acl opt AC 172 ms 13.22 MiB

ケース詳細
Name Status Time Memory
0_00 AC 2 ms 0.68 MiB
1_00 AC 0 ms 0.61 MiB
262143_00 AC 89 ms 7.17 MiB
262144_00 AC 150 ms 9.17 MiB
2_00 AC 2 ms 0.68 MiB
491519_00 AC 172 ms 12.75 MiB
499999_00 AC 172 ms 13.22 MiB
500000_00 AC 169 ms 13.21 MiB
5000_00 AC 3 ms 0.84 MiB
example_00 AC 1 ms 0.70 MiB

#include<bits/stdc++.h> using namespace std; #include<atcoder/all> using namespace atcoder; istream &operator>>(istream &is, modint998244353 &a) { long long v; is >> v; a = v; return is; } ostream &operator<<(ostream &os, const modint998244353 &a) { return os << a.val(); } istream &operator>>(istream &is, modint1000000007 &a) { long long v; is >> v; a = v; return is; } ostream &operator<<(ostream &os, const modint1000000007 &a) { return os << a.val(); } template<int m> istream &operator>>(istream &is, static_modint<m> &a) { long long v; is >> v; a = v; return is; } template<int m> istream &operator>>(istream &is, dynamic_modint<m> &a) { long long v; is >> v; a = v; return is; } template<int m> ostream &operator<<(ostream &os, const static_modint<m> &a) { return os << a.val(); } template<int m> ostream &operator<<(ostream &os, const dynamic_modint<m> &a) { return os << a.val(); } #define rep_(i, a_, b_, a, b, ...) for (int i = (a), lim##i = (b); i < lim##i; ++i) #define rep(i, ...) rep_(i, __VA_ARGS__, __VA_ARGS__, 0, __VA_ARGS__) #define drep_(i, a_, b_, a, b, ...) for (int i = (a)-1, lim##i = (b); i >= lim##i; --i) #define drep(i, ...) drep_(i, __VA_ARGS__, __VA_ARGS__, __VA_ARGS__, 0) using ll = long long; template<class T> istream &operator>>(istream &is, vector<T> &v) { for (auto &e : v) is >> e; return is; } template<class T> ostream &operator<<(ostream &os, const vector<T> &v) { for (auto &e : v) os << e << ' '; return os; } struct fast_ios { fast_ios(){ cin.tie(nullptr); ios::sync_with_stdio(false); cout << fixed << setprecision(20); }; } fast_ios_; // verified by: // https://judge.yosupo.jp/problem/convolution_mod // https://judge.yosupo.jp/problem/inv_of_formal_power_series // https://judge.yosupo.jp/problem/log_of_formal_power_series // https://judge.yosupo.jp/problem/exp_of_formal_power_series // https://judge.yosupo.jp/problem/pow_of_formal_power_series // https://judge.yosupo.jp/problem/polynomial_taylor_shift // https://judge.yosupo.jp/problem/bernoulli_number // https://judge.yosupo.jp/problem/sharp_p_subset_sum using mint = modint998244353; template<typename T> struct Factorial { int MAX; vector<T> fac, finv; Factorial(int m = 0) : MAX(m), fac(m+1, 1), finv(m+1, 1) { rep(i, 2, MAX+1) fac[i] = fac[i-1] * i; finv[MAX] /= fac[MAX]; drep(i, MAX+1, 3) finv[i-1] = finv[i] * i; } T binom(int n, int k) { if (k < 0 || n < k) return 0; return fac[n] * finv[k] * finv[n-k]; } T perm(int n, int k) { if (k < 0 || n < k) return 0; return fac[n] * finv[n-k]; } }; template<class T> struct formal_power_series : vector<T> { using vector<T>::vector; using vector<T>::operator=; using F = formal_power_series; F operator-() const { F res(*this); for (auto &e : res) e = -e; return res; } F &operator*=(const T &g) { for (auto &e : *this) e *= g; return *this; } F &operator/=(const T &g) { assert(g != T(0)); *this *= g.inv(); return *this; } F &operator+=(const F &g) { int n = this->size(), m = g.size(); rep(i, min(n, m)) (*this)[i] += g[i]; return *this; } F &operator-=(const F &g) { int n = this->size(), m = g.size(); rep(i, min(n, m)) (*this)[i] -= g[i]; return *this; } F &operator<<=(const int d) { int n = this->size(); if (d >= n) *this = F(n); this->insert(this->begin(), d, 0); this->resize(n); return *this; } F &operator>>=(const int d) { int n = this->size(); this->erase(this->begin(), this->begin() + min(n, d)); this->resize(n); return *this; } // O(n log n) F inv(int d = -1) const { int n = this->size(); assert(n != 0 && (*this)[0] != 0); if (d == -1) d = n; assert(d >= 0); F res{(*this)[0].inv()}; for (int m = 1; m < d; m *= 2) { F f(this->begin(), this->begin() + min(n, 2*m)); F g(res); f.resize(2*m), internal::butterfly(f); g.resize(2*m), internal::butterfly(g); rep(i, 2*m) f[i] *= g[i]; internal::butterfly_inv(f); f.erase(f.begin(), f.begin() + m); f.resize(2*m), internal::butterfly(f); rep(i, 2*m) f[i] *= g[i]; internal::butterfly_inv(f); T iz = T(2*m).inv(); iz *= -iz; rep(i, m) f[i] *= iz; res.insert(res.end(), f.begin(), f.begin() + m); } res.resize(d); return res; } // fast: FMT-friendly modulus only // O(n log n) F &multiply_inplace(const F &g, int d = -1) { int n = this->size(); if (d == -1) d = n; assert(d >= 0); *this = convolution(move(*this), g); this->resize(d); return *this; } F multiply(const F &g, const int d = -1) const { return F(*this).multiply_inplace(g, d); } // O(n log n) F &divide_inplace(const F &g, int d = -1) { int n = this->size(); if (d == -1) d = n; assert(d >= 0); *this = convolution(move(*this), g.inv(d)); this->resize(d); return *this; } F divide(const F &g, const int d = -1) const { return F(*this).divide_inplace(g, d); } // // naive // // O(n^2) // F &multiply_inplace(const F &g) { // int n = this->size(), m = g.size(); // drep(i, n) { // (*this)[i] *= g[0]; // rep(j, 1, min(i+1, m)) (*this)[i] += (*this)[i-j] * g[j]; // } // return *this; // } // F multiply(const F &g) const { return F(*this).multiply_inplace(g); } // // O(n^2) // F &divide_inplace(const F &g) { // assert(g[0] != T(0)); // T ig0 = g[0].inv(); // int n = this->size(), m = g.size(); // rep(i, n) { // rep(j, 1, min(i+1, m)) (*this)[i] -= (*this)[i-j] * g[j]; // (*this)[i] *= ig0; // } // return *this; // } // F divide(const F &g) const { return F(*this).divide_inplace(g); } // sparse // O(nk) F &multiply_inplace(vector<pair<int, T>> g) { int n = this->size(); auto [d, c] = g.front(); if (d == 0) g.erase(g.begin()); else c = 0; drep(i, n) { (*this)[i] *= c; for (auto &[j, b] : g) { if (j > i) break; (*this)[i] += (*this)[i-j] * b; } } return *this; } F multiply(const vector<pair<int, T>> &g) const { return F(*this).multiply_inplace(g); } // O(nk) F &divide_inplace(vector<pair<int, T>> g) { int n = this->size(); auto [d, c] = g.front(); assert(d == 0 && c != T(0)); T ic = c.inv(); g.erase(g.begin()); rep(i, n) { for (auto &[j, b] : g) { if (j > i) break; (*this)[i] -= (*this)[i-j] * b; } (*this)[i] *= ic; } return *this; } F divide(const vector<pair<int, T>> &g) const { return F(*this).divide_inplace(g); } // multiply and divide (1 + cz^d) // O(n) void multiply_inplace(const int d, const T c) { int n = this->size(); if (c == T(1)) drep(i, n-d) (*this)[i+d] += (*this)[i]; else if (c == T(-1)) drep(i, n-d) (*this)[i+d] -= (*this)[i]; else drep(i, n-d) (*this)[i+d] += (*this)[i] * c; } // O(n) void divide_inplace(const int d, const T c) { int n = this->size(); if (c == T(1)) rep(i, n-d) (*this)[i+d] -= (*this)[i]; else if (c == T(-1)) rep(i, n-d) (*this)[i+d] += (*this)[i]; else rep(i, n-d) (*this)[i+d] -= (*this)[i] * c; } // O(n) T eval(const T &a) const { T x(1), res(0); for (auto e : *this) res += e * x, x *= a; return res; } // O(n) F &integral_inplace() { int n = this->size(); assert(n > 0); if (n == 1) return *this = F{0}; this->insert(this->begin(), 0); this->pop_back(); vector<T> inv(n); inv[1] = 1; int p = T::mod(); rep(i, 2, n) inv[i] = - inv[p%i] * (p/i); rep(i, 2, n) (*this)[i] *= inv[i]; return *this; } F integral() const { return F(*this).integral_inplace(); } // O(n) F &derivative_inplace() { int n = this->size(); assert(n > 0); rep(i, 2, n) (*this)[i] *= i; this->erase(this->begin()); this->push_back(0); return *this; } F derivative() const { return F(*this).derivative_inplace(); } // O(n log n) F &log_inplace(int d = -1) { int n = this->size(); assert(n > 0 && (*this)[0] == 1); if (d == -1) d = n; assert(d >= 0); if (d < n) this->resize(d); F f_inv = this->inv(); this->derivative_inplace(); this->multiply_inplace(f_inv); this->integral_inplace(); return *this; } F log(const int d = -1) const { return F(*this).log_inplace(d); } // O(n log n) // https://arxiv.org/abs/1301.5804 (Figure 1, right) F &exp_inplace(int d = -1) { int n = this->size(); assert(n > 0 && (*this)[0] == 0); if (d == -1) d = n; assert(d >= 0); F g{1}, g_fft{1, 1}; (*this)[0] = 1; this->resize(d); F h_drv(this->derivative()); for (int m = 2; m < d; m *= 2) { // prepare F f_fft(this->begin(), this->begin() + m); f_fft.resize(2*m), internal::butterfly(f_fft); // Step 2.a' { F _g(m); rep(i, m) _g[i] = f_fft[i] * g_fft[i]; internal::butterfly_inv(_g); _g.erase(_g.begin(), _g.begin() + m/2); _g.resize(m), internal::butterfly(_g); rep(i, m) _g[i] *= g_fft[i]; internal::butterfly_inv(_g); _g.resize(m/2); _g /= T(-m) * m; g.insert(g.end(), _g.begin(), _g.begin() + m/2); } // Step 2.b'--d' F t(this->begin(), this->begin() + m); t.derivative_inplace(); { // Step 2.b' F r{h_drv.begin(), h_drv.begin() + m-1}; // Step 2.c' r.resize(m); internal::butterfly(r); rep(i, m) r[i] *= f_fft[i]; internal::butterfly_inv(r); r /= -m; // Step 2.d' t += r; t.insert(t.begin(), t.back()); t.pop_back(); } // Step 2.e' if (2*m < d) { t.resize(2*m); internal::butterfly(t); g_fft = g; g_fft.resize(2*m); internal::butterfly(g_fft); rep(i, 2*m) t[i] *= g_fft[i]; internal::butterfly_inv(t); t.resize(m); t /= 2*m; } else { // この場合分けをしても数パーセントしか速くならない F g1(g.begin() + m/2, g.end()); F s1(t.begin() + m/2, t.end()); t.resize(m/2); g1.resize(m), internal::butterfly(g1); t.resize(m), internal::butterfly(t); s1.resize(m), internal::butterfly(s1); rep(i, m) s1[i] = g_fft[i] * s1[i] + g1[i] * t[i]; rep(i, m) t[i] *= g_fft[i]; internal::butterfly_inv(t); internal::butterfly_inv(s1); rep(i, m/2) t[i+m/2] += s1[i]; t /= m; } // Step 2.f' F v(this->begin() + m, this->begin() + min<int>(d, 2*m)); v.resize(m); t.insert(t.begin(), m-1, 0); t.push_back(0); t.integral_inplace(); rep(i, m) v[i] -= t[m+i]; // Step 2.g' v.resize(2*m); internal::butterfly(v); rep(i, 2*m) v[i] *= f_fft[i]; internal::butterfly_inv(v); v.resize(m); v /= 2*m; // Step 2.h' rep(i, min(d-m, m)) (*this)[m+i] = v[i]; } return *this; } F exp(const int d = -1) const { return F(*this).exp_inplace(d); } // O(n log n) F &pow_inplace(const ll k, int d = -1) { int n = this->size(); if (d == -1) d = n; assert(d >= 0 && k >= 0); if (d == 0) return *this = F(0); if (k == 0) { *this = F(d); (*this)[0] = 1; return *this; } int l = 0; while (l < n && (*this)[l] == 0) ++l; if (l == n || l > (d-1)/k) return *this = F(d); T c{(*this)[l]}; this->erase(this->begin(), this->begin() + l); *this /= c; this->log_inplace(d - l*k); *this *= k; this->exp_inplace(); *this *= c.pow(k); this->insert(this->begin(), l*k, 0); return *this; } F pow(const ll k, const int d = -1) const { return F(*this).pow_inplace(k, d); } // O(n log n) F &shift_inplace(const T c) { int n = this->size(); Factorial<mint> fc(n); rep(i, n) (*this)[i] *= fc.fac[i]; reverse(this->begin(), this->end()); F g(n); T cp = 1; rep(i, n) g[i] = cp * fc.finv[i], cp *= c; this->multiply_inplace(g, n); reverse(this->begin(), this->end()); rep(i, n) (*this)[i] *= fc.finv[i]; return *this; } F shift(const T c) const { return F(*this).shift_inplace(c); } F operator*(const T &g) const { return F(*this) *= g; } F operator/(const T &g) const { return F(*this) /= g; } F operator+(const F &g) const { return F(*this) += g; } F operator-(const F &g) const { return F(*this) -= g; } F operator<<(const int d) const { return F(*this) <<= d; } F operator>>(const int d) const { return F(*this) >>= d; } // for multivariable FPS // F operator*(const F &g) const { return this->multiply(g); } // F operator*=(const F &g) { return this->multiply_inplace(g); } }; using fps = formal_power_series<mint>; // https://judge.yosupo.jp/problem/inv_of_formal_power_series void yosupo_inv() { int n; cin >> n; fps a(n); cin >> a; cout << a.inv() << '\n'; } // https://judge.yosupo.jp/problem/log_of_formal_power_series void yosupo_log() { int n; cin >> n; fps a(n); cin >> a; cout << a.log_inplace() << '\n'; } // https://judge.yosupo.jp/problem/exp_of_formal_power_series void yosupo_exp() { int n; cin >> n; fps a(n); cin >> a; cout << a.exp_inplace() << '\n'; } // https://judge.yosupo.jp/problem/pow_of_formal_power_series void yosupo_pow() { int n, m; cin >> n >> m; fps a(n); cin >> a; cout << a.pow_inplace(m) << '\n'; } // https://judge.yosupo.jp/problem/polynomial_taylor_shift void yosupo_shift() { int n, c; cin >> n >> c; fps a(n); cin >> a; cout << a.shift_inplace(c) << '\n'; } // https://yosupo.jp/problem/stirling_number_of_the_first_kind void yosupo_stirling1() { int n; cin >> n; auto stirling1 = [] (auto self, int n) { if (n < 64) { fps f{1}; rep(i, n) { f.push_back(0); f.multiply_inplace(1, -i); } reverse(f.begin(), f.end()); return f; } int m = n>>1; fps f = self(self, m); f.multiply_inplace(f.shift(-m), n+1); if (n&1) f.multiply_inplace({{0, -n+1}, {1, 1}}); return f; }; cout << stirling1(stirling1, n) << '\n'; } // https://judge.yosupo.jp/problem/bernoulli_number void yosupo_bernoulli() { int n; cin >> n; ++n; Factorial<mint> fc(n); fps f((fc.finv).begin()+1, (fc.finv).end()); f = f.inv(); rep(i, n) f[i] *= fc.fac[i]; cout << f << '\n'; } // https://judge.yosupo.jp/problem/sharp_p_subset_sum void yosupo_count_subset_sum() { int n, t; cin >> n >> t; ++t; vector<int> c(t); rep(i, n) { int s; cin >> s; ++c[s]; } vector<mint> inv(t); inv[1] = 1; int p = mint::mod(); rep(i, 2, t) inv[i] = - inv[p%i] * (p/i); fps f(t); rep(i, t) { if (c[i] == 0) continue; for (int j = 1, d = i; d < t; ++j, d += i) { if (j&1) f[d] += c[i] * inv[j]; else f[d] -= c[i] * inv[j]; } } f.exp_inplace(); f.erase(f.begin()); cout << f << '\n'; } int main() { // yosupo_inv(); // yosupo_log(); // yosupo_exp(); // yosupo_pow(); // yosupo_shift(); yosupo_stirling1(); // yosupo_bernoulli(); // yosupo_count_subset_sum(); }