QOJ.ac
QOJ
ID | 题目 | 提交者 | 结果 | 用时 | 内存 | 语言 | 文件大小 | 提交时间 | 测评时间 |
---|---|---|---|---|---|---|---|---|---|
#304606 | #8008. Fortune Wheel | ucup-team987# | AC ✓ | 150ms | 5968kb | C++20 | 52.3kb | 2024-01-13 21:43:28 | 2024-10-14 08:00:01 |
Judging History
answer
/**
* date : 2024-01-13 22:43:18
* author : Nyaan
*/
#define NDEBUG
using namespace std;
// intrinstic
#include <immintrin.h>
#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <cctype>
#include <cfenv>
#include <cfloat>
#include <chrono>
#include <cinttypes>
#include <climits>
#include <cmath>
#include <complex>
#include <cstdarg>
#include <cstddef>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <deque>
#include <fstream>
#include <functional>
#include <initializer_list>
#include <iomanip>
#include <ios>
#include <iostream>
#include <istream>
#include <iterator>
#include <limits>
#include <list>
#include <map>
#include <memory>
#include <new>
#include <numeric>
#include <ostream>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <stack>
#include <streambuf>
#include <string>
#include <tuple>
#include <type_traits>
#include <typeinfo>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
// utility
namespace Nyaan {
using ll = long long;
using i64 = long long;
using u64 = unsigned long long;
using i128 = __int128_t;
using u128 = __uint128_t;
template <typename T>
using V = vector<T>;
template <typename T>
using VV = vector<vector<T>>;
using vi = vector<int>;
using vl = vector<long long>;
using vd = V<double>;
using vs = V<string>;
using vvi = vector<vector<int>>;
using vvl = vector<vector<long long>>;
template <typename T>
using minpq = priority_queue<T, vector<T>, greater<T>>;
template <typename T, typename U>
struct P : pair<T, U> {
template <typename... Args>
P(Args... args) : pair<T, U>(args...) {}
using pair<T, U>::first;
using pair<T, U>::second;
P &operator+=(const P &r) {
first += r.first;
second += r.second;
return *this;
}
P &operator-=(const P &r) {
first -= r.first;
second -= r.second;
return *this;
}
P &operator*=(const P &r) {
first *= r.first;
second *= r.second;
return *this;
}
template <typename S>
P &operator*=(const S &r) {
first *= r, second *= r;
return *this;
}
P operator+(const P &r) const { return P(*this) += r; }
P operator-(const P &r) const { return P(*this) -= r; }
P operator*(const P &r) const { return P(*this) *= r; }
template <typename S>
P operator*(const S &r) const {
return P(*this) *= r;
}
P operator-() const { return P{-first, -second}; }
};
using pl = P<ll, ll>;
using pi = P<int, int>;
using vp = V<pl>;
constexpr int inf = 1001001001;
constexpr long long infLL = 4004004004004004004LL;
template <typename T>
int sz(const T &t) {
return t.size();
}
template <typename T, typename U>
inline bool amin(T &x, U y) {
return (y < x) ? (x = y, true) : false;
}
template <typename T, typename U>
inline bool amax(T &x, U y) {
return (x < y) ? (x = y, true) : false;
}
template <typename T>
inline T Max(const vector<T> &v) {
return *max_element(begin(v), end(v));
}
template <typename T>
inline T Min(const vector<T> &v) {
return *min_element(begin(v), end(v));
}
template <typename T>
inline long long Sum(const vector<T> &v) {
return accumulate(begin(v), end(v), 0LL);
}
template <typename T>
int lb(const vector<T> &v, const T &a) {
return lower_bound(begin(v), end(v), a) - begin(v);
}
template <typename T>
int ub(const vector<T> &v, const T &a) {
return upper_bound(begin(v), end(v), a) - begin(v);
}
constexpr long long TEN(int n) {
long long ret = 1, x = 10;
for (; n; x *= x, n >>= 1) ret *= (n & 1 ? x : 1);
return ret;
}
template <typename T, typename U>
pair<T, U> mkp(const T &t, const U &u) {
return make_pair(t, u);
}
template <typename T>
vector<T> mkrui(const vector<T> &v, bool rev = false) {
vector<T> ret(v.size() + 1);
if (rev) {
for (int i = int(v.size()) - 1; i >= 0; i--) ret[i] = v[i] + ret[i + 1];
} else {
for (int i = 0; i < int(v.size()); i++) ret[i + 1] = ret[i] + v[i];
}
return ret;
};
template <typename T>
vector<T> mkuni(const vector<T> &v) {
vector<T> ret(v);
sort(ret.begin(), ret.end());
ret.erase(unique(ret.begin(), ret.end()), ret.end());
return ret;
}
template <typename F>
vector<int> mkord(int N, F f) {
vector<int> ord(N);
iota(begin(ord), end(ord), 0);
sort(begin(ord), end(ord), f);
return ord;
}
template <typename T>
vector<int> mkinv(vector<T> &v) {
int max_val = *max_element(begin(v), end(v));
vector<int> inv(max_val + 1, -1);
for (int i = 0; i < (int)v.size(); i++) inv[v[i]] = i;
return inv;
}
vector<int> mkiota(int n) {
vector<int> ret(n);
iota(begin(ret), end(ret), 0);
return ret;
}
template <typename T>
T mkrev(const T &v) {
T w{v};
reverse(begin(w), end(w));
return w;
}
template <typename T>
bool nxp(T &v) {
return next_permutation(begin(v), end(v));
}
// 返り値の型は入力の T に依存
// i 要素目 : [0, a[i])
template <typename T>
vector<vector<T>> product(const vector<T> &a) {
vector<vector<T>> ret;
vector<T> v;
auto dfs = [&](auto rc, int i) -> void {
if (i == (int)a.size()) {
ret.push_back(v);
return;
}
for (int j = 0; j < a[i]; j++) v.push_back(j), rc(rc, i + 1), v.pop_back();
};
dfs(dfs, 0);
return ret;
}
// F : function(void(T&)), mod を取る操作
// T : 整数型のときはオーバーフローに注意する
template <typename T>
T Power(T a, long long n, const T &I, const function<void(T &)> &f) {
T res = I;
for (; n; f(a = a * a), n >>= 1) {
if (n & 1) f(res = res * a);
}
return res;
}
// T : 整数型のときはオーバーフローに注意する
template <typename T>
T Power(T a, long long n, const T &I = T{1}) {
return Power(a, n, I, function<void(T &)>{[](T &) -> void {}});
}
template <typename T>
T Rev(const T &v) {
T res = v;
reverse(begin(res), end(res));
return res;
}
template <typename T>
vector<T> Transpose(const vector<T> &v) {
using U = typename T::value_type;
int H = v.size(), W = v[0].size();
vector res(W, T(H, U{}));
for (int i = 0; i < H; i++) {
for (int j = 0; j < W; j++) {
res[j][i] = v[i][j];
}
}
return res;
}
template <typename T>
vector<T> Rotate(const vector<T> &v, int clockwise = true) {
using U = typename T::value_type;
int H = v.size(), W = v[0].size();
vector res(W, T(H, U{}));
for (int i = 0; i < H; i++) {
for (int j = 0; j < W; j++) {
if (clockwise) {
res[W - 1 - j][i] = v[i][j];
} else {
res[j][H - 1 - i] = v[i][j];
}
}
}
return res;
}
} // namespace Nyaan
// bit operation
namespace Nyaan {
__attribute__((target("popcnt"))) inline int popcnt(const u64 &a) {
return _mm_popcnt_u64(a);
}
inline int lsb(const u64 &a) { return a ? __builtin_ctzll(a) : 64; }
inline int ctz(const u64 &a) { return a ? __builtin_ctzll(a) : 64; }
inline int msb(const u64 &a) { return a ? 63 - __builtin_clzll(a) : -1; }
template <typename T>
inline int gbit(const T &a, int i) {
return (a >> i) & 1;
}
template <typename T>
inline void sbit(T &a, int i, bool b) {
if (gbit(a, i) != b) a ^= T(1) << i;
}
constexpr long long PW(int n) { return 1LL << n; }
constexpr long long MSK(int n) { return (1LL << n) - 1; }
} // namespace Nyaan
// inout
namespace Nyaan {
template <typename T, typename U>
ostream &operator<<(ostream &os, const pair<T, U> &p) {
os << p.first << " " << p.second;
return os;
}
template <typename T, typename U>
istream &operator>>(istream &is, pair<T, U> &p) {
is >> p.first >> p.second;
return is;
}
template <typename T>
ostream &operator<<(ostream &os, const vector<T> &v) {
int s = (int)v.size();
for (int i = 0; i < s; i++) os << (i ? " " : "") << v[i];
return os;
}
template <typename T>
istream &operator>>(istream &is, vector<T> &v) {
for (auto &x : v) is >> x;
return is;
}
istream &operator>>(istream &is, __int128_t &x) {
string S;
is >> S;
x = 0;
int flag = 0;
for (auto &c : S) {
if (c == '-') {
flag = true;
continue;
}
x *= 10;
x += c - '0';
}
if (flag) x = -x;
return is;
}
istream &operator>>(istream &is, __uint128_t &x) {
string S;
is >> S;
x = 0;
for (auto &c : S) {
x *= 10;
x += c - '0';
}
return is;
}
ostream &operator<<(ostream &os, __int128_t x) {
if (x == 0) return os << 0;
if (x < 0) os << '-', x = -x;
string S;
while (x) S.push_back('0' + x % 10), x /= 10;
reverse(begin(S), end(S));
return os << S;
}
ostream &operator<<(ostream &os, __uint128_t x) {
if (x == 0) return os << 0;
string S;
while (x) S.push_back('0' + x % 10), x /= 10;
reverse(begin(S), end(S));
return os << S;
}
void in() {}
template <typename T, class... U>
void in(T &t, U &...u) {
cin >> t;
in(u...);
}
void out() { cout << "\n"; }
template <typename T, class... U, char sep = ' '>
void out(const T &t, const U &...u) {
cout << t;
if (sizeof...(u)) cout << sep;
out(u...);
}
struct IoSetupNya {
IoSetupNya() {
cin.tie(nullptr);
ios::sync_with_stdio(false);
cout << fixed << setprecision(15);
cerr << fixed << setprecision(7);
}
} iosetupnya;
} // namespace Nyaan
// debug
#ifdef NyaanDebug
#define trc(...) (void(0))
#else
#define trc(...) (void(0))
#endif
#ifdef NyaanLocal
#define trc2(...) (void(0))
#else
#define trc2(...) (void(0))
#endif
// macro
#define each(x, v) for (auto&& x : v)
#define each2(x, y, v) for (auto&& [x, y] : v)
#define all(v) (v).begin(), (v).end()
#define rep(i, N) for (long long i = 0; i < (long long)(N); i++)
#define repr(i, N) for (long long i = (long long)(N)-1; i >= 0; i--)
#define rep1(i, N) for (long long i = 1; i <= (long long)(N); i++)
#define repr1(i, N) for (long long i = (N); (long long)(i) > 0; i--)
#define reg(i, a, b) for (long long i = (a); i < (b); i++)
#define regr(i, a, b) for (long long i = (b)-1; i >= (a); i--)
#define fi first
#define se second
#define ini(...) \
int __VA_ARGS__; \
in(__VA_ARGS__)
#define inl(...) \
long long __VA_ARGS__; \
in(__VA_ARGS__)
#define ins(...) \
string __VA_ARGS__; \
in(__VA_ARGS__)
#define in2(s, t) \
for (int i = 0; i < (int)s.size(); i++) { \
in(s[i], t[i]); \
}
#define in3(s, t, u) \
for (int i = 0; i < (int)s.size(); i++) { \
in(s[i], t[i], u[i]); \
}
#define in4(s, t, u, v) \
for (int i = 0; i < (int)s.size(); i++) { \
in(s[i], t[i], u[i], v[i]); \
}
#define die(...) \
do { \
Nyaan::out(__VA_ARGS__); \
return; \
} while (0)
namespace Nyaan {
void solve();
}
int main() { Nyaan::solve(); }
//
using namespace std;
using namespace std;
namespace BinaryGCDImpl {
using u64 = unsigned long long;
using i8 = char;
u64 binary_gcd(u64 a, u64 b) {
if (a == 0 || b == 0) return a + b;
i8 n = __builtin_ctzll(a);
i8 m = __builtin_ctzll(b);
a >>= n;
b >>= m;
n = min(n, m);
while (a != b) {
u64 d = a - b;
i8 s = __builtin_ctzll(d);
bool f = a > b;
b = f ? b : a;
a = (f ? d : -d) >> s;
}
return a << n;
}
using u128 = __uint128_t;
// a > 0
int ctz128(u128 a) {
u64 lo = a & u64(-1);
return lo ? __builtin_ctzll(lo) : 64 + __builtin_ctzll(a >> 64);
}
u128 binary_gcd128(u128 a, u128 b) {
if (a == 0 || b == 0) return a + b;
i8 n = ctz128(a);
i8 m = ctz128(b);
a >>= n;
b >>= m;
n = min(n, m);
while (a != b) {
u128 d = a - b;
i8 s = ctz128(d);
bool f = a > b;
b = f ? b : a;
a = (f ? d : -d) >> s;
}
return a << n;
}
} // namespace BinaryGCDImpl
long long binary_gcd(long long a, long long b) {
return BinaryGCDImpl::binary_gcd(abs(a), abs(b));
}
__int128_t binary_gcd128(__int128_t a, __int128_t b) {
if (a < 0) a = -a;
if (b < 0) b = -b;
return BinaryGCDImpl::binary_gcd128(a, b);
}
/**
* @brief binary GCD
*/
using namespace std;
using namespace std;
namespace internal {
template <typename T>
using is_broadly_integral =
typename conditional_t<is_integral_v<T> || is_same_v<T, __int128_t> ||
is_same_v<T, __uint128_t>,
true_type, false_type>::type;
template <typename T>
using is_broadly_signed =
typename conditional_t<is_signed_v<T> || is_same_v<T, __int128_t>,
true_type, false_type>::type;
template <typename T>
using is_broadly_unsigned =
typename conditional_t<is_unsigned_v<T> || is_same_v<T, __uint128_t>,
true_type, false_type>::type;
#define ENABLE_VALUE(x) \
template <typename T> \
constexpr bool x##_v = x<T>::value;
ENABLE_VALUE(is_broadly_integral);
ENABLE_VALUE(is_broadly_signed);
ENABLE_VALUE(is_broadly_unsigned);
#undef ENABLE_VALUE
#define ENABLE_HAS_TYPE(var) \
template <class, class = void> \
struct has_##var : false_type {}; \
template <class T> \
struct has_##var<T, void_t<typename T::var>> : true_type {}; \
template <class T> \
constexpr auto has_##var##_v = has_##var<T>::value;
#define ENABLE_HAS_VAR(var) \
template <class, class = void> \
struct has_##var : false_type {}; \
template <class T> \
struct has_##var<T, void_t<decltype(T::var)>> : true_type {}; \
template <class T> \
constexpr auto has_##var##_v = has_##var<T>::value;
} // namespace internal
template <uint32_t mod>
struct LazyMontgomeryModInt {
using mint = LazyMontgomeryModInt;
using i32 = int32_t;
using u32 = uint32_t;
using u64 = uint64_t;
static constexpr u32 get_r() {
u32 ret = mod;
for (i32 i = 0; i < 4; ++i) ret *= 2 - mod * ret;
return ret;
}
static constexpr u32 r = get_r();
static constexpr u32 n2 = -u64(mod) % mod;
static_assert(mod < (1 << 30), "invalid, mod >= 2 ^ 30");
static_assert((mod & 1) == 1, "invalid, mod % 2 == 0");
static_assert(r * mod == 1, "this code has bugs.");
u32 a;
constexpr LazyMontgomeryModInt() : a(0) {}
constexpr LazyMontgomeryModInt(const int64_t &b)
: a(reduce(u64(b % mod + mod) * n2)){};
static constexpr u32 reduce(const u64 &b) {
return (b + u64(u32(b) * u32(-r)) * mod) >> 32;
}
constexpr mint &operator+=(const mint &b) {
if (i32(a += b.a - 2 * mod) < 0) a += 2 * mod;
return *this;
}
constexpr mint &operator-=(const mint &b) {
if (i32(a -= b.a) < 0) a += 2 * mod;
return *this;
}
constexpr mint &operator*=(const mint &b) {
a = reduce(u64(a) * b.a);
return *this;
}
constexpr mint &operator/=(const mint &b) {
*this *= b.inverse();
return *this;
}
constexpr mint operator+(const mint &b) const { return mint(*this) += b; }
constexpr mint operator-(const mint &b) const { return mint(*this) -= b; }
constexpr mint operator*(const mint &b) const { return mint(*this) *= b; }
constexpr mint operator/(const mint &b) const { return mint(*this) /= b; }
constexpr bool operator==(const mint &b) const {
return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a);
}
constexpr bool operator!=(const mint &b) const {
return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a);
}
constexpr mint operator-() const { return mint() - mint(*this); }
constexpr mint operator+() const { return mint(*this); }
constexpr mint pow(u64 n) const {
mint ret(1), mul(*this);
while (n > 0) {
if (n & 1) ret *= mul;
mul *= mul;
n >>= 1;
}
return ret;
}
constexpr mint inverse() const {
int x = get(), y = mod, u = 1, v = 0, t = 0, tmp = 0;
while (y > 0) {
t = x / y;
x -= t * y, u -= t * v;
tmp = x, x = y, y = tmp;
tmp = u, u = v, v = tmp;
}
return mint{u};
}
friend ostream &operator<<(ostream &os, const mint &b) {
return os << b.get();
}
friend istream &operator>>(istream &is, mint &b) {
int64_t t;
is >> t;
b = LazyMontgomeryModInt<mod>(t);
return (is);
}
constexpr u32 get() const {
u32 ret = reduce(a);
return ret >= mod ? ret - mod : ret;
}
static constexpr u32 get_mod() { return mod; }
};
template <typename mint>
struct NTT {
static constexpr uint32_t get_pr() {
uint32_t _mod = mint::get_mod();
using u64 = uint64_t;
u64 ds[32] = {};
int idx = 0;
u64 m = _mod - 1;
for (u64 i = 2; i * i <= m; ++i) {
if (m % i == 0) {
ds[idx++] = i;
while (m % i == 0) m /= i;
}
}
if (m != 1) ds[idx++] = m;
uint32_t _pr = 2;
while (1) {
int flg = 1;
for (int i = 0; i < idx; ++i) {
u64 a = _pr, b = (_mod - 1) / ds[i], r = 1;
while (b) {
if (b & 1) r = r * a % _mod;
a = a * a % _mod;
b >>= 1;
}
if (r == 1) {
flg = 0;
break;
}
}
if (flg == 1) break;
++_pr;
}
return _pr;
};
static constexpr uint32_t mod = mint::get_mod();
static constexpr uint32_t pr = get_pr();
static constexpr int level = __builtin_ctzll(mod - 1);
mint dw[level], dy[level];
void setwy(int k) {
mint w[level], y[level];
w[k - 1] = mint(pr).pow((mod - 1) / (1 << k));
y[k - 1] = w[k - 1].inverse();
for (int i = k - 2; i > 0; --i)
w[i] = w[i + 1] * w[i + 1], y[i] = y[i + 1] * y[i + 1];
dw[1] = w[1], dy[1] = y[1], dw[2] = w[2], dy[2] = y[2];
for (int i = 3; i < k; ++i) {
dw[i] = dw[i - 1] * y[i - 2] * w[i];
dy[i] = dy[i - 1] * w[i - 2] * y[i];
}
}
NTT() { setwy(level); }
void fft4(vector<mint> &a, int k) {
if ((int)a.size() <= 1) return;
if (k == 1) {
mint a1 = a[1];
a[1] = a[0] - a[1];
a[0] = a[0] + a1;
return;
}
if (k & 1) {
int v = 1 << (k - 1);
for (int j = 0; j < v; ++j) {
mint ajv = a[j + v];
a[j + v] = a[j] - ajv;
a[j] += ajv;
}
}
int u = 1 << (2 + (k & 1));
int v = 1 << (k - 2 - (k & 1));
mint one = mint(1);
mint imag = dw[1];
while (v) {
// jh = 0
{
int j0 = 0;
int j1 = v;
int j2 = j1 + v;
int j3 = j2 + v;
for (; j0 < v; ++j0, ++j1, ++j2, ++j3) {
mint t0 = a[j0], t1 = a[j1], t2 = a[j2], t3 = a[j3];
mint t0p2 = t0 + t2, t1p3 = t1 + t3;
mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag;
a[j0] = t0p2 + t1p3, a[j1] = t0p2 - t1p3;
a[j2] = t0m2 + t1m3, a[j3] = t0m2 - t1m3;
}
}
// jh >= 1
mint ww = one, xx = one * dw[2], wx = one;
for (int jh = 4; jh < u;) {
ww = xx * xx, wx = ww * xx;
int j0 = jh * v;
int je = j0 + v;
int j2 = je + v;
for (; j0 < je; ++j0, ++j2) {
mint t0 = a[j0], t1 = a[j0 + v] * xx, t2 = a[j2] * ww,
t3 = a[j2 + v] * wx;
mint t0p2 = t0 + t2, t1p3 = t1 + t3;
mint t0m2 = t0 - t2, t1m3 = (t1 - t3) * imag;
a[j0] = t0p2 + t1p3, a[j0 + v] = t0p2 - t1p3;
a[j2] = t0m2 + t1m3, a[j2 + v] = t0m2 - t1m3;
}
xx *= dw[__builtin_ctzll((jh += 4))];
}
u <<= 2;
v >>= 2;
}
}
void ifft4(vector<mint> &a, int k) {
if ((int)a.size() <= 1) return;
if (k == 1) {
mint a1 = a[1];
a[1] = a[0] - a[1];
a[0] = a[0] + a1;
return;
}
int u = 1 << (k - 2);
int v = 1;
mint one = mint(1);
mint imag = dy[1];
while (u) {
// jh = 0
{
int j0 = 0;
int j1 = v;
int j2 = v + v;
int j3 = j2 + v;
for (; j0 < v; ++j0, ++j1, ++j2, ++j3) {
mint t0 = a[j0], t1 = a[j1], t2 = a[j2], t3 = a[j3];
mint t0p1 = t0 + t1, t2p3 = t2 + t3;
mint t0m1 = t0 - t1, t2m3 = (t2 - t3) * imag;
a[j0] = t0p1 + t2p3, a[j2] = t0p1 - t2p3;
a[j1] = t0m1 + t2m3, a[j3] = t0m1 - t2m3;
}
}
// jh >= 1
mint ww = one, xx = one * dy[2], yy = one;
u <<= 2;
for (int jh = 4; jh < u;) {
ww = xx * xx, yy = xx * imag;
int j0 = jh * v;
int je = j0 + v;
int j2 = je + v;
for (; j0 < je; ++j0, ++j2) {
mint t0 = a[j0], t1 = a[j0 + v], t2 = a[j2], t3 = a[j2 + v];
mint t0p1 = t0 + t1, t2p3 = t2 + t3;
mint t0m1 = (t0 - t1) * xx, t2m3 = (t2 - t3) * yy;
a[j0] = t0p1 + t2p3, a[j2] = (t0p1 - t2p3) * ww;
a[j0 + v] = t0m1 + t2m3, a[j2 + v] = (t0m1 - t2m3) * ww;
}
xx *= dy[__builtin_ctzll(jh += 4)];
}
u >>= 4;
v <<= 2;
}
if (k & 1) {
u = 1 << (k - 1);
for (int j = 0; j < u; ++j) {
mint ajv = a[j] - a[j + u];
a[j] += a[j + u];
a[j + u] = ajv;
}
}
}
void ntt(vector<mint> &a) {
if ((int)a.size() <= 1) return;
fft4(a, __builtin_ctz(a.size()));
}
void intt(vector<mint> &a) {
if ((int)a.size() <= 1) return;
ifft4(a, __builtin_ctz(a.size()));
mint iv = mint(a.size()).inverse();
for (auto &x : a) x *= iv;
}
vector<mint> multiply(const vector<mint> &a, const vector<mint> &b) {
int l = a.size() + b.size() - 1;
if (min<int>(a.size(), b.size()) <= 40) {
vector<mint> s(l);
for (int i = 0; i < (int)a.size(); ++i)
for (int j = 0; j < (int)b.size(); ++j) s[i + j] += a[i] * b[j];
return s;
}
int k = 2, M = 4;
while (M < l) M <<= 1, ++k;
setwy(k);
vector<mint> s(M);
for (int i = 0; i < (int)a.size(); ++i) s[i] = a[i];
fft4(s, k);
if (a.size() == b.size() && a == b) {
for (int i = 0; i < M; ++i) s[i] *= s[i];
} else {
vector<mint> t(M);
for (int i = 0; i < (int)b.size(); ++i) t[i] = b[i];
fft4(t, k);
for (int i = 0; i < M; ++i) s[i] *= t[i];
}
ifft4(s, k);
s.resize(l);
mint invm = mint(M).inverse();
for (int i = 0; i < l; ++i) s[i] *= invm;
return s;
}
void ntt_doubling(vector<mint> &a) {
int M = (int)a.size();
auto b = a;
intt(b);
mint r = 1, zeta = mint(pr).pow((mint::get_mod() - 1) / (M << 1));
for (int i = 0; i < M; i++) b[i] *= r, r *= zeta;
ntt(b);
copy(begin(b), end(b), back_inserter(a));
}
};
namespace ArbitraryNTT {
using i64 = int64_t;
using u128 = __uint128_t;
constexpr int32_t m0 = 167772161;
constexpr int32_t m1 = 469762049;
constexpr int32_t m2 = 754974721;
using mint0 = LazyMontgomeryModInt<m0>;
using mint1 = LazyMontgomeryModInt<m1>;
using mint2 = LazyMontgomeryModInt<m2>;
constexpr int r01 = mint1(m0).inverse().get();
constexpr int r02 = mint2(m0).inverse().get();
constexpr int r12 = mint2(m1).inverse().get();
constexpr int r02r12 = i64(r02) * r12 % m2;
constexpr i64 w1 = m0;
constexpr i64 w2 = i64(m0) * m1;
template <typename T, typename submint>
vector<submint> mul(const vector<T> &a, const vector<T> &b) {
static NTT<submint> ntt;
vector<submint> s(a.size()), t(b.size());
for (int i = 0; i < (int)a.size(); ++i) s[i] = i64(a[i] % submint::get_mod());
for (int i = 0; i < (int)b.size(); ++i) t[i] = i64(b[i] % submint::get_mod());
return ntt.multiply(s, t);
}
template <typename T>
vector<int> multiply(const vector<T> &s, const vector<T> &t, int mod) {
auto d0 = mul<T, mint0>(s, t);
auto d1 = mul<T, mint1>(s, t);
auto d2 = mul<T, mint2>(s, t);
int n = d0.size();
vector<int> ret(n);
const int W1 = w1 % mod;
const int W2 = w2 % mod;
for (int i = 0; i < n; i++) {
int n1 = d1[i].get(), n2 = d2[i].get(), a = d0[i].get();
int b = i64(n1 + m1 - a) * r01 % m1;
int c = (i64(n2 + m2 - a) * r02r12 + i64(m2 - b) * r12) % m2;
ret[i] = (i64(a) + i64(b) * W1 + i64(c) * W2) % mod;
}
return ret;
}
template <typename mint>
vector<mint> multiply(const vector<mint> &a, const vector<mint> &b) {
if (a.size() == 0 && b.size() == 0) return {};
if (min<int>(a.size(), b.size()) < 128) {
vector<mint> ret(a.size() + b.size() - 1);
for (int i = 0; i < (int)a.size(); ++i)
for (int j = 0; j < (int)b.size(); ++j) ret[i + j] += a[i] * b[j];
return ret;
}
vector<int> s(a.size()), t(b.size());
for (int i = 0; i < (int)a.size(); ++i) s[i] = a[i].get();
for (int i = 0; i < (int)b.size(); ++i) t[i] = b[i].get();
vector<int> u = multiply<int>(s, t, mint::get_mod());
vector<mint> ret(u.size());
for (int i = 0; i < (int)u.size(); ++i) ret[i] = mint(u[i]);
return ret;
}
template <typename T>
vector<u128> multiply_u128(const vector<T> &s, const vector<T> &t) {
if (s.size() == 0 && t.size() == 0) return {};
if (min<int>(s.size(), t.size()) < 128) {
vector<u128> ret(s.size() + t.size() - 1);
for (int i = 0; i < (int)s.size(); ++i)
for (int j = 0; j < (int)t.size(); ++j) ret[i + j] += i64(s[i]) * t[j];
return ret;
}
auto d0 = mul<T, mint0>(s, t);
auto d1 = mul<T, mint1>(s, t);
auto d2 = mul<T, mint2>(s, t);
int n = d0.size();
vector<u128> ret(n);
for (int i = 0; i < n; i++) {
i64 n1 = d1[i].get(), n2 = d2[i].get();
i64 a = d0[i].get();
i64 b = (n1 + m1 - a) * r01 % m1;
i64 c = ((n2 + m2 - a) * r02r12 + (m2 - b) * r12) % m2;
ret[i] = a + b * w1 + u128(c) * w2;
}
return ret;
}
} // namespace ArbitraryNTT
namespace MultiPrecisionIntegerImpl {
struct TENS {
static constexpr int offset = 30;
constexpr TENS() : _tend() {
_tend[offset] = 1;
for (int i = 1; i <= offset; i++) {
_tend[offset + i] = _tend[offset + i - 1] * 10.0;
_tend[offset - i] = 1.0 / _tend[offset + i];
}
}
long double ten_ld(int n) const {
assert(-offset <= n and n <= offset);
return _tend[n + offset];
}
private:
long double _tend[offset * 2 + 1];
};
} // namespace MultiPrecisionIntegerImpl
// 0 は neg=false, dat={} として扱う
struct MultiPrecisionInteger {
using M = MultiPrecisionInteger;
inline constexpr static MultiPrecisionIntegerImpl::TENS tens = {};
static constexpr int D = 1000000000;
static constexpr int logD = 9;
bool neg;
vector<int> dat;
MultiPrecisionInteger() : neg(false), dat() {}
MultiPrecisionInteger(bool n, const vector<int>& d) : neg(n), dat(d) {}
template <typename I,
enable_if_t<internal::is_broadly_integral_v<I>>* = nullptr>
MultiPrecisionInteger(I x) : neg(false) {
if constexpr (internal::is_broadly_signed_v<I>) {
if (x < 0) neg = true, x = -x;
}
while (x) dat.push_back(x % D), x /= D;
}
MultiPrecisionInteger(const string& S) : neg(false) {
assert(!S.empty());
if (S.size() == 1u && S[0] == '0') return;
int l = 0;
if (S[0] == '-') ++l, neg = true;
for (int ie = S.size(); l < ie; ie -= logD) {
int is = max(l, ie - logD);
long long x = 0;
for (int i = is; i < ie; i++) x = x * 10 + S[i] - '0';
dat.push_back(x);
}
}
friend M operator+(const M& lhs, const M& rhs) {
if (lhs.neg == rhs.neg) return {lhs.neg, _add(lhs.dat, rhs.dat)};
if (_leq(lhs.dat, rhs.dat)) {
// |l| <= |r|
auto c = _sub(rhs.dat, lhs.dat);
bool n = _is_zero(c) ? false : rhs.neg;
return {n, c};
}
auto c = _sub(lhs.dat, rhs.dat);
bool n = _is_zero(c) ? false : lhs.neg;
return {n, c};
}
friend M operator-(const M& lhs, const M& rhs) { return lhs + (-rhs); }
friend M operator*(const M& lhs, const M& rhs) {
auto c = _mul(lhs.dat, rhs.dat);
bool n = _is_zero(c) ? false : (lhs.neg ^ rhs.neg);
return {n, c};
}
friend pair<M, M> divmod(const M& lhs, const M& rhs) {
auto dm = _divmod_newton(lhs.dat, rhs.dat);
bool dn = _is_zero(dm.first) ? false : lhs.neg != rhs.neg;
bool mn = _is_zero(dm.second) ? false : lhs.neg;
return {M{dn, dm.first}, M{mn, dm.second}};
}
friend M operator/(const M& lhs, const M& rhs) {
return divmod(lhs, rhs).first;
}
friend M operator%(const M& lhs, const M& rhs) {
return divmod(lhs, rhs).second;
}
M& operator+=(const M& rhs) { return (*this) = (*this) + rhs; }
M& operator-=(const M& rhs) { return (*this) = (*this) - rhs; }
M& operator*=(const M& rhs) { return (*this) = (*this) * rhs; }
M& operator/=(const M& rhs) { return (*this) = (*this) / rhs; }
M& operator%=(const M& rhs) { return (*this) = (*this) % rhs; }
M operator-() const {
if (is_zero()) return *this;
return {!neg, dat};
}
M operator+() const { return *this; }
friend M abs(const M& m) { return {false, m.dat}; }
bool is_zero() const { return _is_zero(dat); }
friend bool operator==(const M& lhs, const M& rhs) {
return lhs.neg == rhs.neg && lhs.dat == rhs.dat;
}
friend bool operator!=(const M& lhs, const M& rhs) {
return lhs.neg != rhs.neg || lhs.dat != rhs.dat;
}
friend bool operator<(const M& lhs, const M& rhs) {
if (lhs == rhs) return false;
return _neq_lt(lhs, rhs);
}
friend bool operator<=(const M& lhs, const M& rhs) {
if (lhs == rhs) return true;
return _neq_lt(lhs, rhs);
}
friend bool operator>(const M& lhs, const M& rhs) {
if (lhs == rhs) return false;
return _neq_lt(rhs, lhs);
}
friend bool operator>=(const M& lhs, const M& rhs) {
if (lhs == rhs) return true;
return _neq_lt(rhs, lhs);
}
// a * 10^b (1 <= |a| < 10) の形で渡す
// 相対誤差:10^{-16} ~ 10^{-19} 程度 (処理系依存)
pair<long double, int> dfp() const {
if (is_zero()) return {0, 0};
int l = max<int>(0, _size() - 3);
int b = logD * l;
string prefix{};
for (int i = _size() - 1; i >= l; i--) {
prefix += _itos(dat[i], i != _size() - 1);
}
b += prefix.size() - 1;
long double a = 0;
for (auto& c : prefix) a = a * 10.0 + (c - '0');
a *= tens.ten_ld(-((int)prefix.size()) + 1);
a = clamp<long double>(a, 1.0, nextafterl(10.0, 1.0));
if (neg) a = -a;
return {a, b};
}
string to_string() const {
if (is_zero()) return "0";
string res;
if (neg) res.push_back('-');
for (int i = _size() - 1; i >= 0; i--) {
res += _itos(dat[i], i != _size() - 1);
}
return res;
}
long double to_ld() const {
auto [a, b] = dfp();
if (-tens.offset <= b and b <= tens.offset) {
return a * tens.ten_ld(b);
}
return a * powl(10, b);
}
long long to_ll() const {
long long res = _to_ll(dat);
return neg ? -res : res;
}
__int128_t to_i128() const {
__int128_t res = _to_i128(dat);
return neg ? -res : res;
}
friend istream& operator>>(istream& is, M& m) {
string s;
is >> s;
m = M{s};
return is;
}
friend ostream& operator<<(ostream& os, const M& m) {
return os << m.to_string();
}
// 内部の関数をテスト
static void _test_private_function(const M&, const M&);
private:
// size
int _size() const { return dat.size(); }
// a == b
static bool _eq(const vector<int>& a, const vector<int>& b) { return a == b; }
// a < b
static bool _lt(const vector<int>& a, const vector<int>& b) {
if (a.size() != b.size()) return a.size() < b.size();
for (int i = a.size() - 1; i >= 0; i--) {
if (a[i] != b[i]) return a[i] < b[i];
}
return false;
}
// a <= b
static bool _leq(const vector<int>& a, const vector<int>& b) {
return _eq(a, b) || _lt(a, b);
}
// a < b (s.t. a != b)
static bool _neq_lt(const M& lhs, const M& rhs) {
assert(lhs != rhs);
if (lhs.neg != rhs.neg) return lhs.neg;
bool f = _lt(lhs.dat, rhs.dat);
if (f) return !lhs.neg;
return lhs.neg;
}
// a == 0
static bool _is_zero(const vector<int>& a) { return a.empty(); }
// a == 1
static bool _is_one(const vector<int>& a) {
return (int)a.size() == 1 && a[0] == 1;
}
// 末尾 0 を削除
static void _shrink(vector<int>& a) {
while (a.size() && a.back() == 0) a.pop_back();
}
// 末尾 0 を削除
void _shrink() {
while (_size() && dat.back() == 0) dat.pop_back();
}
// a + b
static vector<int> _add(const vector<int>& a, const vector<int>& b) {
vector<int> c(max(a.size(), b.size()) + 1);
for (int i = 0; i < (int)a.size(); i++) c[i] += a[i];
for (int i = 0; i < (int)b.size(); i++) c[i] += b[i];
for (int i = 0; i < (int)c.size() - 1; i++) {
if (c[i] >= D) c[i] -= D, c[i + 1]++;
}
_shrink(c);
return c;
}
// a - b
static vector<int> _sub(const vector<int>& a, const vector<int>& b) {
assert(_leq(b, a));
vector<int> c{a};
int borrow = 0;
for (int i = 0; i < (int)a.size(); i++) {
if (i < (int)b.size()) borrow += b[i];
c[i] -= borrow;
borrow = 0;
if (c[i] < 0) c[i] += D, borrow = 1;
}
assert(borrow == 0);
_shrink(c);
return c;
}
// a * b (fft)
static vector<int> _mul_fft(const vector<int>& a, const vector<int>& b) {
if (a.empty() || b.empty()) return {};
auto m = ArbitraryNTT::multiply_u128(a, b);
vector<int> c;
c.reserve(m.size() + 3);
__uint128_t x = 0;
for (int i = 0;; i++) {
if (i >= (int)m.size() && x == 0) break;
if (i < (int)m.size()) x += m[i];
c.push_back(x % D);
x /= D;
}
_shrink(c);
return c;
}
// a * b (naive)
static vector<int> _mul_naive(const vector<int>& a, const vector<int>& b) {
if (a.empty() || b.empty()) return {};
vector<long long> prod(a.size() + b.size() - 1 + 1);
for (int i = 0; i < (int)a.size(); i++) {
for (int j = 0; j < (int)b.size(); j++) {
long long p = 1LL * a[i] * b[j];
prod[i + j] += p;
if (prod[i + j] >= (4LL * D * D)) {
prod[i + j] -= 4LL * D * D;
prod[i + j + 1] += 4LL * D;
}
}
}
vector<int> c(prod.size() + 1);
long long x = 0;
int i = 0;
for (; i < (int)prod.size(); i++) x += prod[i], c[i] = x % D, x /= D;
while (x) c[i] = x % D, x /= D, i++;
_shrink(c);
return c;
}
// a * b
static vector<int> _mul(const vector<int>& a, const vector<int>& b) {
if (_is_zero(a) || _is_zero(b)) return {};
if (_is_one(a)) return b;
if (_is_one(b)) return a;
if (min<int>(a.size(), b.size()) <= 128) {
return a.size() < b.size() ? _mul_naive(b, a) : _mul_naive(a, b);
}
return _mul_fft(a, b);
}
// 0 <= A < 1e18, 1 <= B < 1e9
static pair<vector<int>, vector<int>> _divmod_li(const vector<int>& a,
const vector<int>& b) {
assert(0 <= (int)a.size() && (int)a.size() <= 2);
assert((int)b.size() == 1);
long long va = _to_ll(a);
int vb = b[0];
return {_integer_to_vec(va / vb), _integer_to_vec(va % vb)};
}
// 0 <= A < 1e18, 1 <= B < 1e18
static pair<vector<int>, vector<int>> _divmod_ll(const vector<int>& a,
const vector<int>& b) {
assert(0 <= (int)a.size() && (int)a.size() <= 2);
assert(1 <= (int)b.size() && (int)b.size() <= 2);
long long va = _to_ll(a), vb = _to_ll(b);
return {_integer_to_vec(va / vb), _integer_to_vec(va % vb)};
}
// 1 <= B < 1e9
static pair<vector<int>, vector<int>> _divmod_1e9(const vector<int>& a,
const vector<int>& b) {
assert((int)b.size() == 1);
if (b[0] == 1) return {a, {}};
if ((int)a.size() <= 2) return _divmod_li(a, b);
vector<int> quo(a.size());
long long d = 0;
int b0 = b[0];
for (int i = a.size() - 1; i >= 0; i--) {
d = d * D + a[i];
assert(d < 1LL * D * b0);
int q = d / b0, r = d % b0;
quo[i] = q, d = r;
}
_shrink(quo);
return {quo, d ? vector<int>{int(d)} : vector<int>{}};
}
// 0 <= A, 1 <= B
static pair<vector<int>, vector<int>> _divmod_naive(const vector<int>& a,
const vector<int>& b) {
if (_is_zero(b)) {
cerr << "Divide by Zero Exception" << endl;
exit(1);
}
assert(1 <= (int)b.size());
if ((int)b.size() == 1) return _divmod_1e9(a, b);
if (max<int>(a.size(), b.size()) <= 2) return _divmod_ll(a, b);
if (_lt(a, b)) return {{}, a};
// B >= 1e9, A >= B
int norm = D / (b.back() + 1);
vector<int> x = _mul(a, {norm});
vector<int> y = _mul(b, {norm});
int yb = y.back();
vector<int> quo(x.size() - y.size() + 1);
vector<int> rem(x.end() - y.size(), x.end());
for (int i = quo.size() - 1; i >= 0; i--) {
if (rem.size() < y.size()) {
// do nothing
} else if (rem.size() == y.size()) {
if (_leq(y, rem)) {
quo[i] = 1, rem = _sub(rem, y);
}
} else {
assert(y.size() + 1 == rem.size());
long long rb = 1LL * rem[rem.size() - 1] * D + rem[rem.size() - 2];
int q = rb / yb;
vector<int> yq = _mul(y, {q});
// 真の商は q-2 以上 q+1 以下だが自信が無いので念のため while を回す
while (_lt(rem, yq)) q--, yq = _sub(yq, y);
rem = _sub(rem, yq);
while (_leq(y, rem)) q++, rem = _sub(rem, y);
quo[i] = q;
}
if (i) rem.insert(begin(rem), x[i - 1]);
}
_shrink(quo), _shrink(rem);
auto [q2, r2] = _divmod_1e9(rem, {norm});
assert(_is_zero(r2));
return {quo, q2};
}
// 0 <= A, 1 <= B
static pair<vector<int>, vector<int>> _divmod_dc(const vector<int>& a,
const vector<int>& b);
// 1 / a を 絶対誤差 B^{-deg} で求める
static vector<int> _calc_inv(const vector<int>& a, int deg) {
assert(!a.empty() && D / 2 <= a.back() and a.back() < D);
int k = deg, c = a.size();
while (k > 64) k = (k + 1) / 2;
vector<int> z(c + k + 1);
z.back() = 1;
z = _divmod_naive(z, a).first;
while (k < deg) {
vector<int> s = _mul(z, z);
s.insert(begin(s), 0);
int d = min(c, 2 * k + 1);
vector<int> t{end(a) - d, end(a)}, u = _mul(s, t);
u.erase(begin(u), begin(u) + d);
vector<int> w(k + 1), w2 = _add(z, z);
copy(begin(w2), end(w2), back_inserter(w));
z = _sub(w, u);
z.erase(begin(z));
k *= 2;
}
z.erase(begin(z), begin(z) + k - deg);
return z;
}
static pair<vector<int>, vector<int>> _divmod_newton(const vector<int>& a,
const vector<int>& b) {
if (_is_zero(b)) {
cerr << "Divide by Zero Exception" << endl;
exit(1);
}
if ((int)b.size() <= 64) return _divmod_naive(a, b);
if ((int)a.size() - (int)b.size() <= 64) return _divmod_naive(a, b);
int norm = D / (b.back() + 1);
vector<int> x = _mul(a, {norm});
vector<int> y = _mul(b, {norm});
int s = x.size(), t = y.size();
int deg = s - t + 2;
vector<int> z = _calc_inv(y, deg);
vector<int> q = _mul(x, z);
q.erase(begin(q), begin(q) + t + deg);
vector<int> yq = _mul(y, {q});
while (_lt(x, yq)) q = _sub(q, {1}), yq = _sub(yq, y);
vector<int> r = _sub(x, yq);
while (_leq(y, r)) q = _add(q, {1}), r = _sub(r, y);
_shrink(q), _shrink(r);
auto [q2, r2] = _divmod_1e9(r, {norm});
assert(_is_zero(r2));
return {q, q2};
}
// int -> string
// 先頭かどうかに応じて zero padding するかを決める
static string _itos(int x, bool zero_padding) {
assert(0 <= x && x < D);
string res;
for (int i = 0; i < logD; i++) {
res.push_back('0' + x % 10), x /= 10;
}
if (!zero_padding) {
while (res.size() && res.back() == '0') res.pop_back();
assert(!res.empty());
}
reverse(begin(res), end(res));
return res;
}
// convert ll to vec
template <typename I,
enable_if_t<internal::is_broadly_integral_v<I>>* = nullptr>
static vector<int> _integer_to_vec(I x) {
if constexpr (internal::is_broadly_signed_v<I>) {
assert(x >= 0);
}
vector<int> res;
while (x) res.push_back(x % D), x /= D;
return res;
}
static long long _to_ll(const vector<int>& a) {
long long res = 0;
for (int i = (int)a.size() - 1; i >= 0; i--) res = res * D + a[i];
return res;
}
static __int128_t _to_i128(const vector<int>& a) {
__int128_t res = 0;
for (int i = (int)a.size() - 1; i >= 0; i--) res = res * D + a[i];
return res;
}
static void _dump(const vector<int>& a, string s = "") {
if (!s.empty()) cerr << s << " : ";
cerr << "{ ";
for (int i = 0; i < (int)a.size(); i++) cerr << a[i] << ", ";
cerr << "}" << endl;
}
};
using bigint = MultiPrecisionInteger;
/**
* @brief 多倍長整数
*/
namespace GCDforBigintImpl {
bigint bigint_pow(bigint a, long long k) {
if (k == 0) return 1;
bigint half = bigint_pow(a, k / 2);
bigint res = half * half;
return k % 2 ? res * a : res;
}
// a = 2^x 5^y の形で表現する
pair<int, int> shrink(int a) {
assert(a > 0);
int x = __builtin_ctz(a);
a >>= x;
int y = a == 1 ? 0
: a == 5 ? 1
: a == 25 ? 2
: a == 125 ? 3
: a == 625 ? 4
: a == 3125 ? 5
: a == 15625 ? 6
: a == 78125 ? 7
: a == 390625 ? 8
: 9;
return {x, y};
}
pair<int, int> shrink(bigint& a) {
assert(a.neg == false);
if (a.dat.empty()) return {0, 0};
pair<int, int> res{0, 0};
while (true) {
int g = gcd(bigint::D, a.dat[0]);
if (g == 1) break;
if (g != bigint::D) a *= bigint::D / g;
a.dat.erase(begin(a.dat));
auto s = shrink(g);
res.first += s.first, res.second += s.second;
}
return res;
}
template <bool FAST = false>
bigint gcd_d_ary(bigint a, bigint b) {
a.neg = b.neg = false;
if constexpr (FAST) {
if (max<int>(a.dat.size(), b.dat.size()) <= 4) {
return __int128_t(BinaryGCDImpl::binary_gcd128(a.to_i128(), b.to_i128()));
}
}
if (a.dat.empty()) return b;
if (b.dat.empty()) return a;
pair<int, int> s = shrink(a), t = shrink(b);
if (a < b) swap(a, b);
while (true) {
if (b.dat.empty()) break;
if constexpr (FAST) {
if ((int)a.dat.size() <= 4) break;
}
a = a - b;
if (!a.dat.empty()) {
while ((a.dat[0] & 1) == 0) {
int m = a.dat[0] ? __builtin_ctz(a.dat[0]) : 32;
m = min(m, bigint::logD);
int mask = (1 << m) - 1, prod = bigint::D >> m;
a.dat[0] >>= m;
for (int i = 1; i < (int)a.dat.size(); i++) {
a.dat[i - 1] += prod * (a.dat[i] & mask);
a.dat[i] >>= m;
}
if (a.dat.back() == 0) a.dat.pop_back();
}
}
if (a < b) swap(a, b);
}
assert(a >= b);
bigint g;
if constexpr (FAST) {
if (b.dat.empty()) {
g = a;
} else {
g = __int128_t(BinaryGCDImpl::binary_gcd128(a.to_i128(), b.to_i128()));
}
} else {
g = a;
}
int e2 = min(s.first, t.first);
int e5 = min(s.second, t.second);
if (e2) g *= bigint_pow(bigint{2}, e2);
if (e5) g *= bigint_pow(bigint{5}, e5);
return g;
}
} // namespace GCDforBigintImpl
MultiPrecisionInteger gcd(const MultiPrecisionInteger& a,
const MultiPrecisionInteger& b) {
return GCDforBigintImpl::gcd_d_ary<true>(a, b);
}
MultiPrecisionInteger lcm(const MultiPrecisionInteger& a,
const MultiPrecisionInteger& b) {
return a / gcd(a, b) * b;
}
using namespace std;
// T : 値, U : 比較用
template <typename T, typename U>
struct RationalBase {
using R = RationalBase;
using Key = T;
T x, y;
RationalBase() : x(0), y(1) {}
template <typename T1>
RationalBase(const T1& _x) : RationalBase<T, U>(_x, T1{1}) {}
template <typename T1, typename T2>
RationalBase(const pair<T1, T2>& _p)
: RationalBase<T, U>(_p.first, _p.second) {}
template <typename T1, typename T2>
RationalBase(const T1& _x, const T2& _y) : x(_x), y(_y) {
assert(y != 0);
if (y == -1) x = -x, y = -y;
if (y != 1) {
T g;
if constexpr (internal::is_broadly_integral_v<T>) {
if constexpr (sizeof(T) == 16) {
g = binary_gcd128(x, y);
} else {
g = binary_gcd(x, y);
}
} else {
g = gcd(x, y);
}
if (g != 0) x /= g, y /= g;
if (y < 0) x = -x, y = -y;
}
}
// y = 0 の代入も認める
static R raw(T _x, T _y) {
R r;
r.x = _x, r.y = _y;
return r;
}
friend R operator+(const R& l, const R& r) {
if (l.y == r.y) return R{l.x + r.x, l.y};
return R{l.x * r.y + l.y * r.x, l.y * r.y};
}
friend R operator-(const R& l, const R& r) {
if (l.y == r.y) return R{l.x - r.x, l.y};
return R{l.x * r.y - l.y * r.x, l.y * r.y};
}
friend R operator*(const R& l, const R& r) { return R{l.x * r.x, l.y * r.y}; }
friend R operator/(const R& l, const R& r) { return R{l.x * r.y, l.y * r.x}; }
R& operator+=(const R& r) { return (*this) = (*this) + r; }
R& operator-=(const R& r) { return (*this) = (*this) - r; }
R& operator*=(const R& r) { return (*this) = (*this) * r; }
R& operator/=(const R& r) { return (*this) = (*this) / r; }
R operator-() const { return raw(-x, y); }
R inverse() const {
assert(x != 0);
R r = raw(y, x);
if (r.y < 0) r.x = -r.x, r.y = -r.y;
return r;
}
R pow(long long p) const {
R res{1}, base{*this};
while (p) {
if (p & 1) res *= base;
base *= base;
p >>= 1;
}
return res;
}
friend bool operator==(const R& l, const R& r) {
return l.x == r.x && l.y == r.y;
};
friend bool operator!=(const R& l, const R& r) {
return l.x != r.x || l.y != r.y;
};
friend bool operator<(const R& l, const R& r) {
return U{l.x} * r.y < U{l.y} * r.x;
};
friend bool operator<=(const R& l, const R& r) { return l < r || l == r; }
friend bool operator>(const R& l, const R& r) {
return U{l.x} * r.y > U{l.y} * r.x;
};
friend bool operator>=(const R& l, const R& r) { return l > r || l == r; }
friend ostream& operator<<(ostream& os, const R& r) {
os << r.x;
if (r.x != 0 && r.y != 1) os << "/" << r.y;
return os;
}
// T にキャストされるので T が bigint の場合は to_ll も要る
T to_mint(T mod) const {
assert(mod != 0);
T a = y, b = mod, u = 1, v = 0, t;
while (b > 0) {
t = a / b;
swap(a -= t * b, b);
swap(u -= t * v, v);
}
return U((u % mod + mod) % mod) * x % mod;
}
};
using Rational = RationalBase<long long, __int128_t>;
using BigRational = RationalBase<bigint, bigint>;
double to_double(const BigRational& r) {
pair<long double, int> a = r.x.dfp();
pair<long double, int> b = r.y.dfp();
return a.first / b.first * powl(10.0, a.second - b.second);
}
using namespace std;
using namespace std;
// x / y (x > 0, y > 0) を管理、デフォルトで 1 / 1
// 入力が互いに素でない場合は gcd を取って格納
// seq : (1, 1) から (x, y) へのパス。右の子が正/左の子が負
template <typename Int>
struct SternBrocotTreeNode {
using Node = SternBrocotTreeNode;
Int lx, ly, x, y, rx, ry;
vector<Int> seq;
SternBrocotTreeNode() : lx(0), ly(1), x(1), y(1), rx(1), ry(0) {}
SternBrocotTreeNode(Int X, Int Y) : SternBrocotTreeNode() {
assert(1 <= X && 1 <= Y);
Int g = gcd(X, Y);
X /= g, Y /= g;
while (min(X, Y) > 0) {
if (X > Y) {
Int d = X / Y;
X -= d * Y;
go_right(d - (X == 0 ? 1 : 0));
} else {
Int d = Y / X;
Y -= d * X;
go_left(d - (Y == 0 ? 1 : 0));
}
}
}
SternBrocotTreeNode(const pair<Int, Int> &xy)
: SternBrocotTreeNode(xy.first, xy.second) {}
SternBrocotTreeNode(const vector<Int> &_seq) : SternBrocotTreeNode() {
for (const Int &d : _seq) {
assert(d != 0);
if (d > 0) go_right(d);
if (d < 0) go_left(d);
}
assert(seq == _seq);
}
// pair<Int, Int> 型で分数を get
pair<Int, Int> get() const { return make_pair(x, y); }
// 区間の下限
pair<Int, Int> lower_bound() const { return make_pair(lx, ly); }
// 区間の上限
pair<Int, Int> upper_bound() const { return make_pair(rx, ry); }
// 根からの深さ
Int depth() const {
Int res = 0;
for (auto &s : seq) res += abs(s);
return res;
}
// 左の子に d 進む
void go_left(Int d = 1) {
if (d <= 0) return;
if (seq.empty() or seq.back() > 0) seq.push_back(0);
seq.back() -= d;
rx += lx * d, ry += ly * d;
x = rx + lx, y = ry + ly;
}
// 右の子に d 進む
void go_right(Int d = 1) {
if (d <= 0) return;
if (seq.empty() or seq.back() < 0) seq.push_back(0);
seq.back() += d;
lx += rx * d, ly += ry * d;
x = rx + lx, y = ry + ly;
}
// 親の方向に d 進む
// d 進めたら true, 進めなかったら false を返す
bool go_parent(Int d = 1) {
if (d <= 0) return true;
while (d != 0) {
if (seq.empty()) return false;
Int d2 = min(d, abs(seq.back()));
if (seq.back() > 0) {
x -= rx * d2, y -= ry * d2;
lx = x - rx, ly = y - ry;
seq.back() -= d2;
} else {
x -= lx * d2, y -= ly * d2;
rx = x - lx, ry = y - ly;
seq.back() += d2;
}
d -= d2;
if (seq.back() == 0) seq.pop_back();
if (d2 == Int{0}) break;
}
return true;
}
// SBT 上の LCA
static Node lca(const Node &lhs, const Node &rhs) {
Node n;
for (int i = 0; i < min<int>(lhs.seq.size(), rhs.seq.size()); i++) {
Int val1 = lhs.seq[i], val2 = rhs.seq[i];
if ((val1 < 0) != (val2 < 0)) break;
if (val1 < 0) n.go_left(min(-val1, -val2));
if (val1 > 0) n.go_right(min(val1, val2));
if (val1 != val2) break;
}
return n;
}
friend ostream &operator<<(ostream &os, const Node &rhs) {
os << "\n";
os << "L : ( " << rhs.lx << ", " << rhs.ly << " )\n";
os << "M : ( " << rhs.x << ", " << rhs.y << " )\n";
os << "R : ( " << rhs.rx << ", " << rhs.ry << " )\n";
os << "seq : {";
for (auto &x : rhs.seq) os << x << ", ";
os << "} \n";
return os;
}
friend bool operator<(const Node &lhs, const Node &rhs) {
return lhs.x * rhs.y < rhs.x * lhs.y;
}
friend bool operator==(const Node &lhs, const Node &rhs) {
return lhs.x == rhs.x and lhs.y == rhs.y;
}
};
/**
* @brief Stern-Brocot Tree
*/
// 分子と分母が INF 以下である非負の既約分数のうち次のものを返す
// first : f(x) が false である最大の既約分数 x
// second : f(x) が true である最小の既約分数 x
// ただし
// - f(0) = true の場合は (0/1, 0/1) を返す
// - true になる分数が存在しない場合は (?, 1/0) を返す
// - INF = 0 の場合は (0/1, 1/0) を返す
template <typename I>
pair<pair<I, I>, pair<I, I>> binary_search_on_stern_brocot_tree(
function<bool(pair<I, I>)> f, const I &INF) {
// INF >= 0
assert(0 <= INF);
SternBrocotTreeNode<I> m;
if (INF == 0) return {m.lower_bound(), m.upper_bound()};
// INF 条件を超える or f(m) = return_value である
auto over = [&](bool return_value) {
return max(m.x, m.y) > INF or f(m.get()) == return_value;
};
if (f(make_pair(0, 1))) return {m.lower_bound(), m.lower_bound()};
int go_left = over(true);
for (; true; go_left ^= 1) {
if (go_left) {
// f(M) = true -> (L, M] に答えがある
// (f(L * b + M) = false) or (INF 超え) になる b の最小は?
I a = 1;
for (; true; a *= 2) {
m.go_left(a);
if (over(false)) {
m.go_parent(a);
break;
}
}
for (a /= 2; a != 0; a /= 2) {
m.go_left(a);
if (over(false)) m.go_parent(a);
}
m.go_left(1);
if (max(m.get().first, m.get().second) > INF)
return {m.lower_bound(), m.upper_bound()};
} else {
// f(M) = false -> (M, R] に答えがある
// (f(M + R * b) = true) or (INF 超え) になる b の最小は?
I a = 1;
for (; true; a *= 2) {
m.go_right(a);
if (over(true)) {
m.go_parent(a);
break;
}
}
for (a /= 2; a != 0; a /= 2) {
m.go_right(a);
if (over(true)) m.go_parent(a);
}
m.go_right(1);
if (max(m.get().first, m.get().second) > INF)
return {m.lower_bound(), m.upper_bound()};
}
}
}
using namespace Nyaan;
void q() {
ini(N, X, K);
vl A(K);
in(A);
vl dp(N, inf);
queue<int> Q;
auto add = [&](int i, int x) {
if (amin(dp[i], x)) Q.push(i);
};
add(0, 0);
while (sz(Q)) {
auto c = Q.front();
Q.pop();
each(a, A) {
int d = (c + N - a) % N;
add(d, dp[c] + 1);
}
}
trc(dp);
auto old_dp = dp;
sort(all(dp));
auto rui = mkrui(dp);
auto f = [&](pair<bigint, bigint> p) {
BigRational r{p};
ll ok = 0, ng = N;
while (ok + 1 < ng) {
int m = (ok + ng) / 2;
(dp[m] <= r + 1 ? ok : ng) = m;
}
// [0, ng) : 振り直さない
// [ng, N) : 振り直す
// E : 期待値の平均
// N E = x + y (1 + E)
// (N - y) E = x + y
// E = (x + y) / (N - y)
ll x = rui[ng];
ll y = N - ng;
BigRational s{x + y, N - y};
return r >= s;
};
auto [L, R] = binary_search_on_stern_brocot_tree<bigint>(f, TEN(18));
trc(L, R);
BigRational E = R;
BigRational ans = min(E + 1, BigRational{old_dp[X], 1});
out(ans.x, ans.y);
}
void Nyaan::solve() {
int t = 1;
// in(t);
while (t--) q();
}
这程序好像有点Bug,我给组数据试试?
详细
Test #1:
score: 100
Accepted
time: 1ms
memory: 3708kb
input:
6 3 2 2 4
output:
8 3
result:
ok 2 number(s): "8 3"
Test #2:
score: 0
Accepted
time: 1ms
memory: 3672kb
input:
5 4 1 1
output:
1 1
result:
ok 2 number(s): "1 1"
Test #3:
score: 0
Accepted
time: 34ms
memory: 5800kb
input:
99999 65238 100 64714 45675 36156 13116 93455 22785 10977 60219 14981 25839 83709 80404 41400 12469 31530 65521 35436 20326 96792 50699 27522 98233 26187 12509 90992 72693 83919 74145 80892 68422 38333 33497 89154 88403 77492 4570 3908 59194 3482 89871 96330 45114 5555 73987 95832 476 949 74649 2084...
output:
3 1
result:
ok 2 number(s): "3 1"
Test #4:
score: 0
Accepted
time: 2ms
memory: 3876kb
input:
10000 23 7 9594 8998 9330 6851 1662 6719 583
output:
42726 4805
result:
ok 2 number(s): "42726 4805"
Test #5:
score: 0
Accepted
time: 1ms
memory: 3640kb
input:
100 3 100 7 68 28 98 19 32 90 79 92 40 96 30 95 91 71 15 33 18 69 1 61 43 5 75 73 64 58 100 88 20 99 37 17 22 82 67 70 55 47 80 66 12 4 24 26 54 74 57 21 77 86 89 83 29 46 31 2 16 49 48 25 93 52 9 85 84 42 39 8 65 10 45 63 87 78 60 23 14 34 59 81 38 41 76 3 13 27 36 35 51 44 62 53 94 6 50 11 97 72 56
output:
1 1
result:
ok 2 number(s): "1 1"
Test #6:
score: 0
Accepted
time: 1ms
memory: 3640kb
input:
100 93 4 63 58 3 89
output:
19 4
result:
ok 2 number(s): "19 4"
Test #7:
score: 0
Accepted
time: 5ms
memory: 5024kb
input:
75057 45721 3 10861 27551 14278
output:
32797 933
result:
ok 2 number(s): "32797 933"
Test #8:
score: 0
Accepted
time: 8ms
memory: 5548kb
input:
97777 94043 1 83579
output:
97619 221
result:
ok 2 number(s): "97619 221"
Test #9:
score: 0
Accepted
time: 2ms
memory: 4156kb
input:
13515 10596 5 11890 9097 4596 13464 13309
output:
6022 489
result:
ok 2 number(s): "6022 489"
Test #10:
score: 0
Accepted
time: 5ms
memory: 4928kb
input:
77777 64477 3 45863 40922 74543
output:
298537 8416
result:
ok 2 number(s): "298537 8416"
Test #11:
score: 0
Accepted
time: 0ms
memory: 3776kb
input:
1 0 1 1
output:
0 1
result:
ok 2 number(s): "0 1"
Test #12:
score: 0
Accepted
time: 3ms
memory: 3924kb
input:
11254 5306 33 4933 97 3341 7991 766 11039 6490 8955 10986 642 421 4570 9198 3221 9106 5076 8660 517 8376 4918 10847 10400 9063 8416 4673 7139 3925 7192 8391 7763 4927 10373 3726
output:
3 1
result:
ok 2 number(s): "3 1"
Test #13:
score: 0
Accepted
time: 1ms
memory: 3704kb
input:
33 32 6 20 26 25 18 5 28
output:
101 32
result:
ok 2 number(s): "101 32"
Test #14:
score: 0
Accepted
time: 33ms
memory: 5464kb
input:
100000 56979 500 6945 45095 52485 23545 72920 30450 64925 31700 11155 65550 80965 77885 97915 26950 86940 50045 46645 74740 35235 13470 18315 72135 25260 88775 4405 25525 61335 97040 13240 2905 37460 51605 2330 54995 17100 30790 68205 53195 63340 85495 48535 94865 87720 52930 30650 47115 68380 24455...
output:
13953 2000
result:
ok 2 number(s): "13953 2000"
Test #15:
score: 0
Accepted
time: 77ms
memory: 5968kb
input:
100000 88341 500 35014 85376 26718 27010 22158 56540 54016 52932 81956 92630 79862 844 37070 30304 9780 50642 56332 91384 20562 17764 94836 28314 23928 46364 59128 5820 21488 60066 5262 8964 53054 28310 79006 92364 20872 34916 10934 56274 1020 23138 70610 85702 56844 99122 56842 8772 61584 5346 4819...
output:
40687 10000
result:
ok 2 number(s): "40687 10000"
Test #16:
score: 0
Accepted
time: 150ms
memory: 5804kb
input:
100000 80536 500 91882 33957 97622 44074 36850 17254 63750 32583 12125 43292 58710 2227 31892 64531 13775 20782 72858 85509 37363 57708 26191 2037 8279 5828 38370 88145 76527 73627 28868 68151 31497 78935 1779 63126 37781 57818 46131 81930 43171 25281 1175 90726 12883 73954 7094 46278 62218 33972 54...
output:
2 1
result:
ok 2 number(s): "2 1"
Test #17:
score: 0
Accepted
time: 3ms
memory: 5528kb
input:
100000 14944 4 76800 80400 34400 91600
output:
101433 250
result:
ok 2 number(s): "101433 250"
Test #18:
score: 0
Accepted
time: 15ms
memory: 5820kb
input:
100000 26729 100 31512 64212 8166 14162 80024 52426 94306 14110 22268 53434 38668 16458 16622 2704 28246 72148 28688 59992 53220 13012 39628 8452 52162 21652 93712 79296 5692 41374 22558 91772 92492 93806 81856 54222 32254 11294 74172 47658 67984 41378 83042 3434 21248 82852 51372 62316 26610 26536 ...
output:
24639 5000
result:
ok 2 number(s): "24639 5000"
Test #19:
score: 0
Accepted
time: 6ms
memory: 5540kb
input:
99999 30274 13 51864 40473 61395 72180 91728 53574 49005 42855 71625 92499 441 95985 99123
output:
308066 33333
result:
ok 2 number(s): "308066 33333"
Test #20:
score: 0
Accepted
time: 0ms
memory: 5100kb
input:
72000 53726 7 70632 15840 58896 57312 46440 38376 64368
output:
9633 125
result:
ok 2 number(s): "9633 125"
Test #21:
score: 0
Accepted
time: 1ms
memory: 3648kb
input:
512 179 5 424 8 128 280 256
output:
707 64
result:
ok 2 number(s): "707 64"
Test #22:
score: 0
Accepted
time: 3ms
memory: 5640kb
input:
99855 34584 8 45675 3780 4725 6930 75285 34335 41895 50400
output:
101037 317
result:
ok 2 number(s): "101037 317"
Test #23:
score: 0
Accepted
time: 0ms
memory: 5388kb
input:
88888 49866 3 88888 44444 22222
output:
22223 1
result:
ok 2 number(s): "22223 1"
Test #24:
score: 0
Accepted
time: 1ms
memory: 3916kb
input:
16 11 1 12
output:
11 2
result:
ok 2 number(s): "11 2"
Extra Test:
score: 0
Extra Test Passed