Linear matroid parity (unweighted) （線形マトロイドパリティ）
(combinatorial_opt/linear_matroid_parity.hpp)

• View this file on GitHub
• Last update: 2023-05-21 18:11:51+09:00
• Include: #include "combinatorial_opt/linear_matroid_parity.hpp"

（重みなし）線形マトロイドパリティ問題とは，$M$ 個のベクトルの組 $(\mathbf{b}_i, \mathbf{c}_i)$ $(i = 1, \dots, M)$ が与えられたとき，$\{1, \dots, M \}$ の部分集合 $S$ であって $2 |S |$ 個のベクトル ${\{b_i\}}_{i \in S} + {\{ c_i \}}_{i \in S}$ が線形独立であるようなもののうち要素数最大のものを求める問題．

本コードは，線形マトロイドパリティ問題を $O(r^2 (r + m))$ で解く乱択アルゴリズムである．試行一回あたりの失敗確率は $O(r / p)$ （$p$ は有限体の位数）である．

使用方法

vector<pair<mint, mint>> vector_pairs;

vector<bool> sol_vec = linear_matroid_parity(vector_pairs);

int answer_size = linear_matroid_parity_size(vector_pairs); // Return size only, a bit faster

問題例

• No.1773 Love Triangle - yukicoder

文献・リンク集

• [1] H. Y. Cheung, L. C. Lau, K. M. Leung, “Algebraic Algorithms for Linear Matroid Parity Problems,” ACM Transactions on Algorithms, 10(3), 1-26, 2014.
• Matroid parity problem - Wikipedia

Depends on

• linear_algebra_matrix/matrix.hpp

Verified with

• combinatorial_opt/test/linear_matroid_parity.yuki1773.test.cpp
• combinatorial_opt/test/linear_matroid_parity_size.yuki1773.test.cpp

Code

#pragma once
#include "../linear_algebra_matrix/matrix.hpp"
#include <cassert>
#include <chrono>
#include <numeric>
#include <random>
#include <utility>
#include <vector>

// Solve linear matroid parity problem and return (especially lexicographically smallest) binary
// vector
// Complexity: O(d^2 (d + m)), d: dimension, m: number of input vector pairs
// Reference:
// [1] H. Y. Cheung, L. C. Lau, K. M. Leung,
//     "Algebraic Algorithms for Linear Matroid Parity Problems,"
//     ACM Transactions on Algorithms, 10(3), 1-26, 2014.
template <class ModInt>
std::vector<bool>
linear_matroid_parity(std::vector<std::pair<std::vector<ModInt>, std::vector<ModInt>>> bcs,
                      long long seed = std::chrono::steady_clock::now().time_since_epoch().count()) {
    if (bcs.empty()) return {};
    const int r = bcs[0].first.size(), r2 = (r + 1) / 2;
    int m = bcs.size();
    for (auto &v : bcs) v.first.resize(r2 * 2), v.second.resize(r2 * 2);

    std::mt19937 mt(seed);
    std::uniform_int_distribution<int> d(0, ModInt::mod() - 1);

    std::vector<ModInt> x(m), xadd(r2);

    matrix<ModInt> Yinv; // (r2 * 2) * (r2 * 2) matrix
    int rankY = -1;
    while (rankY < r2 * 2) {
        Yinv = matrix<ModInt>(r2 * 2, r2 * 2);

        auto add_i = [&](int i) {
            x[i] = d(mt);
            const auto &b = bcs[i].first, &c = bcs[i].second;
            for (int j = 0; j < r2 * 2; j++) {
                for (int k = 0; k < r2 * 2; k++) Yinv[j][k] += x[i] * (b[j] * c[k] - c[j] * b[k]);
            }
        };
        for (int i = 0; i < m; ++i) add_i(i);

        int num_add_vec = (r2 * 2 - Yinv.rank()) / 2;

        bcs.resize(m + num_add_vec,
                   std::make_pair(std::vector<ModInt>(r2 * 2), std::vector<ModInt>(r2 * 2)));
        x.resize(bcs.size());
        for (int i = m; i < int(bcs.size()); ++i) {
            for (auto &x : bcs[i].first) x = d(mt);
            for (auto &x : bcs[i].second) x = d(mt);
        }

        for (int i = m; i < int(bcs.size()); i++) add_i(i);
        rankY = Yinv.inverse();
    }

    std::vector<bool> ret(bcs.size(), 1);

    auto try_erase_i = [&](int i) {
        auto b = bcs[i].first, c = bcs[i].second;
        b.resize(r2 * 2, 0), c.resize(r2 * 2, 0);
        std::vector<ModInt> Yib = Yinv * b, Yic = Yinv * c;
        ModInt bYic = std::inner_product(b.begin(), b.end(), Yic.begin(), ModInt());
        ModInt a00 = bYic * x[i] + 1;
        if (a00 == ModInt()) return;
        ret[i] = 0;
        const ModInt f = x[i] / a00;
        for (int j = 0; j < r2 * 2; ++j) {
            for (int k = 0; k < r2 * 2; ++k) {
                Yinv[j][k] -= f * (Yib[j] * Yic[k] - Yic[j] * Yib[k]);
            }
        }
    };

    for (int i = m; i < int(bcs.size()); ++i) try_erase_i(i);
    for (int i = 0; i < m; i++) try_erase_i(i);
    ret.resize(m);
    return ret;
}

// Solve linear matroid parity problem, size-only (no construction)
template <class ModInt>
int linear_matroid_parity_size(
    const std::vector<std::pair<std::vector<ModInt>, std::vector<ModInt>>> &bcs,
    long long seed = std::chrono::steady_clock::now().time_since_epoch().count()) {

    if (bcs.empty()) return 0;

    std::mt19937 mt(seed);
    std::uniform_int_distribution<int> d(0, ModInt::mod() - 1);

    const int r = bcs[0].first.size();
    matrix<ModInt> mat(r, r);
    for (const auto &bc : bcs) {
        const auto &b = bc.first, &c = bc.second;
        ModInt x = d(mt);
        for (int i = 0; i < r; ++i) {
            for (int j = 0; j < r; ++j) mat[i][j] += x * (b[i] * c[j] - b[j] * c[i]);
        }
    }
    return mat.rank() / 2;
}

#line 2 "linear_algebra_matrix/matrix.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <iterator>
#include <type_traits>
#include <utility>
#include <vector>

namespace matrix_ {
struct has_id_method_impl {
    template <class T_> static auto check(T_ *) -> decltype(T_::id(), std::true_type());
    template <class T_> static auto check(...) -> std::false_type;
};
template <class T_> struct has_id : decltype(has_id_method_impl::check<T_>(nullptr)) {};
} // namespace matrix_

template <typename T> struct matrix {
    int H, W;
    std::vector<T> elem;
    typename std::vector<T>::iterator operator[](int i) { return elem.begin() + i * W; }
    inline T &at(int i, int j) { return elem[i * W + j]; }
    inline T get(int i, int j) const { return elem[i * W + j]; }
    int height() const { return H; }
    int width() const { return W; }
    std::vector<std::vector<T>> vecvec() const {
        std::vector<std::vector<T>> ret(H);
        for (int i = 0; i < H; i++) {
            std::copy(elem.begin() + i * W, elem.begin() + (i + 1) * W, std::back_inserter(ret[i]));
        }
        return ret;
    }
    operator std::vector<std::vector<T>>() const { return vecvec(); }
    matrix() = default;
    matrix(int H, int W) : H(H), W(W), elem(H * W) {}
    matrix(const std::vector<std::vector<T>> &d) : H(d.size()), W(d.size() ? d[0].size() : 0) {
        for (auto &raw : d) std::copy(raw.begin(), raw.end(), std::back_inserter(elem));
    }

    template <typename T2, typename std::enable_if<matrix_::has_id<T2>::value>::type * = nullptr>
    static T2 _T_id() {
        return T2::id();
    }
    template <typename T2, typename std::enable_if<!matrix_::has_id<T2>::value>::type * = nullptr>
    static T2 _T_id() {
        return T2(1);
    }

    static matrix Identity(int N) {
        matrix ret(N, N);
        for (int i = 0; i < N; i++) ret.at(i, i) = _T_id<T>();
        return ret;
    }

    matrix operator-() const {
        matrix ret(H, W);
        for (int i = 0; i < H * W; i++) ret.elem[i] = -elem[i];
        return ret;
    }
    matrix operator*(const T &v) const {
        matrix ret = *this;
        for (auto &x : ret.elem) x *= v;
        return ret;
    }
    matrix operator/(const T &v) const {
        matrix ret = *this;
        const T vinv = _T_id<T>() / v;
        for (auto &x : ret.elem) x *= vinv;
        return ret;
    }
    matrix operator+(const matrix &r) const {
        matrix ret = *this;
        for (int i = 0; i < H * W; i++) ret.elem[i] += r.elem[i];
        return ret;
    }
    matrix operator-(const matrix &r) const {
        matrix ret = *this;
        for (int i = 0; i < H * W; i++) ret.elem[i] -= r.elem[i];
        return ret;
    }
    matrix operator*(const matrix &r) const {
        matrix ret(H, r.W);
        for (int i = 0; i < H; i++) {
            for (int k = 0; k < W; k++) {
                for (int j = 0; j < r.W; j++) ret.at(i, j) += this->get(i, k) * r.get(k, j);
            }
        }
        return ret;
    }
    matrix &operator*=(const T &v) { return *this = *this * v; }
    matrix &operator/=(const T &v) { return *this = *this / v; }
    matrix &operator+=(const matrix &r) { return *this = *this + r; }
    matrix &operator-=(const matrix &r) { return *this = *this - r; }
    matrix &operator*=(const matrix &r) { return *this = *this * r; }
    bool operator==(const matrix &r) const { return H == r.H and W == r.W and elem == r.elem; }
    bool operator!=(const matrix &r) const { return H != r.H or W != r.W or elem != r.elem; }
    bool operator<(const matrix &r) const { return elem < r.elem; }
    matrix pow(int64_t n) const {
        matrix ret = Identity(H);
        bool ret_is_id = true;
        if (n == 0) return ret;
        for (int i = 63 - __builtin_clzll(n); i >= 0; i--) {
            if (!ret_is_id) ret *= ret;
            if ((n >> i) & 1) ret *= (*this), ret_is_id = false;
        }
        return ret;
    }
    std::vector<T> pow_vec(int64_t n, std::vector<T> vec) const {
        matrix x = *this;
        while (n) {
            if (n & 1) vec = x * vec;
            x *= x;
            n >>= 1;
        }
        return vec;
    };
    matrix transpose() const {
        matrix ret(W, H);
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < W; j++) ret.at(j, i) = this->get(i, j);
        }
        return ret;
    }
    // Gauss-Jordan elimination
    // - Require inverse for every non-zero element
    // - Complexity: O(H^2 W)
    template <typename T2, typename std::enable_if<std::is_floating_point<T2>::value>::type * = nullptr>
    static int choose_pivot(const matrix<T2> &mtr, int h, int c) noexcept {
        int piv = -1;
        for (int j = h; j < mtr.H; j++) {
            if (mtr.get(j, c) and (piv < 0 or std::abs(mtr.get(j, c)) > std::abs(mtr.get(piv, c))))
                piv = j;
        }
        return piv;
    }
    template <typename T2, typename std::enable_if<!std::is_floating_point<T2>::value>::type * = nullptr>
    static int choose_pivot(const matrix<T2> &mtr, int h, int c) noexcept {
        for (int j = h; j < mtr.H; j++) {
            if (mtr.get(j, c) != T2()) return j;
        }
        return -1;
    }
    matrix gauss_jordan() const {
        int c = 0;
        matrix mtr(*this);
        std::vector<int> ws;
        ws.reserve(W);
        for (int h = 0; h < H; h++) {
            if (c == W) break;
            int piv = choose_pivot(mtr, h, c);
            if (piv == -1) {
                c++;
                h--;
                continue;
            }
            if (h != piv) {
                for (int w = 0; w < W; w++) {
                    std::swap(mtr[piv][w], mtr[h][w]);
                    mtr.at(piv, w) *= -_T_id<T>(); // To preserve sign of determinant
                }
            }
            ws.clear();
            for (int w = c; w < W; w++) {
                if (mtr.at(h, w) != T()) ws.emplace_back(w);
            }
            const T hcinv = _T_id<T>() / mtr.at(h, c);
            for (int hh = 0; hh < H; hh++)
                if (hh != h) {
                    const T coeff = mtr.at(hh, c) * hcinv;
                    for (auto w : ws) mtr.at(hh, w) -= mtr.at(h, w) * coeff;
                    mtr.at(hh, c) = T();
                }
            c++;
        }
        return mtr;
    }
    int rank_of_gauss_jordan() const {
        for (int i = H * W - 1; i >= 0; i--) {
            if (elem[i] != 0) return i / W + 1;
        }
        return 0;
    }
    int rank() const { return gauss_jordan().rank_of_gauss_jordan(); }

    T determinant_of_upper_triangle() const {
        T ret = _T_id<T>();
        for (int i = 0; i < H; i++) ret *= get(i, i);
        return ret;
    }
    int inverse() {
        assert(H == W);
        std::vector<std::vector<T>> ret = Identity(H), tmp = *this;
        int rank = 0;
        for (int i = 0; i < H; i++) {
            int ti = i;
            while (ti < H and tmp[ti][i] == T()) ti++;
            if (ti == H) {
                continue;
            } else {
                rank++;
            }
            ret[i].swap(ret[ti]), tmp[i].swap(tmp[ti]);
            T inv = _T_id<T>() / tmp[i][i];
            for (int j = 0; j < W; j++) ret[i][j] *= inv;
            for (int j = i + 1; j < W; j++) tmp[i][j] *= inv;
            for (int h = 0; h < H; h++) {
                if (i == h) continue;
                const T c = -tmp[h][i];
                for (int j = 0; j < W; j++) ret[h][j] += ret[i][j] * c;
                for (int j = i + 1; j < W; j++) tmp[h][j] += tmp[i][j] * c;
            }
        }
        *this = ret;
        return rank;
    }
    friend std::vector<T> operator*(const matrix &m, const std::vector<T> &v) {
        assert(m.W == int(v.size()));
        std::vector<T> ret(m.H);
        for (int i = 0; i < m.H; i++) {
            for (int j = 0; j < m.W; j++) ret[i] += m.get(i, j) * v[j];
        }
        return ret;
    }
    friend std::vector<T> operator*(const std::vector<T> &v, const matrix &m) {
        assert(int(v.size()) == m.H);
        std::vector<T> ret(m.W);
        for (int i = 0; i < m.H; i++) {
            for (int j = 0; j < m.W; j++) ret[j] += v[i] * m.get(i, j);
        }
        return ret;
    }
    std::vector<T> prod(const std::vector<T> &v) const { return (*this) * v; }
    std::vector<T> prod_left(const std::vector<T> &v) const { return v * (*this); }
    template <class OStream> friend OStream &operator<<(OStream &os, const matrix &x) {
        os << "[(" << x.H << " * " << x.W << " matrix)";
        os << "\n[column sums: ";
        for (int j = 0; j < x.W; j++) {
            T s = T();
            for (int i = 0; i < x.H; i++) s += x.get(i, j);
            os << s << ",";
        }
        os << "]";
        for (int i = 0; i < x.H; i++) {
            os << "\n[";
            for (int j = 0; j < x.W; j++) os << x.get(i, j) << ",";
            os << "]";
        }
        os << "]\n";
        return os;
    }
    template <class IStream> friend IStream &operator>>(IStream &is, matrix &x) {
        for (auto &v : x.elem) is >> v;
        return is;
    }
};
#line 4 "combinatorial_opt/linear_matroid_parity.hpp"
#include <chrono>
#include <numeric>
#include <random>
#line 9 "combinatorial_opt/linear_matroid_parity.hpp"

// Solve linear matroid parity problem and return (especially lexicographically smallest) binary
// vector
// Complexity: O(d^2 (d + m)), d: dimension, m: number of input vector pairs
// Reference:
// [1] H. Y. Cheung, L. C. Lau, K. M. Leung,
//     "Algebraic Algorithms for Linear Matroid Parity Problems,"
//     ACM Transactions on Algorithms, 10(3), 1-26, 2014.
template <class ModInt>
std::vector<bool>
linear_matroid_parity(std::vector<std::pair<std::vector<ModInt>, std::vector<ModInt>>> bcs,
                      long long seed = std::chrono::steady_clock::now().time_since_epoch().count()) {
    if (bcs.empty()) return {};
    const int r = bcs[0].first.size(), r2 = (r + 1) / 2;
    int m = bcs.size();
    for (auto &v : bcs) v.first.resize(r2 * 2), v.second.resize(r2 * 2);

    std::mt19937 mt(seed);
    std::uniform_int_distribution<int> d(0, ModInt::mod() - 1);

    std::vector<ModInt> x(m), xadd(r2);

    matrix<ModInt> Yinv; // (r2 * 2) * (r2 * 2) matrix
    int rankY = -1;
    while (rankY < r2 * 2) {
        Yinv = matrix<ModInt>(r2 * 2, r2 * 2);

        auto add_i = [&](int i) {
            x[i] = d(mt);
            const auto &b = bcs[i].first, &c = bcs[i].second;
            for (int j = 0; j < r2 * 2; j++) {
                for (int k = 0; k < r2 * 2; k++) Yinv[j][k] += x[i] * (b[j] * c[k] - c[j] * b[k]);
            }
        };
        for (int i = 0; i < m; ++i) add_i(i);

        int num_add_vec = (r2 * 2 - Yinv.rank()) / 2;

        bcs.resize(m + num_add_vec,
                   std::make_pair(std::vector<ModInt>(r2 * 2), std::vector<ModInt>(r2 * 2)));
        x.resize(bcs.size());
        for (int i = m; i < int(bcs.size()); ++i) {
            for (auto &x : bcs[i].first) x = d(mt);
            for (auto &x : bcs[i].second) x = d(mt);
        }

        for (int i = m; i < int(bcs.size()); i++) add_i(i);
        rankY = Yinv.inverse();
    }

    std::vector<bool> ret(bcs.size(), 1);

    auto try_erase_i = [&](int i) {
        auto b = bcs[i].first, c = bcs[i].second;
        b.resize(r2 * 2, 0), c.resize(r2 * 2, 0);
        std::vector<ModInt> Yib = Yinv * b, Yic = Yinv * c;
        ModInt bYic = std::inner_product(b.begin(), b.end(), Yic.begin(), ModInt());
        ModInt a00 = bYic * x[i] + 1;
        if (a00 == ModInt()) return;
        ret[i] = 0;
        const ModInt f = x[i] / a00;
        for (int j = 0; j < r2 * 2; ++j) {
            for (int k = 0; k < r2 * 2; ++k) {
                Yinv[j][k] -= f * (Yib[j] * Yic[k] - Yic[j] * Yib[k]);
            }
        }
    };

    for (int i = m; i < int(bcs.size()); ++i) try_erase_i(i);
    for (int i = 0; i < m; i++) try_erase_i(i);
    ret.resize(m);
    return ret;
}

// Solve linear matroid parity problem, size-only (no construction)
template <class ModInt>
int linear_matroid_parity_size(
    const std::vector<std::pair<std::vector<ModInt>, std::vector<ModInt>>> &bcs,
    long long seed = std::chrono::steady_clock::now().time_since_epoch().count()) {

    if (bcs.empty()) return 0;

    std::mt19937 mt(seed);
    std::uniform_int_distribution<int> d(0, ModInt::mod() - 1);

    const int r = bcs[0].first.size();
    matrix<ModInt> mat(r, r);
    for (const auto &bc : bcs) {
        const auto &b = bc.first, &c = bc.second;
        ModInt x = d(mt);
        for (int i = 0; i < r; ++i) {
            for (int j = 0; j < r; ++j) mat[i][j] += x * (b[i] * c[j] - b[j] * c[i]);
        }
    }
    return mat.rank() / 2;
}

Back to top page