GF2::GaussSolver — Linear system solve function

A standalone function that attempts to solve the system of linear equations \(A \cdot x = b\) is invoked as follows:

std::optional<GF2::Vector> GF2::solve(const GF2::Matrix &A, const GF2::Vector &b) (1)
1 A must be square and b must have the same size as there are rows in A.

If everything goes well the std::optional return value can be safely dereferenced as a bit-vector which will be a solution \(x\) to the system \(A \cdot x = b\). The solution may or may not be unique.

If there is a problem the return value will be a std::nullopt. This happens if the system of equations has no solution. It will also be the case if A is not square, or if the size of b is not the same as the number of rows in A.

The idea here is to get a solution for a system of equations with the least possible fuss.
Over \(\mathbb{F}_2\) any free variable can take on one of the two values 0 and 1. Hence if the system is consistent and has \(f\) free variables it will have \(2^f\) possible solutions. So a consistent system will have a unique solution only in the case \(A\) has full rank. The GaussSolver::operator() method can be used to iterate through potentially non-unique solutions if that is required.
Example
#include <GF2/GF2.h>
int
main()
{
    std::size_t m = 12;

    auto A = GF2::Matrix<>::random(m);
    auto b = GF2::Vector<>::random(m);
    auto x = GF2::solve(A, b);

    if(x) {
        // Check that x is indeed a solution by computing A.x and comparing that to b
        auto Ax = GF2::dot(A, *x);
        std::cout << "Matrix A, solution vector x, product A.x, and right hand side b\n";
        GF2::print(A, *x, Ax, b);
        std::cout << "So A.x == b? " << (Ax == b ? "YES" : "NO") << '\n';
    }
    else {
        std::cout << "System A.x = b has NO solutions for A and b as follows\n";
        GF2::print(A, b);
    }
}
Output for a consistent system (details depends on the values of the random inputs)
Matrix A, solution vector x, product A.x, and right hand side b
011000011001    0       0       0
000100011010    0       0       0
100001011010    0       0       0
111111010000    1       1       1
101011100101    1       1       1
100001111100    1       0       0
111100111110    0       0       0
101111011010    0       0       0
111100010110    1       1       1
011011010000    0       0       0
010011100101    0       1       1
000111101001    1       1       1
So A.x == b? YES
Output for an inconsistent system (details depends on the values of the random inputs)
System A.x = b has NO solutions for A and b as follows
110000000111    0
011111100001    1
011111000101    0
110110111011    1
100111001101    0
000010010010    1
001110011110    0
100010000001    0
110001110110    1
000100100010    0
001101100010    0
000000110000    0