Xor Algebra (Matrix mod 2) - YessineJallouli/Competitive-Programming GitHub Wiki

#include <bits/stdc++.h>
#define XOR_TEMPLATE template<size_t Bits=64>

using namespace std;
/*
 * Linear Algebra mod 2, works for N up to 4096
 */

// don't read M[i][j] directly use  char c; cin >> c; M[i][j] = (c == '1');

XOR_TEMPLATE
struct VecXor
{
    size_t n{};
    vector<bitset<Bits>> content;
    VecXor() = default;
    VecXor(size_t n,bool b={}): content((n+Bits-1)/Bits),n(n)
    {
        if(b) for(auto &c:content) c.flip();
    }

    size_t size() const
    {
        return n;
    }

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

    VecXor& operator^=(const VecXor & O)
    {
        for(int i=0;i< content.size();i++) content[i]^=O.content[i];
        return *this;
    }

    VecXor operator^(const VecXor & O) const
    {
        auto A=*this;
        return A^=O;
    }

    auto operator[](size_t t)
    {
        return content[t/Bits][t%Bits];
    }

    auto operator[](size_t t) const
    {
        return content[t/Bits][t%Bits];
    }

    bool inner_product(const VecXor & O) const
    {
        bool z{};
        for(int i=0;i<content.size();i++)
            z^=(content[i]&O.content[i]).count()&1;
        // Remove left-over
        if (!content.empty()) z^= ((content.back() >> Bits) & (O.content.back() >> Bits)).count()&1;
        return z;
    }
};

XOR_TEMPLATE
struct MatXor : vector<VecXor<Bits>>
{
    using Par = vector<VecXor<Bits>>;
    using Par::Par;
    using Par::size;
    using Par::data;
    MatXor(size_t rows,size_t cols,bool b=false): Par(rows,VecXor<Bits>(cols,b)){}
    static MatXor Id(size_t n) {
        MatXor I(n,n);
        for (int i=0;i<n;i++)
            I[i][i]=true;
        return I;
    }

    static MatXor Diag(vector<bool> v) {
        MatXor D(v.size(),v.size());
        for (int i=0;i<v.size();i++) D[i][i]=v[i];
        return D;
    }

    int rows() const {
        return size();
    }

    int cols() const {
        return size()?data()->size():0;
    }

    MatXor operator*(const MatXor& B) const
    {
        auto& A=*this;
        int n=A.rows(),m=B.cols();
        auto Bt=B.T(); // Transpose needed to speed up computation
        MatXor C(n,m);
        for (int i=0;i<n;i++) for (int j=0;j<m;j++) C[i][j]=A[i].inner_product(Bt[j]);
        return C;
    }

    VecXor<Bits> operator*(const VecXor<Bits>& U) const
    {
        auto& A=*this;
        int n=A.rows();
        VecXor<Bits> V(n);
        for (int i=0;i<n;i++) V[i]=A[i].inner_product(U);
        return V;
    }

    MatXor& operator*=(const MatXor& B) {
        return *this = *this * B;
    }

    // Trace of a matrix. Sum of its eigenvalues
    bool tr() const {
        bool w{};
        for (int i=0;i<min(rows(),cols());i++)
            w^=data()[i][i];
        return w;
    }

    // Transpose of a matrix
    MatXor T() const {
        MatXor Zt(cols(),rows());
        for (int i=0;i<rows();i++) for (int j=0;j<cols();j++)
                Zt[j][i]=data()[i][j];
        return Zt;
    }

    struct MatSolve {
        MatXor E,Y,X; // Results of the row echelon form
        size_t rank; // Rank of E
        bool dir; // Number of swaps modulo 2 in row_echelon is 0. Denotes that E and the original matrix have the same direction
        vector<pair<int,int>> mapper; // Mapping denoting the pivot cells
        MatSolve(MatXor E_,MatXor Y_,size_t rank,bool dir,vector<pair<int,int>> mapper):
                E(move(E_)),Y(move(Y_)),rank(rank),mapper(move(mapper)),dir(dir),X(E.cols(),Y.cols())
        {
        }

        void solve_general() // Solve the system, in general
        {
            int n=E.rows(),l=Y.cols();
            for (int i=rank;i<n;i++) for (int j=0;j<l;j++) if (Y[i][j])
                    {
                        X={};
                        return;
                    }
            for (int k=mapper.size()-1;k>=0;k--)
            {
                auto [r,i] = mapper[k];
                for (int j=0;j<r;j++) if (E[j][i])
                    {
                        Y[j]^=Y[r];
                        E[j]^=E[r];
                    }
            }
            for (auto [r,i]: mapper) X[i] = Y[r];
        }

        // Solve the system, assuming that E is square and det E != 0
        void solve_invertible()
        {
            int n=E.rows();
            X=Y;
            for (int i=n-1;i>=0;i--)
            {
                if (!E[i][i]) throw std::invalid_argument("Matrix is not invertible");
                for (int j=0;j<i;j++) if (E[j][i]) X[j]^=X[i];
            }
        }
        operator MatXor()  {
            return X;
        }
    };

    MatSolve row_echelon(MatXor Y) const
    {
        size_t rnk=0;
        auto A=*this;
        int n=rows(),m=cols();
        bool dir=true;
        vector<pair<int,int>> mapper;
        for (int i=0;i<m && rnk < n;i++)
        {
            int r=rnk;
            while (r<n && !A[r][i]) r++; // Pivot rule
            if (r==n) continue;
            mapper.emplace_back(rnk,i);
            if (r!=rnk)
            {
                swap(A[rnk],A[r]);
                swap(Y[rnk],Y[r]);
                dir=!dir;
            }
            for (int j=rnk+1;j<n;j++) if (A[j][i])
                {
                    A[j]^=A[rnk];
                    Y[j]^=Y[rnk];
                }
            rnk++;
        }
        return {A,Y,rnk,dir,mapper};
    }

    // Solve AX = Y, where A and Y are known
    MatXor solve(const MatXor& Y, bool invertible=false) const {
        invertible = invertible && rows() == cols();
        auto dec=row_echelon(Y);
        if (invertible) dec.solve_invertible();
        else dec.solve_general();
        return dec.X;
    }

    // Solve Ax = y, where A and b are known
    // 1. If A is known to be invertible, set
    VecXor<Bits> solve(const VecXor<Bits> &V , bool invertible =false) const {
        MatXor Y(V.size(),1);
        for (int i=0;i<V.size();i++) Y[i][0]=V[i];
        auto X=solve(Y,invertible);
        if (X.empty()) return {}; // If no solutions, return the empty vector
        VecXor<Bits> U(cols());
        for (int i=0;i<cols();i++) U[i]=X[i][0];
        return U;
    }

    // Basis of vectors (e_1,..,e_r) such that Ae_k=0 for all k
    MatXor null_basis() const {
        int n=rows(),m=cols();
        MatXor Z(*this);
        for (int i=0;i<m;i++) {
            Z.emplace_back(m);
            Z[rows()+i][i]=1;
        }
        auto C = Z.T().row_echelon(MatXor(m,0)).E;
        MatXor B;
        for (int i=0;i<m;i++)
        {
            bool all_zeros=true;
            for (int j=0;j<n && all_zeros;j++)
                all_zeros = !C[i][j];
            if (!all_zeros)
                continue;
            VecXor<Bits> u(m);
            for (int j=0;j<m;j++) u[j] = C[i][n+j];
            B.push_back(u);
        }
        return B;
    }

    // Basis induced by the matrix
    MatXor image_basis() const {
        auto dec = row_echelon(*this);
        dec.E.resize(dec.rank);
        return dec.E;
    }

    // Inverse of a square matrix.
    MatXor inv() const
    {
        return solve(MatXor::Id(rows()),true);
    }

    // Determinant of a square matrix
    bool det() const
    {
        auto dec = row_echelon(MatXor(rows(),0));
        bool w=true;
        for (int i=0;w && i<rows();i++) w= w && dec.E[i][i];
        return w;
    }

    // Rank of a matrix:
    // 1. The number of independent columns/rows in that matrix
    // 2. The dimension of vector space induced by the matrix
    // 3. Size of the basis generated by the matrix
    size_t rank() const {
        return row_echelon(MatXor(rows(),0)).rank;
    }

    // Nullity of the matrix:
    // 1. Size of the basis of elements that annihilate the matrix: Ax = 0
    // 2. S
    size_t nullity() const {
        return cols() - rank();
    }
};

int main() {
    int n,m; cin >> n >> m;
    MatXor A(n,m);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            char c;
            cin >> c;
            A[i][j] = (c == '1');
        }
    }
    MatXor B(n,1);
    for (int i = 0; i < n; i++) {
        char c; cin >> c;
        B[i][0] = (c == '1');
    }
    MatXor X = A.solve(B);
    if (X.empty()) {
        cout << -1 << '\n';
    }
    else {
        auto kernel = A.null_basis();
        cout << kernel.size() << '\n';
        for (int i = 0; i < m; i++) {
            cout << X[i][0];
        }
        cout << '\n';
        for (auto e : kernel) {
            for (int i = 0; i < m; i++) {
                cout << e[i];
            }
            cout << '\n';
        }
    }
}
⚠️ **GitHub.com Fallback** ⚠️