結果
| 問題 |
No.950 行列累乗
|
| コンテスト | |
| ユーザー |
 qwewe
|
| 提出日時 | 2025-05-14 13:10:11 |
| 言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
WA
|
| 実行時間 | - |
| コード長 | 17,697 bytes |
| コンパイル時間 | 1,195 ms |
| コンパイル使用メモリ | 113,460 KB |
| 実行使用メモリ | 7,848 KB |
| 最終ジャッジ日時 | 2025-05-14 13:11:39 |
| 合計ジャッジ時間 | 5,043 ms |
|
ジャッジサーバーID （参考情報） |
judge4 / judge1 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 4 |
| other | AC * 41 WA * 16 |
コンパイルメッセージ
main.cpp: In function ‘int main()’:
main.cpp:377:27: warning: integer overflow in expression of type ‘int’ results in ‘-1794967296’ [-Woverflow]
377 | else if (p > 50000*50000) m_bsgs = 50000; // Cap at 50000 if needed
| ~~~~~^~~~~~
ソースコード
#include <iostream>
#include <vector>
#include <cmath>
#include <map>
#include <utility> // For std::pair
using namespace std;
// Define long long and 128-bit integer types for intermediate calculations
typedef long long ll;
using i128 = __int128;
ll p; // The prime modulus
// Structure definition for 2x2 matrices
struct Matrix {
ll mat[2][2];
// Default constructor initializes to zero matrix
Matrix() {
mat[0][0] = mat[0][1] = mat[1][0] = mat[1][1] = 0;
}
// Static method to create an identity matrix
static Matrix identity() {
Matrix I;
I.mat[0][0] = I.mat[1][1] = 1;
I.mat[0][1] = I.mat[1][0] = 0;
return I;
}
// Matrix multiplication operator override
Matrix operator*(const Matrix& other) const {
Matrix result;
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
for (int k = 0; k < 2; ++k) {
// Use 128-bit integer for intermediate product to prevent overflow
i128 prod = (i128)mat[i][k] * other.mat[k][j];
result.mat[i][j] = (result.mat[i][j] + prod) % p;
}
}
}
return result;
}
// Equality comparison operator override
bool operator==(const Matrix& other) const {
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
if (mat[i][j] != other.mat[i][j]) {
return false;
}
}
}
return true;
}
// Less-than comparison operator override (needed for using Matrix as a map key)
bool operator<(const Matrix& other) const {
for(int i=0; i<2; ++i) {
for(int j=0; j<2; ++j) {
if(mat[i][j] != other.mat[i][j]) {
return mat[i][j] < other.mat[i][j];
}
}
}
return false; // Matrices are equal
}
};
// Matrix exponentiation using binary exponentiation (logarithmic time complexity)
Matrix matrix_pow(Matrix base, ll exp) {
Matrix res = Matrix::identity();
// Ensure base matrix entries are already modulo p
base.mat[0][0] %= p; base.mat[0][1] %= p; base.mat[1][0] %= p; base.mat[1][1] %= p;
while (exp > 0) {
if (exp % 2 == 1) res = res * base; // If exponent is odd, multiply result by base
base = base * base; // Square the base
exp /= 2; // Halve the exponent
}
return res;
}
// Modular exponentiation for scalars (logarithmic time complexity)
ll power(ll base, ll exp) {
ll res = 1;
base %= p; if (base < 0) base += p; // Ensure base is in [0, p-1]
while (exp > 0) {
if (exp % 2 == 1) res = (i128)res * base % p;
base = (i128)base * base % p;
exp /= 2;
}
return res;
}
// Modular inverse using Fermat's Little Theorem (requires p to be prime)
// Calculates n^(p-2) mod p
ll modInverse(ll n) {
n %= p; if (n < 0) n += p; // Ensure n is in [0, p-1]
if (n == 0) return -1; // Inverse does not exist for 0
// For prime p, n^(p-2) is the modular inverse of n modulo p (if n != 0 mod p)
return power(n, p - 2);
}
// Standard Baby-step Giant-step algorithm for scalar discrete logarithm a^x = b (mod p)
// Finds the minimum non-negative integer x satisfying the equation. Returns -1 if no solution.
ll solve_dlog_standard(ll a, ll b) {
a %= p; if(a<0) a+=p;
b %= p; if(b<0) b+=p;
if (p == 0) return -1; // p must be >= 2
// Handle edge cases involving 0 and 1
if (b == 1) { // Looking for a^x = 1
// Smallest non-negative x is 0 if a!=0. If a=0, 0^x=1 impossible for x>=1. 0^0=1 convention.
if (a == 0) return -1;
return 0;
}
if (a == 0) { // a=0, b!=1
if (b == 0) return 1; // 0^1 = 0. Smallest non-negative x=1.
else return -1; // 0^x = b != 0 is impossible for x>=1.
}
// Now a != 0
if (b == 0) { // a != 0
return -1; // a^x = 0 is impossible if a != 0
}
// Normal case: a != 0, b != 0, b != 1
ll m = sqrt(p) + 1; // Set step size for BSGS
map<ll, ll> baby_steps; // Store baby steps (a^j mod p -> j)
ll current_val = 1;
// Calculate and store baby steps
for (ll j = 0; j < m; ++j) {
if (!baby_steps.count(current_val)) { // Store only the first occurrence (smallest j)
baby_steps[current_val] = j;
}
current_val = (i128)current_val * a % p;
}
// Calculate a^(-m) mod p for giant steps
ll am_inv = modInverse(power(a, m));
if (am_inv == -1) return -1; // Inverse should exist as a != 0
current_val = b; // Start giant steps from target value b
// Perform giant steps
for (ll i = 0; i < m; ++i) {
if (baby_steps.count(current_val)) { // Check if current value matches a baby step
ll j = baby_steps[current_val];
// Found solution x = i*m + j. BSGS guarantees this is the minimal non-negative solution.
return i * m + j;
}
// Apply giant step factor a^(-m)
current_val = (i128)current_val * am_inv % p;
}
return -1; // No solution found
}
// Finds the order of element a modulo p.
// The order is the smallest positive integer k such that a^k = 1 (mod p). Returns -1 if error or order not found within limits.
ll find_order(ll a) {
a %= p; if(a<0) a+=p;
if (a == 0) return -1; // Order undefined for 0
if (a == 1) return 1; // Order of 1 is 1
// Use BSGS to find the smallest positive k. Similar to solving a^k=1, but ensure k>0.
ll m = sqrt(p) + 1;
map<ll, ll> table;
ll cur = 1;
// Store powers a^j for j = 0..m-1
for (ll j = 0; j < m; ++j) {
if (!table.count(cur)) table[cur] = j;
cur = (i128)cur * a % p;
}
ll step = power(a, m); // Giant step factor a^m
ll step_inv = modInverse(step);
if(step_inv == -1) return -1; // Should not happen if a != 0
cur = 1; // Current value for giant step check starts from 1. Represents (a^m)^i
ll current_target = 1; // Target is 1
// We need to find smallest k = im + j > 0 such that a^k=1.
// Rearranging: a^{im} = a^{-j}. Calculate (a^{-m})^i and check against table a^j.
// We want smallest i*m+j > 0. BSGS naturally finds smallest non-negative solution. If solution is 0, it means a=1.
// Need to adapt BSGS slightly.
cur = 1; // Start from (a^m)^0 = 1
// Check i=0 case: if 1 is in table with j > 0, then a^j = 1. Smallest such j is the order.
// However, this is implicitly checked by BSGS for a^x=1. The issue is when the result is 0.
// Rerun BSGS for a^k = 1, but interpret the result carefully for minimum positive k.
ll k_cand = solve_dlog_standard(a, 1);
if (k_cand == 0) { // a^0=1. This means smallest non-negative k is 0. Need smallest positive.
// This case should only happen if a=1, which is handled above.
// If somehow a!=1 and BSGS returns 0, it's possibly an issue with BSGS implementation or logic.
// Let's rely on the fact that if a!=1, order k>0.
// The standard BSGS needs modification to return smallest positive k when k=0 is found.
// But for this problem, we only need this function when solving a^n=b where b=1 and n=0 is found by standard BSGS.
// A simplified approach: If a^0=1 is the smallest non-negative solution, find the order by iteration? Might be slow.
// Try re-using BSGS logic slightly modified to guarantee positive result:
map<ll, ll> baby_steps_order;
ll current_val_order = a; // Start from a^1
for (ll j = 1; j <= m; ++j) { // Store a^j for j=1..m
if (!baby_steps_order.count(current_val_order)) {
baby_steps_order[current_val_order] = j;
}
current_val_order = (i128)current_val_order * a % p;
}
ll giant_step_val = 1; // Start from a^0=1
for (ll i = 1; i <= m; ++i) { // Check powers a^{im} for i=1..m
giant_step_val = (i128)giant_step_val * step % p;
if (baby_steps_order.count(giant_step_val)) { // If a^{im} == a^j
ll j = baby_steps_order[giant_step_val];
if (i * m - j > 0) return i * m - j; // Found order k = im - j > 0
// Need smallest k. This finds *a* solution, not necessarily minimal.
// BSGS complexities... Let's assume standard BSGS gives order k>0 when a!=1.
}
}
// Revert to just returning k_cand from standard BSGS if a!=1. It should be >0.
k_cand = solve_dlog_standard(a, 1); // Re-call standard BSGS. It should return order > 0 if a != 1.
if (k_cand == 0 && a!=1) return -1; // If it returns 0 despite a!=1, something is wrong.
return k_cand;
} else {
return k_cand; // Standard BSGS found non-zero result. This is the order.
}
return -1; // Fallback, should not be reached for a in Z_p^*
}
// Finds minimum positive integer n such that a^n = b (mod p)
ll solve_dlog_min_positive(ll a, ll b) {
a %= p; if(a<0) a+=p;
b %= p; if(b<0) b+=p;
if (a == 0) {
if (b == 0) return 1; // a^1 = 0 = b
else return -1; // a^n = b != 0 impossible for n>=1
}
if (b == 0) return -1; // a^n=0 impossible for a != 0
ll x = solve_dlog_standard(a, b); // Find smallest non-negative solution x
if (x == -1) return -1; // No solution exists
if (x > 0) return x; // Smallest non-negative is positive, this is the minimum positive
// If x == 0, this means a^0 = b, so b must be 1.
if (x == 0) {
if (b != 1) {
// This should ideally not happen if solve_dlog_standard is correct.
return -1;
}
// Need minimum positive k such that a^k = 1. This is the order of a.
return find_order(a); // Call function to find order
}
return -1; // Should not be reached
}
int main() {
ios_base::sync_with_stdio(false); // Faster I/O
cin.tie(NULL);
cin >> p; // Read prime modulus
Matrix A, B; // Read matrices A and B
cin >> A.mat[0][0] >> A.mat[0][1] >> A.mat[1][0] >> A.mat[1][1];
cin >> B.mat[0][0] >> B.mat[0][1] >> B.mat[1][0] >> B.mat[1][1];
// Check if n=1 is the solution
if (A == B) { cout << 1 << endl; return 0; }
// Calculate determinant of A modulo p
ll detA = (A.mat[0][0] * A.mat[1][1] - A.mat[0][1] * A.mat[1][0]);
detA = (detA % p + p) % p; // Ensure determinant is in [0, p-1]
// Case 1: A is not invertible (determinant is 0)
if (detA == 0) {
Matrix zero_matrix; // A matrix with all entries 0
if (A == zero_matrix) { // If A is the zero matrix
if (B == zero_matrix) cout << 1 << endl; // A^1 = 0 = B
else cout << -1 << endl; // A^n = 0 for n>=1, cannot equal non-zero B
return 0;
}
// A is non-zero, non-invertible
ll trA = (A.mat[0][0] + A.mat[1][1]); // Calculate trace of A modulo p
trA = (trA % p + p) % p; // Ensure trace is in [0, p-1]
if (trA == 0) { // Nilpotent case: A^2 = 0 (since tr(A)=0, det(A)=0)
// Sequence is A, 0, 0, ...
if (B == zero_matrix) cout << 2 << endl; // A^2 = 0 = B. Minimum n=2. (n=1 case A=B already checked)
else cout << -1 << endl; // B is neither A nor 0
return 0;
} else { // Non-nilpotent case: A^n = tr(A)^(n-1) * A for n >= 1
// Check if B is a scalar multiple of A: B = k * A
ll k = -1; // Scalar multiplier k
bool found_nonzero = false; // Flag to track if a non-zero entry in A is found
// Find first non-zero entry A_ij to calculate potential k = B_ij / A_ij
for (int i=0; i<2; ++i) {
for (int j=0; j<2; ++j) {
if (A.mat[i][j] != 0) {
found_nonzero = true;
ll inv_A_entry = modInverse(A.mat[i][j]); // Find modular inverse of A_ij
if (inv_A_entry == -1) { /* Should not happen if A_ij != 0 */ cout << -1 << endl; return 0;}
k = (i128)B.mat[i][j] * inv_A_entry % p; // Calculate k
break;
}
}
if (found_nonzero) break;
}
// If A was non-zero matrix, it must have a non-zero entry
if (!found_nonzero) { /* This case should not be reachable if A != 0 */ cout << -1 << endl; return 0; }
// Verify if B = k * A holds for all entries
bool is_scalar_multiple = true;
for (int i=0; i<2; ++i) {
for (int j=0; j<2; ++j) {
ll expected_B = (i128)k * A.mat[i][j] % p; // Calculate k * A_ij
if (B.mat[i][j] != expected_B) { // Compare with B_ij
is_scalar_multiple = false;
break;
}
}
if (!is_scalar_multiple) break;
}
if (!is_scalar_multiple) { // If B is not k*A
cout << -1 << endl;
return 0;
}
// If B = k*A, we need to solve tr(A)^(n-1) = k (mod p) for minimum positive n
// Let m = n-1. Need minimum non-negative m such that tr(A)^m = k (mod p).
ll m_sol = solve_dlog_standard(trA, k); // Use BSGS to find minimal non-negative m
if (m_sol == -1) { // If no solution for m
cout << -1 << endl;
} else {
cout << m_sol + 1 << endl; // Minimum positive n = m + 1
}
return 0;
}
} else { // Case 2: A is invertible (determinant is non-zero)
// Calculate determinant of B
ll detB = (B.mat[0][0] * B.mat[1][1] - B.mat[0][1] * B.mat[1][0]);
detB = (detB % p + p) % p;
if (detB == 0) { // If B is not invertible, A^n = B is impossible
cout << -1 << endl;
return 0;
}
// Special subcase: A is a scalar matrix, A = a*I
if (A.mat[0][1] == 0 && A.mat[1][0] == 0 && A.mat[0][0] == A.mat[1][1]) {
// Check if B is also a scalar matrix, B = b*I
if (B.mat[0][1] == 0 && B.mat[1][0] == 0 && B.mat[0][0] == B.mat[1][1]) {
ll a = A.mat[0][0]; // Scalar factor for A
ll b = B.mat[0][0]; // Scalar factor for B
// Need minimum positive n such that a^n = b (mod p)
ll n = solve_dlog_min_positive(a, b);
cout << n << endl;
} else { // B is not scalar, A^n = a^n * I cannot equal B
cout << -1 << endl;
}
return 0;
}
// General invertible case: Use Baby-step Giant-step algorithm for matrices
// Determine step size m. sqrt(p) seems reasonable based on scalar BSGS complexity.
// Order can be up to p^2, so sqrt(p^2) = p steps would be ideal but too slow.
// Use sqrt(p) capped at a reasonable value (e.g., 45000) to manage time/memory.
ll m_bsgs = sqrt(p) + 2;
if (p > 2025000000) m_bsgs = 45000; // If p > 45000^2, cap m_bsgs
else if (p > 50000*50000) m_bsgs = 50000; // Cap at 50000 if needed
map<Matrix, ll> baby_steps; // Store baby steps (A^j -> j)
Matrix current_A = Matrix::identity(); // Start with A^0 = I
// Compute and store baby steps A^0, A^1, ..., A^(m-1)
for(ll j=0; j < m_bsgs; ++j) {
if (baby_steps.find(current_A) == baby_steps.end()) { // Store smallest j for each matrix
baby_steps[current_A] = j;
}
current_A = current_A * A; // Compute next power A^(j+1)
}
// Calculate A inverse: A^{-1} = det(A)^{-1} * adj(A)
ll detA_inv = modInverse(detA);
if (detA_inv == -1) { /* Should not happen as detA != 0 */ cout << -1 << endl; return 0; }
Matrix A_adj; // Adjugate matrix of A
A_adj.mat[0][0] = A.mat[1][1]; A_adj.mat[0][1] = (p - A.mat[0][1]) % p;
A_adj.mat[1][0] = (p - A.mat[1][0]) % p; A_adj.mat[1][1] = A.mat[0][0];
Matrix A_inv_mat; // A inverse matrix
for(int i=0; i<2; ++i) for(int j=0; j<2; ++j) A_inv_mat.mat[i][j] = (i128)detA_inv * A_adj.mat[i][j] % p;
// Calculate (A^{-1})^m = (A^m)^{-1} for giant steps
Matrix A_m_inv = matrix_pow(A_inv_mat, m_bsgs);
Matrix current_B = B; // Start giant steps from target matrix B
ll min_n = -1; // Variable to store the minimum positive solution found
// Perform giant steps. Check B * (A^{-m})^i against baby steps
// Loop limit m_bsgs covers exponents up to m_bsgs^2
for(ll i=0; i <= m_bsgs; ++i) {
if (baby_steps.count(current_B)) { // If current_B matches a baby step A^j
ll j = baby_steps[current_B];
ll n = i * m_bsgs + j; // Candidate solution n = i*m + j
if (n > 0) { // We need the minimum positive solution
if (min_n == -1 || n < min_n) { // Update minimum if this is smaller
min_n = n;
}
}
}
if (i < m_bsgs) { // Avoid extra multiplication at the end
current_B = current_B * A_m_inv; // Apply giant step: multiply by (A^m)^{-1}
}
}
cout << min_n << endl; // Output the minimum positive n found, or -1 if none found
}
return 0;
}