cp-library-cpp

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

View the Project on GitHub suisen-cp/cp-library-cpp

:warning: Gaussian Elimination $\mathbb{F} _ 2$
(library/linear_algebra/gaussian_elimination_f2.hpp)

Gaussian Elimination $\mathbb{F} _ 2$

Depends on

Code

#ifndef SUISEN_GAUSSIAN_ELIMINATION_F2
#define SUISEN_GAUSSIAN_ELIMINATION_F2

#include "library/linear_algebra/matrix_f2.hpp"

namespace suisen {
    struct SolutionSpace {
        SolutionSpace() : _solution(std::nullopt), _basis{} {}
        SolutionSpace(const std::optional<DynamicBitSet> &solution, const std::vector<DynamicBitSet> &basis = {}) : _solution(solution), _basis(basis) {}
        int dim() const { return _basis.size(); }
        int rank() const { return _solution->size() - _basis.size(); }
        bool has_solution() const { return bool(_solution); }
        const DynamicBitSet& solution() const { return *_solution; }
        const std::vector<DynamicBitSet>& basis() const { return _basis; }

        SolutionSpace normalized() {
            SolutionSpace res{ _solution };
            res += *this;
            res._normalized = true;
            return res;
        }
        void normalize() {
            *this = normalized();
        }
        bool add_base(DynamicBitSet v) {
            assert(_normalized);
            for (const auto& e : _basis) if (DynamicBitSet w = e ^ v; w < v) v = std::move(w);
            return v ? (_basis.push_back(v), true) : false;
        }
        friend SolutionSpace& operator+=(SolutionSpace& s, const SolutionSpace& t) {
            for (const DynamicBitSet& base : t._basis) s.add_base(base);
            return s;
        }
        friend SolutionSpace operator+(SolutionSpace s, const SolutionSpace& t) {
            s += t;
            return s;
        }
        friend SolutionSpace& operator&=(SolutionSpace& s, const SolutionSpace& t) {
            return s = s & t;
        }
        friend SolutionSpace operator&(const SolutionSpace& s, const SolutionSpace& t) {
            if (not s._normalized and t._normalized) return t & s;
            assert(s._normalized);
            std::vector<std::pair<DynamicBitSet, DynamicBitSet>> basis;
            for (const DynamicBitSet& e : s._basis) basis.emplace_back(e, e);
            SolutionSpace res;
            for (DynamicBitSet e : t._basis) {
                DynamicBitSet s(e.size());
                for (const auto& [v, t] : basis) if (DynamicBitSet w = e ^ v; w < e) e = std::move(w), s ^= t;
                if (e) basis.emplace_back(e, s);
                else res.add_base(s);
            }
            return res;
        }
    private:
        bool _normalized = false;
        std::optional<DynamicBitSet> _solution;
        std::vector<DynamicBitSet> _basis;
    };

    SolutionSpace gaussian_elimination_f2(MatrixF2 A, const DynamicBitSet &b) {
        assert(A.row_size() == b.size());
        const std::size_t n = A.row_size(), m = A.col_size();
        for (std::size_t i = 0; i < n; ++i) A[i].push_back(b[i]);

        std::vector<std::size_t> ones, zeros;
        std::size_t rank = 0;
        for (std::size_t j = 0; rank < n and j < m; ++j) {
            for (std::size_t i = rank + 1; i < n; ++i) if (A[i][j]) {
                std::swap(A[rank], A[i]);
                if (A[i][j]) A[i] ^= A[rank];
            }
            if (A[rank][j]) {
                for (std::size_t i = 0; i < rank; ++i) if (A[i][j]) A[i] ^= A[rank];
                ++rank;
                ones.push_back(j);
            } else {
                zeros.push_back(j);
            }
        }
        for (std::size_t i = rank; i < n; ++i) {
            if (A[i][m]) return SolutionSpace{};
        }
        DynamicBitSet solution(m);
        for (std::size_t i = 0; i < rank; ++i) {
            solution[ones[i]] = A[i][m];
        }
        std::vector<DynamicBitSet> basis;
        for (std::size_t j : zeros) {
            DynamicBitSet base(m);
            for (std::size_t i = 0; i < rank; ++i) base[ones[i]] = A[i][j];
            base[j] = 1;
            basis.push_back(std::move(base));
        }
        return SolutionSpace{ solution, basis };
    }
} // namespace suisen


#endif // SUISEN_GAUSSIAN_ELIMINATION_F2
#line 1 "library/linear_algebra/gaussian_elimination_f2.hpp"



#line 1 "library/linear_algebra/matrix_f2.hpp"



#include <cassert>
#include <optional>
#include <vector>

#line 1 "library/datastructure/util/dynamic_bitset.hpp"



#line 5 "library/datastructure/util/dynamic_bitset.hpp"
#include <limits>
#include <utility>
#line 8 "library/datastructure/util/dynamic_bitset.hpp"

namespace suisen {
    struct DynamicBitSet {
    private:
        using block = unsigned long long;
        static constexpr std::size_t block_size = std::numeric_limits<block>::digits;
        static constexpr std::size_t log_block_size = __builtin_ctz(block_size);

        struct bitref {
            block& b;
            std::size_t i;
            operator bool() const { return (b >> i) & 1; }
            bool test() const { return (b >> i) & 1; }
            void set() { b |= block(1) << i; }
            void reset() { b &= ~(block(1) << i); }
            void flip() { b ^= block(1) << i; }
            bitref& operator&=(bool val) { b &= block(val) << i; return *this; }
            bitref& operator|=(bool val) { b |= block(val) << i; return *this; }
            bitref& operator^=(bool val) { b ^= block(val) << i; return *this; }
            bitref& operator =(bool val) { val ? set() : reset(); return *this; }
            bitref& operator =(const bitref& v) { return (*this) = bool(v); }
        };

        std::size_t n;
        std::vector<block> blocks;
    public:
        DynamicBitSet(std::size_t n = 0, bool fill_value = false) : n(n), blocks((n + block_size - 1) >> log_block_size, fill_value ? ~block(0) : 0) {}

        bool empty() const { return n == 0; }
        int size() const { return n; }

        void resize(std::size_t new_size, bool fill_value = false) {
            std::size_t new_block_num = (new_size + block_size - 1) >> log_block_size;
            if (new_block_num < block_num()) {
                n = new_size;
                return blocks.resize(new_block_num);
            }
            blocks.resize(new_block_num);
            std::size_t old_size = std::exchange(n, new_size);
            if (old_size <= new_size) range_update(old_size, new_size, fill_value);
        }
        void push_back(bool val) {
            if (n & (block_size - 1)) {
                (*this)[n] = val;
            } else {
                blocks.push_back(val);
            }
            ++n;
        }
        void pop_back() {
            if ((n & (block_size - 1)) == 1) blocks.pop_back();
            --n;
        }

        friend bool operator==(const DynamicBitSet& x, const DynamicBitSet& y) {
            if (x.n != y.n) return false;
            if (x.empty()) return true;
            for (std::size_t i = 0; i < x.block_num() - 1; ++i) {
                if (x.blocks[i] != y.blocks[i]) return false;
            }
            const std::size_t num = x.n - ((x.block_num() - 1) << log_block_size);
            return get_lower_bits(x.blocks.back(), num) == get_lower_bits(y.blocks.back(), num);
        }
        friend bool operator!=(const DynamicBitSet& x, const DynamicBitSet& y) {
            return not (x == y);
        }
        friend bool operator<(const DynamicBitSet& x, const DynamicBitSet& y) {
            assert(x.n == y.n);
            if (x.empty()) return false;
            std::size_t num = x.n - ((x.block_num() - 1) << log_block_size);
            block tx = get_lower_bits(x.blocks.back(), num);
            block ty = get_lower_bits(y.blocks.back(), num);
            if (tx != ty) return tx < ty;
            for (std::size_t i = x.block_num() - 1; i-- > 0;) {
                if (x.blocks[i] != y.blocks[i]) return x.blocks[i] < y.blocks[i];
            }
            return false;
        }
        friend bool operator<=(const DynamicBitSet& x, const DynamicBitSet& y) {
            assert(x.n == y.n);
            if (x.empty()) return true;
            std::size_t num = x.n - ((x.block_num() - 1) << log_block_size);
            block tx = get_lower_bits(x.blocks.back(), num);
            block ty = get_lower_bits(y.blocks.back(), num);
            if (tx != ty) return tx < ty;
            for (std::size_t i = x.block_num() - 1; i-- > 0;) {
                if (x.blocks[i] != y.blocks[i]) return x.blocks[i] < y.blocks[i];
            }
            return true;
        }
        friend bool operator>(const DynamicBitSet& x, const DynamicBitSet& y) {
            return not (x <= y);
        }
        friend bool operator>=(const DynamicBitSet& x, const DynamicBitSet& y) {
            return not (x < y);
        }

        operator bool() const { return any(); }

        friend DynamicBitSet& operator&=(DynamicBitSet& x, const DynamicBitSet& y) {
            assert(x.n == y.n);
            for (std::size_t i = 0; i < y.block_num(); ++i) x.blocks[i] &= y.blocks[i];
            return x;
        }
        friend DynamicBitSet& operator|=(DynamicBitSet& x, const DynamicBitSet& y) {
            assert(x.n == y.n);
            for (std::size_t i = 0; i < y.block_num(); ++i) x.blocks[i] |= y.blocks[i];
            return x;
        }
        friend DynamicBitSet& operator^=(DynamicBitSet& x, const DynamicBitSet& y) {
            assert(x.n == y.n);
            for (std::size_t i = 0; i < y.block_num(); ++i) x.blocks[i] ^= y.blocks[i];
            return x;
        }
        friend DynamicBitSet operator&(DynamicBitSet x, const DynamicBitSet& y) { x &= y; return x; }
        friend DynamicBitSet operator|(DynamicBitSet x, const DynamicBitSet& y) { x |= y; return x; }
        friend DynamicBitSet operator^(DynamicBitSet x, const DynamicBitSet& y) { x ^= y; return x; }

        friend DynamicBitSet& operator<<=(DynamicBitSet &x, std::size_t shamt) {
            return x = x << shamt;
        }
        friend DynamicBitSet& operator>>=(DynamicBitSet &x, std::size_t shamt) {
            return x = x >> shamt;
        }
        friend DynamicBitSet operator<<(const DynamicBitSet &x, std::size_t shamt) {
            if (shamt >= x.n) return DynamicBitSet(x.size());
            DynamicBitSet res(x.size());
            std::size_t block_shamt = shamt >> log_block_size;
            std::size_t bit_shamt = shamt & (block_size - 1);
            for (std::size_t i = 0; i + block_shamt < res.block_num(); ++i) {
                if (bit_shamt == 0) {
                    res.blocks[i + block_shamt] = x.blocks[i];
                } else {
                    res.blocks[i + block_shamt] |= x.blocks[i] << bit_shamt;
                    if (i + block_shamt + 1 != res.block_num()) {
                        res.blocks[i + block_shamt + 1] |= x.blocks[i] >> (block_size - bit_shamt);
                    }
                }
            }
            return res;
        }
        friend DynamicBitSet operator>>(const DynamicBitSet& x, std::size_t shamt) {
            if (shamt >= x.n) return DynamicBitSet(x.size());
            DynamicBitSet res(x.size());
            std::size_t block_shamt = shamt >> log_block_size;
            std::size_t bit_shamt = shamt & (block_size - 1);
            for (std::size_t i = 0; i + block_shamt < x.block_num(); ++i) {
                if (bit_shamt == 0) {
                    res.blocks[i] = x.blocks[i + block_shamt];
                } else {
                    res.blocks[i] |= x.blocks[i + block_shamt] >> bit_shamt;
                    if (i + block_shamt + 1 != x.block_num()) {
                        res.blocks[i] |= x.blocks[i + block_shamt + 1] << (block_size - bit_shamt);
                    }
                }
            }
            res.range_reset(x.n - shamt, x.n);
            return res;
        }

        DynamicBitSet operator~() const {
            DynamicBitSet neg(n);
            for (std::size_t i = 0; i < block_num(); ++i) neg.blocks[i] = ~blocks[i];
            return neg;
        }

        bool operator[](std::size_t i) const {
            return (blocks[block_index(i)] >> bit_index(i)) & 1;
        }
        bitref operator[](std::size_t i) {
            return { blocks[block_index(i)], bit_index(i) };
        }

        void range_set(std::size_t l, std::size_t r) {
            assert(l <= r and r <= n);
            if (l == r) return;
            std::size_t lb = block_index(l), rb = block_index(r - 1);
            std::size_t li = bit_index(l), ri = bit_index(r);
            if (ri == 0) ri = block_size;
            if (lb == rb) {
                blocks[lb] |= mask_range_bits(~block(0), li, ri);
                return;
            }
            blocks[lb] |= mask_upper_bits(~block(0), block_size - li);
            blocks[rb] |= mask_lower_bits(~block(0), ri);
            for (std::size_t i = lb + 1; i < rb; ++i) blocks[i] = ~block(0);
        }
        void range_reset(std::size_t l, std::size_t r) {
            assert(l <= r and r <= n);
            if (l == r) return;
            std::size_t lb = block_index(l), rb = block_index(r - 1);
            std::size_t li = bit_index(l), ri = bit_index(r);
            if (ri == 0) ri = block_size;
            if (lb == rb) {
                blocks[lb] &= ~mask_range_bits(~block(0), li, ri);
                return;
            }
            blocks[lb] &= ~mask_upper_bits(~block(0), block_size - li);
            blocks[rb] &= ~mask_lower_bits(~block(0), ri);
            for (std::size_t i = lb + 1; i < rb; ++i) blocks[i] = block(0);
        }
        void range_flip(std::size_t l, std::size_t r) {
            assert(l <= r and r <= n);
            if (l == r) return;
            std::size_t lb = block_index(l), rb = block_index(r - 1);
            std::size_t li = bit_index(l), ri = bit_index(r);
            if (ri == 0) ri = block_size;
            if (lb == rb) {
                blocks[lb] ^= mask_range_bits(~block(0), li, ri);
                return;
            }
            blocks[lb] ^= mask_upper_bits(~block(0), block_size - li);
            blocks[rb] ^= mask_lower_bits(~block(0), ri);
            for (std::size_t i = lb + 1; i < rb; ++i) blocks[i] ^= ~block(0);
        }
        void range_update(std::size_t l, std::size_t r, bool val) {
            val ? range_set(l, r) : range_reset(l, r);
        }
        int range_count(std::size_t l, std::size_t r) const {
            assert(l <= r and r <= n);
            if (l == r) return 0;
            std::size_t lb = block_index(l), rb = block_index(r - 1);
            std::size_t li = bit_index(l), ri = bit_index(r);
            if (ri == 0) ri = block_size;
            
            if (lb == rb) {
                return __builtin_popcountll(blocks[lb] & mask_range_bits(~block(0), li, ri));
            }
            int res = 0;
            res += __builtin_popcountll(blocks[lb] & mask_upper_bits(~block(0), block_size - li));
            res += __builtin_popcountll(blocks[rb] & mask_lower_bits(~block(0), ri));
            for (std::size_t i = lb + 1; i < rb; ++i) res += __builtin_popcountll(blocks[i]);
            return res;
        }

        void set() {
            for (block& b : blocks) b = ~block(0);
        }
        void reset() {
            for (block& b : blocks) b = 0;
        }
        bool all() const {
            if (empty()) return true;
            for (std::size_t i = 0; i < block_num() - 1; ++i) {
                if (blocks[i] != ~block(0)) return false;
            }
            const std::size_t num = n - ((block_num() - 1) << log_block_size);
            assert(num);
            const block upper = ((block(1) << (block_size - num)) - 1) << num;
            return (upper | blocks.back()) == ~block(0);
        }
        bool none() const {
            if (empty()) return true;
            for (std::size_t i = 0; i < block_num() - 1; ++i) {
                if (blocks[i] != 0) return false;
            }
            const std::size_t num = n - ((block_num() - 1) << log_block_size);
            return get_lower_bits(blocks.back(), num) == 0;
        }
        bool any() const {
            return not none();
        }
        int count() const {
            if (empty()) return 0;
            int res = 0;
            for (std::size_t i = 0; i < block_num() - 1; ++i) {
                res += __builtin_popcountll(blocks[i]);
            }
            const std::size_t num = n - ((block_num() - 1) << log_block_size);
            return res + __builtin_popcountll(get_lower_bits(blocks.back(), num));
        }

        // Returns the position of first set bit. If there is no such positions, then returns size().
        int find_first() const {
            if (empty()) return size();
            for (std::size_t i = 0; i < block_num(); ++i) {
                if (blocks[i] != 0) return std::min(n, __builtin_ctzll(blocks[i]) | (i << log_block_size));
            }
            return n;
        }
        // Returns the position of first set bit after the given position (exclusive). If there is no such positions, then returns size().
        int find_next(std::size_t pos) const {
            std::size_t i = block_index(++pos);
            if (i >= blocks.size()) return n;
            block upper = mask_upper_bits(blocks[i], block_size - bit_index(pos));
            if (upper != 0) return std::min(n, __builtin_ctzll(upper) | (i << log_block_size));
            while (++i < block_num()) {
                if (blocks[i] != 0) return std::min(n, __builtin_ctzll(blocks[i]) | (i << log_block_size));
            }
            return n;
        }

        bool has_intersection(const DynamicBitSet& y) const {
            if (n > y.n) return y.has_intersection(*this);
            if (empty()) return false;
            for (std::size_t i = 0; i < block_num() - 1; ++i) {
                if (blocks[i] & y.blocks[i]) return true;
            }
            const std::size_t num = n - ((block_num() - 1) << log_block_size);
            return get_lower_bits(blocks.back(), num) & y.blocks[block_num() - 1];
        }
        bool is_disjoint(const DynamicBitSet& y) const {
            return not has_intersection(y);
        }

    private:
        static constexpr std::size_t block_index(std::size_t i) {
            return i >> log_block_size;
        }
        static constexpr std::size_t bit_index(std::size_t i) {
            return i & (block_size - 1);
        }
        static constexpr block get_lower_bits(block b, std::size_t num) {
            return num ? (b << (block_size - num) >> (block_size - num)) : block(0);
        }
        static constexpr block get_upper_bits(block b, std::size_t num) {
            return num ? (b >> (block_size - num)) : block(0);
        }
        static constexpr block get_range_bits(block b, std::size_t l, std::size_t r) {
            return l < r ? b << (block_size - r) >> (block_size - r + l) : block(0);
        }
        static constexpr block mask_lower_bits(block b, std::size_t num) {
            return get_lower_bits(b, num);
        }
        static constexpr block mask_upper_bits(block b, std::size_t num) {
            return num ? (b >> (block_size - num) << (block_size - num)) : block(0);
        }
        static constexpr block mask_range_bits(block b, std::size_t l, std::size_t r) {
            return l < r ? b << (block_size - r) >> (block_size - r + l) << l : block(0);
        }

        std::size_t block_num() const {
            return blocks.size();
        }
    };
} // namespace suisen



#line 9 "library/linear_algebra/matrix_f2.hpp"

namespace suisen {
    struct MatrixF2 {
        MatrixF2() : MatrixF2(0, 0) {}
        MatrixF2(int n, int m, bool fill_value = false) : n(n), m(m), dat(n, DynamicBitSet(m, fill_value)) {}

        const DynamicBitSet& operator[](std::size_t i) const { return dat[i]; }
        DynamicBitSet& operator[](std::size_t i) { return dat[i]; }

        operator std::vector<DynamicBitSet>() const { return dat; }

        friend bool operator==(const MatrixF2& x, const MatrixF2& y) { return x.dat == y.dat; }
        friend bool operator!=(const MatrixF2& x, const MatrixF2& y) { return x.dat != y.dat; }

        std::pair<int, int> shape() const { return { n, m }; }
        int row_size() const { return n; }
        int col_size() const { return m; }

        MatrixF2 transposed() const {
            MatrixF2 t(m, n);
            for (std::size_t i = 0; i < n; ++i) for (std::size_t j = 0; j < m; ++j) t[j][i] = dat[i][j];
            return t;
        }

        friend MatrixF2& operator+=(MatrixF2& x, const MatrixF2& y) {
            assert(x.n == y.n and x.m == y.m);
            for (std::size_t i = 0; i < x.n; ++i) x[i] ^= y[i];
            return x;
        }
        friend MatrixF2& operator-=(MatrixF2& x, const MatrixF2& y) { return x += y; }
        friend MatrixF2& operator*=(MatrixF2& x, const MatrixF2& y) { return x = x * y; }
        friend MatrixF2& operator*=(MatrixF2& x, bool val) {
            if (not val) for (auto& row : x.dat) row.reset();
            return x;
        }
        friend MatrixF2& operator/=(MatrixF2& x, const MatrixF2& y) { return x = x * *y.inv(); }
        friend MatrixF2& operator/=(MatrixF2& x, bool val) {
            assert(val);
            return x;
        }
        friend MatrixF2 operator+(MatrixF2 x, const MatrixF2& y) { x += y; return x; }
        friend MatrixF2 operator-(MatrixF2 x, const MatrixF2& y) { x -= y; return x; }
        friend MatrixF2 operator*(const MatrixF2& x, MatrixF2 y) {
            y = y.transposed();
            assert(x.m == y.m);
            MatrixF2 z(x.n, y.n);
            for (std::size_t i = 0; i < x.n; ++i) for (std::size_t j = 0; j < y.n; ++j) {
                z[i][j] = (x[i] & y[j]).count() & 1;
            }
            return z;
        }
        friend MatrixF2 operator*(MatrixF2 x, bool val) { x *= val; return x; }
        friend MatrixF2 operator*(bool val, MatrixF2 x) { x *= val; return x; }
        friend MatrixF2 operator/(const MatrixF2 &x, const MatrixF2& y) { return x * *y.inv(); }
        friend MatrixF2 operator/(MatrixF2 x, bool val) { x /= val; return x; }

        DynamicBitSet operator*(const DynamicBitSet& x) const {
            assert(m == std::size_t(x.size()));
            DynamicBitSet y(n);
            for (std::size_t i = 0; i < n; ++i) y[i] = (dat[i] & x).count() & 1;
            return y;
        }

        MatrixF2 pow(long long b) const {
            assert(n == m);
            MatrixF2 p = *this, res = e1(n);
            for (; b; b >>= 1) {
                if (b & 1) res *= p;
                p *= p;
            }
            return res;
        }

        static MatrixF2 e0(std::size_t n) {
            return MatrixF2(n, n);
        }
        static MatrixF2 e1(std::size_t n) {
            MatrixF2 res(n, n);
            for (std::size_t i = 0; i < n; ++i) res[i][i] = 1;
            return res;
        }

        std::optional<MatrixF2> inv() const {
            assert(n == m);
            MatrixF2 A = *this, B = e1(n);
            for (std::size_t i = 0; i < n; ++i) {
                for (std::size_t j = i + 1; j < n; ++j) if (A[j][i]) {
                    std::swap(A[i], A[j]), std::swap(B[i], B[j]);
                    if (A[j][i]) A[j] ^= A[i], B[j] ^= B[i];
                }
                if (not A[i][i]) return std::nullopt;
            }
            for (std::size_t i = n; i-- > 0;) {
                for (std::size_t j = 0; j < i; ++j) {
                    if (A[j][i]) A[j] ^= A[i], B[j] ^= B[i];
                }
            }
            return B;
        }
        bool det() const {
            MatrixF2 A = *this;
            for (std::size_t i = 0; i < n; ++i) {
                for (std::size_t j = i + 1; j < n; ++j) if (A[j][i]) {
                    std::swap(A[i], A[j]);
                    if (A[j][i]) A[j] ^= A[i];
                }
                if (not A[i][i]) return false;
            }
            return true;
        }
        std::size_t rank() const {
            MatrixF2 A = *this;
            std::size_t r = 0;
            for (std::size_t j = 0; j < m; ++j) {
                for (std::size_t i = r + 1; i < n; ++i) if (A[i][j]) {
                    std::swap(A[r], A[i]);
                    if (A[i][j]) A[i] ^= A[r];
                }
                r += A[r][j];
            }
            return r;
        }

    private:
        std::size_t n, m;
        std::vector<DynamicBitSet> dat;
    };
} // namespace suisen



#line 5 "library/linear_algebra/gaussian_elimination_f2.hpp"

namespace suisen {
    struct SolutionSpace {
        SolutionSpace() : _solution(std::nullopt), _basis{} {}
        SolutionSpace(const std::optional<DynamicBitSet> &solution, const std::vector<DynamicBitSet> &basis = {}) : _solution(solution), _basis(basis) {}
        int dim() const { return _basis.size(); }
        int rank() const { return _solution->size() - _basis.size(); }
        bool has_solution() const { return bool(_solution); }
        const DynamicBitSet& solution() const { return *_solution; }
        const std::vector<DynamicBitSet>& basis() const { return _basis; }

        SolutionSpace normalized() {
            SolutionSpace res{ _solution };
            res += *this;
            res._normalized = true;
            return res;
        }
        void normalize() {
            *this = normalized();
        }
        bool add_base(DynamicBitSet v) {
            assert(_normalized);
            for (const auto& e : _basis) if (DynamicBitSet w = e ^ v; w < v) v = std::move(w);
            return v ? (_basis.push_back(v), true) : false;
        }
        friend SolutionSpace& operator+=(SolutionSpace& s, const SolutionSpace& t) {
            for (const DynamicBitSet& base : t._basis) s.add_base(base);
            return s;
        }
        friend SolutionSpace operator+(SolutionSpace s, const SolutionSpace& t) {
            s += t;
            return s;
        }
        friend SolutionSpace& operator&=(SolutionSpace& s, const SolutionSpace& t) {
            return s = s & t;
        }
        friend SolutionSpace operator&(const SolutionSpace& s, const SolutionSpace& t) {
            if (not s._normalized and t._normalized) return t & s;
            assert(s._normalized);
            std::vector<std::pair<DynamicBitSet, DynamicBitSet>> basis;
            for (const DynamicBitSet& e : s._basis) basis.emplace_back(e, e);
            SolutionSpace res;
            for (DynamicBitSet e : t._basis) {
                DynamicBitSet s(e.size());
                for (const auto& [v, t] : basis) if (DynamicBitSet w = e ^ v; w < e) e = std::move(w), s ^= t;
                if (e) basis.emplace_back(e, s);
                else res.add_base(s);
            }
            return res;
        }
    private:
        bool _normalized = false;
        std::optional<DynamicBitSet> _solution;
        std::vector<DynamicBitSet> _basis;
    };

    SolutionSpace gaussian_elimination_f2(MatrixF2 A, const DynamicBitSet &b) {
        assert(A.row_size() == b.size());
        const std::size_t n = A.row_size(), m = A.col_size();
        for (std::size_t i = 0; i < n; ++i) A[i].push_back(b[i]);

        std::vector<std::size_t> ones, zeros;
        std::size_t rank = 0;
        for (std::size_t j = 0; rank < n and j < m; ++j) {
            for (std::size_t i = rank + 1; i < n; ++i) if (A[i][j]) {
                std::swap(A[rank], A[i]);
                if (A[i][j]) A[i] ^= A[rank];
            }
            if (A[rank][j]) {
                for (std::size_t i = 0; i < rank; ++i) if (A[i][j]) A[i] ^= A[rank];
                ++rank;
                ones.push_back(j);
            } else {
                zeros.push_back(j);
            }
        }
        for (std::size_t i = rank; i < n; ++i) {
            if (A[i][m]) return SolutionSpace{};
        }
        DynamicBitSet solution(m);
        for (std::size_t i = 0; i < rank; ++i) {
            solution[ones[i]] = A[i][m];
        }
        std::vector<DynamicBitSet> basis;
        for (std::size_t j : zeros) {
            DynamicBitSet base(m);
            for (std::size_t i = 0; i < rank; ++i) base[ones[i]] = A[i][j];
            base[j] = 1;
            basis.push_back(std::move(base));
        }
        return SolutionSpace{ solution, basis };
    }
} // namespace suisen
Back to top page