Adjugate Matrix
(math/matrix/adjugate.hpp)

View this file on GitHub
Last update: 2026-07-18 19:05:08+09:00
Include: #include "math/matrix/adjugate.hpp"

Overview

For a square matrix $A$, adjugate returns the transpose of its cofactor matrix. The result satisfies

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

The implementation performs one rank-revealing Gauss–Jordan elimination. It uses the inverse formula at full rank, returns zero below rank $N-1$, and reconstructs the rank-one adjugate from left and right null vectors at rank $N-1$. It is deterministic and does not compute individual minors.

Requirements

T must be a field type supporting construction from 0 and 1, equality, addition, subtraction, multiplication, and division by nonzero values. The input matrix must be square.

Public Interface

template <class T>
Matrix<T> adjugate(Matrix<T> matrix);

Function	Description	Complexity
`adjugate(matrix)`	Returns the adjugate. The argument is copied and the caller’s matrix is unchanged.	$O(N^3)$ time and $O(N^2)$ memory

The empty matrix produces an empty matrix. The adjugate of every 1 x 1 matrix, including the zero matrix, is the matrix whose only entry is 1.

Example

#include "math/matrix/adjugate.hpp"
#include "math/modint.hpp"

int main() {
    using Mint = m1une::math::modint998244353;
    m1une::matrix::Matrix<Mint> matrix(2, 2);
    matrix[0][0] = 1;
    matrix[0][1] = 2;
    matrix[1][0] = 3;
    matrix[1][1] = 4;
    auto result = m1une::matrix::adjugate(matrix);
    return result[0][0] == Mint(4) && result[0][1] == Mint(0) - Mint(2) &&
                   result[1][0] == Mint(0) - Mint(3) && result[1][1] == Mint(1)
               ? 0
               : 1;
}

Depends on

Dense Matrix (math/matrix/matrix.hpp)

Required by

Verified with

Code

#ifndef M1UNE_MATRIX_ADJUGATE_HPP
#define M1UNE_MATRIX_ADJUGATE_HPP 1

#include <cassert>
#include <vector>

#include "matrix.hpp"

namespace m1une {
namespace matrix {

template <class T>
Matrix<T> adjugate(Matrix<T> matrix) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    Matrix<T> augmented(size, size * 2);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            augmented[row][col] = matrix[row][col];
        }
        augmented[row][size + row] = T(1);
    }

    std::vector<int> pivot_columns;
    T pivot_product = T(1);
    bool negate = false;
    for (int col = 0; col < size && int(pivot_columns.size()) < size; col++) {
        const int pivot_row = int(pivot_columns.size());
        int pivot = pivot_row;
        while (pivot < size && augmented[pivot][col] == T()) pivot++;
        if (pivot == size) continue;
        if (pivot != pivot_row) {
            augmented.swap_rows(pivot, pivot_row);
            negate = !negate;
        }

        const T pivot_value = augmented[pivot_row][col];
        pivot_product *= pivot_value;
        const T inverse_pivot = T(1) / pivot_value;
        for (int index = col; index < size; index++) {
            augmented[pivot_row][index] *= inverse_pivot;
        }
        for (int index = size; index < size * 2; index++) {
            augmented[pivot_row][index] *= inverse_pivot;
        }

        for (int row = 0; row < size; row++) {
            if (row == pivot_row || augmented[row][col] == T()) continue;
            const T factor = augmented[row][col];
            augmented[row][col] = T();
            for (int index = col + 1; index < size; index++) {
                augmented[row][index] -= factor * augmented[pivot_row][index];
            }
            for (int index = size; index < size * 2; index++) {
                augmented[row][index] -= factor * augmented[pivot_row][index];
            }
        }
        pivot_columns.push_back(col);
    }

    const int rank = int(pivot_columns.size());
    Matrix<T> result(size, size);
    if (rank + 1 < size) return result;

    if (rank == size) {
        const T determinant = negate ? T() - pivot_product : pivot_product;
        for (int row = 0; row < size; row++) {
            for (int col = 0; col < size; col++) {
                result[row][col] = determinant * augmented[row][size + col];
            }
        }
        return result;
    }

    int free_column = 0;
    while (free_column < rank && pivot_columns[free_column] == free_column) {
        free_column++;
    }
    std::vector<T> right_null(size);
    right_null[free_column] = T(1);
    for (int row = 0; row < rank; row++) {
        right_null[pivot_columns[row]] = T() - augmented[row][free_column];
    }

    T scale = pivot_product;
    if (negate != bool((size - 1 + free_column) & 1)) scale = T() - scale;
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            result[row][col] =
                scale * right_null[row] * augmented[size - 1][size + col];
        }
    }
    return result;
}

}  // namespace matrix
}  // namespace m1une

#endif  // M1UNE_MATRIX_ADJUGATE_HPP

#line 1 "math/matrix/adjugate.hpp"



#include <cassert>
#include <vector>

#line 1 "math/matrix/matrix.hpp"



#line 5 "math/matrix/matrix.hpp"
#include <cstddef>
#include <cstdint>
#include <utility>
#line 9 "math/matrix/matrix.hpp"

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 8 "math/matrix/adjugate.hpp"

namespace m1une {
namespace matrix {

template <class T>
Matrix<T> adjugate(Matrix<T> matrix) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    Matrix<T> augmented(size, size * 2);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            augmented[row][col] = matrix[row][col];
        }
        augmented[row][size + row] = T(1);
    }

    std::vector<int> pivot_columns;
    T pivot_product = T(1);
    bool negate = false;
    for (int col = 0; col < size && int(pivot_columns.size()) < size; col++) {
        const int pivot_row = int(pivot_columns.size());
        int pivot = pivot_row;
        while (pivot < size && augmented[pivot][col] == T()) pivot++;
        if (pivot == size) continue;
        if (pivot != pivot_row) {
            augmented.swap_rows(pivot, pivot_row);
            negate = !negate;
        }

        const T pivot_value = augmented[pivot_row][col];
        pivot_product *= pivot_value;
        const T inverse_pivot = T(1) / pivot_value;
        for (int index = col; index < size; index++) {
            augmented[pivot_row][index] *= inverse_pivot;
        }
        for (int index = size; index < size * 2; index++) {
            augmented[pivot_row][index] *= inverse_pivot;
        }

        for (int row = 0; row < size; row++) {
            if (row == pivot_row || augmented[row][col] == T()) continue;
            const T factor = augmented[row][col];
            augmented[row][col] = T();
            for (int index = col + 1; index < size; index++) {
                augmented[row][index] -= factor * augmented[pivot_row][index];
            }
            for (int index = size; index < size * 2; index++) {
                augmented[row][index] -= factor * augmented[pivot_row][index];
            }
        }
        pivot_columns.push_back(col);
    }

    const int rank = int(pivot_columns.size());
    Matrix<T> result(size, size);
    if (rank + 1 < size) return result;

    if (rank == size) {
        const T determinant = negate ? T() - pivot_product : pivot_product;
        for (int row = 0; row < size; row++) {
            for (int col = 0; col < size; col++) {
                result[row][col] = determinant * augmented[row][size + col];
            }
        }
        return result;
    }

    int free_column = 0;
    while (free_column < rank && pivot_columns[free_column] == free_column) {
        free_column++;
    }
    std::vector<T> right_null(size);
    right_null[free_column] = T(1);
    for (int row = 0; row < rank; row++) {
        right_null[pivot_columns[row]] = T() - augmented[row][free_column];
    }

    T scale = pivot_product;
    if (negate != bool((size - 1 + free_column) & 1)) scale = T() - scale;
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            result[row][col] =
                scale * right_null[row] * augmented[size - 1][size + col];
        }
    }
    return result;
}

}  // namespace matrix
}  // namespace m1une

m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

Adjugate Matrix
(math/matrix/adjugate.hpp)

Overview

Requirements

Public Interface

Example

Depends on

Required by

Verified with

Code

Adjugate Matrix (math/matrix/adjugate.hpp)

Overview

Requirements

Public Interface

Example

Depends on

Required by

Verified with

Code

Adjugate Matrix
(math/matrix/adjugate.hpp)