結果
| 問題 |
No.1907 DETERMINATION
|
| ユーザー |
 qwewe
|
| 提出日時 | 2025-05-14 13:10:53 |
| 言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
| 結果 |
TLE
|
| 実行時間 | - |
| コード長 | 16,431 bytes |
| コンパイル時間 | 1,499 ms |
| コンパイル使用メモリ | 121,696 KB |
| 実行使用メモリ | 22,620 KB |
| 最終ジャッジ日時 | 2025-05-14 13:12:49 |
| 合計ジャッジ時間 | 36,419 ms |
|
ジャッジサーバーID （参考情報） |
judge4 / judge5 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 4 |
| other | AC * 48 TLE * 2 -- * 13 |
ソースコード
#include <iostream>
#include <vector>
#include <random>
#include <chrono>
#include <numeric> // Not used but potentially useful
// Modular arithmetic constants and functions
const int MOD = 998244353;
// Computes (base^exp) % MOD
long long power(long long base, long long exp) {
long long res = 1;
base %= MOD;
while (exp > 0) {
if (exp % 2 == 1) res = (res * base) % MOD;
base = (base * base) % MOD;
exp /= 2;
}
return res;
}
// Computes modular inverse using Fermat's Little Theorem
// Assumes MOD is prime and n is not a multiple of MOD
long long modInverse(long long n) {
// Ensure n is in [0, MOD-1]
n %= MOD;
if (n < 0) n += MOD;
if (n == 0) return -1; // Error: division by zero
// Use Fermat's Little Theorem: a^(p-2) = a^-1 (mod p) for prime p
long long res = power(n, MOD - 2);
return res;
}
// Type definitions
using matrix = std::vector<std::vector<long long>>;
// Polynomial represented as vector of coefficients, p[k] is coefficient of x^k
using polynomial = std::vector<long long>;
// Matrix operations
// Adds two NxN matrices A and B modulo MOD
matrix add(const matrix& A, const matrix& B, int N) {
matrix C(N, std::vector<long long>(N));
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
C[i][j] = (A[i][j] + B[i][j]) % MOD;
// Ensure result is non-negative
if (C[i][j] < 0) C[i][j] += MOD;
}
}
return C;
}
// Multiplies two NxN matrices A and B modulo MOD
matrix multiply(const matrix& A, const matrix& B, int N) {
matrix C(N, std::vector<long long>(N, 0));
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
for (int k = 0; k < N; ++k) {
C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
}
// Ensure result is non-negative
if (C[i][j] < 0) C[i][j] += MOD;
}
}
return C;
}
// Scales matrix A by scalar s modulo MOD
matrix scale(const matrix& A, long long s, int N) {
matrix C(N, std::vector<long long>(N));
s %= MOD;
if (s < 0) s += MOD; // Ensure scalar is non-negative
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
C[i][j] = (A[i][j] * s) % MOD;
}
}
return C;
}
// Computes the inverse of matrix A and its determinant using Gaussian elimination
// Returns {inverse matrix, determinant}. If singular, determinant is 0.
// Modifies the input matrix A. Pass a copy if original needed.
std::pair<matrix, long long> inverse_and_det(matrix A, int N) {
long long det = 1;
matrix inv(N, std::vector<long long>(N, 0));
for(int i=0; i<N; ++i) inv[i][i] = 1; // Initialize inverse as identity matrix
for (int i = 0; i < N; ++i) {
// Find pivot element in column i at or below row i
int pivot = i;
while (pivot < N && A[pivot][i] == 0) {
pivot++;
}
// If no non-zero element found in column i below row i-1, matrix is singular
if (pivot == N) {
return {inv, 0};
}
// Swap rows i and pivot to bring pivot element to A[i][i]
if (pivot != i) {
std::swap(A[i], A[pivot]);
std::swap(inv[i], inv[pivot]); // Apply same swap to inverse matrix
det = (MOD - det) % MOD; // Swapping rows changes determinant sign
}
long long pivot_val = A[i][i];
// Should not be 0 after the search, unless matrix is singular in complex way? Should be caught by pivot == N check.
if (pivot_val == 0) {
return {inv, 0}; // Matrix is singular
}
// Compute modular inverse of the pivot element
long long inv_pivot_val = modInverse(pivot_val);
if (inv_pivot_val == -1) { // Check if modInverse failed
return {inv, 0}; // Error or singular matrix property
}
// Update determinant
det = (det * pivot_val) % MOD;
// Normalize row i: Divide by pivot element (multiply by modular inverse)
for (int j = i; j < N; ++j) {
A[i][j] = (A[i][j] * inv_pivot_val) % MOD;
}
// Apply same normalization to the corresponding row of inverse matrix
for (int j = 0; j < N; ++j) {
inv[i][j] = (inv[i][j] * inv_pivot_val) % MOD;
}
// Eliminate non-zero elements in column i for all other rows k != i
for (int k = 0; k < N; ++k) {
if (k != i) {
long long factor = A[k][i];
if (factor == 0) continue; // Skip if element is already zero
// Update row k of matrix A: R_k <- R_k - factor * R_i
for (int j = i; j < N; ++j) {
A[k][j] = (A[k][j] - factor * A[i][j]) % MOD;
if (A[k][j] < 0) A[k][j] += MOD; // Ensure positive result
}
// Update row k of inverse matrix: R_k <- R_k - factor * R_i
for (int j = 0; j < N; ++j) {
inv[k][j] = (inv[k][j] - factor * inv[i][j]) % MOD;
if (inv[k][j] < 0) inv[k][j] += MOD; // Ensure positive result
}
}
}
}
// Ensure determinant is non-negative
if (det < 0) det += MOD;
return {inv, det};
}
// Reduces matrix A to upper Hessenberg form using similarity transformations.
// Modifies A in place.
void hessenberg(matrix& A, int N) {
for (int k = 0; k < N - 1; ++k) { // Iterate through columns
// Find pivot: first non-zero element below A[k+1, k]
int pivot_row = k + 1;
while(pivot_row < N && A[pivot_row][k] == 0) {
pivot_row++;
}
// If entire column below A[k+1,k] is zero, continue to next column
if (pivot_row == N) continue;
// Swap rows k+1 and pivot_row if needed
if (pivot_row != k + 1) {
std::swap(A[k+1], A[pivot_row]);
// Apply similarity transformation: swap columns k+1 and pivot_row
for(int i=0; i<N; ++i) std::swap(A[i][k+1], A[i][pivot_row]);
}
long long pivot_val = A[k+1][k];
// If pivot value is zero after swap, this column is already reduced. Should not happen if pivot_row found non-zero element.
if (pivot_val == 0) continue;
long long pivot_inv = modInverse(pivot_val);
if (pivot_inv == -1) {
// Should not happen if pivot_val != 0. Indicates potential issue.
continue;
}
// Eliminate elements below A[k+1, k]
for (int i = k + 2; i < N; ++i) {
long long factor = A[i][k]; // Store value before modification
if (factor == 0) continue; // Skip if already zero
long long m = (factor * pivot_inv) % MOD;
// Row operation: R_i <- R_i - m * R_{k+1}
// Updates elements A[i, j] for j >= k
for (int j = k; j < N; ++j) {
A[i][j] = (A[i][j] - m * A[k+1][j]) % MOD;
if (A[i][j] < 0) A[i][j] += MOD;
}
// Column operation: C_{k+1} <- C_{k+1} + m * C_i
// Updates elements A[j, k+1] for all rows j
for (int j = 0; j < N; ++j) {
A[j][k+1] = (A[j][k+1] + m * A[j][i]) % MOD;
if (A[j][k+1] < 0) A[j][k+1] += MOD;
}
// Set A[i][k] explicitly to zero to avoid potential numerical issues (although not needed with exact arithmetic)
A[i][k] = 0;
}
}
}
// Computes the characteristic polynomial det(lambda*I - H) for an upper Hessenberg matrix H
// Uses O(N^3) recurrence relation method.
polynomial char_poly_hessenberg(const matrix& H, int N) {
// p[k] stores the characteristic polynomial of the leading k x k submatrix H_k
std::vector<polynomial> p(N + 1);
p[0] = {1}; // Base case: p_0(lambda) = 1
for (int k = 1; k <= N; ++k) {
// Initialize p_k(lambda) coefficients to zero. Degree k polynomial needs k+1 coeffs.
p[k].resize(k + 1, 0);
// Get H[k-1][k-1] (diagonal element H_kk in 1-based indexing)
long long Hkk = H[k - 1][k - 1];
// Calculate first term: lambda * p_{k-1}(lambda)
// This shifts coefficients of p_{k-1} up by one degree
for (int m = 0; m < k; ++m) { // p[k-1] has size k (coeffs for degree 0 to k-1)
p[k][m + 1] = (p[k][m + 1] + p[k - 1][m]) % MOD;
}
// Calculate second term: -Hkk * p_{k-1}(lambda)
for (int m = 0; m < k; ++m) {
p[k][m] = (p[k][m] - Hkk * p[k - 1][m]) % MOD;
if (p[k][m] < 0) p[k][m] += MOD; // Ensure non-negative
}
// Calculate third term: - sum_{i=1}^{k-1} C_{ik} p_{i-1}(lambda)
// Where C_{ik} = H_{ik} * product_{j=i}^{k-1} H_{j+1, j}
long long current_prod = 1; // Stores product H_{i+1,i}*H_{i+2, i+1}*...*H_{k, k-1}
for (int i = k - 1; i >= 1; --i) { // Iterate i down from k-1 to 1
// Get H_{i+1, i} (0-based H[i][i-1])
long long H_ip1_i_val = (i < N && i-1 >= 0) ? H[i][i - 1] : 0;
current_prod = (current_prod * H_ip1_i_val) % MOD;
// If any subdiagonal H_{j+1, j} is zero, the product becomes zero for all smaller i
if (current_prod == 0) break;
// Get H_{i, k} (0-based H[i-1][k-1])
long long H_i_k_val = (i-1 >= 0 && k-1 < N) ? H[i - 1][k - 1] : 0;
long long C_ik = (H_i_k_val * current_prod) % MOD;
if (C_ik == 0) continue; // Skip if coefficient is zero
// Subtract C_ik * p_{i-1}(lambda) from p_k(lambda)
for (int m = 0; m < i; ++m) { // p[i-1] has size i (coeffs for degree 0 to i-1)
p[k][m] = (p[k][m] - C_ik * p[i - 1][m]) % MOD;
if (p[k][m] < 0) p[k][m] += MOD; // Ensure non-negative
}
}
}
// Return the characteristic polynomial of the full N x N matrix H
return p[N];
}
// Combination helper functions using precomputed factorials and their inverses
std::vector<long long> fact, invFact;
// Precomputes factorials and inverse factorials modulo MOD up to N
void precompute_combinations(int N) {
fact.resize(N + 1);
invFact.resize(N + 1);
fact[0] = 1;
invFact[0] = 1;
for (int i = 1; i <= N; ++i) {
fact[i] = (fact[i - 1] * i) % MOD;
}
// Compute inverse factorial of N! using modular inverse
invFact[N] = modInverse(fact[N]);
// Check for error in modular inverse
if (invFact[N] == -1 && fact[N] != 0) { /* Error: Modular inverse failed */ }
// Compute remaining inverse factorials using relation invFact[i] = invFact[i+1] * (i+1)
for(int i = N - 1; i >= 0; --i) {
invFact[i] = (invFact[i+1] * (i+1)) % MOD;
}
}
// Computes combinations nCr mod MOD using precomputed values
long long nCr_mod(int n, int r) {
if (r < 0 || r > n) return 0; // Base cases for invalid r
// Compute nCr = n! / (r! * (n-r)!) using modular inverses
return (((fact[n] * invFact[r]) % MOD) * invFact[n - r]) % MOD;
}
int main() {
// Faster I/O
std::ios_base::sync_with_stdio(false);
std::cin.tie(NULL);
int N;
std::cin >> N; // Read matrix size
// Read input matrices M0 and M1
matrix M0(N, std::vector<long long>(N));
matrix M1(N, std::vector<long long>(N));
for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) std::cin >> M0[i][j];
for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) std::cin >> M1[i][j];
// Setup random number generation for choosing random shift 'r'
std::mt19937 rng(std::chrono::steady_clock::now().time_since_epoch().count());
std::uniform_int_distribution<long long> distrib(1, MOD - 1); // Choose r from [1, MOD-1]
matrix M0_prime(N, std::vector<long long>(N));
long long det_M0_prime = 0;
long long r = 0; // Initialize r
matrix M0_prime_inv(N, std::vector<long long>(N));
// Retry loop to find a random 'r' such that M0' = M0 + r*M1 is invertible
int retries = 0;
while (true) {
// Fallback strategy: Try r=0 if M0 is invertible after several failed random attempts
if (retries > 10 && retries % 2 == 1) {
auto res0 = inverse_and_det(M0, N); // Pass copy of M0
if (res0.second != 0) { // Check if M0 is invertible
r = 0;
M0_prime_inv = res0.first;
det_M0_prime = res0.second;
break; // Found invertible M0'=M0, exit loop
}
}
// Choose a random non-zero r
r = distrib(rng);
// Compute M0' = M0 + r*M1
matrix rM1 = scale(M1, r, N);
M0_prime = add(M0, rM1, N);
// Check if M0' is invertible
auto result = inverse_and_det(M0_prime, N); // Pass copy of M0_prime
if (result.second != 0) { // If determinant is non-zero
M0_prime_inv = result.first; // Store inverse
det_M0_prime = result.second; // Store determinant
break; // Found invertible M0', exit loop
}
retries++;
// Failsafe after many retries: Check if determinant polynomial is likely zero.
if (retries > 20) {
// Heuristic check: If det(M0), det(M1), det(M0+M1) are all zero, assume poly is zero.
bool M0_singular = (inverse_and_det(M0, N).second == 0);
bool M1_singular = (inverse_and_det(M1, N).second == 0);
matrix M_test = add(M0, M1, N);
bool M0pM1_singular = (inverse_and_det(M_test, N).second == 0);
if (M0_singular && M1_singular && M0pM1_singular) {
// Output all zero coefficients
for (int i = 0; i <= N; ++i) std::cout << 0 << "\n";
return 0; // Exit successfully
}
// If still failing after checks, output zeros as fallback guess (might be wrong)
std::cerr << "Failed to find invertible matrix after many retries. Outputting zeros as fallback." << std::endl;
for (int i = 0; i <= N; ++i) std::cout << 0 << "\n";
return 0; // Exit
}
}
// Compute B = - (M0')^{-1} * M1
matrix B_neg = multiply(M0_prime_inv, M1, N);
matrix B = scale(B_neg, -1, N);
// Reduce B to Hessenberg form H (modifies B in-place)
hessenberg(B, N);
// Compute characteristic polynomial coefficients c_k for det(lambda*I - B)
polynomial c = char_poly_hessenberg(B, N);
// Compute coefficients q_k of Q(x) = det(M0') * det(I - xB)
// Q(x) = det(M0') * x^N * det( (1/x)I - B ) = det(M0') * x^N * chi_B(1/x)
// If chi_B(lambda) = sum c_k lambda^k, then x^N chi_B(1/x) = sum c_k x^{N-k}
// The coefficient of x^k in Q(x) is det(M0') * c_{N-k}
std::vector<long long> q(N + 1);
for (int k = 0; k <= N; ++k) {
// Ensure N-k index is valid for c (which has size N+1 for degrees 0 to N)
if (N - k >= 0 && N-k < c.size()) {
q[k] = (det_M0_prime * c[N - k]) % MOD;
} else {
q[k] = 0; // Should not happen if c has size N+1
}
if (q[k] < 0) q[k] += MOD; // Ensure non-negative
}
// Precompute combinations and powers needed for polynomial shift
precompute_combinations(N);
std::vector<long long> neg_r_powers(N + 1);
neg_r_powers[0] = 1;
long long neg_r = (MOD - r) % MOD; // Calculate -r mod MOD
for (int i = 1; i <= N; ++i) {
neg_r_powers[i] = (neg_r_powers[i - 1] * neg_r) % MOD;
}
// Compute final coefficients a_j of P(x) using relation P(x) = Q(x-r)
// a_j = sum_{k=j}^N q_k * nCr(k, j) * (-r)^(k-j)
std::vector<long long> a(N + 1);
for (int j = 0; j <= N; ++j) { // Compute coefficient a_j for x^j
long long sum = 0;
for (int k = j; k <= N; ++k) { // Sum over terms from Q(x)
long long term = (q[k] * nCr_mod(k, j)) % MOD; // q_k * nCr(k,j)
term = (term * neg_r_powers[k - j]) % MOD; // * (-r)^(k-j)
sum = (sum + term) % MOD;
}
a[j] = sum;
if (a[j] < 0) a[j] += MOD; // Ensure non-negative
}
// Output the coefficients a_0, ..., a_N
for (int i = 0; i <= N; ++i) {
std::cout << a[i] << "\n";
}
return 0; // Success
}