matrix.hpp¶
提供稠密矩阵乘法、单位矩阵和快速幂；适合线性转移重复作用。
\[
\displaystyle C=AB
\]
Complexity: Time: O(nmk) multiplication and O(n^3) square inversion. Space: O(nm) result plus O(n^2) inversion workspace.
AC 记录：inverse_matrix, matrix_product, pow_of_matrix。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(nmk) multiplication and O(n^3) square inversion.
/// Space: O(nm) result plus O(n^2) inversion workspace.
#include "noya/linear_algebra.hpp"
#include <cassert>
#include <cstdint>
#include <optional>
#include <utility>
#include <vector>
namespace noya {
template <class T> using matrix = std::vector<std::vector<T>>;
template <class T> matrix<T> identity_matrix(int n) {
assert(n >= 0);
matrix<T> res(n, std::vector<T>(n));
for (int idx = 0; idx < n; idx++) {
res[idx][idx] = T(1);
}
return res;
}
/// @brief Dense matrix multiplication in O(nmk).
template <class T>
matrix<T> matrix_multiply(const matrix<T> &a, const matrix<T> &b) {
int rs = int(a.size());
int mid = rs == 0 ? int(b.size()) : int(a[0].size());
for (const auto &row : a) {
assert(int(row.size()) == mid);
}
assert(int(b.size()) == mid);
int cs = mid == 0 ? 0 : int(b[0].size());
for (const auto &row : b) {
assert(int(row.size()) == cs);
}
matrix<T> res(rs, std::vector<T>(cs));
for (int row = 0; row < rs; row++) {
for (int idx = 0; idx < mid; idx++) {
for (int col = 0; col < cs; col++) {
res[row][col] += a[row][idx] * b[idx][col];
}
}
}
return res;
}
template <class T> matrix<T> matrix_power(matrix<T> vl, std::uint64_t exp) {
int n = int(vl.size());
for (const auto &row : vl) {
assert(int(row.size()) == n);
}
matrix<T> res = identity_matrix<T>(n);
while (exp > 0) {
if (exp & 1) {
res = matrix_multiply(res, vl);
}
vl = matrix_multiply(vl, vl);
exp >>= 1;
}
return res;
}
/// @brief Multiply matrices of static modular integers, transposing the right
/// operand and reducing one wide accumulator per output entry.
template <class Mint>
matrix<Mint> matrix_multiply_mod(const matrix<Mint> &a, const matrix<Mint> &b) {
int rs = int(a.size());
int mid = rs == 0 ? int(b.size()) : int(a[0].size());
int cs = mid == 0 ? 0 : int(b[0].size());
matrix<unsigned int> tra(cs, std::vector<unsigned int>(mid));
for (int idx = 0; idx < mid; idx++) {
for (int col = 0; col < cs; col++) {
tra[col][idx] = b[idx][col].val();
}
}
matrix<Mint> res(rs, std::vector<Mint>(cs));
for (int row = 0; row < rs; row++) {
for (int col = 0; col < cs; col++) {
unsigned __int128 sum = 0;
for (int idx = 0; idx < mid; idx++) {
sum += static_cast<std::uint64_t>(a[row][idx].val()) * tra[col][idx];
}
res[row][col] = Mint::raw(unsigned(sum % Mint::mod()));
}
}
return res;
}
/// @brief Binary exponentiation specialized for dense modular matrices.
template <class Mint>
matrix<Mint> matrix_power_mod(matrix<Mint> vl, std::uint64_t exp) {
int n = int(vl.size());
matrix<Mint> res = identity_matrix<Mint>(n);
while (exp > 0) {
if (exp & 1) {
res = matrix_multiply_mod(res, vl);
}
exp >>= 1;
if (exp != 0) {
vl = matrix_multiply_mod(vl, vl);
}
}
return res;
}
/// @brief Invert a square matrix over a field, returning nullopt when singular.
template <class T, class IsZero = exact_zero<T>>
std::optional<matrix<T>> matrix_inverse(matrix<T> vl, IsZero iz = {}) {
int n = int(vl.size());
for (const auto &row : vl) {
assert(int(row.size()) == n);
}
matrix<T> inv = identity_matrix<T>(n);
for (int col = 0; col < n; col++) {
int piv = col;
while (piv < n && iz(vl[piv][col])) {
piv++;
}
if (piv == n) {
return std::nullopt;
}
std::swap(vl[piv], vl[col]);
std::swap(inv[piv], inv[col]);
T scl = T(1) / vl[col][col];
for (int idx = 0; idx < n; idx++) {
vl[col][idx] *= scl;
inv[col][idx] *= scl;
}
for (int row = 0; row < n; row++) {
if (row == col || iz(vl[row][col])) {
continue;
}
T r = vl[row][col];
for (int idx = 0; idx < n; idx++) {
vl[row][idx] -= r * vl[col][idx];
inv[row][idx] -= r * inv[col][idx];
}
}
}
return inv;
}
} // namespace noya
#ifndef NOYA_MATRIX_HPP
#define NOYA_MATRIX_HPP 1
/// @complexity Time: O(nmk) multiplication and O(n^3) square inversion.
/// Space: O(nm) result plus O(n^2) inversion workspace.
#include "noya/linear_algebra.hpp"
#include <cassert>
#include <cstdint>
#include <optional>
#include <utility>
#include <vector>
namespace noya {
template <class T> using matrix = std::vector<std::vector<T>>;
template <class T> matrix<T> identity_matrix(int n) {
assert(n >= 0);
matrix<T> res(n, std::vector<T>(n));
for (int idx = 0; idx < n; idx++) {
res[idx][idx] = T(1);
}
return res;
}
/// @brief Dense matrix multiplication in O(nmk).
template <class T>
matrix<T> matrix_multiply(const matrix<T> &a, const matrix<T> &b) {
int rs = int(a.size());
int mid = rs == 0 ? int(b.size()) : int(a[0].size());
for (const auto &row : a) {
assert(int(row.size()) == mid);
}
assert(int(b.size()) == mid);
int cs = mid == 0 ? 0 : int(b[0].size());
for (const auto &row : b) {
assert(int(row.size()) == cs);
}
matrix<T> res(rs, std::vector<T>(cs));
for (int row = 0; row < rs; row++) {
for (int idx = 0; idx < mid; idx++) {
for (int col = 0; col < cs; col++) {
res[row][col] += a[row][idx] * b[idx][col];
}
}
}
return res;
}
template <class T> matrix<T> matrix_power(matrix<T> vl, std::uint64_t exp) {
int n = int(vl.size());
for (const auto &row : vl) {
assert(int(row.size()) == n);
}
matrix<T> res = identity_matrix<T>(n);
while (exp > 0) {
if (exp & 1) {
res = matrix_multiply(res, vl);
}
vl = matrix_multiply(vl, vl);
exp >>= 1;
}
return res;
}
/// @brief Multiply matrices of static modular integers, transposing the right
/// operand and reducing one wide accumulator per output entry.
template <class Mint>
matrix<Mint> matrix_multiply_mod(const matrix<Mint> &a, const matrix<Mint> &b) {
int rs = int(a.size());
int mid = rs == 0 ? int(b.size()) : int(a[0].size());
int cs = mid == 0 ? 0 : int(b[0].size());
matrix<unsigned int> tra(cs, std::vector<unsigned int>(mid));
for (int idx = 0; idx < mid; idx++) {
for (int col = 0; col < cs; col++) {
tra[col][idx] = b[idx][col].val();
}
}
matrix<Mint> res(rs, std::vector<Mint>(cs));
for (int row = 0; row < rs; row++) {
for (int col = 0; col < cs; col++) {
unsigned __int128 sum = 0;
for (int idx = 0; idx < mid; idx++) {
sum += static_cast<std::uint64_t>(a[row][idx].val()) * tra[col][idx];
}
res[row][col] = Mint::raw(unsigned(sum % Mint::mod()));
}
}
return res;
}
/// @brief Binary exponentiation specialized for dense modular matrices.
template <class Mint>
matrix<Mint> matrix_power_mod(matrix<Mint> vl, std::uint64_t exp) {
int n = int(vl.size());
matrix<Mint> res = identity_matrix<Mint>(n);
while (exp > 0) {
if (exp & 1) {
res = matrix_multiply_mod(res, vl);
}
exp >>= 1;
if (exp != 0) {
vl = matrix_multiply_mod(vl, vl);
}
}
return res;
}
/// @brief Invert a square matrix over a field, returning nullopt when singular.
template <class T, class IsZero = exact_zero<T>>
std::optional<matrix<T>> matrix_inverse(matrix<T> vl, IsZero iz = {}) {
int n = int(vl.size());
for (const auto &row : vl) {
assert(int(row.size()) == n);
}
matrix<T> inv = identity_matrix<T>(n);
for (int col = 0; col < n; col++) {
int piv = col;
while (piv < n && iz(vl[piv][col])) {
piv++;
}
if (piv == n) {
return std::nullopt;
}
std::swap(vl[piv], vl[col]);
std::swap(inv[piv], inv[col]);
T scl = T(1) / vl[col][col];
for (int idx = 0; idx < n; idx++) {
vl[col][idx] *= scl;
inv[col][idx] *= scl;
}
for (int row = 0; row < n; row++) {
if (row == col || iz(vl[row][col])) {
continue;
}
T r = vl[row][col];
for (int idx = 0; idx < n; idx++) {
vl[row][idx] -= r * vl[col][idx];
inv[row][idx] -= r * inv[col][idx];
}
}
}
return inv;
}
} // namespace noya
#endif // NOYA_MATRIX_HPP
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <optional>
#include <utility>
#include <vector>
/// @complexity Time: O(nmk) multiplication and O(n^3) square inversion.
/// Space: O(nm) result plus O(n^2) inversion workspace.
/// @complexity Time: O(rows * columns * min(rows,columns)) elimination; O(n^3) square determinant/inverse.
/// Space: O(rows * columns).
namespace noya {
/// @brief Exact zero predicate used by elimination routines by default.
template <class T> struct exact_zero {
bool operator()(const T &val) const { return val == T{}; }
};
/// @brief Consistency flag, one solution, nullspace basis, and pivot columns.
template <class T> struct linear_system_solution {
bool ok = false;
std::vector<T> sol;
std::vector<std::vector<T>> ker;
std::vector<int> pc;
};
namespace linear_algebra_internal {
template <class T> int column_count(const std::vector<std::vector<T>> &mat) {
if (mat.empty()) {
return 0;
}
int cs = int(mat[0].size());
for (const auto &row : mat) {
assert(int(row.size()) == cs);
}
return cs;
}
} // namespace linear_algebra_internal
/// @brief Compute matrix rank over a field.
template <class T, class IsZero = exact_zero<T>>
int matrix_rank(std::vector<std::vector<T>> mat, IsZero iz = {}) {
int rs = int(mat.size());
int cs = linear_algebra_internal::column_count(mat);
int ran = 0;
for (int col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && iz(mat[piv][col])) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
for (int row = ran + 1; row < rs; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = mat[row][col] / mat[ran][col];
for (int j = col; j < cs; j++) {
mat[row][j] -= r * mat[ran][j];
}
}
ran++;
}
return ran;
}
/// @brief Compute the determinant of a square matrix over a field.
template <class T, class IsZero = exact_zero<T>>
T determinant(std::vector<std::vector<T>> mat, IsZero iz = {}) {
int n = int(mat.size());
assert(linear_algebra_internal::column_count(mat) == n);
T res = T(1);
for (int col = 0; col < n; col++) {
int piv = col;
while (piv < n && iz(mat[piv][col])) {
piv++;
}
if (piv == n) {
return T{};
}
if (piv != col) {
std::swap(mat[piv], mat[col]);
res = -res;
}
T pv = mat[col][col];
res *= pv;
for (int row = col + 1; row < n; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = mat[row][col] / pv;
for (int j = col; j < n; j++) {
mat[row][j] -= r * mat[col][j];
}
}
}
return res;
}
/// @brief Solve A*x=b and return one solution plus a basis of the nullspace.
template <class T, class IsZero = exact_zero<T>>
linear_system_solution<T> solve_linear(std::vector<std::vector<T>> mat,
std::vector<T> rhs, IsZero iz = {}) {
int rs = int(mat.size());
assert(int(rhs.size()) == rs);
int cs = linear_algebra_internal::column_count(mat);
std::vector<int> pc;
int ran = 0;
for (int col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && iz(mat[piv][col])) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
std::swap(rhs[piv], rhs[ran]);
T inv = T(1) / mat[ran][col];
for (int j = col; j < cs; j++) {
mat[ran][j] *= inv;
}
rhs[ran] *= inv;
for (int row = 0; row < rs; row++) {
if (row == ran || iz(mat[row][col])) {
continue;
}
T r = mat[row][col];
for (int j = col; j < cs; j++) {
mat[row][j] -= r * mat[ran][j];
}
rhs[row] -= r * rhs[ran];
}
pc.push_back(col);
ran++;
}
for (int row = ran; row < rs; row++) {
bool az = true;
for (int col = 0; col < cs; col++) {
az &= iz(mat[row][col]);
}
if (az && !iz(rhs[row])) {
return {};
}
}
linear_system_solution<T> res;
res.ok = true;
res.sol.assign(cs, T{});
res.pc = pc;
std::vector<bool> ip(cs);
for (int row = 0; row < ran; row++) {
int col = pc[row];
ip[col] = true;
res.sol[col] = rhs[row];
}
for (int fc = 0; fc < cs; fc++) {
if (ip[fc]) {
continue;
}
std::vector<T> bv(cs, T{});
bv[fc] = T(1);
for (int row = 0; row < ran; row++) {
bv[pc[row]] = -mat[row][fc];
}
res.ker.push_back(std::move(bv));
}
return res;
}
} // namespace noya
namespace noya {
template <class T> using matrix = std::vector<std::vector<T>>;
template <class T> matrix<T> identity_matrix(int n) {
assert(n >= 0);
matrix<T> res(n, std::vector<T>(n));
for (int idx = 0; idx < n; idx++) {
res[idx][idx] = T(1);
}
return res;
}
/// @brief Dense matrix multiplication in O(nmk).
template <class T>
matrix<T> matrix_multiply(const matrix<T> &a, const matrix<T> &b) {
int rs = int(a.size());
int mid = rs == 0 ? int(b.size()) : int(a[0].size());
for (const auto &row : a) {
assert(int(row.size()) == mid);
}
assert(int(b.size()) == mid);
int cs = mid == 0 ? 0 : int(b[0].size());
for (const auto &row : b) {
assert(int(row.size()) == cs);
}
matrix<T> res(rs, std::vector<T>(cs));
for (int row = 0; row < rs; row++) {
for (int idx = 0; idx < mid; idx++) {
for (int col = 0; col < cs; col++) {
res[row][col] += a[row][idx] * b[idx][col];
}
}
}
return res;
}
template <class T> matrix<T> matrix_power(matrix<T> vl, std::uint64_t exp) {
int n = int(vl.size());
for (const auto &row : vl) {
assert(int(row.size()) == n);
}
matrix<T> res = identity_matrix<T>(n);
while (exp > 0) {
if (exp & 1) {
res = matrix_multiply(res, vl);
}
vl = matrix_multiply(vl, vl);
exp >>= 1;
}
return res;
}
/// @brief Multiply matrices of static modular integers, transposing the right
/// operand and reducing one wide accumulator per output entry.
template <class Mint>
matrix<Mint> matrix_multiply_mod(const matrix<Mint> &a, const matrix<Mint> &b) {
int rs = int(a.size());
int mid = rs == 0 ? int(b.size()) : int(a[0].size());
int cs = mid == 0 ? 0 : int(b[0].size());
matrix<unsigned int> tra(cs, std::vector<unsigned int>(mid));
for (int idx = 0; idx < mid; idx++) {
for (int col = 0; col < cs; col++) {
tra[col][idx] = b[idx][col].val();
}
}
matrix<Mint> res(rs, std::vector<Mint>(cs));
for (int row = 0; row < rs; row++) {
for (int col = 0; col < cs; col++) {
unsigned __int128 sum = 0;
for (int idx = 0; idx < mid; idx++) {
sum += static_cast<std::uint64_t>(a[row][idx].val()) * tra[col][idx];
}
res[row][col] = Mint::raw(unsigned(sum % Mint::mod()));
}
}
return res;
}
/// @brief Binary exponentiation specialized for dense modular matrices.
template <class Mint>
matrix<Mint> matrix_power_mod(matrix<Mint> vl, std::uint64_t exp) {
int n = int(vl.size());
matrix<Mint> res = identity_matrix<Mint>(n);
while (exp > 0) {
if (exp & 1) {
res = matrix_multiply_mod(res, vl);
}
exp >>= 1;
if (exp != 0) {
vl = matrix_multiply_mod(vl, vl);
}
}
return res;
}
/// @brief Invert a square matrix over a field, returning nullopt when singular.
template <class T, class IsZero = exact_zero<T>>
std::optional<matrix<T>> matrix_inverse(matrix<T> vl, IsZero iz = {}) {
int n = int(vl.size());
for (const auto &row : vl) {
assert(int(row.size()) == n);
}
matrix<T> inv = identity_matrix<T>(n);
for (int col = 0; col < n; col++) {
int piv = col;
while (piv < n && iz(vl[piv][col])) {
piv++;
}
if (piv == n) {
return std::nullopt;
}
std::swap(vl[piv], vl[col]);
std::swap(inv[piv], inv[col]);
T scl = T(1) / vl[col][col];
for (int idx = 0; idx < n; idx++) {
vl[col][idx] *= scl;
inv[col][idx] *= scl;
}
for (int row = 0; row < n; row++) {
if (row == col || iz(vl[row][col])) {
continue;
}
T r = vl[row][col];
for (int idx = 0; idx < n; idx++) {
vl[row][idx] -= r * vl[col][idx];
inv[row][idx] -= r * inv[col][idx];
}
}
}
return inv;
}
} // namespace noya