Submit Info #46653

Problem Lang User Status Time Memory
$\sum_{i=0}^{\infty} r^i i^d$ cpp14 hitonanode AC 1626 ms 232.66 MiB

ケース詳細
Name Status Time Memory
0_00 AC 2 ms 0.68 MiB
0_01 AC 2 ms 0.68 MiB
0_02 AC 1621 ms 232.66 MiB
2_00 AC 1 ms 0.66 MiB
2_01 AC 1 ms 0.60 MiB
2_02 AC 132 ms 51.52 MiB
2_03 AC 166 ms 80.87 MiB
2_04 AC 1626 ms 232.60 MiB
2_05 AC 1613 ms 232.59 MiB
example_00 AC 1 ms 0.71 MiB

#line 1 "formal_power_series/test/sum_of_exponential_times_polynomial_limit.test.cpp" #define PROBLEM "https://judge.yosupo.jp/problem/sum_of_exponential_times_polynomial_limit" #line 2 "modint.hpp" #include <iostream> #include <set> #include <vector> // CUT begin template <int mod> struct ModInt { #if __cplusplus >= 201402L #define MDCONST constexpr #else #define MDCONST #endif using lint = long long; MDCONST static int get_mod() { return mod; } static int get_primitive_root() { static int primitive_root = 0; if (!primitive_root) { primitive_root = [&]() { std::set<int> fac; int v = mod - 1; for (lint i = 2; i * i <= v; i++) while (v % i == 0) fac.insert(i), v /= i; if (v > 1) fac.insert(v); for (int g = 1; g < mod; g++) { bool ok = true; for (auto i : fac) if (ModInt(g).pow((mod - 1) / i) == 1) { ok = false; break; } if (ok) return g; } return -1; }(); } return primitive_root; } int val; MDCONST ModInt() : val(0) {} MDCONST ModInt &_setval(lint v) { return val = (v >= mod ? v - mod : v), *this; } MDCONST ModInt(lint v) { _setval(v % mod + mod); } MDCONST explicit operator bool() const { return val != 0; } MDCONST ModInt operator+(const ModInt &x) const { return ModInt()._setval((lint)val + x.val); } MDCONST ModInt operator-(const ModInt &x) const { return ModInt()._setval((lint)val - x.val + mod); } MDCONST ModInt operator*(const ModInt &x) const { return ModInt()._setval((lint)val * x.val % mod); } MDCONST ModInt operator/(const ModInt &x) const { return ModInt()._setval((lint)val * x.inv() % mod); } MDCONST ModInt operator-() const { return ModInt()._setval(mod - val); } MDCONST ModInt &operator+=(const ModInt &x) { return *this = *this + x; } MDCONST ModInt &operator-=(const ModInt &x) { return *this = *this - x; } MDCONST ModInt &operator*=(const ModInt &x) { return *this = *this * x; } MDCONST ModInt &operator/=(const ModInt &x) { return *this = *this / x; } friend MDCONST ModInt operator+(lint a, const ModInt &x) { return ModInt()._setval(a % mod + x.val); } friend MDCONST ModInt operator-(lint a, const ModInt &x) { return ModInt()._setval(a % mod - x.val + mod); } friend MDCONST ModInt operator*(lint a, const ModInt &x) { return ModInt()._setval(a % mod * x.val % mod); } friend MDCONST ModInt operator/(lint a, const ModInt &x) { return ModInt()._setval(a % mod * x.inv() % mod); } MDCONST bool operator==(const ModInt &x) const { return val == x.val; } MDCONST bool operator!=(const ModInt &x) const { return val != x.val; } MDCONST bool operator<(const ModInt &x) const { return val < x.val; } // To use std::map<ModInt, T> friend std::istream &operator>>(std::istream &is, ModInt &x) { lint t; return is >> t, x = ModInt(t), is; } MDCONST friend std::ostream &operator<<(std::ostream &os, const ModInt &x) { return os << x.val; } MDCONST ModInt pow(lint n) const { lint ans = 1, tmp = this->val; while (n) { if (n & 1) ans = ans * tmp % mod; tmp = tmp * tmp % mod, n /= 2; } return ans; } static std::vector<long long> facs, invs; MDCONST static void _precalculation(int N) { int l0 = facs.size(); if (N <= l0) return; facs.resize(N), invs.resize(N); for (int i = l0; i < N; i++) facs[i] = facs[i - 1] * i % mod; long long facinv = ModInt(facs.back()).pow(mod - 2).val; for (int i = N - 1; i >= l0; i--) { invs[i] = facinv * facs[i - 1] % mod; facinv = facinv * i % mod; } } MDCONST lint inv() const { if (this->val < std::min(mod >> 1, 1 << 21)) { while (this->val >= int(facs.size())) _precalculation(facs.size() * 2); return invs[this->val]; } else { return this->pow(mod - 2).val; } } MDCONST ModInt fac() const { while (this->val >= int(facs.size())) _precalculation(facs.size() * 2); return facs[this->val]; } MDCONST ModInt doublefac() const { lint k = (this->val + 1) / 2; return (this->val & 1) ? ModInt(k * 2).fac() / (ModInt(2).pow(k) * ModInt(k).fac()) : ModInt(k).fac() * ModInt(2).pow(k); } MDCONST ModInt nCr(const ModInt &r) const { return (this->val < r.val) ? 0 : this->fac() / ((*this - r).fac() * r.fac()); } ModInt sqrt() const { if (val == 0) return 0; if (mod == 2) return val; if (pow((mod - 1) / 2) != 1) return 0; ModInt b = 1; while (b.pow((mod - 1) / 2) == 1) b += 1; int e = 0, m = mod - 1; while (m % 2 == 0) m >>= 1, e++; ModInt x = pow((m - 1) / 2), y = (*this) * x * x; x *= (*this); ModInt z = b.pow(m); while (y != 1) { int j = 0; ModInt t = y; while (t != 1) j++, t *= t; z = z.pow(1LL << (e - j - 1)); x *= z, z *= z, y *= z; e = j; } return ModInt(std::min(x.val, mod - x.val)); } }; template <int mod> std::vector<long long> ModInt<mod>::facs = {1}; template <int mod> std::vector<long long> ModInt<mod>::invs = {0}; // using mint = ModInt<998244353>; // using mint = ModInt<1000000007>; #line 2 "number/sieve.hpp" #include <cassert> #include <map> #line 5 "number/sieve.hpp" // CUT begin // Linear sieve algorithm for fast prime factorization // Complexity: O(N) time, O(N) space: // - MAXN = 10^7: ~44 MB, 80~100 ms (Codeforces / AtCoder GCC, C++17) // - MAXN = 10^8: ~435 MB, 810~980 ms (Codeforces / AtCoder GCC, C++17) // Reference: // [1] D. Gries, J. Misra, "A Linear Sieve Algorithm for Finding Prime Numbers," // Communications of the ACM, 21(12), 999-1003, 1978. // - https://cp-algorithms.com/algebra/prime-sieve-linear.html // - https://37zigen.com/linear-sieve/ struct Sieve { std::vector<int> min_factor; std::vector<int> primes; Sieve(int MAXN) : min_factor(MAXN + 1) { for (int d = 2; d <= MAXN; d++) { if (!min_factor[d]) { min_factor[d] = d; primes.emplace_back(d); } for (const auto &p : primes) { if (p > min_factor[d] or d * p > MAXN) break; min_factor[d * p] = p; } } } // Prime factorization for 1 <= x <= MAXN^2 // Complexity: O(log x) (x <= MAXN) // O(MAXN / log MAXN) (MAXN < x <= MAXN^2) template <typename T> std::map<T, int> factorize(T x) { std::map<T, int> ret; assert(x > 0 and x <= ((long long)min_factor.size() - 1) * ((long long)min_factor.size() - 1)); for (const auto &p : primes) { if (x < T(min_factor.size())) break; while (!(x % p)) x /= p, ret[p]++; } if (x >= T(min_factor.size())) ret[x]++, x = 1; while (x > 1) ret[min_factor[x]]++, x /= min_factor[x]; return ret; } // Enumerate divisors of 1 <= x <= MAXN^2 // Be careful of highly composite numbers https://oeis.org/A002182/list https://gist.github.com/dario2994/fb4713f252ca86c1254d#file-list-txt // (n, (# of div. of n)): 45360->100, 735134400(<1e9)->1344, 963761198400(<1e12)->6720 template <typename T> std::vector<T> divisors(T x) { std::vector<T> ret{1}; for (const auto p : factorize(x)) { int n = ret.size(); for (int i = 0; i < n; i++) { for (T a = 1, d = 1; d <= p.second; d++) { a *= p.first; ret.push_back(ret[i] * a); } } } return ret; // NOT sorted } // Moebius function Table, (-1)^{# of different prime factors} for square-free x // return: [0=>0, 1=>1, 2=>-1, 3=>-1, 4=>0, 5=>-1, 6=>1, 7=>-1, 8=>0, ...] https://oeis.org/A008683 std::vector<int> GenerateMoebiusFunctionTable() { std::vector<int> ret(min_factor.size()); for (unsigned i = 1; i < min_factor.size(); i++) { if (i == 1) ret[i] = 1; else if ((i / min_factor[i]) % min_factor[i] == 0) ret[i] = 0; else ret[i] = -ret[i / min_factor[i]]; } return ret; } // Calculate [0^K, 1^K, ..., nmax^K] in O(nmax) // Note: **0^0 == 1** template <typename MODINT> std::vector<MODINT> enumerate_kth_pows(long long K, int nmax) { assert(nmax < int(min_factor.size())); assert(K >= 0); if (K == 0) return std::vector<MODINT>(nmax + 1, 1); std::vector<MODINT> ret(nmax + 1); ret[0] = 0, ret[1] = 1; for (int n = 2; n <= nmax; n++) { if (min_factor[n] == n) { ret[n] = MODINT(n).pow(K); } else { ret[n] = ret[n / min_factor[n]] * ret[min_factor[n]]; } } return ret; } }; // Sieve sieve(1 << 15); // (can factorize n <= 10^9) #line 3 "convolution/ntt.hpp" #include <algorithm> #include <array> #line 7 "convolution/ntt.hpp" #include <tuple> #line 9 "convolution/ntt.hpp" // CUT begin // Integer convolution for arbitrary mod // with NTT (and Garner's algorithm) for ModInt / ModIntRuntime class. // We skip Garner's algorithm if `skip_garner` is true or mod is in `nttprimes`. // input: a (size: n), b (size: m) // return: vector (size: n + m - 1) template <typename MODINT> std::vector<MODINT> nttconv(std::vector<MODINT> a, std::vector<MODINT> b, bool skip_garner = false); constexpr int nttprimes[3] = {998244353, 167772161, 469762049}; // Integer FFT (Fast Fourier Transform) for ModInt class // (Also known as Number Theoretic Transform, NTT) // is_inverse: inverse transform // ** Input size must be 2^n ** template <typename MODINT> void ntt(std::vector<MODINT> &a, bool is_inverse = false) { int n = a.size(); if (n == 1) return; static const int mod = MODINT::get_mod(); static const MODINT root = MODINT::get_primitive_root(); assert(__builtin_popcount(n) == 1 and (mod - 1) % n == 0); static std::vector<MODINT> w{1}, iw{1}; for (int m = w.size(); m < n / 2; m *= 2) { MODINT dw = root.pow((mod - 1) / (4 * m)), dwinv = 1 / dw; w.resize(m * 2), iw.resize(m * 2); for (int i = 0; i < m; i++) w[m + i] = w[i] * dw, iw[m + i] = iw[i] * dwinv; } if (!is_inverse) { for (int m = n; m >>= 1;) { for (int s = 0, k = 0; s < n; s += 2 * m, k++) { for (int i = s; i < s + m; i++) { MODINT x = a[i], y = a[i + m] * w[k]; a[i] = x + y, a[i + m] = x - y; } } } } else { for (int m = 1; m < n; m *= 2) { for (int s = 0, k = 0; s < n; s += 2 * m, k++) { for (int i = s; i < s + m; i++) { MODINT x = a[i], y = a[i + m]; a[i] = x + y, a[i + m] = (x - y) * iw[k]; } } } int n_inv = MODINT(n).inv(); for (auto &v : a) v *= n_inv; } } template <int MOD> std::vector<ModInt<MOD>> nttconv_(const std::vector<int> &a, const std::vector<int> &b) { int sz = a.size(); assert(a.size() == b.size() and __builtin_popcount(sz) == 1); std::vector<ModInt<MOD>> ap(sz), bp(sz); for (int i = 0; i < sz; i++) ap[i] = a[i], bp[i] = b[i]; ntt(ap, false); if (a == b) bp = ap; else ntt(bp, false); for (int i = 0; i < sz; i++) ap[i] *= bp[i]; ntt(ap, true); return ap; } long long garner_ntt_(int r0, int r1, int r2, int mod) { using mint2 = ModInt<nttprimes[2]>; static const long long m01 = 1LL * nttprimes[0] * nttprimes[1]; static const long long m0_inv_m1 = ModInt<nttprimes[1]>(nttprimes[0]).inv(); static const long long m01_inv_m2 = mint2(m01).inv(); int v1 = (m0_inv_m1 * (r1 + nttprimes[1] - r0)) % nttprimes[1]; auto v2 = (mint2(r2) - r0 - mint2(nttprimes[0]) * v1) * m01_inv_m2; return (r0 + 1LL * nttprimes[0] * v1 + m01 % mod * v2.val) % mod; } template <typename MODINT> std::vector<MODINT> nttconv(std::vector<MODINT> a, std::vector<MODINT> b, bool skip_garner) { int sz = 1, n = a.size(), m = b.size(); while (sz < n + m) sz <<= 1; if (sz <= 16) { std::vector<MODINT> ret(n + m - 1); for (int i = 0; i < n; i++) { for (int j = 0; j < m; j++) ret[i + j] += a[i] * b[j]; } return ret; } int mod = MODINT::get_mod(); if (skip_garner or std::find(std::begin(nttprimes), std::end(nttprimes), mod) != std::end(nttprimes)) { a.resize(sz), b.resize(sz); if (a == b) { ntt(a, false); b = a; } else ntt(a, false), ntt(b, false); for (int i = 0; i < sz; i++) a[i] *= b[i]; ntt(a, true); a.resize(n + m - 1); } else { std::vector<int> ai(sz), bi(sz); for (int i = 0; i < n; i++) ai[i] = a[i].val; for (int i = 0; i < m; i++) bi[i] = b[i].val; auto ntt0 = nttconv_<nttprimes[0]>(ai, bi); auto ntt1 = nttconv_<nttprimes[1]>(ai, bi); auto ntt2 = nttconv_<nttprimes[2]>(ai, bi); a.resize(n + m - 1); for (int i = 0; i < n + m - 1; i++) { a[i] = garner_ntt_(ntt0[i].val, ntt1[i].val, ntt2[i].val, mod); } } return a; } #line 6 "formal_power_series/formal_power_series.hpp" using namespace std; // CUT begin // Formal Power Series (形式的冪級数) based on ModInt<mod> / ModIntRuntime // Reference: https://ei1333.github.io/luzhiled/snippets/math/formal-power-series.html template <typename T> struct FormalPowerSeries : vector<T> { using vector<T>::vector; using P = FormalPowerSeries; void shrink() { while (this->size() and this->back() == T(0)) this->pop_back(); } P operator+(const P &r) const { return P(*this) += r; } P operator+(const T &v) const { return P(*this) += v; } P operator-(const P &r) const { return P(*this) -= r; } P operator-(const T &v) const { return P(*this) -= v; } P operator*(const P &r) const { return P(*this) *= r; } P operator*(const T &v) const { return P(*this) *= v; } P operator/(const P &r) const { return P(*this) /= r; } P operator/(const T &v) const { return P(*this) /= v; } P operator%(const P &r) const { return P(*this) %= r; } P &operator+=(const P &r) { if (r.size() > this->size()) this->resize(r.size()); for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i]; shrink(); return *this; } P &operator+=(const T &v) { if (this->empty()) this->resize(1); (*this)[0] += v; shrink(); return *this; } P &operator-=(const P &r) { if (r.size() > this->size()) this->resize(r.size()); for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i]; shrink(); return *this; } P &operator-=(const T &v) { if (this->empty()) this->resize(1); (*this)[0] -= v; shrink(); return *this; } P &operator*=(const T &v) { for (auto &x : (*this)) x *= v; shrink(); return *this; } P &operator*=(const P &r) { if (this->empty() || r.empty()) this->clear(); else { auto ret = nttconv(*this, r); *this = P(ret.begin(), ret.end()); } return *this; } P &operator%=(const P &r) { *this -= *this / r * r; shrink(); return *this; } P operator-() const { P ret = *this; for (auto &v : ret) v = -v; return ret; } P &operator/=(const T &v) { assert(v != T(0)); for (auto &x : (*this)) x /= v; return *this; } P &operator/=(const P &r) { if (this->size() < r.size()) { this->clear(); return *this; } int n = (int)this->size() - r.size() + 1; return *this = (reversed().pre(n) * r.reversed().inv(n)).pre(n).reversed(n); } P pre(int sz) const { P ret(this->begin(), this->begin() + min((int)this->size(), sz)); ret.shrink(); return ret; } P operator>>(int sz) const { if ((int)this->size() <= sz) return {}; return P(this->begin() + sz, this->end()); } P operator<<(int sz) const { if (this->empty()) return {}; P ret(*this); ret.insert(ret.begin(), sz, T(0)); return ret; } P reversed(int deg = -1) const { assert(deg >= -1); P ret(*this); if (deg != -1) ret.resize(deg, T(0)); reverse(ret.begin(), ret.end()); ret.shrink(); return ret; } P differential() const { // formal derivative (differential) of f.p.s. const int n = (int)this->size(); P ret(max(0, n - 1)); for (int i = 1; i < n; i++) ret[i - 1] = (*this)[i] * T(i); return ret; } P integral() const { const int n = (int)this->size(); P ret(n + 1); ret[0] = T(0); for (int i = 0; i < n; i++) ret[i + 1] = (*this)[i] / T(i + 1); return ret; } P inv(int deg) const { assert(deg >= -1); assert(this->size() and ((*this)[0]) != T(0)); // Requirement: F(0) != 0 const int n = this->size(); if (deg == -1) deg = n; P ret({T(1) / (*this)[0]}); for (int i = 1; i < deg; i <<= 1) { ret = (ret + ret - ret * ret * pre(i << 1)).pre(i << 1); } ret = ret.pre(deg); ret.shrink(); return ret; } P log(int deg = -1) const { assert(deg >= -1); assert(this->size() and ((*this)[0]) == T(1)); // Requirement: F(0) = 1 const int n = (int)this->size(); if (deg == 0) return {}; if (deg == -1) deg = n; return (this->differential() * this->inv(deg)).pre(deg - 1).integral(); } P sqrt(int deg = -1) const { assert(deg >= -1); const int n = (int)this->size(); if (deg == -1) deg = n; if (this->empty()) return {}; if ((*this)[0] == T(0)) { for (int i = 1; i < n; i++) if ((*this)[i] != T(0)) { if ((i & 1) or deg - i / 2 <= 0) return {}; return (*this >> i).sqrt(deg - i / 2) << (i / 2); } return {}; } T sqrtf0 = (*this)[0].sqrt(); if (sqrtf0 == T(0)) return {}; P y = (*this) / (*this)[0], ret({T(1)}); T inv2 = T(1) / T(2); for (int i = 1; i < deg; i <<= 1) { ret = (ret + y.pre(i << 1) * ret.inv(i << 1)) * inv2; } return ret.pre(deg) * sqrtf0; } P exp(int deg = -1) const { assert(deg >= -1); assert(this->empty() or ((*this)[0]) == T(0)); // Requirement: F(0) = 0 const int n = (int)this->size(); if (deg == -1) deg = n; P ret({T(1)}); for (int i = 1; i < deg; i <<= 1) { ret = (ret * (pre(i << 1) + T(1) - ret.log(i << 1))).pre(i << 1); } return ret.pre(deg); } P pow(long long k, int deg = -1) const { assert(deg >= -1); const int n = (int)this->size(); if (deg == -1) deg = n; for (int i = 0; i < n; i++) { if ((*this)[i] != T(0)) { T rev = T(1) / (*this)[i]; P C = (*this) * rev, D(n - i); for (int j = i; j < n; j++) D[j - i] = C.coeff(j); D = (D.log(deg) * T(k)).exp(deg) * (*this)[i].pow(k); if (k * (i > 0) > deg or k * i > deg) return {}; P E(deg); long long S = i * k; for (int j = 0; j + S < deg and j < (int)D.size(); j++) E[j + S] = D[j]; E.shrink(); return E; } } return *this; } // Calculate f(X + c) from f(X), O(NlogN) P shift(T c) const { const int n = (int)this->size(); P ret = *this; for (int i = 0; i < n; i++) { ret[i] *= T(i).fac(); } reverse(ret.begin(), ret.end()); P exp_cx(n, 1); for (int i = 1; i < n; i++) { exp_cx[i] = exp_cx[i - 1] * c / i; } ret = (ret * exp_cx), ret.resize(n); reverse(ret.begin(), ret.end()); for (int i = 0; i < n; i++) { ret[i] /= T(i).fac(); } return ret; } T coeff(int i) const { if ((int)this->size() <= i or i < 0) return T(0); return (*this)[i]; } T eval(T x) const { T ret = 0, w = 1; for (auto &v : *this) ret += w * v, w *= x; return ret; } }; #line 5 "formal_power_series/multipoint_evaluation.hpp" // CUT begin // multipoint polynomial evaluation // input: xs = [x_0, ..., x_{N - 1}]: points to evaluate // f = \sum_i^M f_i x^i // Complexity: O(N (lgN)^2) building, O(N (lgN)^2 + M lg M) evaluation template <typename Tfield> struct MultipointEvaluation { int nx; using polynomial = FormalPowerSeries<Tfield>; std::vector<polynomial> segtree; MultipointEvaluation(const std::vector<Tfield> &xs) : nx(xs.size()) { segtree.resize(nx * 2 - 1); for (int i = 0; i < nx; i++) { segtree[nx - 1 + i] = {-xs[i], 1}; } for (int i = nx - 2; i >= 0; i--) { segtree[i] = segtree[2 * i + 1] * segtree[2 * i + 2]; } } std::vector<Tfield> ret; void _eval_rec(polynomial f, int now) { f %= segtree[now]; if (now - (nx - 1) >= 0) { ret[now - (nx - 1)] = f.coeff(0); return; } _eval_rec(f, 2 * now + 1); _eval_rec(f, 2 * now + 2); } std::vector<Tfield> evaluate_polynomial(const polynomial &f) { ret.resize(nx); _eval_rec(f, 0); return ret; } std::vector<Tfield> evaluate_polynomial(const std::vector<Tfield> &f) { return evaluate_polynomial(polynomial(f.begin(), f.end())); } std::vector<Tfield> _interpolate_coeffs; polynomial _rec_interpolation(int now, const std::vector<Tfield> &ys) const { int i = now - (nx - 1); if (i >= 0) return {ys[i]}; auto retl = _rec_interpolation(2 * now + 1, ys); auto retr = _rec_interpolation(2 * now + 2, ys); return retl * segtree[2 * now + 2] + retr * segtree[2 * now + 1]; } std::vector<Tfield> polynomial_interpolation(std::vector<Tfield> ys) { assert(nx == int(ys.size())); if (_interpolate_coeffs.empty()) { _interpolate_coeffs = evaluate_polynomial(segtree[0].differential()); for (auto &x : _interpolate_coeffs) x = x.inv(); } for (int i = 0; i < nx; i++) ys[i] *= _interpolate_coeffs[i]; return _rec_interpolation(0, ys); } }; #line 4 "formal_power_series/sum_of_exponential_times_polynomial_limit.hpp" // CUT begin // $d$ 次以下の多項式 $f(x)$ と定数 $r$ について, // $\sum_{i=0}^\infty r^i f(i)$ の値を $[f(0), ..., f(d - 1), f(d)]$ の値から $O(d)$ で計算. // Requirement: r != 1 // https://judge.yosupo.jp/problem/sum_of_exponential_times_polynomial_limit // Document: https://hitonanode.github.io/cplib-cpp/formal_power_series/sum_of_exponential_times_polynomial_limit.hpp template <typename MODINT> MODINT sum_of_exponential_times_polynomial_limit(MODINT r, std::vector<MODINT> init) { assert(r != 1); if (init.empty()) return 0; if (init.size() == 1) return init[0] / (1 - r); auto &bs = init; const int d = int(bs.size()) - 1; MODINT rp = 1; for (int i = 1; i <= d; i++) rp *= r, bs[i] = bs[i] * rp + bs[i - 1]; MODINT ret = 0; rp = 1; for (int i = 0; i <= d; i++) { ret += bs[d - i] * MODINT(d + 1).nCr(i) * rp; rp *= -r; } return ret / MODINT(1 - r).pow(d + 1); }; #line 8 "formal_power_series/test/sum_of_exponential_times_polynomial_limit.test.cpp" using namespace std; using mint = ModInt<998244353>; int main() { int r, d; cin >> r >> d; mint::_precalculation(d + 10); auto initial_terms = Sieve(d).enumerate_kth_pows<mint>(d, d); if (d > 0 and d <= 100000) { vector<mint> xs(d + 1); for (int i = 0; i <= d; i++) xs[i] = i; vector<mint> f(d + 1); f[d] = 1; assert(MultipointEvaluation<mint>(xs).evaluate_polynomial(f) == initial_terms); } cout << sum_of_exponential_times_polynomial_limit<mint>(r, initial_terms) << '\n'; }