結果
| 問題 | No.1889 K Consecutive Ks (Hard) |
| ユーザー |
 miscalc
|
| 提出日時 | 2023-06-06 09:16:27 |
| 言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 71 ms / 6,000 ms |
| コード長 | 36,676 bytes |
| コンパイル時間 | 3,961 ms |
| コンパイル使用メモリ | 243,652 KB |
| 最終ジャッジ日時 | 2025-02-13 23:02:03 |
|
ジャッジサーバーID （参考情報） |
judge1 / judge2 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 3 |
| other | AC * 22 |
ソースコード
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using ld = long double;
using pll = pair<ll, ll>;
using tlll = tuple<ll, ll, ll>;
constexpr ll INF = 1LL << 60;
template<class T> bool chmin(T& a, T b) {if (a > b) {a = b; return true;} return false;}
template<class T> bool chmax(T& a, T b) {if (a < b) {a = b; return true;} return false;}
ll safemod(ll A, ll M) {ll res = A % M; if (res < 0) res += M; return res;}
ll divfloor(ll A, ll B) {if (B < 0) A = -A, B = -B; return (A - safemod(A, B)) / B;}
ll divceil(ll A, ll B) {if (B < 0) A = -A, B = -B; return divfloor(A + B - 1, B);}
ll pow_ll(ll A, ll B) {if (A == 0 || A == 1) {return A;} if (A == -1) {return B & 1 ? -1 : 1;} ll res = 1; for (int i = 0; i < B; i++) {res *= A;} return res;}
ll mul_limited(ll A, ll B, ll M = INF) { return B == 0 ? 0 : A > M / B ? M : A * B; }
ll pow_limited(ll A, ll B, ll M = INF) { if (A == 0 || A == 1) {return A;} ll res = 1; for (int i = 0; i < B; i++) {if (res > M / A) return M; res *= A;} return res;}
ll logfloor(ll A, ll B) {assert(A >= 2); ll res = 0; for (ll tmp = 1; tmp <= B / A; tmp *= A) {res++;} return res;}
ll logceil(ll A, ll B) {assert(A >= 2); ll res = 0; for (ll tmp = 1; tmp < B; tmp *= A) {res++;} return res;}
ll arisum_ll(ll a, ll d, ll n) { return n * a + (n & 1 ? ((n - 1) >> 1) * n : (n >> 1) * (n - 1)) * d; }
ll arisum2_ll(ll a, ll l, ll n) { return n & 1 ? ((a + l) >> 1) * n : (n >> 1) * (a + l); }
ll arisum3_ll(ll a, ll l, ll d) { assert((l - a) % d == 0); return arisum2_ll(a, l, (l - a) / d + 1); }
template<class T> void unique(vector<T> &V) {V.erase(unique(V.begin(), V.end()), V.end());}
template<class T> void sortunique(vector<T> &V) {sort(V.begin(), V.end()); V.erase(unique(V.begin(), V.end()), V.end());}
#define FINALANS(A) do {cout << (A) << '\n'; exit(0);} while (false)
template<class T> void printvec(const vector<T> &V) {int _n = V.size(); for (int i = 0; i < _n; i++) cout << V[i] << (i == _n - 1 ? "" : " ");cout << '\n';}
template<class T> void printvect(const vector<T> &V) {for (auto v : V) cout << v << '\n';}
template<class T> void printvec2(const vector<vector<T>> &V) {for (auto &v : V) printvec(v);}
//*
#include <atcoder/modint>
#include <atcoder/math>
#include <atcoder/convolution>
#include <atcoder/internal_math>
using namespace atcoder;
//*/
template<const int MOD = 1000000007, 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) <= 300)
{
using mint = static_modint<MOD>;
vector<mint> A2(N), B2(M);
for (int i = 0; i < N; i++)
A2[i] = A[i];
for (int j = 0; j < M; j++)
B2[j] = B[j];
vector<mint> C2(N + M - 1, 0);
for (int i = 0; i < N; i++)
for (int j = 0; j < M; j++)
C2[i + j] += A2[i] * B2[j];
vector<T> C(N + M - 1);
for (int i = 0; i < N + M - 1; i++)
C[i] = C2[i].val();
return C;
}
constexpr ll MOD1 = 167772161, MOD2 = 469762049, MOD3 = 1224736769;
using mint2 = static_modint<MOD2>;
using mint3 = static_modint<MOD3>;
using mint4 = static_modint<MOD>;
constexpr int i1_2 = internal::inv_gcd(MOD1, MOD2).second;
constexpr int i12_3 = internal::inv_gcd(MOD1 * MOD2, MOD3).second;
constexpr int m12_4 = MOD1 * MOD2 % MOD;
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();
mint3 x2_m3 = mint3::raw(c1) + mint3::raw(t1) * mint3::raw(MOD1);
mint4 x2_m = mint4::raw(c1) + mint4::raw(t1) * mint4::raw(MOD1);
int t2 = ((mint3::raw(c3) - x2_m3) * mint3::raw(i12_3)).val();
C[i] = (x2_m + mint4::raw(t2) * mint4::raw(m12_4)).val();
}
return C;
}
template<class T>
vector<T> convolution_anymod(const vector<T> &A, const vector<T> &B, const int MOD)
{
int N = A.size(), M = B.size();
if (min(N, M) <= 300)
{
using mint = dynamic_modint<100>;
mint::set_mod(MOD);
vector<mint> A2(N), B2(M);
for (int i = 0; i < N; i++)
A2[i] = A[i];
for (int j = 0; j < M; j++)
B2[j] = B[j];
vector<mint> C2(N + M - 1, 0);
for (int i = 0; i < N; i++)
for (int j = 0; j < M; j++)
C2[i + j] += A2[i] * B2[j];
vector<T> C(N + M - 1);
for (int i = 0; i < N + M - 1; i++)
C[i] = C2[i].val();
return C;
}
constexpr ll MOD1 = 167772161, MOD2 = 469762049, MOD3 = 1224736769;
using mint2 = static_modint<MOD2>;
using mint3 = static_modint<MOD3>;
using mint4 = dynamic_modint<100>;
mint4::set_mod(MOD);
constexpr int i1_2 = internal::inv_gcd(MOD1, MOD2).second;
constexpr int i12_3 = internal::inv_gcd(MOD1 * MOD2, MOD3).second;
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();
mint3 x2_m3 = mint3::raw(c1) + mint3::raw(t1) * mint3::raw(MOD1);
mint4 x2_m = mint4::raw(c1) + mint4::raw(t1) * mint4::raw(MOD1);
int t2 = ((mint3::raw(c3) - x2_m3) * mint3::raw(i12_3)).val();
C[i] = (x2_m + mint4::raw(t2) * mint4::raw(MOD1) * mint4::raw(MOD2)).val();
}
return C;
}
template<const int MOD>
vector<static_modint<MOD>> convolution_anymod(const vector<static_modint<MOD>> &A, const vector<static_modint<MOD>> &B)
{
int N = A.size(), M = B.size();
vector<int> A2(N), B2(M);
for (int i = 0; i < N; i++)
A2[i] = A[i].val();
for (int i = 0; i < M; i++)
B2[i] = B[i].val();
vector<int> C2 = convolution_anymod<MOD>(A2, B2);
vector<static_modint<MOD>> C(N + M - 1);
for (int i = 0; i < N + M - 1; i++)
C[i] = static_modint<MOD>::raw(C2[i]);
return C;
}
template<const int id>
vector<dynamic_modint<id>> convolution_anymod(const vector<dynamic_modint<id>> &A, const vector<dynamic_modint<id>> &B)
{
int N = A.size(), M = B.size();
vector<int> A2(N), B2(M);
for (int i = 0; i < N; i++)
A2[i] = A[i].val();
for (int i = 0; i < M; i++)
B2[i] = B[i].val();
vector<int> C2 = convolution_anymod(A2, B2, dynamic_modint<id>::mod());
vector<dynamic_modint<id>> C(N + M - 1);
for (int i = 0; i < N + M - 1; i++)
C[i] = dynamic_modint<id>::raw(C2[i]);
return C;
}
// https://opt-cp.com/fps-implementation/
// https://qiita.com/hotman78/items/f0e6d2265badd84d429a
// https://opt-cp.com/fps-fast-algorithms/
// https://maspypy.com/%E5%A4%9A%E9%A0%85%E5%BC%8F%E3%83%BB%E5%BD%A2%E5%BC%8F%E7%9A%84%E3%81%B9%E3%81%8D%E7%B4%9A%E6%95%B0-%E9%AB%98%E9%80%9F%E3%81%AB%E8%A8%88%E7%AE%97%E3%81%A7%E3%81%8D%E3%82%8B%E3%82%82%E3%81%AE
template<class T, bool is_ntt_friendly>
struct FormalPowerSeries : vector<T>
{
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; }
// (1 + cx^d) を掛ける
F multiply(const int d, const T c)
{
int n = (*this).size();
if (c == T(1))
{
for (int i = n - 1 - d; i >= 0; i--)
(*this)[i + d] += (*this)[i];
}
else if (c == T(-1))
{
for (int i = n - 1 - d; i >= 0; i--)
(*this)[i + d] -= (*this)[i];
}
else
{
for (int i = n - 1 - d; i >= 0; i--)
(*this)[i + d] += (*this)[i] * c;
}
return *this;
}
F multiplication(const int d, const T c) const { return multiply(F(*this)); }
// (1 + cx^d) で割る
F divide(const int d, const T c)
{
int n = (*this).size();
if (c == T(1))
{
for (int i = 0; i < n - d; i++)
(*this)[i + d] -= (*this)[i];
}
else if (c == T(-1))
{
for (int i = 0; i < n - d; i++)
(*this)[i + d] += (*this)[i];
}
else
{
for (int i = 0; i < n - d; i++)
(*this)[i + d] -= (*this)[i] * c;
}
return *this;
}
F division(const int d, const T c) const { return divide(F(*this)); }
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); }
// mod (x^n - 1) をとったものを返す
F circular_mod(int n) const
{
F res(n, T(0));
for (int i = 0; i < (*this).size(); i++)
res[i % n] += (*this)[i];
return res;
}
F inv(int d = -1) const
{
int n = (*this).size();
assert(n != 0 && (*this).front() != 0);
if (d == -1)
d = n;
assert(d > 0);
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);
// a
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);
// b, c
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;
// d, e
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);
// f, g
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);
// h
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 powmod(ll k, const F &g) const
{
F res(2 * g.size(), 0);
res.front() = 1;
F tmp = (*this) % g;
tmp.resize(g.size());
while (k > 0)
{
if (k & 1)
{
res *= tmp;
res %= g;
res.resize(2 * g.size());
}
tmp = tmp.convolution2(tmp, tmp);
tmp %= g;
tmp.resize(g.size());
k >>= 1;
}
return res;
}
// 素数 mod を要求
// 存在しないなら空配列を返す
F sqrt() const
{
int n = (*this).size(), d;
for (d = 0; d < n; d += 2)
{
if ((*this)[d] != 0)
break;
if (d + 1 < n && (*this)[d + 1] != 0)
return F(0);
}
if (d >= n)
return F(n, 0);
T a = (*this)[d];
int p = T::mod();
if (a.pow((p - 1) / 2) == p - 1)
return F(0);
T r;
if (p % 4 == 3)
r = a.pow((p + 1) / 4);
else
{
int q = p - 1, s = 0;
while (q % 2 == 0)
q /= 2, s++;
T z = 2;
while (z.pow((p - 1) / 2) != p - 1)
z++;
int m = s;
T c = z.pow(q);
T t = a.pow(q);
r = a.pow((q + 1) / 2);
while (t != 1)
{
int m2 = 1;
for (T tmp = t * t; tmp != 1; tmp = tmp * tmp, m2++);
T b = c.pow(1 << (m - m2 - 1));
m = m2, c = b * b, t *= c, r *= b;
}
}
T inv_2 = T(2).inv();
F f = F(*this) >> d, res = F{r};
while (res.size() < f.size())
{
res.resize(min(f.size(), 2 * res.size()), T(0));
res = (res + res.inv() * f) * inv_2;
}
res <<= d / 2;
return res;
}
F div_poly(const F &g) const
{
int n = (*this).size(), m = g.size();
int k = n - m + 1;
if (k <= 0)
return F{};
F f2 = F(*this), g2 = F(g);
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; }
T eval(const T &x)
{
T res(0);
for (int i = (int)(*this).size() - 1; i >= 0; i--)
{
res *= x;
res += (*this)[i];
}
return res;
}
F taylor_shift(const T &c)
{
int n = (*this).size();
F fac(n), finv(n);
fac[0] = 1;
for (int i = 1; i < n; i++)
fac[i] = fac[i - 1] * i;
finv[n - 1] = fac[n - 1].inv();
for (int i = n - 2; i >= 0; i--)
finv[i] = finv[i + 1] * (i + 1);
F f = F(*this), g = F(n);
for (int i = 0; i < n; i++)
f[i] *= fac[i];
g[0] = 1;
for (int i = 1; i < n; i++)
g[i] = c * g[i - 1];
for (int i = 0; i < n; i++)
g[i] *= finv[i];
reverse(f.begin(), f.end());
F h = f * g;
reverse(h.begin(), h.end());
for (int i = 0; i < n; i++)
h[i] *= finv[i];
return h;
}
vector<T> eval_multipoint(const vector<T> &xs)
{
int m0 = xs.size(), m = 1;
while (m < m0)
m <<= 1;
vector<F> node(2 * m, F{1});
for (int i = 0; i < m0; i++)
node[m + i] = {-xs[i], T(1)};
for (int i = m - 1; i > 0; i--)
node[i] = convolution2(node[i << 1], node[(i << 1) | 1]);
node[1] = (*this).divmod(node[1]).second;
for (int i = 2; i < m + m0; i++)
node[i] = node[i >> 1].divmod(node[i]).second;
vector<T> res(m0);
for (int i = 0; i < m0; i++)
res[i] = node[m + i].empty() ? T(0) : node[m + i][0];
return res;
}
};
// (次数, 係数) を昇順に並べたもの
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 (*this) *= k; }
S operator/(const T &k) const { return (*this) /= 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] : (*this))
res.emplace_back(make_pair(d + e, a * b));
sort(res.begin(), res.end());
S res2;
for (auto da : res)
{
auto [d, a] = da;
if (res2.empty() || res2.back() != d)
res2.emplace_back(da);
else
res2.back() += a;
}
return res;
}
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 exp(int n) const
{
vector<T> minv(n);
minv[1] = T(1);
for (int i = 2; i < n; i++)
minv[i] = -minv[T::mod() % i] * (T::mod() / i);
S fd = (*this).differential();
F g(n, T(0));
g[0] = T(1);
for (int i = 0; i < n - 1; i++)
{
for (auto [d, a] : fd)
{
if (i - d < 0)
break;
g[i + 1] += a * g[i - d];
}
g[i + 1] *= minv[i + 1];
}
return g;
}
// バグっています
F pow(ll m, int n) const
{
if (m == 0)
{
F res(n, T(0));
res.front() = T(1);
return res;
}
if ((*this).empty())
return F(n, T(0));
vector<T> minv(n);
minv[1] = T(1);
for (int i = 2; i < n; i++)
minv[i] = -minv[T::mod() % i] * (T::mod() / i);
S f = (*this) >> (*this).front().first;
S fd = f.differential();
F g(n, T(0)), gd(n, T(0));
g[0] = f.front().second.pow(m);
int len = m > n ? n - 1 : min(f.back().first * m, ll(n - 1));
for (int i = 0; i < len; i++)
{
for (auto [d, a] : fd)
{
if (i - d < 0)
break;
gd[i] += a * g[i - d];
}
gd[i] *= m;
for (auto [d, a] : f)
{
if (d == 0)
continue;
if (i - d < 0)
break;
gd[i] -= a * gd[i - d];
}
g[i + 1] = gd[i] * minv[i + 1];
}
return g << ((*this).front().first != 0 && m > n ? n : (*this).front().first * m);
}
};
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>
FormalPowerSeries<T, is_ntt_friendly> convolution_many(const vector<FormalPowerSeries<T, is_ntt_friendly>> &fs, int d = -1)
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
if ((int)fs.size() == 0)
return F{1};
deque<F> deq;
for (auto f : fs)
deq.push_back(f);
while ((int)deq.size() > 1)
{
F f = deq.front();
deq.pop_front();
F g = deq.front();
deq.pop_front();
f = f.convolution2(f, g, d);
deq.push_back(f);
}
if (d != -1)
deq.front().resize(d);
return deq.front();
}
template <class T, bool is_ntt_friendly>
RationalFormalPowerSeries<T, is_ntt_friendly> rational_sum(const vector<RationalFormalPowerSeries<T, is_ntt_friendly>> &rs, int d = -1)
{
using R = RationalFormalPowerSeries<T, is_ntt_friendly>;
if (rs.size() == 0)
return R{{1}, {1}};
vector<R> res = vector<R>(rs);
while (res.size() > 1)
{
vector<R> nxt;
for (int i = 0; i < (int)res.size(); i += 2)
{
if (i + 1 < (int)res.size())
nxt.emplace_back(res[i] + res[i + 1]);
else
nxt.emplace_back(res[i]);
if (d != -1)
{
if ((int)nxt.back().num.size() > d)
nxt.back().num.resize(d);
if ((int)nxt.back().den.size() > d)
nxt.back().den.resize(d);
}
}
res = nxt;
}
if (d != -1)
res.front().num.resize(d), res.front().den.resize(d);
return res.front();
}
template <class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> interpolation(const vector<T> &xs, const vector<T> &ys)
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
using R = RationalFormalPowerSeries<T, is_ntt_friendly>;
int n = xs.size();
assert(n == ys.size());
vector<F> fs(n);
for (int i = 0; i < n; i++)
fs[i] = F{-xs[i], T(1)};
F g = convolution_many(fs);
F h = g.differential();
vector<T> a = h.eval_multipoint(xs);
vector<R> rs(n);
for (int i = 0; i < n; i++)
rs[i] = R{F{ys[i] / a[i]}, fs[i]};
R q = rational_sum(rs, n);
return q.num;
}
// prod[d in D](1 + cx^d) を M 次の項まで求める
template <class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> multiply_many(const int &M, const T &c, const vector<int> &D)
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
vector<int> cnt(M + 1, 0);
for (auto d : D)
{
if (d < 0 || M < d)
continue;
cnt[d]++;
}
vector<T> inv(M + 1);
inv[1] = T(1);
for (int i = 2; i <= M; i++)
inv[i] = -inv[T::mod() % i] * (T::mod() / i);
F f(M + 1, 0);
for (int k = 1; k <= M; k++)
{
T pw = 1;
for (int i = 1; k * i <= M; i++)
{
pw *= c;
if (i & 1)
f[k * i] += T::raw(cnt[k]) * pw * inv[i];
else
f[k * i] -= T::raw(cnt[k]) * pw * inv[i];
}
}
return f.exp();
}
// 多重集合 S の要素から何個か選んで総和を 0, 1, …, M にする方法の数
template <class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> subset_sum(const int &M, const vector<int> &S)
{
return multiply_many<T, is_ntt_friendly>(M, T(1), S);
}
// 集合 S の各要素が無限個ある集合 T から何個か選んで総和を 0, 1, …, M にする方法の数
template <class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> partition(const int &M, const vector<int> &S)
{
return multiply_many<T, is_ntt_friendly>(M, T(-1), S).inv();
}
template<class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> stirling1(const int &N)
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
using S = vector<pair<int, T>>;
if (N == 0)
return {1};
if (N == 1)
return {0, 1};
if (N & 1)
{
F f = stirling1<T, is_ntt_friendly>(N - 1);
f.resize(N + 1, T(0));
return f * S{{0, 1 - N}, {1, 1}};
}
else
{
F f = stirling1<T, is_ntt_friendly>(N / 2);
f.resize(N + 1, T(0));
F g = f.taylor_shift(-(N / 2));
return f * g;
}
}
template<class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> stirling2(const int &N)
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
vector<T> fac(N + 1, T(0)), finv(N + 1, T(0));
fac[0] = T(1);
for (int i = 1; i <= N; i++)
fac[i] = fac[i - 1] * i;
finv[N] = fac[N].inv();
for (int i = N - 1; i >= 0; i--)
finv[i] = finv[i + 1] * (i + 1);
vector<int> minfactor(N + 1, -1);
for (int i = 2; i <= N; i++)
{
if (minfactor[i] != -1)
continue;
for (int k = 2 * i; k <= N; k += i)
minfactor[k] = i;
}
vector<T> power(N + 1);
for (int i = 0; i <= N; i++)
{
if (minfactor[i] == -1)
power[i] = T(i).pow(N);
else
power[i] = power[minfactor[i]] * power[i / minfactor[i]];
}
F A(N + 1), B(N + 1);
for (int i = 0; i <= N; i++)
{
A[i] = power[i] * finv[i];
B[i] = (i & 1) ? -finv[i] : finv[i];
}
return A * B;
}
template<class T, bool is_ntt_friendly>
FormalPowerSeries<T, is_ntt_friendly> bernoulli_number(const int &N)
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
F fac(N + 2, T(0)), finv(N + 2, T(0));
fac[0] = T(1);
for (int i = 1; i <= N + 1; i++)
fac[i] = fac[i - 1] * i;
finv[N + 1] = fac[N + 1].inv();
for (int i = N; i >= 0; i--)
finv[i] = finv[i + 1] * (i + 1);
F f = (finv >> 1).inv();
for (int i = 0; i <= N; i++)
f[i] *= fac[i];
f.pop_back();
return f;
}
// [x^N] P(x)/Q(x) を求める（P の次数は Q の次数より小さい）
template<class T, bool is_ntt_friendly>
T bostan_mori(const FormalPowerSeries<T, is_ntt_friendly> &P, const FormalPowerSeries<T, is_ntt_friendly> &Q, ll N)
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
int d = (int)Q.size() - 1;
assert((int)P.size() <= d);
if (is_ntt_friendly)
{
int z = 1;
while (z < 2 * d + 1)
z <<= 1;
T iz = T(z).inv();
F U = F(P), V = F(Q);
U.resize(z), V.resize(z);
while (N > 0)
{
U.butterfly2(U), V.butterfly2(V);
for (int i = 0; i < z; i += 2)
{
T x = V[i + 1], y = V[i];
U[i] *= x, V[i] *= x;
U[i + 1] *= y, V[i + 1] *= y;
}
U.butterfly_inv2(U), V.butterfly_inv2(V);
for (int i = 0; i < (z >> 1); i++)
{
U[i] = U[2 * i + (N & 1)] * iz;
V[i] = V[2 * i] * iz;
}
for (int i = (z >> 1); i < z; i++)
U[i] = 0, V[i] = 0;
N >>= 1;
}
return U.front() / V.front();
}
else
{
F U = F(P), V = F(Q);
U.resize(d), V.resize(d + 1);
while (N > 0)
{
F U2 = F(U), V2 = F(V), V3 = F(V);
for (int i = 1; i <= d; i += 2)
V3[i] = -V3[i];
U2 *= V3, V2 *= V3;
for (int i = 0; i <= d; i++)
{
U[i] = U2[2 * i + (N & 1)];
V[i] = V2[2 * i];
}
N >>= 1;
}
return U.front() / V.front();
}
}
// a_n = sum[i = 1..d] c_i a_{n-i}（n ≥ d）を満たすとき、a_N を求める（A は 0-indexed で C は 1-indexed）
template<class T, bool is_ntt_friendly>
T linear_recurrence(const vector<T> &A, const vector<T> &C, ll N)
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
int d = C.size();
assert((int)A.size() >= d);
F Ga(d), Q(d + 1);
Q[0] = 1;
for (int i = 0; i < d; i++)
Ga[i] = A[i], Q[i + 1] = -C[i];
F P = Ga * Q;
return bostan_mori(P, Q, N);
}
// https://37zigen.com/multipoint-evaluation/#i-2
template<class T, bool is_ntt_friendly>
T factorial_fast(ll N)
{
using F = FormalPowerSeries<T, is_ntt_friendly>;
if (N >= T::mod())
return 0;
int M = sqrt(N);
vector<F> fs(M);
for (int i = 0; i < M; i++)
fs[i] = {i + 1, 1};
F f = convolution_many(fs);
vector<T> xs(M);
for (int i = 0; i < M; i++)
xs[i] = i * M;
vector<T> ys = f.eval_multipoint(xs);
T res = 1;
for (auto y : ys)
res *= y;
for (int i = M * M + 1; i <= N; i++)
res *= i;
return res;
}
//*
using mint = modint998244353;
const bool ntt = true;
//*/
/*
using mint = modint1000000007;
const bool ntt = false;
//*/
/*
using mint = modint;
const bool ntt = false;
//*/
using fps = FormalPowerSeries<mint, ntt>;
using sfps = SparseFormalPowerSeries<mint, ntt>;
using rfps = RationalFormalPowerSeries<mint, ntt>;
int main()
{
ll N, M;
cin >> N >> M;
fps g(N + 1, 0);
for (ll i = 1; i <= M; i++)
{
for (ll j = 0; i * j <= N; j++)
{
if (1 + i * j <= N)
g.at(1 + i * j) += 1;
if (i + i * j <= N)
g.at(i + i * j) -= 1;
}
}
fps f = (fps{1} - g).inv();
mint ans = mint(M).pow(N) - f.at(N);
cout << ans.val() << endl;
}