cp-library-cpp
This documentation is automatically generated by online-judge-tools/verification-helper
test/src/sequence/eulerian_number/yuki2005.test.cpp
- View this file on GitHub
- Last update: 2026-06-28 22:26:07+09:00
- Problem: https://yukicoder.me/problems/no/2005
Depends on
- SWAG を Deque に拡張したやつ
(library/datastructure/deque_aggregation.hpp)
- 階乗テーブル
(library/math/factorial.hpp)
- 線形篩
(library/number/linear_sieve.hpp)
- Eulerian Number
(library/sequence/eulerian_number.hpp)
- Powers
(library/sequence/powers.hpp)
Code
#define PROBLEM "https://yukicoder.me/problems/no/2005"
#include <iostream>
#include <atcoder/modint>
using mint = atcoder::modint998244353;
std::istream& operator>>(std::istream& in, mint &a) {
long long e; in >> e; a = e;
return in;
}
std::ostream& operator<<(std::ostream& out, const mint &a) {
out << a.val();
return out;
}
#include "library/sequence/eulerian_number.hpp"
#include "library/datastructure/deque_aggregation.hpp"
mint op(mint x, mint y) {
return x * y;
}
mint e() {
return 1;
}
constexpr uint32_t K_MAX = 5000;
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
uint32_t n;
uint64_t m;
std::cin >> n >> m;
std::vector<mint> c(K_MAX + 1);
for (uint32_t i = 0; i < n; ++i) {
uint32_t k;
std::cin >> k;
++c[k];
}
suisen::factorial<mint> fac(n + K_MAX);
mint ans = 0;
auto en = suisen::eulerian_number_table<mint>(K_MAX);
suisen::DequeAggregation<mint, op, e> dq;
for (uint32_t d = 0; d < n; ++d) dq.push_front(m + d);
for (uint32_t k = 1; k <= K_MAX; ++k) {
dq.push_front(m + n + k - 1);
mint sum = 0;
const uint32_t p = std::min(uint64_t(k), m);
for (uint32_t i = 0; i < p; ++i) {
sum += en[k][i] * dq.prod();
dq.pop_front();
dq.push_back(m - i - 1);
}
ans += c[k] * sum * fac.fac_inv(n + k);
for (uint32_t i = p; i --> 0;) {
dq.push_front(m - i + n + k - 1);
dq.pop_back();
}
}
std::cout << ans.val() << std::endl;
return 0;
}#line 1 "test/src/sequence/eulerian_number/yuki2005.test.cpp"
#define PROBLEM "https://yukicoder.me/problems/no/2005"
#include <iostream>
#include <atcoder/modint>
using mint = atcoder::modint998244353;
std::istream& operator>>(std::istream& in, mint &a) {
long long e; in >> e; a = e;
return in;
}
std::ostream& operator<<(std::ostream& out, const mint &a) {
out << a.val();
return out;
}
#line 1 "library/sequence/eulerian_number.hpp"
#line 1 "library/math/factorial.hpp"
#include <cassert>
#include <vector>
namespace suisen {
// 引数として与える値に対して、法が十分大きいことを仮定する
template <typename T, typename U = T>
struct factorial {
factorial() = default;
factorial(int n) { ensure(n); }
static void ensure(const int n) {
int sz = _fac.size();
if (n + 1 <= sz) return;
int new_size = std::max(n + 1, sz * 2);
_fac.resize(new_size), _fac_inv.resize(new_size);
for (int i = sz; i < new_size; ++i) _fac[i] = _fac[i - 1] * i;
_fac_inv[new_size - 1] = U(1) / _fac[new_size - 1];
for (int i = new_size - 1; i > sz; --i) _fac_inv[i - 1] = _fac_inv[i] * i;
}
T fac(const int i) {
ensure(i);
return _fac[i];
}
T operator()(int i) {
return fac(i);
}
U fac_inv(const int i) {
ensure(i);
return _fac_inv[i];
}
// i の逆数
// i = 0 の場合は assert 違反となる
U inv(const int i) {
assert(i > 0);
ensure(i);
return _fac_inv[i] * _fac[i - 1];
}
U binom(const int n, const int r) {
if (n < 0 or r < 0 or n < r) return 0;
ensure(n);
return _fac[n] * _fac_inv[r] * _fac_inv[n - r];
}
// binom(n, r) の逆数
// binom(n, r) = 0 の場合は assert 違反となる
U binom_inv(const int n, const int r) {
assert(r >= 0 and n >= r);
ensure(n);
return _fac_inv[n] * _fac[r] * _fac[n - r];
}
// n 種類から重複を許して r 個選ぶ場合の数
// x_1+x_2+...+x_n=r（x_i は非負整数）となる x の個数でもある
// multichoose(n, r) = binom(n + r - 1, r)
U multichoose(const int n, const int r) {
if (n < 0 or r < 0) return 0;
return r > 0 ? binom(n + r - 1, r) : U(1);
}
// n 種類から重複を許して r 個選ぶ場合の数 multichoose(n, r) の逆数
// x_1+x_2+...+x_n=r（x_i は非負整数）となる x の個数の逆数でもある
// multichoose(n, r) = binom(n + r - 1, r)
// multichoose(n, r) = 0 の場合は assert 違反となる
U multichoose_inv(const int n, const int r) {
assert(n >= 0 and r >= 0);
return r > 0 ? binom_inv(n + r - 1, r) : U(1);
}
template <typename ...Ds, std::enable_if_t<std::conjunction_v<std::is_integral<Ds>...>, std::nullptr_t> = nullptr>
U polynom(const int n, const Ds& ...ds) {
if (n < 0) return 0;
ensure(n);
int sumd = 0;
U res = _fac[n];
for (int d : { ds... }) {
if (d < 0 or d > n) return 0;
sumd += d;
res *= _fac_inv[d];
}
if (sumd > n) return 0;
res *= _fac_inv[n - sumd];
return res;
}
U perm(const int n, const int r) {
if (n < 0 or r < 0 or n < r) return 0;
ensure(n);
return _fac[n] * _fac_inv[n - r];
}
// perm(n, r) の逆数
// perm(n, r) = 0 の場合は assert 違反となる
U perm_inv(const int n, const int r) {
assert(r >= 0 and n >= r);
ensure(n);
return _fac_inv[n] * _fac[n - r];
}
private:
static std::vector<T> _fac;
static std::vector<U> _fac_inv;
};
template <typename T, typename U>
std::vector<T> factorial<T, U>::_fac{ 1 };
template <typename T, typename U>
std::vector<U> factorial<T, U>::_fac_inv{ 1 };
} // namespace suisen
#line 1 "library/sequence/powers.hpp"
#include <cstdint>
#line 1 "library/number/linear_sieve.hpp"
#line 5 "library/number/linear_sieve.hpp"
#include <numeric>
#line 7 "library/number/linear_sieve.hpp"
namespace suisen {
// reference: https://37zigen.com/linear-sieve/
class LinearSieve {
public:
LinearSieve(const int n) : _n(n), min_prime_factor(std::vector<int>(n + 1)) {
std::iota(min_prime_factor.begin(), min_prime_factor.end(), 0);
prime_list.reserve(_n / 20);
for (int d = 2; d <= _n; ++d) {
if (min_prime_factor[d] == d) prime_list.push_back(d);
const int prime_max = std::min(min_prime_factor[d], _n / d);
for (int prime : prime_list) {
if (prime > prime_max) break;
min_prime_factor[prime * d] = prime;
}
}
}
int prime_num() const noexcept { return prime_list.size(); }
/**
* Returns a vector of primes in [0, n].
* It is guaranteed that the returned vector is sorted in ascending order.
*/
const std::vector<int>& get_prime_list() const noexcept {
return prime_list;
}
const std::vector<int>& get_min_prime_factor() const noexcept { return min_prime_factor; }
/**
* Returns a vector of `{ prime, index }`.
* It is guaranteed that the returned vector is sorted in ascending order.
*/
std::vector<std::pair<int, int>> factorize(int n) const noexcept {
assert(0 < n and n <= _n);
std::vector<std::pair<int, int>> prime_powers;
while (n > 1) {
int p = min_prime_factor[n], c = 0;
do { n /= p, ++c; } while (n % p == 0);
prime_powers.emplace_back(p, c);
}
return prime_powers;
}
private:
const int _n;
std::vector<int> min_prime_factor;
std::vector<int> prime_list;
};
} // namespace suisen
#line 6 "library/sequence/powers.hpp"
namespace suisen {
// returns { 0^k, 1^k, ..., n^k }
template <typename mint>
std::vector<mint> powers(uint32_t n, uint64_t k) {
const auto mpf = LinearSieve(n).get_min_prime_factor();
std::vector<mint> res(n + 1);
res[0] = k == 0;
for (uint32_t i = 1; i <= n; ++i) res[i] = i == 1 ? 1 : uint32_t(mpf[i]) == i ? mint(i).pow(k) : res[mpf[i]] * res[i / mpf[i]];
return res;
}
} // namespace suisen
#line 6 "library/sequence/eulerian_number.hpp"
// reference: https://en.wikipedia.org/wiki/Eulerian_number
namespace suisen {
template <typename FPSType>
std::vector<typename FPSType::value_type> eulerian_number(uint32_t n) {
using mint = typename FPSType::value_type;
if (n == 0) return {};
factorial<mint> fac(n + 1);
const uint32_t h = (n + 1) >> 1;
FPSType f = powers<mint>(h, n);
f.erase(f.begin());
FPSType g(h);
for (uint32_t i = 0; i < h; ++i) {
mint v = fac.binom(n + 1, i);
g[i] = i & 1 ? -v : v;
}
FPSType res = f * g;
res.resize(n);
for (uint32_t i = h; i < n; ++i) res[i] = res[n - 1 - i];
return res;
}
template <typename mint>
std::vector<std::vector<mint>> eulerian_number_table(uint32_t n) {
if (n == 0) return {};
std::vector dp(n + 1, std::vector<mint>{});
for (uint32_t i = 1; i <= n; ++i) {
dp[i].resize(i);
dp[i][0] = dp[i][i - 1] = 1;
for (uint32_t j = 1; j < i - 1; ++j) dp[i][j] = (i - j) * dp[i - 1][j - 1] + (j + 1) * dp[i - 1][j];
}
return dp;
}
} // namespace suisen
#line 1 "library/datastructure/deque_aggregation.hpp"
#line 6 "library/datastructure/deque_aggregation.hpp"
/**
* [Idea] reference : https://motsu-xe.hatenablog.com/entry/2021/05/13/224016
*
* SWAG + simulate a deque with 2 stacks
*
* [Operations] reference : https://www.slideshare.net/catupper/amortize-analysis-of-deque-with-2-stack
*
* `l`, `r` is a stack of { value, sum }
*
* accumulate
* <---------- ------> fold values from inside
* ( l ][ r )
*
* pop_front:
* 1. `l` is not empty
* ( l ][ r ) -> ( l ][ r ) # pop from `l`. O(1)
* 2. `l` is empty
* (][ r ) -> ( l ][ r ) # split `r` at its middle point. amortized O(1)
* ( l ][ r ) -> ( l ][ r ) # pop from `l`. O(1)
*
* pop_back:
* 1. `r` is not empty
* ( l ][ r ) -> ( l ][ r ) # pop from `r`. O(1)
* 2. `r` is empty
* ( l ][) -> ( l ][ r ) # split `l` at its middle point. amortized O(1)
* ( l ][ r ) -> ( l ][ r ) # pop from `r`. O(1)
*
* push_front:
* ( l ][ r ) -> ( l ][ r ) # push to `l`. O(1)
*
* push_back:
* ( l ][ r ) -> ( l ][ r ) # push to `r`. O(1)
*/
namespace suisen {
template <typename T, T(*op)(T, T), T(*e)()>
struct DequeAggregation {
struct DequeAggregationIterator {
using difference_type = int;
using value_type = T;
using pointer = value_type*;
using reference = value_type&;
using iterator_category = std::random_access_iterator_tag;
using fi_iterator_type = typename std::vector<std::pair<value_type, value_type>>::const_reverse_iterator;
using se_iterator_type = typename std::vector<std::pair<value_type, value_type>>::const_iterator;
fi_iterator_type it_l;
fi_iterator_type it_l_end;
se_iterator_type it_r_begin;
se_iterator_type it_r;
DequeAggregationIterator& operator++() {
if (it_l == it_l_end) ++it_r;
else ++it_l;
return *this;
}
DequeAggregationIterator operator++(int) { DequeAggregationIterator ret = *this; ++(*this); return ret; }
DequeAggregationIterator& operator--() {
if (it_r == it_r_begin) --it_l;
else --it_r;
return *this;
}
DequeAggregationIterator operator--(int) { DequeAggregationIterator ret = *this; --(*this); return ret; }
DequeAggregationIterator& operator+=(difference_type dif) {
if (dif < 0) return *this -= -dif;
if (int d = it_l_end - it_l; d < dif) it_l = it_l_end, it_r += dif - d;
else it_l += dif;
return *this;
}
friend DequeAggregationIterator operator+(DequeAggregationIterator it, difference_type dif) { it += dif; return it; }
friend DequeAggregationIterator operator+(difference_type dif, DequeAggregationIterator it) { it += dif; return it; }
DequeAggregationIterator& operator-=(difference_type dif) {
if (dif < 0) return *this += -dif;
if (int d = it_r - it_r_begin; d < dif) it_r = it_r_begin, it_l -= dif - d;
else it_r -= dif;
return *this;
}
friend DequeAggregationIterator operator-(DequeAggregationIterator it, difference_type dif) { it -= dif; return it; }
difference_type operator-(const DequeAggregationIterator &rhs) const {
difference_type d1 = it_l == it_l_end ? it_r - it_r_begin : it_l - it_l_end;
difference_type d2 = rhs.it_l == rhs.it_l_end ? rhs.it_r - rhs.it_r_begin : rhs.it_l - rhs.it_l_end;
return d1 - d2;
}
const value_type& operator[](difference_type i) const { return *((*this) + i); }
const value_type& operator*() const { return it_l == it_l_end ? it_r->first : it_l->first; }
bool operator!=(const DequeAggregationIterator &rhs) const { return it_l != rhs.it_l or it_r != rhs.it_r; }
bool operator==(const DequeAggregationIterator &rhs) const { return not (*this != rhs); }
bool operator< (const DequeAggregationIterator &rhs) const { return (*this) - rhs < 0; }
bool operator<=(const DequeAggregationIterator &rhs) const { return (*this) - rhs <= 0; }
bool operator> (const DequeAggregationIterator &rhs) const { return (*this) - rhs > 0; }
bool operator>=(const DequeAggregationIterator &rhs) const { return (*this) - rhs >= 0; }
};
using iterator = DequeAggregationIterator;
using difference_type = typename iterator::difference_type;
using value_type = typename iterator::value_type;
using pointer = typename iterator::pointer;
using reference = typename iterator::reference;
DequeAggregation() = default;
template <typename InputIterator, std::enable_if_t<std::is_constructible_v<value_type, typename InputIterator::value_type>, std::nullptr_t> = nullptr>
DequeAggregation(InputIterator first, InputIterator last) {
for (; first != last; ++first) push_back(*first);
}
template <typename Container, std::enable_if_t<std::is_constructible_v<value_type, typename Container::value_type>, std::nullptr_t> = nullptr>
DequeAggregation(const Container &c) : DequeAggregation(std::begin(c), std::end(c)) {}
value_type prod() const {
return op(prod(_st_l), prod(_st_r));
}
void push_back(const value_type &val) { _st_r.emplace_back(val, op(prod(_st_r), val)); }
void push_front(const value_type &val) { _st_l.emplace_back(val, op(val, prod(_st_l))); }
void pop_back() {
if (_st_r.size()) return _st_r.pop_back();
const int size = _st_l.size();
const int l = size >> 1, r = size - l;
assert(r); // <=> size > 0
for (int i = r - 1; i > 0; --i) push_back(std::move(_st_l[i].first));
_st_l.erase(_st_l.begin(), _st_l.begin() + r);
if (l == 0) return;
_st_l[0].second = _st_l[0].first;
for (int i = 1; i < l; ++i) _st_l[i].second = op(_st_l[i].first, _st_l[i - 1].second);
}
void pop_front() {
if (_st_l.size()) return _st_l.pop_back();
const int size = _st_r.size();
const int r = size >> 1, l = size - r;
assert(l); // <=> size > 0
for (int i = l - 1; i > 0; --i) push_front(std::move(_st_r[i].first));
_st_r.erase(_st_r.begin(), _st_r.begin() + l);
if (r == 0) return;
_st_r[0].second = _st_r[0].first;
for (int i = 1; i < r; ++i) _st_r[i].second = op(_st_r[i - 1].second, _st_r[i].first);
}
const value_type& front() const { return _st_l.size() ? _st_l.back().first : _st_r.front().first; }
const value_type& back() const { return _st_r.size() ? _st_r.back().first : _st_l.front().first; }
const value_type& operator[](int i) const {
const int k = i - _st_l.size();
return k < 0 ? _st_l[~k].first : _st_r[k].first;
}
int size() const { return _st_l.size() + _st_r.size(); }
void clear() { _st_l.clear(), _st_r.clear(); }
void shrink_to_fit() { _st_l.shrink_to_fit(), _st_r.shrink_to_fit(); }
iterator begin() const { return iterator { _st_l.rbegin(), _st_l.rend(), _st_r.begin(), _st_r.begin() }; }
iterator end() const { return iterator { _st_l.rend(), _st_l.rend(), _st_r.begin(), _st_r.end() }; }
iterator cbegin() const { return begin(); }
iterator cend() const { return end(); }
private:
std::vector<std::pair<value_type, value_type>> _st_l, _st_r;
value_type prod(const std::vector<std::pair<value_type, value_type>> &st) const {
return st.empty() ? e() : st.back().second;
}
};
} // namespace suisen
#line 20 "test/src/sequence/eulerian_number/yuki2005.test.cpp"
mint op(mint x, mint y) {
return x * y;
}
mint e() {
return 1;
}
constexpr uint32_t K_MAX = 5000;
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
uint32_t n;
uint64_t m;
std::cin >> n >> m;
std::vector<mint> c(K_MAX + 1);
for (uint32_t i = 0; i < n; ++i) {
uint32_t k;
std::cin >> k;
++c[k];
}
suisen::factorial<mint> fac(n + K_MAX);
mint ans = 0;
auto en = suisen::eulerian_number_table<mint>(K_MAX);
suisen::DequeAggregation<mint, op, e> dq;
for (uint32_t d = 0; d < n; ++d) dq.push_front(m + d);
for (uint32_t k = 1; k <= K_MAX; ++k) {
dq.push_front(m + n + k - 1);
mint sum = 0;
const uint32_t p = std::min(uint64_t(k), m);
for (uint32_t i = 0; i < p; ++i) {
sum += en[k][i] * dq.prod();
dq.pop_front();
dq.push_back(m - i - 1);
}
ans += c[k] * sum * fac.fac_inv(n + k);
for (uint32_t i = p; i --> 0;) {
dq.push_front(m - i + n + k - 1);
dq.pop_back();
}
}
std::cout << ans.val() << std::endl;
return 0;
}