結果
| 問題 |
No.802 だいたい等差数列
|
| ユーザー |
 miscalc
|
| 提出日時 | 2025-04-11 09:45:51 |
| 言語 | C++23 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 1,541 ms / 2,000 ms |
| コード長 | 37,240 bytes |
| コンパイル時間 | 6,417 ms |
| コンパイル使用メモリ | 337,248 KB |
| 実行使用メモリ | 70,464 KB |
| 最終ジャッジ日時 | 2025-04-11 09:46:32 |
| 合計ジャッジ時間 | 36,228 ms |
|
ジャッジサーバーID （参考情報） |
judge2 / judge4 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 4 |
| other | AC * 30 |
ソースコード
#include <bits/stdc++.h>
using namespace std;
namespace {
//*
#define INCLUDE_MODINT
//*/
#include <atcoder/modint>
using namespace atcoder;
//*
#define FAST_IO
//*/
#define SINGLE_TESTCASE
/*
using mint = modint998244353;
const bool ntt = true;
//*/
//*
using mint = modint1000000007;
const bool ntt = false;
//*/
/*
using mint = modint;
const bool ntt = false;
//*/
namespace mytemplate
{
using ll = long long;
template <class T> using vc = vector<T>;
template <class T> using vvc = vector<vector<T>>;
template <class T> using vvvc = vector<vector<vector<T>>>;
template <class T> using pql = priority_queue<T, vc<T>, greater<T>>;
template <class T> using pqg = priority_queue<T>;
#ifdef __SIZEOF_INT128__
using i128 = __int128_t;
i128 stoi128(string s) { i128 res = 0; if (s.front() == '-') { for (int i = 1; i < (int)s.size(); i++) res = 10 * res + s[i] - '0'; res = -res; } else { for (auto c : s) res = 10 * res + c - '0'; } return res; }
string i128tos(i128 x) { string sign = "", res = ""; if (x < 0) x = -x, sign = "-"; while (x > 0) { res += '0' + x % 10; x /= 10; } reverse(res.begin(), res.end()); if (res == "") return "0"; return sign + res; }
istream &operator>>(istream &is, i128 &a) { string s; is >> s; a = stoi128(s); return is; }
ostream &operator<<(ostream &os, const i128 &a) { os << i128tos(a); return os; }
#endif
#define cauto const auto
#define overload4(_1, _2, _3, _4, name, ...) name
#define rep1(i, n) for (ll i = 0, nnnnn = ll(n); i < nnnnn; i++)
#define rep2(i, l, r) for (ll i = ll(l), rrrrr = ll(r); i < rrrrr; i++)
#define rep3(i, l, r, d) for (ll i = ll(l), rrrrr = ll(r), ddddd = ll(d); ddddd > 0 ? i < rrrrr : i > rrrrr; i += d)
#define rep(...) overload4(__VA_ARGS__, rep3, rep2, rep1)(__VA_ARGS__)
#define repi1(i, n) for (int i = 0, nnnnn = int(n); i < nnnnn; i++)
#define repi2(i, l, r) for (int i = int(l), rrrrr = int(r); i < rrrrr; i++)
#define repi3(i, l, r, d) for (int i = int(l), rrrrr = int(r), ddddd = int(d); ddddd > 0 ? i < rrrrr : i > rrrrr; i += d)
#define repi(...) overload4(__VA_ARGS__, repi3, repi2, repi1)(__VA_ARGS__)
#define fe(...) for (auto __VA_ARGS__)
#define fec(...) for (cauto &__VA_ARGS__)
#define fem(...) for (auto &__VA_ARGS__)
#define ALL(a) (a).begin(), (a).end()
#define SZ(a) (ll)((a).size())
#define SZI(a) (int)((a).size())
const ll INF = 4'000'000'000'000'000'037;
template <class T = ll>
inline T divfloor(cauto &a, cauto &b) { return T(a) / T(b) - (T(a) % T(b) && (T(a) ^ T(b)) < 0); }
template <class T = ll>
inline T divround(cauto &a, cauto &b) { return divfloor<T>(2 * a + b, 2 * b); }
template <class T1 = ll> T1 mul_limited(auto a, auto b, T1 m = INF) { return b == 0 ? 0 : a > m / b ? m : a * b; }
template <class T = ll> struct max_op { T operator()(const T &a, const T &b) const { return max(a, b); } };
template <class T = ll> struct min_op { T operator()(const T &a, const T &b) const { return min(a, b); } };
template <class T, const T val> struct const_fn { T operator()() const { return val; } };
template <class T, size_t d, size_t i = 0>
auto dvec(cauto (&sz)[d], const T &init)
{
if constexpr (i < d)
return vc(sz[i], dvec<T, d, i + 1>(sz, init));
else
return init;
}
template <class T> void unique(vector<T> &v) { v.erase(unique(v.begin(), v.end()), v.end()); }
template <class T> void rotate(vector<T> &v, int k) { rotate(v.begin(), v.begin() + k, v.end()); }
void unique(string &s) { s.erase(unique(s.begin(), s.end()), s.end()); }
void rotate(string &s, int k) { rotate(s.begin(), s.begin() + k, s.end()); }
#if __cplusplus < 202002L
#define LB(v, ...) (ll)(lower_bound(ALL(v), __VA_ARGS__) - begin(v))
#define UB(v, ...) (ll)(upper_bound(ALL(v), __VA_ARGS__) - begin(v))
#define LBI(v, ...) (int)(lower_bound(ALL(v), __VA_ARGS__) - begin(v))
#define UBI(v, ...) (int)(upper_bound(ALL(v), __VA_ARGS__) - begin(v))
#else
#define LB(v, ...) (ll)(ranges::lower_bound((v), __VA_ARGS__) - begin(v))
#define UB(v, ...) (ll)(ranges::upper_bound((v), __VA_ARGS__) - begin(v))
#define LBI(v, ...) (int)(ranges::lower_bound((v), __VA_ARGS__) - begin(v))
#define UBI(v, ...) (int)(ranges::upper_bound((v), __VA_ARGS__) - begin(v))
#endif
#if __cplusplus < 202002L
#endif
template <class T>
struct bsubsets
{
private:
T x;
public:
bsubsets(T x) : x(x) {}
struct Iterator
{
private:
T y;
bool is_end;
const bsubsets &bs;
public:
Iterator(T y, bool is_end, const bsubsets &bs) : y(y), is_end(is_end), bs(bs) {}
T operator*() const { return y; }
Iterator& operator++()
{
if (y == 0)
is_end = true;
y = (y - 1) & bs.x;
return *this;
}
bool operator!=(const Iterator &other) const { return y != other.y || is_end != other.is_end; }
};
Iterator begin() const { return Iterator(x, false, *this); }
Iterator end() const { return Iterator(x, true, *this); }
};
template <class T>
struct bsupsets
{
private:
int n;
T x;
public:
bsupsets(int n, T x) : n(n), x(x) {}
struct Iterator
{
private:
T y;
const bsupsets &bs;
public:
Iterator(T y, const bsupsets &bs) : y(y), bs(bs) {}
T operator*() const { return y; }
Iterator& operator++()
{
y = (y + 1) | bs.x;
return *this;
}
bool operator!=(const Iterator &other) const { return y != other.y; }
};
Iterator begin() const { return Iterator(x, *this); }
Iterator end() const { return Iterator((T(1) << n) | x, *this); }
};
template <class Tuple, size_t... I> Tuple tuple_add(Tuple &a, const Tuple &b, const index_sequence<I...>) { ((get<I>(a) += get<I>(b)), ...); return a; }
template <class Tuple> Tuple operator+=(Tuple &a, const Tuple &b) { return tuple_add(a, b, make_index_sequence<tuple_size_v<Tuple>>{}); }
template <class Tuple> Tuple operator+(Tuple a, const Tuple &b) { return a += b; }
using namespace atcoder;
template <class T, internal::is_modint_t<T> * = nullptr> istream &operator>>(istream &is, T &a) { ll v; is >> v; a = v; return is; }
template <class T, internal::is_modint_t<T> * = nullptr> ostream &operator<<(ostream &os, const T &a) { os << a.val(); return os; }
#define MINT(...) mint __VA_ARGS__; INPUT(__VA_ARGS__)
template <class Tuple, enable_if_t<__is_tuple_like<Tuple>::value == true> * = nullptr> istream &operator>>(istream &is, Tuple &t) { apply([&](auto&... a){ (is >> ... >> a); }, t); return is; }
template <class... T> void INPUT(T&... a) { (cin >> ... >> a); }
template <class T> void INPUTVEC(int n, vector<T> &v) { v.resize(n); rep(i, n) cin >> v[i]; }
template <class T, class... Ts> void INPUTVEC(int n, vector<T>& v, vector<Ts>&... vs) { INPUTVEC(n, v); INPUTVEC(n, vs...); }
template <class T> void INPUTVEC2(int n, int m, vector<vector<T>> &v) { v.assign(n, vector<T>(m)); rep(i, n) rep(j, m) cin >> v[i][j]; }
template <class T, class... Ts> void INPUTVEC2(int n, int m, vector<T>& v, vector<Ts>&... vs) { INPUTVEC2(n, m, v); INPUTVEC2(n, m, vs...); }
#define CHAR(...) char __VA_ARGS__; INPUT(__VA_ARGS__)
#define INT(...) int __VA_ARGS__; INPUT(__VA_ARGS__)
#define LL(...) ll __VA_ARGS__; INPUT(__VA_ARGS__)
#define STR(...) string __VA_ARGS__; INPUT(__VA_ARGS__)
#define ARR(T, n, ...) array<T, n> __VA_ARGS__; INPUT(__VA_ARGS__)
#define VEC(T, n, ...) vector<T> __VA_ARGS__; INPUTVEC(n, __VA_ARGS__)
#define VEC2(T, n, m, ...) vector<vector<T>> __VA_ARGS__; INPUTVEC2(n, m, __VA_ARGS__)
#define ENDL endl
template <class T> void PRINT(const T &a) { cout << a << ENDL; }
template <class T, class... Ts> void PRINT(const T& a, const Ts&... b) { cout << a; (cout << ... << (cout << ' ', b)); cout << ENDL; }
#define PRINTEXIT(...) do { PRINT(__VA_ARGS__); exit(0); } while (false)
#define PRINTRETURN(...) do { PRINT(__VA_ARGS__); return; } while (false)
}
using namespace mytemplate;
#define dump(...)
//*
#include <atcoder/modint>
#include <atcoder/math>
#include <atcoder/convolution>
#include <atcoder/internal_math>
using namespace atcoder;
//*/
template<class T>
vector<T> convolution_anymod(const vector<T> &A, const vector<T> &B)
{
int N = A.size(), M = B.size();
if (min(N, M) <= 250)
{
vector<T> C(N + M - 1, 0);
for (int i = 0; i < N; i++)
for (int j = 0; j < M; j++)
C[i + j] += A[i] * B[j];
return C;
}
constexpr ll MOD1 = 167772161, MOD2 = 469762049, MOD3 = 1224736769;
using mint2 = static_modint<MOD2>;
using mint3 = static_modint<MOD3>;
constexpr int i1_2 = internal::inv_gcd(MOD1, MOD2).second;
constexpr int i12_3 = internal::inv_gcd(MOD1 * MOD2, MOD3).second;
T m12 = T(MOD1) * T(MOD2);
vector<int> A_(N), B_(M);
for (int i = 0; i < N; i++)
A_[i] = A[i].val();
for (int i = 0; i < M; i++)
B_[i] = B[i].val();
auto C1 = convolution<MOD1>(A_, B_);
auto C2 = convolution<MOD2>(A_, B_);
auto C3 = convolution<MOD3>(A_, B_);
vector<T> C(N + M - 1);
for (ll i = 0; i < N + M - 1; i++)
{
int c1 = C1[i], c2 = C2[i], c3 = C3[i];
int t1 = (mint2(c2 - c1) * mint2::raw(i1_2)).val();
int t2 = ((mint3(c3 - c1) - mint3::raw(t1) * mint3::raw(MOD1)) * mint3::raw(i12_3)).val();
C[i] = T(c1) + T(t1) * T(MOD1) + T(t2) * m12;
}
return C;
}
template<class T1>
struct LagrangeInterpolation
{
int D;
vector<T1> Y, fac, finv, prodl, prodr;
template<class T2>
LagrangeInterpolation(const vector<T2> &y)
{
D = (int)y.size() - 1;
Y.resize(D + 1);
for (int i = 0; i <= D; i++)
{
Y[i] = y[i];
}
fac.resize(D + 1), finv.resize(D + 1);
fac[0] = 1;
for (int i = 1; i <= D; i++)
fac[i] = fac[i - 1] * i;
finv[D] = fac[D].inv();
for (int i = D - 1; i >= 0; i--)
finv[i] = finv[i + 1] * (i + 1);
prodl.resize(D + 2), prodr.resize(D + 2);
}
};
template<class T, bool is_ntt_friendly>
struct FormalPowerSeries : vector<T>
{
private:
static vector<T> fac, finv, invmint;
void calc(int n)
{
while ((int)fac.size() <= n)
{
int i = fac.size();
fac.emplace_back(fac[i - 1] * i);
invmint.emplace_back(-invmint[T::mod() % i] * (T::mod() / i));
finv.emplace_back(finv[i - 1] * invmint[i]);
}
}
public:
T get_fac(int n) { calc(n); return fac[n]; }
T get_finv(int n) { calc(n); return finv[n]; }
using vector<T>::vector;
using vector<T>::operator=;
using F = FormalPowerSeries;
using S = vector<pair<ll, T>>;
FormalPowerSeries(const S &f, int n = -1)
{
if (n == -1)
n = f.back().first + 1;
(*this).assign(n, T(0));
for (auto [d, a] : f)
(*this)[d] += a;
}
F operator-() const
{
F res(*this);
for (auto &a : res)
a = -a;
return res;
}
F &operator*=(const T &k)
{
for (auto &a : *this)
a *= k;
return *this;
}
F operator*(const T &k) const { return F(*this) *= k; }
friend F operator*(const T k, const F &f) { return f * k; }
F &operator/=(const T &k)
{
*this *= k.inv();
return *this;
}
F operator/(const T &k) const { return F(*this) /= k; }
F &operator+=(const F &g)
{
int n = (*this).size(), m = g.size();
(*this).resize(max(n, m), T(0));
for (int i = 0; i < m; i++)
(*this)[i] += g[i];
return *this;
}
F operator+(const F &g) const { return F(*this) += g; }
F &operator-=(const F &g)
{
int n = (*this).size(), m = g.size();
(*this).resize(max(n, m), T(0));
for (int i = 0; i < m; i++)
(*this)[i] -= g[i];
return *this;
}
F operator-(const F &g) const { return F(*this) -= g; }
F &operator<<=(const ll d)
{
int n = (*this).size();
(*this).insert((*this).begin(), min(ll(n), d), T(0));
(*this).resize(n);
return *this;
}
F operator<<(const ll d) const { return F(*this) <<= d; }
F &operator>>=(const ll d)
{
int n = (*this).size();
(*this).erase((*this).begin(), (*this).begin() + min(ll(n), d));
(*this).resize(n, T(0));
return *this;
}
F operator>>(const ll d) const { return F(*this) >>= d; }
F &operator*=(const S &g)
{
int n = (*this).size();
auto [d, c] = g.front();
if (d != 0)
c = 0;
for (int i = n - 1; i >= 0; i--)
{
(*this)[i] *= c;
for (auto &[j, b] : g)
{
if (j == 0)
continue;
if (j > i)
break;
(*this)[i] += (*this)[i - j] * b;
}
}
return *this;
}
F operator*(const S &g) const { return F(*this) *= g; }
F &operator/=(const S &g)
{
int n = (*this).size();
auto [d, c] = g.front();
assert(d == 0 && c != T(0));
T inv_c = c.inv();
for (int i = 0; i < n; i++)
{
for (auto &[j, b] : g)
{
if (j == 0)
continue;
if (j > i)
break;
(*this)[i] -= (*this)[i - j] * b;
}
(*this)[i] *= inv_c;
}
return *this;
}
F operator/(const S &g) const { return F(*this) /= g; }
template<const int MOD>
F convolution2(const vector<static_modint<MOD>> &A, const vector<static_modint<MOD>> &B, const int d = -1) const
{
F res;
if (is_ntt_friendly)
res = convolution(A, B);
else
res = convolution_anymod(A, B);
if (d != -1 && (int)res.size() > d)
res.resize(d);
return res;
}
template<const int id>
F convolution2(const vector<dynamic_modint<id>> &A, const vector<dynamic_modint<id>> &B, const int d = -1) const
{
F res;
res = convolution_anymod(A, B);
if (d != -1 && (int)res.size() > d)
res.resize(d);
return res;
}
F &operator*=(const F &g)
{
int n = (*this).size();
if (n == 0)
return *this;
*this = convolution2(*this, g, n);
return *this;
}
F operator*(const F &g) const { return F(*this) *= g; }
template <const int MOD>
void butterfly2(FormalPowerSeries<static_modint<MOD>, true> &A) const { internal::butterfly(A); }
template <const int MOD>
void butterfly2(FormalPowerSeries<static_modint<MOD>, false> &A) const { assert(false); }
template <const int id>
void butterfly2(FormalPowerSeries<dynamic_modint<id>, false> &A) const { assert(false); }
template <const int MOD>
void butterfly_inv2(FormalPowerSeries<static_modint<MOD>, true> &A) const { internal::butterfly_inv(A); }
template <const int MOD>
void butterfly_inv2(FormalPowerSeries<static_modint<MOD>, false> &A) const { assert(false); }
template <const int id>
void butterfly_inv2(FormalPowerSeries<dynamic_modint<id>, false> &A) const { assert(false); }
F circular_mod(int n) const
{
F res(n, T(0));
for (int i = 0; i < (int)(*this).size(); i++)
res[i % n] += (*this)[i];
return res;
}
F inv(int d = -1) const
{
int n = (*this).size();
assert(!(*this).empty() && (*this).at(0) != T(0));
if (d == -1)
d = n;
F f, g2;
F g{(*this).front().inv()};
while ((int)g.size() < d)
{
if (is_ntt_friendly)
{
int m = g.size();
f = F{(*this).begin(), (*this).begin() + min(n, 2 * m)};
g2 = F(g);
f.resize(2 * m, T(0)), butterfly2(f);
g2.resize(2 * m, T(0)), butterfly2(g2);
for (int i = 0; i < 2 * m; i++)
f[i] *= g2[i];
butterfly_inv2(f);
f.erase(f.begin(), f.begin() + m);
f.resize(2 * m, T(0)), butterfly2(f);
for (int i = 0; i < 2 * m; i++)
f[i] *= g2[i];
butterfly_inv2(f);
T iz = T(2 * m).inv();
iz *= -iz;
for (int i = 0; i < m; i++)
f[i] *= iz;
g.insert(g.end(), f.begin(), f.begin() + m);
}
else
{
g.resize(2 * g.size(), T(0));
g *= F{T(2)} - g * (*this);
}
}
return {g.begin(), g.begin() + d};
}
F &operator/=(const F &g)
{
*this *= g.inv((*this).size());
return *this;
}
F operator/(const F &g) const { return F(*this) *= g.inv((*this).size()); }
F differentiate()
{
*this >>= 1;
for (int i = 0; i < int((*this).size()) - 1; i++)
(*this)[i] *= i + 1;
return *this;
}
F differential() const { return F(*this).differentiate(); }
F integrate()
{
int n = (*this).size();
vector<T> minv(n);
minv[1] = T(1);
*this <<= 1;
for (int i = 2; i < n; i++)
{
minv[i] = -minv[T::mod() % i] * (T::mod() / i);
(*this)[i] *= minv[i];
}
return *this;
}
F integral() const { return F(*this).integrate(); }
F log() const
{
assert((*this).front() == T(1));
return ((*this).differential() / (*this)).integral();
}
F exp() const // https://arxiv.org/pdf/1301.5804.pdf
{
int n = (*this).size();
assert(n != 0 && (*this).front() == T(0));
//*
if (is_ntt_friendly)
{
F f{T(1)}, g{T(1)};
F dh = (*this).differential();
F f2, g2, f3, q, s, h, u;
g2 = {T(0)};
while ((int)f.size() < n)
{
int m = f.size();
T im = T(m).inv(), i2m = T(2 * m).inv();
f2 = F(f);
f2.resize(2 * m), butterfly2(f2);
F f3(f);
butterfly2(f3);
for (int i = 0; i < m; i++)
f3[i] *= g2[i];
butterfly_inv2(f3);
f3.erase(f3.begin(), f3.begin() + m / 2);
f3.resize(m, T(0)), butterfly2(f3);
for (int i = 0; i < m; i++)
f3[i] *= g2[i];
butterfly_inv2(f3);
for (int i = 0; i < m / 2; i++)
f3[i] *= -im * im;
g.insert(g.end(), f3.begin(), f3.begin() + m / 2);
g2 = F(g), g2.resize(2 * m), butterfly2(g2);
q = F(dh);
q.resize(2 * m);
for (int i = m - 1; i < 2 * m; i++)
q[i] = T(0);
butterfly2(q);
for (int i = 0; i < 2 * m; i++)
q[i] *= f2[i];
butterfly_inv2(q);
q = q.circular_mod(m);
for (int i = 0; i < m; i++)
q[i] *= i2m;
q.resize(m + 1);
s = ((f.differential() - q) << 1).circular_mod(m);
s.resize(2 * m);
butterfly2(s);
for (int i = 0; i < 2 * m; i++)
s[i] *= g2[i];
butterfly_inv2(s);
for (int i = 0; i < m; i++)
s[i] *= i2m;
s.resize(m);
h = (*this);
h.resize(2 * m), s.resize(2 * m);
u = (h - (s << (m - 1)).integral()) >> m;
butterfly2(u);
for (int i = 0; i < 2 * m; i++)
u[i] *= f2[i];
butterfly_inv2(u);
for (int i = 0; i < m; i++)
u[i] *= i2m;
u.resize(m);
f.insert(f.end(), u.begin(), u.end());
}
return {f.begin(), f.begin() + n};
}
else
//*/
{
F f{T(1)}, g{T(1)};
while ((int)f.size() < n)
{
int m = f.size();
g = convolution2(g, F{T(2)} - f * g, m);
F q = (*this).differential();
q.resize(m - 1);
F r = f.convolution2(f, q).circular_mod(m);
r.resize(m + 1);
F s = ((f.differential() - r) << 1).circular_mod(m);
F t = g * s;
F h = (*this);
h.resize(2 * m), t.resize(2 * m);
F u = (h - (t << (m - 1)).integral()) >> m;
F v = f * u;
f.insert(f.end(), v.begin(), v.end());
}
return {f.begin(), f.begin() + n};
/*
F f{T(1)};
while ((int)f.size() < n)
{
int m = f.size();
f.resize(min(n, 2 * m), T(0));
f *= (*this) + F{T(1)} - f.log();
}
return f;
//*/
}
}
F pow(const ll k) const
{
if (k == 0)
{
F res((*this).size(), T(0));
res[0] = T(1);
return res;
}
int n = (*this).size(), d;
for (d = 0; d < n; d++)
{
if ((*this)[d] != T(0))
break;
}
if (d == n)
return F(n, 0);
F res = F(*this) >> d;
T c = res[0];
res /= c;
res = (res.log() * T(k)).exp();
res *= c.pow(k), res <<= (d != 0 && k > n ? n : d * k);
return res;
}
F div_poly(const F &g) const
{
F f2 = F(*this), g2 = F(g);
while (!f2.empty() && f2.back() == T(0))
f2.pop_back();
while (!g2.empty() && g2.back() == T(0))
g2.pop_back();
int n = f2.size(), m = g2.size();
int k = n - m + 1;
if (k <= 0)
return F{};
reverse(f2.begin(), f2.end());
reverse(g2.begin(), g2.end());
f2.resize(k, T(0)), g2.resize(k, T(0));
F q = f2 / g2;
reverse(q.begin(), q.end());
while (!q.empty() && q.back() == T(0))
q.pop_back();
return q;
}
pair<F, F> divmod(const F &g) const
{
int m = g.size();
assert(m != 0);
F q = (*this).div_poly(g);
F f3 = F(*this), g3 = F(g), q3 = F(q);
f3.resize(m - 1, T(0)), g3.resize(m - 1, T(0)), q3.resize(m - 1, T(0));
F r = f3 - q3 * g3;
while (!r.empty() && r.back() == T(0))
r.pop_back();
return make_pair(q, r);
}
F operator%(const F &g) const { return (*this).divmod(g).second; }
F &operator%=(const F &g) { return (*this) = (*this) % g; }
F div_poly(const S &g) const
{
F f2 = F(*this);
while (!f2.empty() && f2.back() == T(0))
f2.pop_back();
assert(!g.empty());
int n = f2.size(), m = g.back().first + 1;
int k = n - m + 1;
if (k <= 0)
return F{};
reverse(f2.begin(), f2.end());
S g2(g.size());
for (int i = 0; i < (int)g.size(); i++)
g2[(int)g.size() - 1 - i] = make_pair(m - 1 - g[i].first, g[i].second);
f2.resize(k, T(0));
F q = f2 / g2;
reverse(q.begin(), q.end());
while (!q.empty() && q.back() == T(0))
q.pop_back();
return q;
}
pair<F, F> divmod(const S &g) const
{
assert(!g.empty());
int m = g.back().first + 1;
F q = (*this).div_poly(g);
F f3 = F(*this), q3 = F(q);
f3.resize(m - 1, T(0)), q3.resize(m - 1, T(0));
F r = f3 - q3 * g;
while (!r.empty() && r.back() == T(0))
r.pop_back();
return make_pair(q, r);
}
F operator%(const S &g) const { return (*this).divmod(g).second; }
F &operator%=(const S &g) { return (*this) = (*this) % g; }
F to_egf()
{
for (int i = 0; i < (int)(*this).size(); i++)
(*this)[i] *= get_finv(i);
return (*this);
}
F to_ogf()
{
for (int i = 0; i < (int)(*this).size(); i++)
(*this)[i] *= get_fac(i);
return (*this);
}
F get_ogf() const { return F(*this).to_ogf(); }
F taylor_shift(const T &c) const
{
int n = (*this).size();
F f = F(*this).get_ogf();
reverse(f.begin(), f.end());
F g = F(n);
g[0] = 1;
for (int i = 1; i < n; i++)
g[i] = c * g[i - 1];
g.to_egf();
F h = f * g;
reverse(h.begin(), h.end());
return h.to_egf();
}
};
template <class T, bool is_ntt_friendly>
struct SparseFormalPowerSeries : vector<pair<ll, T>>
{
using vector<pair<ll, T>>::vector;
using vector<pair<ll, T>>::operator=;
using F = FormalPowerSeries<T, is_ntt_friendly>;
using S = SparseFormalPowerSeries;
F to_fps(int n) const
{
F res(n, T(0));
for (auto [d, a] : (*this))
res[d] += a;
return res;
}
SparseFormalPowerSeries(const F &f)
{
(*this).clear();
for (int i = 0; i < (int)f.size(); i++)
{
if (f[i] != T(0))
(*this).emplace_back(make_pair(i, f[i]));
}
}
S operator-() const
{
S res(*this);
for (auto &[d, a] : res)
a = -a;
return res;
}
S operator*=(const T &k)
{
for (auto &[d, a] : (*this))
a *= k;
return (*this);
}
S operator/=(const T &k)
{
(*this) *= k.inv();
return (*this);
}
S operator*(const T &k) const { return S(*this) *= k; }
S operator/(const T &k) const { return S(*this) /= k; }
friend S operator*(const T k, const S &f) { return f * k; }
S operator+(const S &g) const
{
S res;
int n = (*this).size(), m = g.size(), i = 0, j = 0;
while (i < n || j < m)
{
pair<ll, T> tmp;
if (j == m || (i != n && (*this)[i].first <= g[j].first))
tmp = (*this)[i++];
else
tmp = g[j++];
if (!res.empty() && res.back().first == tmp.first)
res.back().second += tmp.second;
else
res.emplace_back(tmp);
}
return res;
}
S operator-(const S &g) const
{
S res;
int n = (*this).size(), m = g.size(), i = 0, j = 0;
while (i < n || j < m)
{
pair<ll, T> tmp;
if (j == m || (i != n && (*this)[i].first <= g[j].first))
tmp = (*this)[i++];
else
{
tmp = g[j++];
tmp.second = -tmp.second;
}
if (!res.empty() && res.back().first == tmp.first)
res.back().second += tmp.second;
else
res.emplace_back(tmp);
}
return res;
}
S operator*(const S &g) const
{
S res;
for (auto [d, a] : (*this))
for (auto [e, b] : g)
res.emplace_back(make_pair(d + e, a * b));
sort(res.begin(), res.end(), [&](pair<ll, T> p1, pair<ll, T> p2)
{ return p1.first < p2.first; });
S res2;
for (auto da : res)
{
auto [d, a] = da;
if (res2.empty() || res2.back().first != d)
res2.emplace_back(da);
else
res2.back().second += a;
}
return res2;
}
S operator+=(const S &g) { return (*this) = (*this) + g; }
S operator-=(const S &g) { return (*this) = (*this) - g; }
S operator*=(const S &g) { return (*this) = (*this) * g; }
S operator<<=(ll k)
{
for (auto &[d, a] : (*this))
d += k;
return (*this);
}
S operator<<(ll k) const { return (*this) <<= k; }
S operator>>(ll k) const
{
S res;
for (auto [d, a] : (*this))
{
d -= k;
if (d >= 0)
res.emplace_back(make_pair(d, a));
}
return res;
}
S operator>>=(ll k) { return (*this) = (*this) >> k; }
F inv(int n) const
{
F f(n, T(0));
f.front() = T(1);
return f / (*this);
}
S differentiate()
{
for (auto &[d, a] : (*this))
a *= d--;
if (!(*this).empty() && (*this).front().first == -1)
(*this).erase((*this).begin());
return (*this);
}
S differential() const { return S(*this).differentiate(); }
S integrate()
{
for (auto &[d, a] : (*this))
a /= T(++d);
return (*this);
}
S integral() const { return S(*this).integrate(); }
F log(int n) const
{
F f = (*this).to_fps(n);
return (f.differential() / (*this)).integral();
}
F diffeq(const S &a, const S &b, int n) const
{
assert(a.front().first == 0 && a.front().second == 1);
vector<T> minv(n);
minv[1] = T(1);
for (int i = 2; i < n; i++)
minv[i] = -minv[T::mod() % i] * (T::mod() / i);
F f(n, T(0));
f[0] = T(1);
for (int k = 0; k < n - 1; k++)
{
for (auto [i, ai] : a)
{
if (0 <= k - i + 1 && k - i + 1 < k + 1)
f[k + 1] -= ai * (k - i + 1) * f[k - i + 1];
}
for (auto [j, bj] : b)
{
if (0 <= k - j && k - j < k + 1)
f[k + 1] -= bj * f[k - j];
}
f[k + 1] *= minv[k + 1];
}
return f;
}
F exp(int n) const
{
return diffeq(S{{0, 1}}, -((*this).differential()), n);
}
F pow(ll m, int n) const
{
S f(*this);
if (f.empty())
{
F res(n, T(0));
if (m == 0)
res.front() = T(1);
return res;
}
auto [d0, a0] = f.front();
T a0_inv = a0.inv();
for (auto &[d, a] : f)
d -= d0, a *= a0_inv;
if (m >= 0)
{
F g = diffeq(f, -m * f.differential(), n);
return (g * a0.pow(m)) << mul_limited(d0, m);
}
else
{
F g = diffeq(f, -m * f.differential(), n + (d0 * (-m)));
F h = (g * a0_inv.pow(-m)) >> (d0 * (-m));
h.resize(n);
return h;
}
}
};
template <class T, bool is_ntt_friendly>
vector<T> FormalPowerSeries<T, is_ntt_friendly>::fac{1, 1};
template <class T, bool is_ntt_friendly>
vector<T> FormalPowerSeries<T, is_ntt_friendly>::finv{1, 1};
template<class T, bool is_ntt_friendly>
vector<T> FormalPowerSeries<T, is_ntt_friendly>::invmint{0, 1};
template<class T, bool is_ntt_friendly>
struct RationalFormalPowerSeries
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
using R = RationalFormalPowerSeries;
F num, den;
R operator-() const
{
R res(*this);
res.num = -res.num;
return res;
}
R operator*=(const T &k)
{
(*this).num *= k;
return *this;
}
R operator*(const T &k) const { return R(*this) *= k; }
friend R operator*(const T k, const R &r) { return r * k; }
R operator/=(const T &k)
{
(*this).den *= k;
return k;
}
R operator/(const T &k) const { return R(*this) /= k; }
R &operator+=(const R &r)
{
F f, g;
f = f.convolution2((*this).num, r.den);
g = g.convolution2((*this).den, r.num);
(*this).num = f + g;
(*this).den = (*this).den.convolution2((*this).den, r.den);
return *this;
}
R operator+(const R &r) const { return R(*this) += r; }
R &operator-=(const R &r)
{
F f, g;
f = f.convolution2((*this).num, r.den);
g = g.convolution2((*this).den, r.num);
(*this).num = f - g;
(*this).den = (*this).den.convolution2((*this).den, r.den);
return *this;
}
R operator-(const R &r) const { return R(*this) -= r; }
R operator*=(const R &r)
{
(*this).num = (*this).num.convolution2((*this).num, r.num);
(*this).den = (*this).den.convolution2((*this).den, r.den);
return *this;
}
R operator*(const R &r) const { return R(*this) *= r; }
R operator/=(const R &r)
{
(*this).num = (*this).num.convolution2((*this).num, r.den);
(*this).den = (*this).den.convolution2((*this).den, r.num);
return *this;
}
R operator/(const R &r) const { return R(*this) /= r; }
R inv()
{
R res(*this);
swap(res.num, res.den);
return res;
}
};
template<class T, bool is_ntt_friendly>
vector<T> sample_points_shift(const vector<T> &ys, int M, T c)
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
F f;
int N = ys.size();
vector<T> a;
{
vector<T> p(N), q(N);
for (int i = 0; i < N; i++)
{
p[i] = ys[i] * f.get_finv(i);
q[i] = i % 2 == 0 ? f.get_finv(i) : -f.get_finv(i);
}
a = f.convolution2(p, q);
a.resize(N);
}
vector<T> b;
{
vector<T> p(N), q(N);
T tmp = 1;
for (int i = 0; i < N; i++)
{
p[i] = a[i] * f.get_fac(i);
q[i] = tmp * f.get_finv(i);
tmp *= c - i;
}
reverse(q.begin(), q.end());
b = f.convolution2(p, q);
b.erase(b.begin(), b.begin() + N - 1);
for (int i = 0; i < N; i++)
b[i] *= f.get_finv(i);
}
vector<T> res;
{
vector<T> p(M);
for (int i = 0; i < M; i++)
p[i] = f.get_finv(i);
res = f.convolution2(b, p);
res.resize(M);
for (int i = 0; i < M; i++)
res[i] *= f.get_fac(i);
}
return res;
}
template<class T, bool is_ntt_friendly>
struct FactorialFast
{
private:
const int P, K;
vector<T> Y, Z, fac;
public:
FactorialFast(const int K = 9) : P(T::mod()), K(K)
{
Y = {1};
for (int i = 0; i < K; i++)
{
Z = sample_points_shift<T, is_ntt_friendly>(Y, (1 << (i + 2)) - (1 << i), 1 << i);
Z.insert(Z.begin(), Y.begin(), Y.end());
Y.resize(1 << (i + 1));
for (int j = 0; j < (1 << (i + 1)); j++)
Y[j] = Z[2 * j] * Z[2 * j + 1] * T::raw((1 << i) * (2 * j + 1));
}
if ((1 << K) <= P / (1 << K))
{
Z = sample_points_shift<T, is_ntt_friendly>(Y, P / (1 << K), 1 << K);
Y.insert(Y.end(), Z.begin(), Z.end());
}
fac.resize(P / (1 << K) + 1);
fac.at(0) = 1;
for (int i = 0; i < P / (1 << K); i++)
fac[i + 1] = fac[i] * Y[i] * T::raw((1 + i) * (1 << K));
}
};
template<class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> stirling1_fixed_n(const int &N)
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
using S = SparseFormalPowerSeries<T, is_ntt_friendly>;
if (N == 0)
return {1};
if (N == 1)
return {0, 1};
if (N & 1)
{
F f = stirling1_fixed_n<T, is_ntt_friendly>(N - 1);
f.resize(N + 1, T(0));
return f * S{{0, 1 - N}, {1, 1}};
}
else
{
F f = stirling1_fixed_n<T, is_ntt_friendly>(N / 2);
f.resize(N + 1, T(0));
F g = f.taylor_shift(-(N / 2));
return f * g;
}
}
using fps = FormalPowerSeries<mint, ntt>;
using sfps = SparseFormalPowerSeries<mint, ntt>;
//*
#include <atcoder/all>
//*/
namespace fast_prime
{
bool is_prime(ll n)
{
if (n <= 1) return false;
static const ll as[7] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
static const ll ps[9] = {2, 3, 5, 13, 19, 73, 193, 407521, 299210837};
for (auto &&p : ps)
{
if (n == p) return true;
if (n % p == 0) return false;
}
ll d = n - 1;
int s = 0;
while (d % 2 == 0) d /= 2, s++;
for (auto &&a : as)
{
ll a0 = 1;
for (ll d2 = d, tmp = a; d2 > 0; d2 /= 2, tmp = __int128_t(tmp) * tmp % n)
{
if (d2 % 2 != 0)
a0 = __int128_t(a0) * tmp % n;
}
if (a0 == 1 || a0 == n - 1)
continue;
for (int r = 1; r <= s; r++)
{
if (r == s)
return false;
a0 = __int128_t(a0) * a0 % n;
if (a0 == n - 1)
break;
}
}
return true;
}
ll get_prime_factor(ll n)
{
int m = pow(n, .125);
for (int c = 1; ; c++)
{
auto f = [&](ll a) -> ll
{ return (__int128_t(a) * a + c) % n; };
ll x = 2, y = 2, prod = 1, g = 1;
for (int t = 1; g == 1; t = min(2 * t, m))
{
for (int i = 0; i < t; i++)
{
x = f(x), y = f(f(y));
prod = __int128_t(prod) * (x - y) % n;
}
g = gcd(prod, n);
}
if (g == n) continue;
return is_prime(g) ? g : is_prime(n / g) ? n / g : get_prime_factor(g);
}
}
vector<ll> factorize(ll n)
{
vector<ll> res;
for (int p = 2; p < 100; p++)
{
while (n % p == 0)
{
n /= p;
res.emplace_back(p);
}
}
while (n > 1)
{
if (is_prime(n))
{
res.emplace_back(n);
break;
}
ll p = get_prime_factor(n);
n /= p;
res.emplace_back(p);
}
sort(res.begin(), res.end());
return res;
}
vector<tuple<ll, int, ll>> ord_pow(const vector<ll> &ps)
{
vector<tuple<ll, int, ll>> res;
for (auto &&p : ps)
{
if (res.empty() || get<0>(res.back()) != p)
res.emplace_back(make_tuple(p, 1, p));
else
get<1>(res.back())++, get<2>(res.back()) *= p;
}
return res;
}
vector<ll> divisors(const vector<tuple<ll, int, ll>> &peqs)
{
vector<ll> ds;
auto dfs = [&](auto self, ll d, int i) -> void
{
if (i == (int)peqs.size())
{
ds.emplace_back(d);
return;
}
auto &&[p, e, q] = peqs[i];
for (ll r = 1, j = 0; j <= e; r *= p, j++) self(self, d * r, i + 1);
};
dfs(dfs, 1, 0);
sort(ds.begin(), ds.end());
return ds;
}
vector<ll> divisors(const vector<ll> &ps) { return divisors(ord_pow(ps)); }
vector<ll> divisors(ll n) { return divisors(factorize(n)); }
}
using namespace fast_prime;
template <class T>
struct Binomial
{
private:
static vector<T> _fac, _finv, _inv;
public:
static void calc(int n)
{
int i = _fac.size();
if (n < i)
return;
_fac.resize(n + 1), _finv.resize(n + 1), _inv.resize(n + 1);
for (; i <= n; i++)
{
_fac[i] = _fac[i - 1] * i;
_inv[i] = -_inv[T::mod() % i] * (T::mod() / i);
_finv[i] = _finv[i - 1] * _inv[i];
}
}
static T inv(int n) { assert(n > 0); calc(n); return _inv[n]; }
};
template <class T> vector<T> Binomial<T>::_fac{1, 1};
template <class T> vector<T> Binomial<T>::_finv{1, 1};
template <class T> vector<T> Binomial<T>::_inv{0, 1};
template <class T = ll, class U = i128>
pair<T, T> svp2(const pair<T, T> &a, const pair<T, T> &b)
{
assert(a != make_pair(0, 0) && b != make_pair(0, 0));
auto [a1, a2] = a;
auto [b1, b2] = b;
if ((U)a1 * a1 + (U)a2 * a2 < (U)b1 * b1 + (U)b2 * b2)
swap(a1, b1), swap(a2, b2);
while ((U)a1 * a1 + (U)a2 * a2 > (U)b1 * b1 + (U)b2 * b2)
{
swap(a1, b1), swap(a2, b2);
T k = divround<U>((U)a1 * b1 + (U)a2 * b2, (U)a1 * a1 + (U)a2 * a2);
b1 -= k * a1, b2 -= k * a2;
}
return make_pair(a1, a2);
}
template <class T = int, class U = ll>
pair<T, T> mint_to_rat(const T &r, const T &m)
{
auto [p, q] = svp2<T, U>(make_pair(r, T(1)), make_pair(m, T(0)));
if (q < 0)
p = -p, q = -q;
return make_pair(p, q);
}
template <class M, class T = int, class U = ll>
pair<T, T> mint_to_rat(const M &x) { return mint_to_rat<T, U>((T)x.val(), (T)M::mod()); }
void init() {}
void main2()
{
LL(N, M, L, R);
fps f(M + 1);
rep(i, 1, M + 1) f.at(i) = 1;
sfps gnum = {{L, 1}, {R + 1, -1}};
sfps gden = {{0, 1}, {1, -1}};
fps h = f * gnum.pow(N - 1, M + 1) * gden.pow(-(N - 1), M + 1);
mint ans = 0;
rep(i, M + 1) ans += h.at(i);
PRINT(ans);
}
void test()
{
/*
//*/
}
}
int main()
{
cauto CERR = [](string val, string color)
{
string s = "\033[" + color + "m" + val + "\033[m";
/* コードテストで確認する際にコメントアウトを外す
cerr << val;
//*/
};
CERR("\n[FAST_IO]\n\n", "34");
cin.tie(0);
ios::sync_with_stdio(false);
CERR("\n[FAST_IO]\n\n", "32");
cout << fixed << setprecision(20);
test();
init();
CERR("\n[SINGLE_TESTCASE]\n\n", "36");
main2();
}