結果
| 問題 | No.1648 Sum of Powers |
| ユーザー |
 NyaanNyaan
|
| 提出日時 | 2021-08-13 21:50:08 |
| 言語 | C++17 (gcc 12.3.0 + boost 1.83.0) |
| 結果 |
RE
|
| 実行時間 | - |
| コード長 | 22,336 bytes |
| コンパイル時間 | 3,517 ms |
| コンパイル使用メモリ | 276,040 KB |
| 実行使用メモリ | 21,156 KB |
| 最終ジャッジ日時 | 2024-10-03 18:26:29 |
| 合計ジャッジ時間 | 10,045 ms |
|
ジャッジサーバーID （参考情報） |
judge5 / judge2 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 119 ms
21,156 KB |
| testcase_01 | AC | 40 ms
6,820 KB |
| testcase_02 | RE | - |
| testcase_03 | RE | - |
| testcase_04 | RE | - |
| testcase_05 | RE | - |
| testcase_06 | RE | - |
| testcase_07 | RE | - |
| testcase_08 | RE | - |
| testcase_09 | RE | - |
| testcase_10 | RE | - |
| testcase_11 | RE | - |
| testcase_12 | RE | - |
| testcase_13 | RE | - |
| testcase_14 | RE | - |
| testcase_15 | RE | - |
| testcase_16 | RE | - |
| testcase_17 | RE | - |
| testcase_18 | RE | - |
| testcase_19 | RE | - |
| testcase_20 | RE | - |
| testcase_21 | RE | - |
| testcase_22 | TLE | - |
| testcase_23 | -- | - |
| testcase_24 | -- | - |
| testcase_25 | -- | - |
| testcase_26 | -- | - |
| testcase_27 | -- | - |
| testcase_28 | -- | - |
| testcase_29 | -- | - |
| testcase_30 | -- | - |
| testcase_31 | -- | - |
| testcase_32 | -- | - |
| testcase_33 | -- | - |
| testcase_34 | -- | - |
| testcase_35 | -- | - |
| testcase_36 | -- | - |
| testcase_37 | -- | - |
| testcase_38 | -- | - |
| testcase_39 | -- | - |
| testcase_40 | -- | - |
| testcase_41 | -- | - |
| testcase_42 | -- | - |
| testcase_43 | -- | - |
| testcase_44 | -- | - |
| testcase_45 | -- | - |
| testcase_46 | -- | - |
| testcase_47 | -- | - |
| testcase_48 | -- | - |
| testcase_49 | -- | - |
| testcase_50 | -- | - |
| testcase_51 | -- | - |
| testcase_52 | -- | - |
| testcase_53 | -- | - |
| testcase_54 | -- | - |
| testcase_55 | -- | - |
| testcase_56 | -- | - |
| testcase_57 | -- | - |
ソースコード
/**
* date : 2021-08-13 21:50:05
*/
#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, 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;
T &x() { return first; }
const T &x() const { return first; }
U &y() { return second; }
const U &y() const { return 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;
}
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; }
};
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;
}
} // 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;
}
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...);
}
void outr() {}
template <typename T, class... U, char sep = ' '>
void outr(const T &t, const U &... u) {
cout << t;
outr(u...);
}
struct IoSetupNya {
IoSetupNya() {
cin.tie(nullptr);
ios::sync_with_stdio(false);
cout << fixed << setprecision(15);
cerr << fixed << setprecision(7);
}
} iosetupnya;
} // namespace Nyaan
// debug
namespace DebugImpl {
template <typename U, typename = void>
struct is_specialize : false_type {};
template <typename U>
struct is_specialize<
U, typename conditional<false, typename U::iterator, void>::type>
: true_type {};
template <typename U>
struct is_specialize<
U, typename conditional<false, decltype(U::first), void>::type>
: true_type {};
template <typename U>
struct is_specialize<U, enable_if_t<is_integral<U>::value, void>> : true_type {
};
void dump(const char& t) { cerr << t; }
void dump(const string& t) { cerr << t; }
void dump(const bool& t) { cerr << (t ? "true" : "false"); }
template <typename U,
enable_if_t<!is_specialize<U>::value, nullptr_t> = nullptr>
void dump(const U& t) {
cerr << t;
}
template <typename T>
void dump(const T& t, enable_if_t<is_integral<T>::value>* = nullptr) {
string res;
if (t == Nyaan::inf) res = "inf";
if constexpr (is_signed<T>::value) {
if (t == -Nyaan::inf) res = "-inf";
}
if constexpr (sizeof(T) == 8) {
if (t == Nyaan::infLL) res = "inf";
if constexpr (is_signed<T>::value) {
if (t == -Nyaan::infLL) res = "-inf";
}
}
if (res.empty()) res = to_string(t);
cerr << res;
}
template <typename T, typename U>
void dump(const pair<T, U>&);
template <typename T>
void dump(const pair<T*, int>&);
template <typename T>
void dump(const T& t,
enable_if_t<!is_void<typename T::iterator>::value>* = nullptr) {
cerr << "[ ";
for (auto it = t.begin(); it != t.end();) {
dump(*it);
cerr << (++it == t.end() ? "" : ", ");
}
cerr << " ]";
}
template <typename T, typename U>
void dump(const pair<T, U>& t) {
cerr << "( ";
dump(t.first);
cerr << ", ";
dump(t.second);
cerr << " )";
}
template <typename T>
void dump(const pair<T*, int>& t) {
cerr << "[ ";
for (int i = 0; i < t.second; i++) {
dump(t.first[i]);
cerr << (i == t.second - 1 ? "" : ", ");
}
cerr << " ]";
}
void trace() { cerr << endl; }
template <typename Head, typename... Tail>
void trace(Head&& head, Tail&&... tail) {
cerr << " ";
dump(head);
if (sizeof...(tail) != 0) cerr << ",";
trace(forward<Tail>(tail)...);
}
} // namespace DebugImpl
#ifdef NyaanDebug
#define trc(...) \
do { \
cerr << "## " << #__VA_ARGS__ << " = "; \
DebugImpl::trace(__VA_ARGS__); \
} while (0)
#else
#define trc(...) (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(); }
//
template <typename mint>
std::pair<int, mint> GaussElimination(vector<vector<mint>> &a,
int pivot_end = -1,
bool diagonalize = false) {
int H = a.size(), W = a[0].size();
int rank = 0, je = pivot_end;
if (je == -1) je = W;
mint det = 1;
for (int j = 0; j < je; j++) {
int idx = -1;
for (int i = rank; i < H; i++) {
if (a[i][j] != mint(0)) {
idx = i;
break;
}
}
if (idx == -1) {
det = 0;
continue;
}
if (rank != idx) {
det = -det;
swap(a[rank], a[idx]);
}
det *= a[rank][j];
if (diagonalize && a[rank][j] != mint(1)) {
mint coeff = a[rank][j].inverse();
for (int k = j; k < W; k++) a[rank][k] *= coeff;
}
int is = diagonalize ? 0 : rank + 1;
for (int i = is; i < H; i++) {
if (i == rank) continue;
if (a[i][j] != mint(0)) {
mint coeff = a[i][j] / a[rank][j];
for (int k = j; k < W; k++) a[i][k] -= a[rank][k] * coeff;
}
}
rank++;
}
return make_pair(rank, det);
}
template <typename mint>
vector<vector<mint>> inverse_matrix(const vector<vector<mint>>& a) {
int N = a.size();
assert(N > 0);
assert(N == (int)a[0].size());
vector<vector<mint>> m(N, vector<mint>(2 * N));
for (int i = 0; i < N; i++) {
copy(begin(a[i]), end(a[i]), begin(m[i]));
m[i][N + i] = 1;
}
auto [rank, det] = GaussElimination(m, N, true);
if (rank != N) return {};
vector<vector<mint>> b(N);
for (int i = 0; i < N; i++) {
copy(begin(m[i]) + N, end(m[i]), back_inserter(b[i]));
}
return b;
}
template <class T>
struct Matrix {
vector<vector<T> > A;
Matrix() = default;
Matrix(int n, int m) : A(n, vector<T>(m, T())) {}
Matrix(int n) : A(n, vector<T>(n, T())){};
int H() const { return A.size(); }
int W() const { return A[0].size(); }
int size() const { return A.size(); }
inline const vector<T> &operator[](int k) const { return A[k]; }
inline vector<T> &operator[](int k) { return A[k]; }
static Matrix I(int n) {
Matrix mat(n);
for (int i = 0; i < n; i++) mat[i][i] = 1;
return (mat);
}
Matrix &operator+=(const Matrix &B) {
int n = H(), m = W();
assert(n == B.H() && m == B.W());
for (int i = 0; i < n; i++)
for (int j = 0; j < m; j++) (*this)[i][j] += B[i][j];
return (*this);
}
Matrix &operator-=(const Matrix &B) {
int n = H(), m = W();
assert(n == B.H() && m == B.W());
for (int i = 0; i < n; i++)
for (int j = 0; j < m; j++) (*this)[i][j] -= B[i][j];
return (*this);
}
Matrix &operator*=(const Matrix &B) {
int n = H(), m = B.W(), p = W();
assert(p == B.H());
vector<vector<T> > C(n, vector<T>(m, T{}));
for (int i = 0; i < n; i++)
for (int k = 0; k < p; k++)
for (int j = 0; j < m; j++) C[i][j] += (*this)[i][k] * B[k][j];
A.swap(C);
return (*this);
}
Matrix &operator^=(long long k) {
Matrix B = Matrix::I(H());
while (k > 0) {
if (k & 1) B *= *this;
*this *= *this;
k >>= 1LL;
}
A.swap(B.A);
return (*this);
}
Matrix operator+(const Matrix &B) const { return (Matrix(*this) += B); }
Matrix operator-(const Matrix &B) const { return (Matrix(*this) -= B); }
Matrix operator*(const Matrix &B) const { return (Matrix(*this) *= B); }
Matrix operator^(const long long k) const { return (Matrix(*this) ^= k); }
bool operator==(const Matrix &B) const {
assert(H() == B.H() && W() == B.W());
for (int i = 0; i < H(); i++)
for (int j = 0; j < W(); j++)
if (A[i][j] != B[i][j]) return false;
return true;
}
bool operator!=(const Matrix &B) const {
assert(H() == B.H() && W() == B.W());
for (int i = 0; i < H(); i++)
for (int j = 0; j < W(); j++)
if (A[i][j] != B[i][j]) return true;
return false;
}
friend ostream &operator<<(ostream &os, const Matrix &p) {
int n = p.H(), m = p.W();
for (int i = 0; i < n; i++) {
os << (i ? " " : "") << "[";
for (int j = 0; j < m; j++) {
os << p[i][j] << (j + 1 == m ? "]\n" : ",");
}
}
return (os);
}
T determinant() const {
Matrix B(*this);
assert(H() == W());
T ret = 1;
for (int i = 0; i < H(); i++) {
int idx = -1;
for (int j = i; j < W(); j++) {
if (B[j][i] != 0) {
idx = j;
break;
}
}
if (idx == -1) return 0;
if (i != idx) {
ret *= T(-1);
swap(B[i], B[idx]);
}
ret *= B[i][i];
T inv = T(1) / B[i][i];
for (int j = 0; j < W(); j++) {
B[i][j] *= inv;
}
for (int j = i + 1; j < H(); j++) {
T a = B[j][i];
if (a == 0) continue;
for (int k = i; k < W(); k++) {
B[j][k] -= B[i][k] * a;
}
}
}
return ret;
}
};
/**
* @brief 行列ライブラリ
*/
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(r * mod == 1, "invalid, r * mod != 1");
static_assert(mod < (1 << 30), "invalid, mod >= 2 ^ 30");
static_assert((mod & 1) == 1, "invalid, mod % 2 == 0");
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 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 { return pow(mod - 2); }
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; }
};
struct ArbitraryLazyMontgomeryModInt {
using mint = ArbitraryLazyMontgomeryModInt;
using i32 = int32_t;
using u32 = uint32_t;
using u64 = uint64_t;
static u32 mod;
static u32 r;
static u32 n2;
static u32 get_r() {
u32 ret = mod;
for (i32 i = 0; i < 4; ++i) ret *= 2 - mod * ret;
return ret;
}
static void set_mod(u32 m) {
assert(m < (1 << 30));
assert((m & 1) == 1);
mod = m;
n2 = -u64(m) % m;
r = get_r();
assert(r * mod == 1);
}
u32 a;
ArbitraryLazyMontgomeryModInt() : a(0) {}
ArbitraryLazyMontgomeryModInt(const int64_t &b)
: a(reduce(u64(b % mod + mod) * n2)){};
static u32 reduce(const u64 &b) {
return (b + u64(u32(b) * u32(-r)) * mod) >> 32;
}
mint &operator+=(const mint &b) {
if (i32(a += b.a - 2 * mod) < 0) a += 2 * mod;
return *this;
}
mint &operator-=(const mint &b) {
if (i32(a -= b.a) < 0) a += 2 * mod;
return *this;
}
mint &operator*=(const mint &b) {
a = reduce(u64(a) * b.a);
return *this;
}
mint &operator/=(const mint &b) {
*this *= b.inverse();
return *this;
}
mint operator+(const mint &b) const { return mint(*this) += b; }
mint operator-(const mint &b) const { return mint(*this) -= b; }
mint operator*(const mint &b) const { return mint(*this) *= b; }
mint operator/(const mint &b) const { return mint(*this) /= b; }
bool operator==(const mint &b) const {
return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a);
}
bool operator!=(const mint &b) const {
return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a);
}
mint operator-() const { return mint() - mint(*this); }
mint pow(u64 n) const {
mint ret(1), mul(*this);
while (n > 0) {
if (n & 1) ret *= mul;
mul *= mul;
n >>= 1;
}
return ret;
}
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 = ArbitraryLazyMontgomeryModInt(t);
return (is);
}
mint inverse() const { return pow(mod - 2); }
u32 get() const {
u32 ret = reduce(a);
return ret >= mod ? ret - mod : ret;
}
static u32 get_mod() { return mod; }
};
typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::mod;
typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::r;
typename ArbitraryLazyMontgomeryModInt::u32 ArbitraryLazyMontgomeryModInt::n2;
int64_t mod_sqrt(const int64_t &a, const int64_t &p) {
assert(0 <= a && a < p);
if (a < 2) return a;
using Mint = ArbitraryLazyMontgomeryModInt;
Mint::set_mod(p);
if (Mint(a).pow((p - 1) >> 1) != 1) return -1;
Mint b = 1, one = 1;
while (b.pow((p - 1) >> 1) == 1) b += one;
int64_t m = p - 1, e = 0;
while (m % 2 == 0) m >>= 1, e += 1;
Mint x = Mint(a).pow((m - 1) >> 1);
Mint y = Mint(a) * x * x;
x *= a;
Mint z = Mint(b).pow(m);
while (y != 1) {
int64_t j = 0;
Mint t = y;
while (t != one) {
j += 1;
t *= t;
}
z = z.pow(int64_t(1) << (e - j - 1));
x *= z;
z *= z;
y *= z;
e = j;
}
return x.get();
}
/**
* @brief mod sqrt(Tonelli-Shanks algorithm)
* @docs docs/modulo/mod-sqrt.md
*/
template <typename mint>
vector<mint> QuadraticEquation(mint a, mint b, mint c) {
assert(mint::get_mod() % 2 != 0);
if (a == mint()) {
if(b == mint()) {
assert(c != mint());
return {};
}
return vector<mint>{-c * b.inverse()};
}
mint ia = a.inverse(), inv2 = mint(2).inverse();
b *= ia, c *= ia;
auto D = mod_sqrt(((b * inv2).pow(2) - c).get(), mint::get_mod());
if(D == -1) return {};
if(D <= 1) return vector<mint>{-b * inv2 + D};
return vector<mint>{-b * inv2 + D, -b * inv2 - D};
}
//
using namespace Nyaan;
using mint = LazyMontgomeryModInt<998244353>;
// using mint = LazyMontgomeryModInt<1000000007>;
using vm = vector<mint>;
using vvm = vector<vm>;
using namespace Nyaan;
void Nyaan::solve() {
inl(A, B, P, Q);
if(B == 0) exit(1);
using mat = Matrix<mint>;
mat M(2);
M[0][0] = A, M[0][1] = -B, M[1][0] = 1;
map<vi, int> memo;
mat MM;
mat S1(2, 1);
S1[0][0] = A, S1[1][0] = 2;
mat PQ(2, 1);
PQ[0][0] = P, PQ[1][0] = Q;
{
mat a = mat::I(2);
rep(i, TEN(5)) {
auto b = a * S1;
vi v(2);
v[0] = b[0][0].get(), v[1] = b[1][0].get();
memo[v] = i;
a *= M;
}
MM = M;
}
mat base = mat::I(2);
for (int i = 0;; i++) {
if (base.determinant() != 0) {
mat inv(2);
inv.A = inverse_matrix(base.A);
mat mt = inv * PQ;
vi v(2);
v[0] = mt[0][0].get(), v[1] = mt[1][0].get();
if (memo.find(v) != memo.end()) {
out(i * TEN(5) + memo[v] + 1);
exit(0);
}
}
base *= MM;
}
}