determinant_mod.hpp

SECTIONMath INCLUDEnoya/determinant_mod.hpp

Determinant of a square integer matrix modulo any positive modulus in O(n^3 log modulus), without requiring modular inverses.

Verified by matrix_det_arbitrary_mod.

\[ \displaystyle \det(A)\bmod m \]

Implementation

View on GitHub

SourceExpanded

#ifndef NOYA_DETERMINANT_MOD_HPP
#define NOYA_DETERMINANT_MOD_HPP 1

/// @complexity Time: O(n^3 log m) arithmetic bit operations.
/// Space: O(n^2).

#include <cassert>
#include <cstdint>
#include <utility>
#include <vector>

namespace noya {

/// @brief Determinant of a square integer matrix modulo any positive modulus
/// in O(n^3 log modulus), without requiring modular inverses.
inline std::int64_t
determinant_mod(std::vector<std::vector<std::int64_t>> matrix,
                std::int64_t modulus) {
  assert(modulus > 0);
  const int n = int(matrix.size());
  for (const auto &row : matrix) {
    assert(int(row.size()) == n);
  }
  for (auto &row : matrix) {
    for (std::int64_t &value : row) {
      value %= modulus;
      if (value < 0) {
        value += modulus;
      }
    }
  }

  std::int64_t result = 1 % modulus;
  for (int column = 0; column < n; column++) {
    int pivot = column;
    while (pivot < n && matrix[pivot][column] == 0) {
      pivot++;
    }
    if (pivot == n) {
      return 0;
    }
    if (pivot != column) {
      std::swap(matrix[pivot], matrix[column]);
      result = (modulus - result) % modulus;
    }

    for (int row = column + 1; row < n; row++) {
      while (matrix[row][column] != 0) {
        std::int64_t quotient = matrix[column][column] / matrix[row][column];
        for (int next = column; next < n; next++) {
          __int128 value =
              matrix[column][next] - __int128(quotient) * matrix[row][next];
          matrix[column][next] = std::int64_t(value % modulus);
          if (matrix[column][next] < 0) {
            matrix[column][next] += modulus;
          }
        }
        std::swap(matrix[column], matrix[row]);
        result = (modulus - result) % modulus;
      }
    }
    result = std::int64_t(__int128(result) * matrix[column][column] % modulus);
  }
  return result;
}

} // namespace noya

#endif // NOYA_DETERMINANT_MOD_HPP

#include <cassert>
#include <cstdint>
#include <utility>
#include <vector>

/// @complexity Time: O(n^3 log m) arithmetic bit operations.
/// Space: O(n^2).

namespace noya {

/// @brief Determinant of a square integer matrix modulo any positive modulus
/// in O(n^3 log modulus), without requiring modular inverses.
inline std::int64_t
determinant_mod(std::vector<std::vector<std::int64_t>> matrix,
                std::int64_t modulus) {
  assert(modulus > 0);
  const int n = int(matrix.size());
  for (const auto &row : matrix) {
    assert(int(row.size()) == n);
  }
  for (auto &row : matrix) {
    for (std::int64_t &value : row) {
      value %= modulus;
      if (value < 0) {
        value += modulus;
      }
    }
  }

  std::int64_t result = 1 % modulus;
  for (int column = 0; column < n; column++) {
    int pivot = column;
    while (pivot < n && matrix[pivot][column] == 0) {
      pivot++;
    }
    if (pivot == n) {
      return 0;
    }
    if (pivot != column) {
      std::swap(matrix[pivot], matrix[column]);
      result = (modulus - result) % modulus;
    }

    for (int row = column + 1; row < n; row++) {
      while (matrix[row][column] != 0) {
        std::int64_t quotient = matrix[column][column] / matrix[row][column];
        for (int next = column; next < n; next++) {
          __int128 value =
              matrix[column][next] - __int128(quotient) * matrix[row][next];
          matrix[column][next] = std::int64_t(value % modulus);
          if (matrix[column][next] < 0) {
            matrix[column][next] += modulus;
          }
        }
        std::swap(matrix[column], matrix[row]);
        result = (modulus - result) % modulus;
      }
    }
    result = std::int64_t(__int128(result) * matrix[column][column] % modulus);
  }
  return result;
}

} // namespace noya