adjugate_matrix.hpp

SECTIONMath INCLUDEnoya/adjugate_matrix.hpp

计算矩阵的伴随矩阵，即使矩阵奇异也能得到所有代数余子式。

\[ \displaystyle \operatorname{adj}(A)A=\det(A)I \]

Complexity: Time: O(n^3). Space: O(n^2).

AC 记录：adjugate_matrix。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

当前头文件，省略 include guard；依赖见 #include。

/// @complexity Time: O(n^3).
/// Space: O(n^2).

#include <algorithm>
#include <cassert>
#include <vector>

namespace noya {

/// @brief Compute the adjugate of a square matrix over a field.  Gauss--Jordan
/// elimination tracks the row-operation matrix E with `E A = R`.  For full
/// rank, `E = A^{-1}`.  For rank n-1, `adj(R)` is an outer product of the
/// unique right-null vector and the final zero-row selector; multiplying by E
/// recovers `adj(A)`.  Rank at most n-2 makes every (n-1)-minor zero.
template <class T>
std::vector<std::vector<T>> adjugate_matrix(std::vector<std::vector<T>> mat) {
  int sz = int(mat.size());
  for (const auto &row : mat) {
    assert(int(row.size()) == sz);
  }
  std::vector<std::vector<T>> ops(sz, std::vector<T>(sz, T{}));
  for (int row = 0; row < sz; row++) {
    ops[row][row] = T(1);
  }

  T det = T(1);
  std::vector<int> pc;
  int ran = 0;
  for (int col = 0; col < sz && ran < sz; col++) {
    int piv = ran;
    while (piv < sz && mat[piv][col] == T{}) {
      piv++;
    }
    if (piv == sz) {
      continue;
    }
    if (piv != ran) {
      std::swap(mat[piv], mat[ran]);
      std::swap(ops[piv], ops[ran]);
      det = -det;
    }
    T pv = mat[ran][col];
    det *= pv;
    T inv = T(1) / pv;
    for (int nxt = col; nxt < sz; nxt++) {
      mat[ran][nxt] *= inv;
    }
    for (int nxt = 0; nxt < sz; nxt++) {
      ops[ran][nxt] *= inv;
    }
    for (int row = 0; row < sz; row++) {
      if (row == ran || mat[row][col] == T{}) {
        continue;
      }
      T r = mat[row][col];
      for (int nxt = col; nxt < sz; nxt++) {
        mat[row][nxt] -= r * mat[ran][nxt];
      }
      for (int nxt = 0; nxt < sz; nxt++) {
        ops[row][nxt] -= r * ops[ran][nxt];
      }
    }
    pc.push_back(col);
    ran++;
  }

  std::vector<std::vector<T>> res(sz, std::vector<T>(sz, T{}));
  if (ran == sz) {
    for (int row = 0; row < sz; row++) {
      for (int col = 0; col < sz; col++) {
        res[row][col] = det * ops[row][col];
      }
    }
    return res;
  }
  if (ran + 1 != sz) {
    return res;
  }

  std::vector<bool> ip(sz);
  for (int col : pc) {
    ip[col] = true;
  }
  int fc = 0;
  while (ip[fc]) {
    fc++;
  }
  std::vector<T> rv(sz, T{});
  rv[fc] = T(1);
  for (int row = 0; row < ran; row++) {
    rv[pc[row]] = -mat[row][fc];
  }
  T scl = det;
  if ((ran + fc) % 2 != 0) {
    scl = -scl;
  }
  for (int row = 0; row < sz; row++) {
    for (int col = 0; col < sz; col++) {
      res[row][col] = scl * rv[row] * ops[ran][col];
    }
  }
  return res;
}

} // namespace noya

#ifndef NOYA_ADJUGATE_MATRIX_HPP
#define NOYA_ADJUGATE_MATRIX_HPP 1

/// @complexity Time: O(n^3).
/// Space: O(n^2).

#include <algorithm>
#include <cassert>
#include <vector>

namespace noya {

/// @brief Compute the adjugate of a square matrix over a field.  Gauss--Jordan
/// elimination tracks the row-operation matrix E with `E A = R`.  For full
/// rank, `E = A^{-1}`.  For rank n-1, `adj(R)` is an outer product of the
/// unique right-null vector and the final zero-row selector; multiplying by E
/// recovers `adj(A)`.  Rank at most n-2 makes every (n-1)-minor zero.
template <class T>
std::vector<std::vector<T>> adjugate_matrix(std::vector<std::vector<T>> mat) {
  int sz = int(mat.size());
  for (const auto &row : mat) {
    assert(int(row.size()) == sz);
  }
  std::vector<std::vector<T>> ops(sz, std::vector<T>(sz, T{}));
  for (int row = 0; row < sz; row++) {
    ops[row][row] = T(1);
  }

  T det = T(1);
  std::vector<int> pc;
  int ran = 0;
  for (int col = 0; col < sz && ran < sz; col++) {
    int piv = ran;
    while (piv < sz && mat[piv][col] == T{}) {
      piv++;
    }
    if (piv == sz) {
      continue;
    }
    if (piv != ran) {
      std::swap(mat[piv], mat[ran]);
      std::swap(ops[piv], ops[ran]);
      det = -det;
    }
    T pv = mat[ran][col];
    det *= pv;
    T inv = T(1) / pv;
    for (int nxt = col; nxt < sz; nxt++) {
      mat[ran][nxt] *= inv;
    }
    for (int nxt = 0; nxt < sz; nxt++) {
      ops[ran][nxt] *= inv;
    }
    for (int row = 0; row < sz; row++) {
      if (row == ran || mat[row][col] == T{}) {
        continue;
      }
      T r = mat[row][col];
      for (int nxt = col; nxt < sz; nxt++) {
        mat[row][nxt] -= r * mat[ran][nxt];
      }
      for (int nxt = 0; nxt < sz; nxt++) {
        ops[row][nxt] -= r * ops[ran][nxt];
      }
    }
    pc.push_back(col);
    ran++;
  }

  std::vector<std::vector<T>> res(sz, std::vector<T>(sz, T{}));
  if (ran == sz) {
    for (int row = 0; row < sz; row++) {
      for (int col = 0; col < sz; col++) {
        res[row][col] = det * ops[row][col];
      }
    }
    return res;
  }
  if (ran + 1 != sz) {
    return res;
  }

  std::vector<bool> ip(sz);
  for (int col : pc) {
    ip[col] = true;
  }
  int fc = 0;
  while (ip[fc]) {
    fc++;
  }
  std::vector<T> rv(sz, T{});
  rv[fc] = T(1);
  for (int row = 0; row < ran; row++) {
    rv[pc[row]] = -mat[row][fc];
  }
  T scl = det;
  if ((ran + fc) % 2 != 0) {
    scl = -scl;
  }
  for (int row = 0; row < sz; row++) {
    for (int col = 0; col < sz; col++) {
      res[row][col] = scl * rv[row] * ops[ran][col];
    }
  }
  return res;
}

} // namespace noya

#endif // NOYA_ADJUGATE_MATRIX_HPP

#include <algorithm>
#include <cassert>
#include <vector>

/// @complexity Time: O(n^3).
/// Space: O(n^2).

namespace noya {

/// @brief Compute the adjugate of a square matrix over a field.  Gauss--Jordan
/// elimination tracks the row-operation matrix E with `E A = R`.  For full
/// rank, `E = A^{-1}`.  For rank n-1, `adj(R)` is an outer product of the
/// unique right-null vector and the final zero-row selector; multiplying by E
/// recovers `adj(A)`.  Rank at most n-2 makes every (n-1)-minor zero.
template <class T>
std::vector<std::vector<T>> adjugate_matrix(std::vector<std::vector<T>> mat) {
  int sz = int(mat.size());
  for (const auto &row : mat) {
    assert(int(row.size()) == sz);
  }
  std::vector<std::vector<T>> ops(sz, std::vector<T>(sz, T{}));
  for (int row = 0; row < sz; row++) {
    ops[row][row] = T(1);
  }

  T det = T(1);
  std::vector<int> pc;
  int ran = 0;
  for (int col = 0; col < sz && ran < sz; col++) {
    int piv = ran;
    while (piv < sz && mat[piv][col] == T{}) {
      piv++;
    }
    if (piv == sz) {
      continue;
    }
    if (piv != ran) {
      std::swap(mat[piv], mat[ran]);
      std::swap(ops[piv], ops[ran]);
      det = -det;
    }
    T pv = mat[ran][col];
    det *= pv;
    T inv = T(1) / pv;
    for (int nxt = col; nxt < sz; nxt++) {
      mat[ran][nxt] *= inv;
    }
    for (int nxt = 0; nxt < sz; nxt++) {
      ops[ran][nxt] *= inv;
    }
    for (int row = 0; row < sz; row++) {
      if (row == ran || mat[row][col] == T{}) {
        continue;
      }
      T r = mat[row][col];
      for (int nxt = col; nxt < sz; nxt++) {
        mat[row][nxt] -= r * mat[ran][nxt];
      }
      for (int nxt = 0; nxt < sz; nxt++) {
        ops[row][nxt] -= r * ops[ran][nxt];
      }
    }
    pc.push_back(col);
    ran++;
  }

  std::vector<std::vector<T>> res(sz, std::vector<T>(sz, T{}));
  if (ran == sz) {
    for (int row = 0; row < sz; row++) {
      for (int col = 0; col < sz; col++) {
        res[row][col] = det * ops[row][col];
      }
    }
    return res;
  }
  if (ran + 1 != sz) {
    return res;
  }

  std::vector<bool> ip(sz);
  for (int col : pc) {
    ip[col] = true;
  }
  int fc = 0;
  while (ip[fc]) {
    fc++;
  }
  std::vector<T> rv(sz, T{});
  rv[fc] = T(1);
  for (int row = 0; row < ran; row++) {
    rv[pc[row]] = -mat[row][fc];
  }
  T scl = det;
  if ((ran + fc) % 2 != 0) {
    scl = -scl;
  }
  for (int row = 0; row < sz; row++) {
    for (int col = 0; col < sz; col++) {
      res[row][col] = scl * rv[row] * ops[ran][col];
    }
  }
  return res;
}

} // namespace noya