Sparse Determinant
(math/matrix/sparse_determinant.hpp)

• View this file on GitHub
• Last update: 2026-07-14 02:42:28+09:00
• Include: #include "math/matrix/sparse_determinant.hpp"

Overview

sparse_determinant computes the determinant of a sparse square matrix over a large finite field without materializing a dense matrix. It combines random diagonal preconditioning, a scalar Wiedemann sequence, and Berlekamp–Massey.

The algorithm is Las Vegas: the returned value is exact, while the number of preconditioning attempts is random and has constant expectation over a large field.

Requirements

T must be a field type supporting construction from integers, equality, addition, subtraction, multiplication, and division.

The convenience overload also requires T::mod() and assumes the nonzero residues modulo T::mod() are suitable random field elements. This includes the repository’s static and dynamic modular integer types when their modulus is prime.

Entries use zero-based row and column indices. Duplicate positions are allowed and their values are added. Explicit zero entries are allowed.

API

template <class T>
struct SparseMatrixEntry {
    int row;
    int col;
    T value;
};

template <class T, class RandomValue>
T sparse_determinant_with_randomizer(
    int size,
    const std::vector<SparseMatrixEntry<T>>& entries,
    RandomValue random_nonzero
);

template <class T>
T sparse_determinant(
    int size,
    const std::vector<SparseMatrixEntry<T>>& entries,
    std::uint64_t seed = 0x243f6a8885a308d3ULL
);

random_nonzero() must return independent nonzero values of T. The default overload uses a deterministic SplitMix64 stream; supplying the same seed makes the sequence of attempts reproducible.

Function	Description	Complexity
sparse_determinant_with_randomizer(size, entries, random_nonzero)	Computes the determinant using the supplied field sampler.	Expected $O(N(M+N))$ time, $O(N+M)$ memory
sparse_determinant(size, entries, seed)	Uses the built-in deterministic pseudorandom stream.	Expected $O(N(M+N))$ time, $O(N+M)$ memory

Here $M$ is entries.size(). The determinant of the empty matrix is T(1). The functions do not mutate entries.

Example

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

#include <iostream>
#include <vector>

int main() {
    using mint = m1une::math::modint998244353;
    using Entry = m1une::matrix::SparseMatrixEntry<mint>;
    std::vector<Entry> entries;
    entries.push_back(Entry{0, 0, mint(2)});
    entries.push_back(Entry{1, 1, mint(3)});
    std::cout << m1une::matrix::sparse_determinant(2, entries) << "\n";  // 6
}

Required by

• Math All (math/all.hpp)
• Matrix Bundle (math/matrix/all.hpp)

Verified with

• verify/math/math_algorithms.test.cpp
• verify/math/matrix/matrix.test.cpp
• verify/math/matrix/sparse_determinant.test.cpp

Code

#ifndef M1UNE_MATRIX_SPARSE_DETERMINANT_HPP
#define M1UNE_MATRIX_SPARSE_DETERMINANT_HPP 1

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

namespace m1une {
namespace matrix {

template <class T>
struct SparseMatrixEntry {
    int row;
    int col;
    T value;
};

namespace internal {

struct SparseDeterminantRandom {
    std::uint64_t state;

    explicit SparseDeterminantRandom(std::uint64_t seed) : state(seed) {}

    std::uint64_t operator()() {
        std::uint64_t value = (state += 0x9e3779b97f4a7c15ULL);
        value = (value ^ (value >> 30)) * 0xbf58476d1ce4e5b9ULL;
        value = (value ^ (value >> 27)) * 0x94d049bb133111ebULL;
        return value ^ (value >> 31);
    }
};

template <class T>
std::vector<T> berlekamp_massey(const std::vector<T>& sequence) {
    std::vector<T> recurrence(1, T(1));
    std::vector<T> previous(1, T(1));
    int degree = 0;
    int shift = 1;
    T previous_discrepancy = T(1);

    for (int index = 0; index < int(sequence.size()); index++) {
        T discrepancy = sequence[index];
        for (int i = 1; i <= degree; i++) {
            discrepancy += recurrence[i] * sequence[index - i];
        }
        if (discrepancy == T()) {
            shift++;
            continue;
        }

        const T factor = discrepancy / previous_discrepancy;
        std::vector<T> old_recurrence = recurrence;
        if (int(recurrence.size()) < int(previous.size()) + shift) {
            recurrence.resize(previous.size() + std::size_t(shift), T());
        }
        for (int i = 0; i < int(previous.size()); i++) {
            recurrence[i + shift] -= factor * previous[i];
        }

        if (2 * degree <= index) {
            degree = index + 1 - degree;
            previous = std::move(old_recurrence);
            previous_discrepancy = discrepancy;
            shift = 1;
        } else {
            shift++;
        }
    }
    recurrence.resize(std::size_t(degree + 1));
    return recurrence;
}

}  // namespace internal

// Randomized black-box determinant over a finite field. random_nonzero must
// return independent nonzero field elements.
template <class T, class RandomValue>
T sparse_determinant_with_randomizer(
    int size, const std::vector<SparseMatrixEntry<T>>& entries,
    RandomValue random_nonzero
) {
    assert(size >= 0);
    for (const SparseMatrixEntry<T>& entry : entries) {
        assert(0 <= entry.row && entry.row < size);
        assert(0 <= entry.col && entry.col < size);
    }
    if (size == 0) return T(1);

    auto random_vector = [&]() {
        std::vector<T> result(size);
        for (T& value : result) {
            value = random_nonzero();
            assert(value != T());
        }
        return result;
    };

    while (true) {
        std::vector<T> diagonal = random_vector();
        std::vector<T> left = random_vector();
        std::vector<T> state = random_vector();
        std::vector<T> sequence(std::size_t(2 * size));

        for (int step = 0; step < 2 * size; step++) {
            for (int i = 0; i < size; i++) sequence[step] += left[i] * state[i];
            for (int i = 0; i < size; i++) state[i] *= diagonal[i];

            std::vector<T> next(size);
            for (const SparseMatrixEntry<T>& entry : entries) {
                next[entry.row] += entry.value * state[entry.col];
            }
            state = std::move(next);
        }

        std::vector<T> recurrence = internal::berlekamp_massey(sequence);
        if (recurrence.back() == T()) return T();
        if (int(recurrence.size()) != size + 1) continue;

        T determinant = recurrence.back();
        if (size % 2 == 1) determinant = T() - determinant;
        for (const T& value : diagonal) determinant /= value;
        return determinant;
    }
}

template <class T>
T sparse_determinant(
    int size, const std::vector<SparseMatrixEntry<T>>& entries,
    std::uint64_t seed = 0x243f6a8885a308d3ULL
) {
    const std::uint64_t modulus = T::mod();
    assert(modulus > 1);
    internal::SparseDeterminantRandom random(seed);
    auto random_nonzero = [&]() {
        return T(1 + random() % (modulus - 1));
    };
    return sparse_determinant_with_randomizer<T>(size, entries, random_nonzero);
}

}  // namespace matrix
}  // namespace m1une

#endif  // M1UNE_MATRIX_SPARSE_DETERMINANT_HPP

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



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

namespace m1une {
namespace matrix {

template <class T>
struct SparseMatrixEntry {
    int row;
    int col;
    T value;
};

namespace internal {

struct SparseDeterminantRandom {
    std::uint64_t state;

    explicit SparseDeterminantRandom(std::uint64_t seed) : state(seed) {}

    std::uint64_t operator()() {
        std::uint64_t value = (state += 0x9e3779b97f4a7c15ULL);
        value = (value ^ (value >> 30)) * 0xbf58476d1ce4e5b9ULL;
        value = (value ^ (value >> 27)) * 0x94d049bb133111ebULL;
        return value ^ (value >> 31);
    }
};

template <class T>
std::vector<T> berlekamp_massey(const std::vector<T>& sequence) {
    std::vector<T> recurrence(1, T(1));
    std::vector<T> previous(1, T(1));
    int degree = 0;
    int shift = 1;
    T previous_discrepancy = T(1);

    for (int index = 0; index < int(sequence.size()); index++) {
        T discrepancy = sequence[index];
        for (int i = 1; i <= degree; i++) {
            discrepancy += recurrence[i] * sequence[index - i];
        }
        if (discrepancy == T()) {
            shift++;
            continue;
        }

        const T factor = discrepancy / previous_discrepancy;
        std::vector<T> old_recurrence = recurrence;
        if (int(recurrence.size()) < int(previous.size()) + shift) {
            recurrence.resize(previous.size() + std::size_t(shift), T());
        }
        for (int i = 0; i < int(previous.size()); i++) {
            recurrence[i + shift] -= factor * previous[i];
        }

        if (2 * degree <= index) {
            degree = index + 1 - degree;
            previous = std::move(old_recurrence);
            previous_discrepancy = discrepancy;
            shift = 1;
        } else {
            shift++;
        }
    }
    recurrence.resize(std::size_t(degree + 1));
    return recurrence;
}

}  // namespace internal

// Randomized black-box determinant over a finite field. random_nonzero must
// return independent nonzero field elements.
template <class T, class RandomValue>
T sparse_determinant_with_randomizer(
    int size, const std::vector<SparseMatrixEntry<T>>& entries,
    RandomValue random_nonzero
) {
    assert(size >= 0);
    for (const SparseMatrixEntry<T>& entry : entries) {
        assert(0 <= entry.row && entry.row < size);
        assert(0 <= entry.col && entry.col < size);
    }
    if (size == 0) return T(1);

    auto random_vector = [&]() {
        std::vector<T> result(size);
        for (T& value : result) {
            value = random_nonzero();
            assert(value != T());
        }
        return result;
    };

    while (true) {
        std::vector<T> diagonal = random_vector();
        std::vector<T> left = random_vector();
        std::vector<T> state = random_vector();
        std::vector<T> sequence(std::size_t(2 * size));

        for (int step = 0; step < 2 * size; step++) {
            for (int i = 0; i < size; i++) sequence[step] += left[i] * state[i];
            for (int i = 0; i < size; i++) state[i] *= diagonal[i];

            std::vector<T> next(size);
            for (const SparseMatrixEntry<T>& entry : entries) {
                next[entry.row] += entry.value * state[entry.col];
            }
            state = std::move(next);
        }

        std::vector<T> recurrence = internal::berlekamp_massey(sequence);
        if (recurrence.back() == T()) return T();
        if (int(recurrence.size()) != size + 1) continue;

        T determinant = recurrence.back();
        if (size % 2 == 1) determinant = T() - determinant;
        for (const T& value : diagonal) determinant /= value;
        return determinant;
    }
}

template <class T>
T sparse_determinant(
    int size, const std::vector<SparseMatrixEntry<T>>& entries,
    std::uint64_t seed = 0x243f6a8885a308d3ULL
) {
    const std::uint64_t modulus = T::mod();
    assert(modulus > 1);
    internal::SparseDeterminantRandom random(seed);
    auto random_nonzero = [&]() {
        return T(1 + random() % (modulus - 1));
    };
    return sparse_determinant_with_randomizer<T>(size, entries, random_nonzero);
}

}  // namespace matrix
}  // namespace m1une

Back to top page