Determinant Modulo a Composite Modulus
(math/matrix/determinant_mod.hpp)
- View this file on GitHub
- Last update: 2026-07-18 19:05:08+09:00
- Include:
#include "math/matrix/determinant_mod.hpp"
Overview
determinant_mod computes the determinant of an integral square matrix modulo
an arbitrary positive modulus. Unlike field Gaussian elimination, it never
divides by a matrix entry, so the modulus may be composite and pivot values do
not need to be invertible.
Pairs of rows are reduced with the extended Euclidean algorithm. Every row transformation is unimodular, which preserves the determinant up to a tracked sign.
Requirements
The matrix entry type must be an integral type no wider than 64 bits. Entries
may be signed or unsigned and are normalized modulo modulus. The modulus must
be positive; modulus == 1 is supported.
Public Interface
template <class Integer>
std::uint64_t determinant_mod(
const Matrix<Integer>& matrix,
std::uint64_t modulus
);
| Function | Description | Complexity |
|---|---|---|
determinant_mod(matrix, modulus) |
Returns det(matrix) in [0, modulus). The matrix must be square and is not mutated. |
$O(N^3 + N^2 \log M)$ time and $O(N^2)$ memory |
Here $N$ is the matrix size and $M$ is the modulus. The determinant of the
empty matrix is 1 % modulus.
Example
#include "math/matrix/determinant_mod.hpp"
#include <cstdint>
int main() {
m1une::matrix::Matrix<long long> matrix(2, 2);
matrix[0][0] = 2;
matrix[0][1] = 3;
matrix[1][0] = 4;
matrix[1][1] = 5;
std::uint64_t determinant = m1une::matrix::determinant_mod(matrix, 12);
return determinant == 10 ? 0 : 1;
}
Depends on
Required by
Verified with
- verify/math/math_algorithms.test.cpp
- verify/math/matrix/matrix.test.cpp
- verify/math/matrix/matrix_det_arbitrary_mod.test.cpp
Code
#ifndef M1UNE_MATRIX_DETERMINANT_MOD_HPP
#define M1UNE_MATRIX_DETERMINANT_MOD_HPP 1
#include <cassert>
#include <cstdint>
#include <type_traits>
#include "matrix.hpp"
namespace m1une {
namespace matrix {
namespace detail {
inline std::uint64_t determinant_multiply_mod(std::uint64_t lhs,
std::uint64_t rhs,
std::uint64_t modulus) {
return std::uint64_t(static_cast<unsigned __int128>(lhs) * rhs % modulus);
}
inline std::uint64_t determinant_subtract_product_mod(
std::uint64_t value, std::uint64_t lhs, std::uint64_t rhs,
std::uint64_t modulus) {
const std::uint64_t product = determinant_multiply_mod(lhs, rhs, modulus);
return std::uint64_t((static_cast<unsigned __int128>(value) + modulus - product) %
modulus);
}
inline std::uint64_t determinant_add_products_mod(
std::uint64_t first_lhs, std::uint64_t first_rhs,
std::uint64_t second_lhs, std::uint64_t second_rhs,
std::uint64_t modulus) {
const std::uint64_t first =
determinant_multiply_mod(first_lhs, first_rhs, modulus);
const std::uint64_t second =
determinant_multiply_mod(second_lhs, second_rhs, modulus);
return std::uint64_t((static_cast<unsigned __int128>(first) + second) % modulus);
}
template <class Integer>
std::uint64_t determinant_normalize(Integer value, std::uint64_t modulus) {
static_assert(std::is_integral_v<Integer>);
static_assert(sizeof(Integer) <= sizeof(std::uint64_t));
if constexpr (std::is_signed_v<Integer>) {
__int128 residue = static_cast<__int128>(value) % static_cast<__int128>(modulus);
if (residue < 0) residue += modulus;
return std::uint64_t(residue);
} else {
return std::uint64_t(static_cast<unsigned __int128>(value) % modulus);
}
}
} // namespace detail
template <class Integer>
std::uint64_t determinant_mod(const Matrix<Integer>& matrix,
std::uint64_t modulus) {
static_assert(std::is_integral_v<Integer>);
assert(matrix.rows() == matrix.cols());
assert(modulus > 0);
const int size = matrix.rows();
if (size == 0) return std::uint64_t(1) % modulus;
Matrix<std::uint64_t> reduced(size, size);
for (int row = 0; row < size; row++) {
for (int col = 0; col < size; col++) {
reduced[row][col] =
detail::determinant_normalize(matrix[row][col], modulus);
}
}
std::uint64_t result = std::uint64_t(1) % modulus;
bool negate = false;
for (int col = 0; col < size; col++) {
int pivot = col;
while (pivot < size && reduced[pivot][col] == 0) pivot++;
if (pivot == size) return 0;
if (pivot != col) {
reduced.swap_rows(pivot, col);
negate = !negate;
}
for (int row = col + 1; row < size; row++) {
std::uint64_t upper = reduced[col][col];
std::uint64_t lower = reduced[row][col];
if (lower == 0) continue;
std::uint64_t upper_upper = 1 % modulus;
std::uint64_t upper_lower = 0;
std::uint64_t lower_upper = 0;
std::uint64_t lower_lower = 1 % modulus;
while (upper != 0 && lower != 0) {
if (upper < lower) {
const std::uint64_t quotient = lower / upper;
lower -= quotient * upper;
lower_upper = detail::determinant_subtract_product_mod(
lower_upper, quotient, upper_upper, modulus);
lower_lower = detail::determinant_subtract_product_mod(
lower_lower, quotient, upper_lower, modulus);
} else {
const std::uint64_t quotient = upper / lower;
upper -= quotient * lower;
upper_upper = detail::determinant_subtract_product_mod(
upper_upper, quotient, lower_upper, modulus);
upper_lower = detail::determinant_subtract_product_mod(
upper_lower, quotient, lower_lower, modulus);
}
}
for (int index = col; index < size; index++) {
const std::uint64_t old_upper = reduced[col][index];
const std::uint64_t old_lower = reduced[row][index];
reduced[col][index] = detail::determinant_add_products_mod(
upper_upper, old_upper, upper_lower, old_lower, modulus);
reduced[row][index] = detail::determinant_add_products_mod(
lower_upper, old_upper, lower_lower, old_lower, modulus);
}
if (upper == 0) {
reduced.swap_rows(col, row);
negate = !negate;
}
}
result = detail::determinant_multiply_mod(
result, reduced[col][col], modulus);
if (result == 0) return 0;
}
return negate ? modulus - result : result;
}
} // namespace matrix
} // namespace m1une
#endif // M1UNE_MATRIX_DETERMINANT_MOD_HPP#line 1 "math/matrix/determinant_mod.hpp"
#include <cassert>
#include <cstdint>
#include <type_traits>
#line 1 "math/matrix/matrix.hpp"
#line 5 "math/matrix/matrix.hpp"
#include <cstddef>
#line 7 "math/matrix/matrix.hpp"
#include <utility>
#include <vector>
namespace m1une {
namespace matrix {
template <class T>
class Matrix {
private:
int _rows;
int _cols;
std::vector<T> _data;
static std::size_t storage_size(int rows, int cols) {
assert(rows >= 0);
assert(cols >= 0);
return std::size_t(rows) * std::size_t(cols);
}
public:
using value_type = T;
Matrix() : _rows(0), _cols(0) {}
Matrix(int rows, int cols, const T& value = T())
: _rows(rows), _cols(cols), _data(storage_size(rows, cols), value) {}
Matrix(int rows, int cols, std::vector<T> values)
: _rows(rows), _cols(cols), _data(std::move(values)) {
assert(rows >= 0);
assert(cols >= 0);
assert(_data.size() == std::size_t(rows) * std::size_t(cols));
}
explicit Matrix(const std::vector<std::vector<T>>& values)
: _rows(int(values.size())), _cols(values.empty() ? 0 : int(values[0].size())),
_data(storage_size(_rows, _cols)) {
for (int row = 0; row < _rows; row++) {
assert(int(values[std::size_t(row)].size()) == _cols);
for (int col = 0; col < _cols; col++) {
(*this)[row][col] = values[std::size_t(row)][std::size_t(col)];
}
}
}
int rows() const {
return _rows;
}
int cols() const {
return _cols;
}
bool empty() const {
return _rows == 0 || _cols == 0;
}
std::vector<T>& data() {
return _data;
}
const std::vector<T>& data() const {
return _data;
}
T* operator[](int row) {
assert(0 <= row && row < _rows);
return _data.data() + std::size_t(row) * std::size_t(_cols);
}
const T* operator[](int row) const {
assert(0 <= row && row < _rows);
return _data.data() + std::size_t(row) * std::size_t(_cols);
}
T& operator()(int row, int col) {
assert(0 <= col && col < _cols);
return (*this)[row][col];
}
const T& operator()(int row, int col) const {
assert(0 <= col && col < _cols);
return (*this)[row][col];
}
static Matrix identity(int size) {
assert(size >= 0);
Matrix result(size, size);
for (int i = 0; i < size; i++) result[i][i] = T(1);
return result;
}
Matrix transposed() const {
Matrix result(_cols, _rows);
for (int row = 0; row < _rows; row++) {
for (int col = 0; col < _cols; col++) {
result[col][row] = (*this)[row][col];
}
}
return result;
}
void swap_rows(int first, int second) {
assert(0 <= first && first < _rows);
assert(0 <= second && second < _rows);
if (first == second) return;
for (int col = 0; col < _cols; col++) {
std::swap((*this)[first][col], (*this)[second][col]);
}
}
Matrix& operator+=(const Matrix& rhs) {
assert(_rows == rhs._rows && _cols == rhs._cols);
for (std::size_t i = 0; i < _data.size(); i++) _data[i] += rhs._data[i];
return *this;
}
Matrix& operator-=(const Matrix& rhs) {
assert(_rows == rhs._rows && _cols == rhs._cols);
for (std::size_t i = 0; i < _data.size(); i++) _data[i] -= rhs._data[i];
return *this;
}
Matrix& operator*=(const T& scalar) {
for (T& value : _data) value *= scalar;
return *this;
}
Matrix& operator/=(const T& scalar) {
for (T& value : _data) value /= scalar;
return *this;
}
Matrix& operator*=(const Matrix& rhs) {
return *this = *this * rhs;
}
Matrix operator+() const {
return *this;
}
Matrix operator-() const {
Matrix result = *this;
for (T& value : result._data) value = T() - value;
return result;
}
friend Matrix operator+(Matrix lhs, const Matrix& rhs) {
return lhs += rhs;
}
friend Matrix operator-(Matrix lhs, const Matrix& rhs) {
return lhs -= rhs;
}
friend Matrix operator*(Matrix lhs, const T& rhs) {
return lhs *= rhs;
}
friend Matrix operator*(const T& lhs, Matrix rhs) {
return rhs *= lhs;
}
friend Matrix operator/(Matrix lhs, const T& rhs) {
return lhs /= rhs;
}
friend Matrix operator*(const Matrix& lhs, const Matrix& rhs) {
assert(lhs._cols == rhs._rows);
Matrix result(lhs._rows, rhs._cols);
for (int row = 0; row < lhs._rows; row++) {
T* output = result[row];
for (int middle = 0; middle < lhs._cols; middle++) {
const T coefficient = lhs[row][middle];
if (coefficient == T()) continue;
const T* input = rhs[middle];
for (int col = 0; col < rhs._cols; col++) {
output[col] += coefficient * input[col];
}
}
}
return result;
}
friend std::vector<T> operator*(const Matrix& lhs, const std::vector<T>& rhs) {
assert(lhs._cols == int(rhs.size()));
std::vector<T> result(std::size_t(lhs._rows));
for (int row = 0; row < lhs._rows; row++) {
T value = T();
for (int col = 0; col < lhs._cols; col++) {
value += lhs[row][col] * rhs[std::size_t(col)];
}
result[std::size_t(row)] = value;
}
return result;
}
friend std::vector<T> operator*(const std::vector<T>& lhs, const Matrix& rhs) {
assert(int(lhs.size()) == rhs._rows);
std::vector<T> result(std::size_t(rhs._cols));
for (int row = 0; row < rhs._rows; row++) {
if (lhs[std::size_t(row)] == T()) continue;
for (int col = 0; col < rhs._cols; col++) {
result[std::size_t(col)] += lhs[std::size_t(row)] * rhs[row][col];
}
}
return result;
}
bool operator==(const Matrix& rhs) const {
return _rows == rhs._rows && _cols == rhs._cols && _data == rhs._data;
}
bool operator!=(const Matrix& rhs) const {
return !(*this == rhs);
}
Matrix pow(std::uint64_t exponent) const {
assert(_rows == _cols);
Matrix result = identity(_rows);
Matrix base = *this;
while (exponent > 0) {
if (exponent & 1) result *= base;
exponent >>= 1;
if (exponent > 0) base *= base;
}
return result;
}
};
} // namespace matrix
} // namespace m1une
#line 9 "math/matrix/determinant_mod.hpp"
namespace m1une {
namespace matrix {
namespace detail {
inline std::uint64_t determinant_multiply_mod(std::uint64_t lhs,
std::uint64_t rhs,
std::uint64_t modulus) {
return std::uint64_t(static_cast<unsigned __int128>(lhs) * rhs % modulus);
}
inline std::uint64_t determinant_subtract_product_mod(
std::uint64_t value, std::uint64_t lhs, std::uint64_t rhs,
std::uint64_t modulus) {
const std::uint64_t product = determinant_multiply_mod(lhs, rhs, modulus);
return std::uint64_t((static_cast<unsigned __int128>(value) + modulus - product) %
modulus);
}
inline std::uint64_t determinant_add_products_mod(
std::uint64_t first_lhs, std::uint64_t first_rhs,
std::uint64_t second_lhs, std::uint64_t second_rhs,
std::uint64_t modulus) {
const std::uint64_t first =
determinant_multiply_mod(first_lhs, first_rhs, modulus);
const std::uint64_t second =
determinant_multiply_mod(second_lhs, second_rhs, modulus);
return std::uint64_t((static_cast<unsigned __int128>(first) + second) % modulus);
}
template <class Integer>
std::uint64_t determinant_normalize(Integer value, std::uint64_t modulus) {
static_assert(std::is_integral_v<Integer>);
static_assert(sizeof(Integer) <= sizeof(std::uint64_t));
if constexpr (std::is_signed_v<Integer>) {
__int128 residue = static_cast<__int128>(value) % static_cast<__int128>(modulus);
if (residue < 0) residue += modulus;
return std::uint64_t(residue);
} else {
return std::uint64_t(static_cast<unsigned __int128>(value) % modulus);
}
}
} // namespace detail
template <class Integer>
std::uint64_t determinant_mod(const Matrix<Integer>& matrix,
std::uint64_t modulus) {
static_assert(std::is_integral_v<Integer>);
assert(matrix.rows() == matrix.cols());
assert(modulus > 0);
const int size = matrix.rows();
if (size == 0) return std::uint64_t(1) % modulus;
Matrix<std::uint64_t> reduced(size, size);
for (int row = 0; row < size; row++) {
for (int col = 0; col < size; col++) {
reduced[row][col] =
detail::determinant_normalize(matrix[row][col], modulus);
}
}
std::uint64_t result = std::uint64_t(1) % modulus;
bool negate = false;
for (int col = 0; col < size; col++) {
int pivot = col;
while (pivot < size && reduced[pivot][col] == 0) pivot++;
if (pivot == size) return 0;
if (pivot != col) {
reduced.swap_rows(pivot, col);
negate = !negate;
}
for (int row = col + 1; row < size; row++) {
std::uint64_t upper = reduced[col][col];
std::uint64_t lower = reduced[row][col];
if (lower == 0) continue;
std::uint64_t upper_upper = 1 % modulus;
std::uint64_t upper_lower = 0;
std::uint64_t lower_upper = 0;
std::uint64_t lower_lower = 1 % modulus;
while (upper != 0 && lower != 0) {
if (upper < lower) {
const std::uint64_t quotient = lower / upper;
lower -= quotient * upper;
lower_upper = detail::determinant_subtract_product_mod(
lower_upper, quotient, upper_upper, modulus);
lower_lower = detail::determinant_subtract_product_mod(
lower_lower, quotient, upper_lower, modulus);
} else {
const std::uint64_t quotient = upper / lower;
upper -= quotient * lower;
upper_upper = detail::determinant_subtract_product_mod(
upper_upper, quotient, lower_upper, modulus);
upper_lower = detail::determinant_subtract_product_mod(
upper_lower, quotient, lower_lower, modulus);
}
}
for (int index = col; index < size; index++) {
const std::uint64_t old_upper = reduced[col][index];
const std::uint64_t old_lower = reduced[row][index];
reduced[col][index] = detail::determinant_add_products_mod(
upper_upper, old_upper, upper_lower, old_lower, modulus);
reduced[row][index] = detail::determinant_add_products_mod(
lower_upper, old_upper, lower_lower, old_lower, modulus);
}
if (upper == 0) {
reduced.swap_rows(col, row);
negate = !negate;
}
}
result = detail::determinant_multiply_mod(
result, reduced[col][col], modulus);
if (result == 0) return 0;
}
return negate ? modulus - result : result;
}
} // namespace matrix
} // namespace m1une