結果
| 問題 |
No.1648 Sum of Powers
|
| コンテスト | |
| ユーザー |
 Kude
|
| 提出日時 | 2021-08-18 13:02:28 |
| 言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 30 ms / 2,000 ms |
| コード長 | 7,026 bytes |
| コンパイル時間 | 4,560 ms |
| コンパイル使用メモリ | 265,516 KB |
| 最終ジャッジ日時 | 2025-01-23 22:53:38 |
|
ジャッジサーバーID （参考情報） |
judge5 / judge4 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 2 |
| other | AC * 56 |
ソースコード
#include<bits/stdc++.h>
#include<atcoder/all>
using namespace std;
using namespace atcoder;
#define rep(i,n)for (int i = 0; i < int(n); ++i)
#define rrep(i,n)for (int i = int(n)-1; i >= 0; --i)
#define all(x) (x).begin(), (x).end()
#define rall(x) (x).rbegin(), (x).rend()
template<class T> void chmax(T& a, const T& b) {a = max(a, b);}
template<class T> void chmin(T& a, const T& b) {a = min(a, b);}
using ll = long long;
using P = pair<int,int>;
using VI = vector<int>;
using VVI = vector<VI>;
using VL = vector<ll>;
using VVL = vector<VL>;
using mint = modint998244353;
template<class T, int N>
struct Mat: std::array<std::array<T, N>, N> {
friend Mat<T, N> operator*(const Mat& A, const Mat& B) {
Mat<T, N> C = {};
for(int i = 0; i < N; i++) {
for(int k = 0; k < N; k++) {
for(int j = 0; j < N; j++) {
C[i][j] += A[i][k] * B[k][j];
}
}
}
return C;
}
friend std::array<T, N> operator*(const Mat& A, const std::array<T, N>& v) {
std::array<T, N> x = {};
for(int i = 0; i < N; i++) for(int j = 0; j < N; j++) x[i] += A[i][j] * v[j];
return x;
}
friend Mat<T, N> operator+(const Mat& A, const Mat& B) {
Mat<T, N> C;
for(int i = 0; i < N; i++) {
for(int j = 0; j < N; j++) {
C[i][j] = A[i][j] + B[i][j];
}
}
return C;
}
Mat<T, N>& operator*=(const Mat& A) { return *this = *this * A; }
Mat<T, N>& operator+=(const Mat& A) {
for(int i = 0; i < N; i++) {
for(int j = 0; j < N; j++) {
(*this)[i][j] += A[i][j];
}
}
}
static Mat<T, N> I() {
Mat<T, N> X = {};
for(int i = 0; i < N; i++) X[i][i] = 1;
return X;
}
Mat<T, N> pow(long long k) const {
assert(k >= 0);
auto X = *this;
auto Y = I();
while(k) {
if (k & 1) Y *= X;
k >>= 1;
if (k) X *= X;
}
return Y;
}
Mat<T, N> inv() const {
auto X = *this;
auto Y = I();
for(int p = 0; p < N; p++) {
bool ok = false;
for(int i = p; i < N; i++) {
if (X[i][p] != T()) {
ok = true;
if (i != p) {
// std::swap(X[i], X[p]);
for(int j = p; j < N; j++) std::swap(X[i][j], X[p][j]);
std::swap(Y[i], Y[p]);
}
break;
}
}
assert(ok);
T c = 1 / X[p][p];
for(int j = p; j < N; j++) X[p][j] *= c;
for(int j = 0; j < N; j++) Y[p][j] *= c;
for(int i = 0; i < N; i++) if (i != p) {
T c = X[i][p];
for(int j = p; j < N; j++) X[i][j] -= c * X[p][j];
for(int j = 0; j < N; j++) Y[i][j] -= c * Y[p][j];
}
}
return Y;
}
T det() const {
if (N == 0) return T();
auto X = *this;
bool flag = false;
T res = 1;
bool ok = false;
for(int p = 0; p < N; p++) {
for(int i = p; i < N; i++) {
if (X[i][p] != T()) {
ok = true;
if (i != p) {
flag = !flag;
for(int j = p; j < N; j++) std::swap(X[i][j], X[p][j]);
}
break;
}
}
if (!ok) return T();
// T = modint
for(int i = p + 1; i < N; i++) {
T c0 = 1, c1 = 0, d0 = 0, d1 = 1;
int a = X[p][p].val(), b = X[i][p].val();
while(true) {
if (b == 0) {
break;
}
int q = a / b;
a -= q * b;
c0 -= q * d0;
c1 -= q * d1;
if (a == 0) {
flag = !flag;
for(int j = p; j < N; j++) std::swap(X[i][j], X[p][j]);
std::swap(c0, d1);
std::swap(c1, d0);
break;
}
q = b / a;
b -= q * a;
d0 -= q * c0;
d1 -= q * c1;
}
for(int j = p; j < N; j++) {
T a = X[p][j], b = X[i][j];
T x = c0 * a + c1 * b, y = d0 * a + d1 * b;
X[p][j] = x;
X[i][j] = y;
}
}
res *= X[p][p];
}
if (flag) res = -res;
return res;
}
friend std::ostream& operator<<(std::ostream& os, const Mat<T, N>& A) {
for(int i = 0; i < N; i++) for(int j = 0; j < N; j++) os << A[i][j] << " \n"[j + 1 == N];
return os;
}
};
std::ostream& operator<<(std::ostream& os, const mint& x) {
return os << x.val();
}
template<class mint>
int BSGS(mint a, mint b) {
// constraint: a.inv() is defined if a != 0
// returns x s.t. a ^ x = b
const int p = mint::mod();
if (p == 1) return 0;
if (a == 0) {
if (b == 0) return 1;
else if (b == 1) return 0;
else return -1;
}
const int sq_ceil = [p]{
int l = 0, r = 46340;
while(r - l > 1) {
int c = (r + l) >> 1;
(c * c >= p ? r : l) = c;
}
return r;
}();
std::unordered_map<int, int> log_a;
mint aj = 1;
for(int j = 0; j < sq_ceil; j++, aj *= a) {
if (!log_a.count(aj.val())) {
log_a[aj.val()] = j;
}
}
mint a_inv = a.inv();
mint step = a_inv.pow(sq_ceil);
mint t = b;
for(int i = 0; i < p; i += sq_ceil, t *= step) {
auto it = log_a.find(t.val());
if (it != log_a.end()) return i + it->second;
}
return -1;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
int a, b, p, q;
cin >> a >> b >> p >> q;
if (b == 0) {
if (a == 0) {
cout << 2 << '\n';
} else {
// a^(n-1) = q
int n = BSGS<mint>(a, mint(q) / a) + 2;
cout << n << '\n';
}
return 0;
}
Mat<mint, 2> X{array<mint, 2>{a, -b}, array<mint, 2>{1, 0}};
array<mint, 2> v{a, 2};
v = X * v;
unordered_map<ll, int> log_a;
for(int j = 0; j < 100000; j++) {
ll t = v[0].val() | ll(v[1].val()) << 32;
if (!log_a.count(t)) log_a[t] = j;
v = X * v;
}
auto Y = X.inv().pow(100000);
v = array<mint, 2>{p, q};
for(int i = 0; i < 100000; i++) {
ll t = v[0].val() | ll(v[1].val()) << 32;
auto it = log_a.find(t);
if (it != log_a.end()) {
cout << 100000LL * i + it->second + 2 << '\n';
break;
}
v = Y * v;
}
}