algorithm – Linear Algebra Library in C++;


This is actually an extension of the already written Matrix Library from this post.
This Matrix class is the result of the changes made thanks to this answer by Toby Speight, and having added a few other functionalities.

The library is composed of a few classes namely: a Fraction which contains the numbers that will be used in library, the Matrix class and the newly LA Vector class which contains functions such as:

bool is_linearly_dependent(std::initializer_list<Vector> vec_set);
bool is_linear_combination(std::initializer_list<Vector> vec_set, Vector vec);
bool spans_space(std::initializer_list<Vector> vec_set);
std::vector<Vector> row_space_basis(Matrix mx);
std::vector<Vector> null_space(Matrix mx);

The library is compiled in GCC 10.2.0, using boost format from boost 1.74.0, in Codeblocks on Windows 10.
While using boost format I ran into an unknown compiler error which I think I solved applying the changes suggested by this answer in this boostorg/format issue.

Fraction.h

#ifndef FRACTION_H_INCLUDED
#define FRACTION_H_INCLUDED

#include <iostream>
#include <ostream>
#include <cstring>
#include <assert.h>

class Fraction
{
    long long gcf(long long a, long long b);
    void simplify();

public:
    long long num;
    long long den;

    Fraction(long long _num = 0, long long _den = 1) : num{std::move(_num)}, den{std::move(_den)}
    {
        assert(_den != 0);
        simplify();
    }

    Fraction (Fraction n, Fraction d) : num(n.num * d.den), den(n.den * d.num)
    {
        assert(den != 0);
        simplify();
    }

    friend std::ostream& operator<< (std::ostream& os, const Fraction& fr);

    std::string to_string() const;

    bool is_integer()
    {
        return den == 1;
    }

    explicit operator bool() const
    {
        return num != 0;
    }

    bool operator== (const Fraction& fr) const
    {
        return num == fr.num && den == fr.den;
    }

    bool operator!= (const Fraction& fr) const
    {
        return !(*this == fr);
    }

    bool operator== (int n) const
    {
        return (n * den) == num;
    }

    bool operator!= (int n) const
    {
        return !(*this == n);
    }

    Fraction operator-() const
    {
        return { -num, den };
    }

    Fraction operator+() const
    {
        return *this;
    }

    long double to_double()
    {
        return double(num) / den;
    }

    float to_float()
    {
        return double(num) / den;
    }

    Fraction operator++()
    {
        num += den;
        return *this;
    }

    Fraction operator++(int)
    {
        Fraction fr = *this;
        ++(*this);
        return fr;
    }

    Fraction operator--()
    {
        num -= den;
        return *this;
    }

    Fraction operator--(int)
    {
        Fraction fr = *this;
        --(*this);
        return fr;
    }

    Fraction operator+(const Fraction& fr) const;
    Fraction operator/(const Fraction& fr) const;
    Fraction operator-(const Fraction& fr) const;
    Fraction operator*(const Fraction& fr) const;

    friend Fraction operator+(const Fraction& fr, long long n);
    friend Fraction operator+(long long n, const Fraction& fr);
    friend Fraction operator-(const Fraction& fr, long long n);
    friend Fraction operator-(long long n, const Fraction& fr);
    friend Fraction operator/(const Fraction& fr, long long n);
    friend Fraction operator/(long long n, const Fraction& fr);
    friend Fraction operator*(const Fraction& fr, long long n);
    friend Fraction operator*(long long n, const Fraction& fr);

    void operator+= (const Fraction& fr);
    void operator-= (const Fraction& fr);
    void operator*= ( const Fraction& fr);
    void operator/= (const Fraction& fr);

    void operator+=(long long n);
    void operator-=(long long n);
    void operator*=(long long n);
    void operator/=(long long n);
};

Fraction pow_fract(const Fraction& fr, int x);

#endif // FRACTION_H_INCLUDED

Fraction.cpp

#include "Fraction.h"

#include <iostream>
#include <ostream>
#include <sstream>

using namespace std;

std::ostream& operator<< (std::ostream& os, const Fraction& fr)
{
    if(fr.den == 1)
        os << fr.num;
    else
        os << fr.num << "/" << fr.den;

    return os;
}

string Fraction::to_string() const
{
    ostringstream os;
    os << *this;
    return os.str();
}

long long Fraction::gcf(long long a, long long b)
{
    if( b == 0)
        return abs(a);
    else
        return gcf(b, a%b);
}

void Fraction::simplify()
{
    if (den == 0 || num == 0)
    {
        num = 0;
        den = 1;
    }
    // Put neg. sign in numerator only.
    if (den < 0)
    {
        num *= -1;
        den *= -1;
    }

    // Factor out GCF from numerator and denominator.
    long long n = gcf(num, den);
    num = num / n;
    den = den / n;
}

Fraction Fraction::operator- (const Fraction& fr) const
{
    Fraction sub( (num * fr.den) - (fr.num * den), den * fr.den );

    int nu = sub.num;
    int de = sub.den;

    sub.simplify();

    return sub;
}

Fraction Fraction::operator+(const Fraction& fr) const
{
    Fraction addition ((num * fr.den) + (fr.num * den), den * fr.den );

    addition.simplify();

    return addition;
}

Fraction Fraction::operator*(const Fraction& fr) const
{
    Fraction multiplication(num * fr.num, den * fr.den);

    multiplication.simplify();

    return multiplication;
}

Fraction Fraction::operator / (const Fraction& fr) const
{
    Fraction sub(num * fr.den, den * fr.num);

    sub.simplify();

    return sub;
}

Fraction operator+(const Fraction& fr, long long n)
{
    return (Fraction(n) + fr);
}

Fraction operator+(long long n, const Fraction& fr)
{
    return (Fraction(n) + fr);
}

Fraction operator-(const Fraction& fr, long long n)
{
    return (fr - Fraction(n));
}

Fraction operator-(long long n, const Fraction& fr)
{
    return (Fraction(n) - fr);
}

Fraction operator/(const Fraction& fr, long long n)
{
    return (fr / Fraction(n));
}

Fraction operator/(long long n, const Fraction& fr)
{
    return (Fraction(n) / fr);
}

Fraction operator*(const Fraction& fr, long long n)
{
    return (Fraction(n) * fr);
}

Fraction operator*(long long n, const Fraction& fr)
{
    return (Fraction(n) * fr);
}

void Fraction::operator+=(const Fraction& fr)
{
   *this = *this + fr;
}

void Fraction::operator-=(const Fraction& fr)
{
   *this = *this - fr;
}

void Fraction::operator/=(const Fraction& fr)
{
  *this = *this / fr;
}

void Fraction::operator*=(const Fraction& fr)
{
    *this = *this * fr;
}

void Fraction::operator+=(long long n)
{
    *this = *this + n;
}

void Fraction::operator-=(long long n)
{
    *this = *this - n;
}

void Fraction::operator*=(long long n)
{
   *this = *this * n;
}

void Fraction::operator/=(long long n)
{
    *this = *this / n;
}

Fraction pow_fract(const Fraction& fr, int x)
{
    Fraction p(fr);

    for(int i = 0; i < x - 1; ++i)
        p *= fr;

    return p;
}

Matrix.h

#ifndef MATRIX_H_INCLUDED
#define MATRIX_H_INCLUDED

#include <vector>
#include <ostream>
#include <assert.h>
#include "Fraction.h"

namespace L_Algebra
{

class Matrix
{
private:
    std::size_t rows_num;
    std::size_t cols_num;

    std::vector<Fraction> data;

    Fraction& at(std::size_t r, std::size_t c)
    {
        return data.at( r * cols_num + c );
    }

    const Fraction& at(std::size_t r, std::size_t c) const
    {
        return data.at(r * cols_num + c);
    }

public:
    Matrix () = default;

    Matrix(std::size_t r, std::size_t c, Fraction n = 0 ) : rows_num(r), cols_num(c), data(r * c, n)
    {
        assert(r > 0 && c > 0);
    }

    Matrix(std::size_t r, std::size_t c, std::initializer_list<Fraction> values ) : rows_num(r), cols_num(c), data(values)
    {
        assert(r > 0 && c > 0);
        assert(values.size() == size());
    }

    Matrix(std::initializer_list<std::initializer_list<Fraction>> values );

    friend std::ostream& operator<<(std::ostream& out, const Matrix& mx);
    //friend std::vector<Fraction> operator<<(std::ostream& os, std::vector<Fraction> diag);

    explicit operator bool() const
    {
        return ! is_zero();
    }

    bool operator== (const Matrix& mx) const
    {
        return data == mx.data;
    }
    bool operator!= (const Matrix& mx) const
    {
        return !(*this == mx);
    }

    Matrix operator-()
    {
        return ( (*this) * (-1) );
    }

    Matrix operator+()
    {
        return (*this);
    }

    Matrix operator+(const Matrix& mx) const;
    Matrix operator-(const Matrix& mx) const;
    Matrix operator*(const Matrix& mx) const;

    Matrix& operator+=(const Matrix& mx);
    Matrix& operator-=(const Matrix& mx);
    Matrix& operator*=(const Matrix& mx);
    Matrix& operator*=(const Fraction& n);

    friend Matrix operator*(const Matrix& mx, Fraction n);
    friend Matrix operator*(Fraction n, const Matrix& mx);

    Matrix operator/(const Fraction& n) const;

    Fraction& operator()(std::size_t r, std::size_t c)
    {
        return at(r,c);
    }

    const Fraction& operator()(std::size_t r, std::size_t c) const
    {
        return at(r,c);
    }

    constexpr std::size_t size() const
    {
        return rows_num * cols_num;
    }

    void clear()
    {
        data.clear();
    }

    void resize(int r, int c, long long n = 0)
    {
        data.clear();

        data.resize( r * c, n );

        rows_num = r;
        cols_num = c;
    }

    size_t rows() const
    {
        return rows_num;
    }

    size_t cols() const
    {
        return cols_num;
    }

    static Matrix Identity(int n);
    static Matrix Constant(int r, int c, long long n);

    bool is_square() const
    {
        return rows_num == cols_num;
    }

    bool is_identity() const;
    bool is_symmetric() const;
    bool is_skewSymmetric() const;
    bool is_diagonal() const;
    bool is_zero() const;
    bool is_constant() const;
    bool is_orthogonal() const;
    bool is_invertible() const;
    bool is_linearly_dependent() const;
    bool is_linearly_independent() const;
    bool is_upperTriangular() const;
    bool is_lowerTriangular() const;
    bool is_consistent() const;

    Matrix transpose() const;
    Fraction determinant() const;
    Matrix inverse() const;
    Matrix adjoint() const;
    Matrix gaussElimination() const;
    Matrix gaussJordanElimination() const;
    Fraction trace() const;
    std::size_t rank() const;
    std::vector<Fraction> main_diagonal();
    std::vector<Fraction> secondary_diagonal();

    friend Matrix transitionMatrix(Matrix from, Matrix to);

private:
    void swapRows(int row1, int row2);
    bool pivotEqualTo_one_Found(int pivot_row, int pivot_col, int& alternative_pivot_row );
    bool pivotNot_zero_Found(int pivot_row, int pivot_col, int& col_dif_zero );
    bool firstNumberNot_zero(int row_num, int& num_coluna_num_dif_zero);
    void changePivotTo_one(int row_num, Fraction constant);
    void zeroOutTheColumn(int row_num, int num_pivot_row, Fraction constant);

    bool has_one_row_zero() const;
};

extern std::ostream& operator << (std::ostream& os,  const std::vector<Fraction>& v);

} // L_Algebra namespace



#endif // MATRIX_H_INCLUDED

Matrix.cpp

#include "Matrix.h"

#include <iostream>
#include <assert.h>
#include <algorithm>
#include <numeric>
#include <iomanip>
#include <boost/format.hpp>

using namespace std;

namespace L_Algebra
{

Matrix::Matrix(std::initializer_list<std::initializer_list<Fraction>> values )
{
    size_t len = 0;
    for(auto iter = values.begin(); iter != values.end(); ++iter)
        if(iter->size() != 0)
        {
            len = iter->size();
            break;
        }

    assert(len > 0);

    for(auto iter = values.begin(); iter != values.end(); ++iter)
    {
        if(iter->size() != 0)
            assert(iter->size() == len);

        if(iter->size() == 0)
            for(size_t i = 0; i < len; ++i)
                data.push_back(0);
        else
            for(auto iterj = iter->begin(); iterj != iter->end(); ++iterj)
                data.push_back(*iterj);
    }

    rows_num = values.size();
    cols_num = len;
}

bool Matrix::has_one_row_zero() const
{
    bool has;

    for(int i = 0; i < rows_num; ++i)
    {
        has = true;
        for(int j = 0; j < cols_num; ++j)
            if(at(i,j) != 0)
            {
                has = false;
                break;
            }

        if(has)
            return true;
    }

    return false;
}

ostream& operator<<(ostream& os, const Matrix& mx)
{
    size_t width = 1;
    for(const auto element : mx.data)
    {
        auto w = element.to_string().size();
        if(width < w)
            width = w;
    }

    string w = "%" + to_string(width + 4) + "d";

    for (int i = 0; i < mx.rows(); i++)
    {
        for (int j = 0; j < mx.cols(); j++)
            os << boost::format(w.c_str()) %  mx.at(i, j);

        os << 'n';
    }

    return os;
}

// to print the diagonal
std::ostream& operator<<(std::ostream& os,  const std::vector<Fraction>& v)
{
    for (auto e: v)
        os << e << " ";

    return os;
}

Matrix Matrix::operator+(const Matrix& mx) const
{
    assert(rows_num == mx.rows_num && cols_num == mx.cols_num);

    Matrix addition(rows_num, cols_num);

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            addition.at(i, j)= at(i, j) + mx.at(i, j);

    return addition;
}

Matrix Matrix::operator-(const Matrix& mx) const
{
    assert(rows_num == mx.rows_num && cols_num == mx.cols_num);

    Matrix sub(rows_num, cols_num);

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            sub.at(i, j) = at(i, j) - mx.at(i, j);

    return sub;
}

Matrix Matrix::operator*(const Matrix& mx) const
{
    assert(cols_num == mx.rows_num);

    Matrix multiplication(rows_num, mx.cols_num);

    for(int i = 0; i < rows_num; ++i)
        for (int j = 0; j < mx.cols_num; ++j)
            for(int x = 0; x < cols_num; ++x)
                multiplication.at(i,j) += at(i, x) * mx.at(x, j);

    return multiplication;
}

Matrix& Matrix::operator*=(const Matrix& mx)
{
    assert(cols_num == mx.rows_num);

    return *this = (*this * mx);
}

Matrix& Matrix::operator-=(const Matrix& mx)
{
    assert(rows_num == mx.rows_num && cols_num == mx.cols_num);

    transform(data.begin(), data.end(), mx.data.begin(), data.end(), minus{});

    return *this;
}

Matrix& Matrix::operator+=(const Matrix& mx)
{
    assert(rows_num == mx.rows_num && cols_num == mx.cols_num);

    transform(data.begin(), data.end(), mx.data.begin(), data.end(), plus{});

    return *this;
}

Matrix operator*(const Matrix& mx, Fraction n)
{
    Matrix multiplication(mx.rows_num, mx.cols_num);

    for(int i = 0; i < mx.rows_num; ++i)
        for(int j = 0; j < mx.cols_num; ++j)
            multiplication.at(i, j) = mx.at(i, j) * n;

    return multiplication;
}

Matrix operator*(Fraction n, const Matrix& mx)
{
    Matrix multiplication(mx.rows_num, mx.cols_num);

    for(int i = 0; i < mx.rows_num; ++i)
        for(int j = 0; j < mx.cols_num; ++j)
            multiplication.at(i, j) = mx.at(i, j) * n;

    return multiplication;
}

Matrix& Matrix::operator*=(const Fraction& n)
{
    return *this = *this * n;
}

Matrix Matrix::operator/(const Fraction& n) const
{
    assert(n != 0);

    Matrix division(rows_num, cols_num);

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            division.at(i, j) = at(i, j) / n;

    return division;
}

Matrix Matrix::Identity(int n)
{
    assert(n > 0);

    Matrix mx(n,n);

    for(int i = 0; i < n; ++i)
        mx.at(i, i) = 1;

    return mx;
}

Matrix Matrix::Constant(int r, int c, long long n)
{
    Matrix mx(r,c, n);

    return mx;
}

bool Matrix::is_identity() const
{
    if(! is_square())
        return false;

    return *this == Identity(cols_num);
}

bool Matrix::is_symmetric() const
{
    if(! is_square())
        return false;

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            if(at(i,j) != at(j,i))
                return false;

    return true;
}

bool Matrix::is_skewSymmetric() const
{
    if(! is_square())
        return false;

    for(int i = 0; i < rows_num; ++i)
        for(int j = i+1; j < cols_num; ++j)
            if(at(i,j) != -at(j,i))
                return false;

    return true;
}

bool Matrix::is_diagonal() const
{
    if(! is_square())
        return false;

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            if(i != j)
                if( at(i, j) != 0 )
                    return false;

    return true;
}

bool Matrix::is_zero() const
{
    return all_of( data.begin(), data.end(), ( ) (const auto& x)
    {
        return x == 0;
    } );
}

bool Matrix::is_constant() const
{
    return adjacent_find( data.begin(), data.end(), not_equal_to{} ) == data.end();
}

bool Matrix::is_orthogonal() const
{
    if(! is_square())
        return false;

    return (*this * transpose() == Identity(cols_num));
}

bool Matrix::is_invertible() const
{
    return this->determinant() != 0;
}

bool Matrix::is_linearly_dependent() const
{
    return this->determinant() == 0;
}

bool Matrix::is_linearly_independent() const
{
    return ! this->is_linearly_dependent();
}

bool Matrix::is_lowerTriangular() const
{
    if(! is_square())
        return false;

    for(int i = 0; i < rows_num; ++i)
        for(int j = i + 1; j < cols_num; ++j)
            if( at(i,j) )
                return false;

    return true;
}

bool Matrix::is_upperTriangular() const
{
    if(! is_square())
        return false;

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < i; ++j)
            if( at(i,j) )
                return false;

    return true;
}

bool Matrix::is_consistent( ) const
{
    Matrix mx1 = gaussJordanElimination();

    bool square = is_square();

    int num_non_zero_numbers = 0;
    for(int i = 0; i < rows_num; ++i)
    {
        if (square)
            for(int j = 0; j < cols_num; ++j)
            {
                if(mx1(i, j) != 0)
                    ++num_non_zero_numbers;
            }
        else
            for(int j = 0; j < cols_num - 1; ++j)
            {
                if(mx1(i, j) != 0)
                    ++num_non_zero_numbers;
            }

        if( ! square && num_non_zero_numbers == 0 && mx1(i, cols_num - 1) != 0)
            return false;

        if(num_non_zero_numbers > 1)
            return false;

        num_non_zero_numbers = 0;
    }

    return true;
}

Matrix Matrix::transpose() const
{
    Matrix trans(cols_num, rows_num);

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            trans.at(j, i) = at(i, j);

    return trans;
}

Fraction Matrix::trace() const
{
    assert(is_square());

    Fraction tr;
    for(int i = 0; i < rows_num; ++i)
        tr += at(i,i);

    return tr;
}

size_t Matrix::rank() const
{
    Matrix mx = this->gaussJordanElimination();

    int rank = 0;

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            if(mx(i, j) != 0)
            {
                ++rank;

                break;
            }

    return rank;
}

Fraction Matrix::determinant() const
{
    assert(is_square());

    if(is_zero())
        return {0};

    if(has_one_row_zero())
        return {0};

    if(rows_num == 1)
        return at(0,0);

    if(is_identity())
        return {1};

    if(is_constant())
        return {0};

    if(cols_num == 2)
        return at(0,0) * at(1,1) - at(0,1) * at(1,0);

    bool alternative_pivot_1_found;

    bool pivot_not_zero_found;

    bool number_not_zero_found;

    int row_with_alternative_pivot;

    int row_with_pivot_not_zero;

    int pivot_row = 0;
    int pivot_col = 0;

    Matrix mx(*this);
    vector<Fraction> row_mults;
    int sign = 1;

    while (pivot_row < (rows_num - 1))
    {
        alternative_pivot_1_found = mx.pivotEqualTo_one_Found ( pivot_row, pivot_col, row_with_alternative_pivot);

        pivot_not_zero_found = mx.pivotNot_zero_Found(pivot_row, pivot_col, row_with_pivot_not_zero);

        if (mx.at(pivot_row, pivot_col) != 1 && alternative_pivot_1_found )
        {
            mx.swapRows(pivot_row, row_with_alternative_pivot);

            sign *= (-1);
        }
        else if (mx.at(pivot_row, pivot_col) == 0 && pivot_not_zero_found )
        {
            mx.swapRows(pivot_row, row_with_pivot_not_zero);

            sign *= (-1);
        }

        int col_dif_zero;

        number_not_zero_found = mx.firstNumberNot_zero(pivot_row, col_dif_zero);

        if(number_not_zero_found)
        {
            if (mx.at(pivot_row, col_dif_zero) != 1)
            {
                row_mults.push_back(mx.at(pivot_row, col_dif_zero));

                mx.changePivotTo_one(pivot_row, mx.at(pivot_row, col_dif_zero));
            }
        }

        for (int i = pivot_row + 1; i < rows_num; ++i)
            mx.zeroOutTheColumn(i, pivot_row, mx.at(i, col_dif_zero));

        ++pivot_row;
        ++pivot_col;
    }

    Fraction det(sign);

    for(int i = 0; i < rows_num; ++i)
        det  *= mx.at(i,i);

    return accumulate(row_mults.begin(), row_mults.end(), det, multiplies());
}

Matrix Matrix::inverse() const
{
    assert(is_square());

    if( ! is_invertible())
        throw runtime_error("aNOT INVERTIBLEn");

    Matrix mx = *this;
    Matrix inverse = Matrix::Identity(rows_num);

    bool alternative_pivot_1_found;

    bool pivot_not_zero_found;

    bool number_not_zero_found;

    int row_with_alternative_pivot;

    int row_with_pivot_not_zero;

    int pivot_row = 0;
    int pivot_col = 0;

    //Gauss Elimination
    while (pivot_row < (rows_num - 1))
    {
        alternative_pivot_1_found = mx.pivotEqualTo_one_Found (pivot_row, pivot_col, row_with_alternative_pivot);

        pivot_not_zero_found = mx.pivotNot_zero_Found(pivot_row, pivot_col, row_with_pivot_not_zero);

        if (mx.at(pivot_row, pivot_col) != 1 && alternative_pivot_1_found )
        {
            inverse.swapRows(pivot_row, row_with_alternative_pivot);
            mx.swapRows(pivot_row, row_with_alternative_pivot);
        }
        else if (mx.at(pivot_row, pivot_col) == 0 && pivot_not_zero_found )
        {
            inverse.swapRows(pivot_row, row_with_pivot_not_zero);
            mx.swapRows(pivot_row, row_with_pivot_not_zero );
        }

        int col_dif_zero;

        number_not_zero_found = mx.firstNumberNot_zero(pivot_row, col_dif_zero);

        if(number_not_zero_found)
        {
            if ( mx.at(pivot_row, col_dif_zero) != 1)
            {
                inverse.changePivotTo_one(pivot_row, mx.at(pivot_row, col_dif_zero));
                mx.changePivotTo_one(pivot_row, mx.at(pivot_row, col_dif_zero));
            }
        }

        if(number_not_zero_found)
        {
            for (int i = pivot_row + 1; i < cols_num; ++i)
            {
                inverse.zeroOutTheColumn(i, pivot_row, mx.at(i, col_dif_zero));
                mx.zeroOutTheColumn(i, pivot_row, mx.at(i, col_dif_zero));
            }
        }

        ++pivot_row;
        ++pivot_col;
    }

    //Jordan Elimination
    while(pivot_row > 0)
    {
        int col_dif_zero;

        number_not_zero_found = mx.firstNumberNot_zero(pivot_row, col_dif_zero);

        if(number_not_zero_found)
        {
            if ( mx.at(pivot_row, col_dif_zero) != 1)
            {
                inverse.changePivotTo_one(pivot_row, mx.at(pivot_row, col_dif_zero));
                mx.changePivotTo_one(pivot_row, mx.at(pivot_row, col_dif_zero));
            }
        }

        if(number_not_zero_found)
        {
            for (int i = pivot_row - 1; i >= 0; --i)
            {
                inverse.zeroOutTheColumn(i, pivot_row, mx.at(i, col_dif_zero));
                mx.zeroOutTheColumn(i, pivot_row, mx.at(i, col_dif_zero));
            }
        }

        --pivot_row;
    }

    return inverse;
}

Matrix Matrix::adjoint() const
{
    assert(is_square());
    assert(cols_num > 1);

    if(is_zero())
        return Matrix(rows_num, cols_num);

    if(is_constant())
        return Matrix(rows_num, cols_num);

    if(is_identity())
        return *this;

    Matrix cofact(rows_num, cols_num);

    int r = 0, c = 0;

    Matrix temp(rows_num - 1, cols_num - 1);

    for(int i = 0; i < rows_num; ++i)
    {
        for(int j = 0; j < cols_num; ++j)
        {
            for(int k = 0; k < rows_num; ++k)
            {
                for(int h = 0; h < cols_num; ++h)
                {
                    if (k != i && h != j)
                    {
                        temp(r, c++) = at(k, h);

                        if(c == cols_num - 1)
                        {
                            c = 0;
                            ++r;
                        }
                    }
                }
            }

            c = 0;
            r = 0;

            int sign;

            sign = ( ( i + j ) % 2 == 0 ) ? 1 : -1;

            cofact.at(i, j) = sign * temp.determinant();
        }
    }

    return cofact.transpose();
}

Matrix Matrix::gaussJordanElimination() const
{
    Matrix mx = *this;

    bool alternative_pivot_1_found;

    bool pivot_not_zero_found;

    bool number_not_zero_found;

    int row_with_alternative_pivot;

    int row_with_pivot_not_zero;

    int pivot_row = 0;
    int pivot_col = 0;

    ///Gauss Elimination
    while (pivot_row < (rows_num - 1) && pivot_row < (cols_num))
    {
        alternative_pivot_1_found = mx.pivotEqualTo_one_Found ( pivot_row, pivot_col,
                                    row_with_alternative_pivot);

        pivot_not_zero_found = mx.pivotNot_zero_Found(
                                   pivot_row, pivot_col, row_with_pivot_not_zero);

        if (mx.at( pivot_row, pivot_col) != 1 && alternative_pivot_1_found )
        {
            mx.swapRows(pivot_row, row_with_alternative_pivot);
        }
        else if (mx.at( pivot_row, pivot_col) == 0 && pivot_not_zero_found )
        {
            mx.swapRows( pivot_row, row_with_pivot_not_zero );
        }

        int col_dif_zero;

        number_not_zero_found = mx.firstNumberNot_zero(pivot_row, col_dif_zero);

        if(number_not_zero_found)
        {
            if (( mx.at(pivot_row, col_dif_zero) ) != 1)
            {
                mx.changePivotTo_one(pivot_row,
                                     mx.at(pivot_row, col_dif_zero) );
            }
        }

        if(number_not_zero_found)
        {
            for(int i = pivot_row + 1; i < rows_num; ++i)
            {
                mx.zeroOutTheColumn( i, pivot_row, mx.at(i, col_dif_zero));
            }
        }

        ++pivot_row;
        ++pivot_col;
    }

    //Jordan Elimination
    while(pivot_row > 0)
    {
        int col_dif_zero;

        number_not_zero_found = mx.firstNumberNot_zero(pivot_row, col_dif_zero);

        if(number_not_zero_found)
            if ( mx.at(pivot_row, col_dif_zero) != 1)
            {
                mx.changePivotTo_one(pivot_row, mx.at(pivot_row, col_dif_zero));
            }

        if(number_not_zero_found)
            for (int i = pivot_row - 1; i >= 0; --i)
                mx.zeroOutTheColumn(i, pivot_row, mx.at(i, col_dif_zero));

        --pivot_row;
    }

    return mx;
}

Matrix Matrix::gaussElimination() const
{
    Matrix mx = *this;

    bool alternative_pivot_1_found;

    bool pivot_not_zero_found;

    bool number_not_zero_found;

    int row_with_alternative_pivot;

    int row_with_pivot_not_zero;

    int pivot_row = 0;
    int pivot_col = 0;

    ///Gauss Elimination
    while (pivot_row < (rows_num - 1) && pivot_row < (cols_num) )
    {
        alternative_pivot_1_found = mx.pivotEqualTo_one_Found ( pivot_row, pivot_col,
                                    row_with_alternative_pivot);

        pivot_not_zero_found = mx.pivotNot_zero_Found(
                                   pivot_row, pivot_col, row_with_pivot_not_zero);

        if (mx.at( pivot_row, pivot_col) != 1 && alternative_pivot_1_found )
        {
            mx.swapRows(pivot_row, row_with_alternative_pivot);
        }
        else if (mx.at( pivot_row, pivot_col) == 0 && pivot_not_zero_found )
        {
            mx.swapRows( pivot_row, row_with_pivot_not_zero );
        }

        int col_dif_zero;

        number_not_zero_found = mx.firstNumberNot_zero(pivot_row, col_dif_zero);

        if(number_not_zero_found)
        {
            if (( mx.at(pivot_row, col_dif_zero) ) != 1)
            {
                mx.changePivotTo_one(pivot_row,
                                     mx.at(pivot_row, col_dif_zero) );
            }
        }

        if(number_not_zero_found)
        {
            for(int i = pivot_row + 1; i < rows_num; ++i)
            {
                mx.zeroOutTheColumn( i, pivot_row, mx.at(i, col_dif_zero));
            }
        }

        ++pivot_row;
        ++pivot_col;
    }

    int col_dif_zero;

    number_not_zero_found = mx.firstNumberNot_zero(pivot_row, col_dif_zero);

    if(number_not_zero_found)
        if ( mx.at(pivot_row, col_dif_zero) != 1)
        {
            mx.changePivotTo_one(pivot_row, mx.at(pivot_row, col_dif_zero));
        }

    return mx;
}

vector<Fraction> Matrix::main_diagonal()
{
    assert(is_square());

    vector<Fraction> diag;

    for(int i = 0; i < rows_num; ++i)
        diag.push_back(at(i,i));

    return diag;
}

vector<Fraction> Matrix::secondary_diagonal()
{
    assert(is_square());

    vector<Fraction> diag;

    for(int i = 0, j = rows_num - 1; i < rows_num; ++i, --j)
        diag.push_back(at(i,j));

    return diag;
}

void Matrix::swapRows( int row1, int row2)
{
    for (int i = 0; i < cols_num; i++ )
        std::swap( at(row1,i ), at(row2, i) );
}

bool Matrix::pivotEqualTo_one_Found( int pivot_row, int pivot_col, int& alternative_pivot_row )
{
    for (int i = pivot_row + 1; i < rows_num; ++i)
    {
        if(at(i, pivot_col) == 1)
        {
            alternative_pivot_row = i;

            return true;
        }
    }

    return false;
}

bool Matrix::pivotNot_zero_Found(int pivot_row, int pivot_col,int& col_dif_zero )
{
    for (int i = pivot_row + 1; i < rows_num; ++i)
        if(at(i, pivot_col) != 0)
        {
            col_dif_zero = i;

            return true;
        }

    return false;
}

bool Matrix::firstNumberNot_zero(int row_num, int& num_coluna_num_dif_zero)
{
    for (int i = 0; i < cols_num; ++i)
        if (at(row_num, i) != 0)
        {
            num_coluna_num_dif_zero = i;

            return true;
        }

    return false;
}

void Matrix::changePivotTo_one( int row_num, Fraction constant)
{
    for(int i = 0; i < cols_num; ++i)
        if (at(row_num, i).num != 0)
            at(row_num, i) = (at(row_num, i) / constant);
}

void Matrix::zeroOutTheColumn( int row_num, int num_pivot_row, Fraction constant)
{
    for(int i = 0; i < cols_num; ++i)
        at(row_num, i) = at(row_num, i) -  (constant * at(num_pivot_row, i));
}

}// L_Algebra namespace

LA_Vector.h

#ifndef LA_VECTOR_H
#define LA_VECTOR_H

#include "Fraction.h"
#include "Matrix.h"
#include <initializer_list>
#include <deque>
#include <ostream>

namespace L_Algebra
{

class Vector
{
    std::deque<Fraction> data;

    Fraction& at(std::size_t i)
    {
        return data.at(i);
    }

    const Fraction& at(std::size_t i) const
    {
        return data.at(i);
    }

    void push_back(Fraction n)
    {
        data.push_back(n);
    }

    friend std::vector<Vector> null_space(Matrix mx);
    friend std::vector<Vector> null_space_(Matrix mx);

public:
    Vector() = default;

    Vector(std::vector<int> d)
    {
        assert(d.size() > 0);

        for(auto const &e: d)
            data.push_back(e);
    }

    Vector(std::deque<int> d)
    {
        assert(d.size() > 0);

        for(auto const &e: d)
            data.push_back(e);
    }

    Vector(std::vector<Fraction> d)
    {
        assert(d.size() > 0);

        for(auto const &e: d)
            data.push_back(e);
    }

    Vector(std::deque<Fraction> d) : data(d)
    {
        assert(data.size() > 0);
    }

    Vector(int d) : data(d, 0)
    {
        assert(data.size() > 0);
    }

    Vector(int d, long long int n) : data(d, n)
    {
        assert(data.size() > 0);
    }

    Vector(std::initializer_list<Fraction> values) : data(values)
    {
        assert(data.size() > 0);
    }

    friend std::ostream& operator<< (std::ostream& os, const Vector& lav);

    explicit operator bool() const
    {
        return dimension() != 0;
    }

    bool operator==(const Vector& lav) const
    {
        return data == lav.data;
    }

    bool operator!=(const Vector& lav) const
    {
        return data != lav.data;
    }

    Fraction& operator()(size_t i)
    {
        return data.at(i);
    }

    const Fraction& operator()(size_t i) const
    {
        return data.at(i);
    }

    Vector operator+(const Vector& lav) const;
    Vector operator-(const Vector& lav) const;
    Vector operator->*(const Vector& lav) const; // vectorial product
    Fraction operator*(const Vector& lav) const; // dot product

    Vector& operator+=(const Vector& lav);
    Vector& operator-=(const Vector& lav);

    friend Vector operator*(const Vector& mx, Fraction n);
    friend Vector operator*(Fraction n, const Vector& mx);

    std::size_t dimension() const
    {
        return data.size();
    }

    Fraction norm_Power2() const;
    double norm() const;
};

Vector proj(Vector u, Vector a);
Vector proj_orthogonal(Vector u, Vector a);

bool is_orthogonal(std::initializer_list<Vector> vec_set);

bool is_linearly_dependent(std::initializer_list<Vector> vec_set);
bool is_linearly_dependent(std::initializer_list<Matrix> matrices_set);
bool is_linearly_independent(std::initializer_list<Vector> vec_set);
bool is_linearly_independent(std::initializer_list<Matrix> matrices_set);

bool is_linear_combination(std::initializer_list<Vector> vec_set, Vector vec);
bool is_linear_combination(std::initializer_list<Matrix> matrices_set, Matrix mx);

bool spans_space(std::initializer_list<Vector> vec_set);
bool spans_space(std::initializer_list<Matrix> matrix_set);
bool is_in_span(Vector vec, std::initializer_list<Vector> span);

bool is_basis(std::initializer_list<Vector> vec_set);
bool is_basis(std::initializer_list<Matrix> matrices_set);

Vector change_basis(Vector vec, std::initializer_list<Vector> basis_from, std::initializer_list<Vector> basis_to);
Vector change_basis(Vector vec_in_standard_basis, std::initializer_list<Vector> destination_basis);

std::vector<Vector> row_space_basis(Matrix mx);
std::vector<Vector> column_space_basis(Matrix mx);
std::vector<Vector> null_space(Matrix mx);

std::size_t row_space_dim(Matrix mx);
std::size_t column_space_dim(Matrix mx);
std::size_t nullity(Matrix mx);

Vector coordinate_vector_relative_to_basis(std::initializer_list<Vector> basis, Vector vec);
Vector vector_with_coordinate_relative_to_basis(std::initializer_list<Vector> basis, Vector coordinate_vec);

Matrix vectorsToMatrix(std::vector<Vector>vec_set);

Matrix turnMatricesIntoLinearCombination(std::vector<Matrix>matrix_set);

/*
Vector rowOfMatrixToVector(Matrix mx, int row);
Vector columnOfMatrixToVector(Matrix mx, int column);
*/

} // L_Algebra namespace

#endif // LA_VECTOR_H

LA_Vector.cpp

#include "LA_Vector.h"

#include <iostream>
#include <math.h>
#include <assert.h>
#include <set>
#include <deque>
#include <algorithm>

using namespace std;

namespace L_Algebra
{

Matrix transitionMatrix(Matrix from, Matrix to)
{
    assert(from.size() == to.size());

    int rows_num = to.rows();
    int cols_num = to.cols();

    bool alternative_pivot_1_found;

    bool pivot_not_zero_found;

    bool number_not_zero_found;

    int row_with_alternative_pivot;

    int row_with_pivot_not_zero;

    int pivot_row = 0;
    int pivot_col = 0;

    //Gauss Elimination
    while (pivot_row < (rows_num - 1))
    {
        alternative_pivot_1_found = to.pivotEqualTo_one_Found (pivot_row, pivot_col, row_with_alternative_pivot);

        pivot_not_zero_found = to.pivotNot_zero_Found(pivot_row, pivot_col, row_with_pivot_not_zero);

        if (to.at(pivot_row, pivot_col) != 1 && alternative_pivot_1_found )
        {
            from.swapRows(pivot_row, row_with_alternative_pivot);
            to.swapRows(pivot_row, row_with_alternative_pivot);
        }
        else if (to.at(pivot_row, pivot_col) == 0 && pivot_not_zero_found )
        {
            from.swapRows(pivot_row, row_with_pivot_not_zero);
            to.swapRows(pivot_row, row_with_pivot_not_zero );
        }

        int col_dif_zero;

        number_not_zero_found = to.firstNumberNot_zero(pivot_row, col_dif_zero);

        if(number_not_zero_found)
        {
            if ( to.at(pivot_row, col_dif_zero) != 1)
            {
                from.changePivotTo_one(pivot_row, to.at(pivot_row, col_dif_zero));
                to.changePivotTo_one(pivot_row, to.at(pivot_row, col_dif_zero));
            }
        }

        if(number_not_zero_found)
        {
            for (int i = pivot_row + 1; i < cols_num; ++i)
            {
                from.zeroOutTheColumn(i, pivot_row, to.at(i, col_dif_zero));
                to.zeroOutTheColumn(i, pivot_row, to.at(i, col_dif_zero));
            }
        }

        ++pivot_row;
        ++pivot_col;
    }

    //Jordan Elimination
    while(pivot_row > 0)
    {
        int col_dif_zero;

        number_not_zero_found = to.firstNumberNot_zero(pivot_row, col_dif_zero);

        if(number_not_zero_found)
        {
            if ( to.at(pivot_row, col_dif_zero) != 1)
            {
                from.changePivotTo_one(pivot_row, to.at(pivot_row, col_dif_zero));
                to.changePivotTo_one(pivot_row, to.at(pivot_row, col_dif_zero));
            }
        }

        if(number_not_zero_found)
        {
            for (int i = pivot_row - 1; i >= 0; --i)
            {
                from.zeroOutTheColumn(i, pivot_row, to.at(i, col_dif_zero));
                to.zeroOutTheColumn(i, pivot_row, to.at(i, col_dif_zero));
            }
        }

        --pivot_row;
    }

    return from;
}

bool is_consistent(const Matrix& mx)
{
    int rows_num = mx.rows();
    int cols_num = mx.cols();

    Matrix mx1 = mx.gaussJordanElimination();

    bool square = mx.is_square();

    int num_non_zero_numbers = 0;
    for(int i = 0; i < rows_num; ++i)
    {
        if (square)
            for(int j = 0; j < cols_num; ++j)
            {
                if(mx1(i, j) != 0)
                    ++num_non_zero_numbers;
            }
        else
            for(int j = 0; j < cols_num - 1; ++j)
            {
                if(mx1(i, j) != 0)
                    ++num_non_zero_numbers;
            }

        if(num_non_zero_numbers > 1)
            return false;

        if( ! square && num_non_zero_numbers == 0 && mx1(i, cols_num - 1) != 0)
            return false;

        num_non_zero_numbers = 0;
    }

    return true;
}

Matrix vectorsToMatrix(std::vector<Vector>vec_set)
{
    assert(vec_set.size() > 0);

    int len = vec_set.size();
    for(int i = 0; i < len; ++i)
        assert(vec_set(i).dimension() == vec_set(0).dimension());

    int rows_num = vec_set(0).dimension();
    int cols_num = len;

    Matrix mx(rows_num, cols_num);

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
        {
            mx(i, j) = vec_set.at(j)(i);
        }

    return mx;
}

Matrix turnMatricesIntoLinearCombination(std::vector<Matrix>matrix_set)
{
    assert(matrix_set.size() > 0);

    int len = matrix_set.size();
    for(int i = 0; i < len; ++i)
        assert(matrix_set(i).size() == matrix_set(0).size());
    /*
        int rows_num = matrix_set(0).size();
        int cols_num = len;

        int r = matrix_set(0).rows();
        int c = matrix_set(0).cols();

        Matrix m(rows_num, cols_num);

        Vector lav(r * c);

        size_t vec_lav_size = cols_num;
        vector<Vector> vec_lav(vec_lav_size, r * c);

        // pass the values from the set of matrices to a set of la_vectors
        int ind = 0;
        for(size_t h = 0; h < vec_lav_size; ++h)
        {
            for(int i = 0; i < r; ++i)
                for(int j = 0; j < c; ++j)
                    vec_lav.at(h)(ind++) = matrix_set.at(h)(i, j);

            ind = 0;
        }

         transform the values from the set of the matrices into a new matrix;
        for(int i = 0; i < rows_num; ++i)
            for(int j = 0; j < cols_num; ++j)
                m(i, j) = vec_lav.at(j)(i);
    */

    int rows_num = matrix_set(0).size();
    int cols_num = len;

    int r = matrix_set(0).rows();
    int c = matrix_set(0).cols();

    Matrix m(rows_num, cols_num);

    for(int i = 0; i < cols_num; ++i)
    {
        int id = 0;

        for(int x = 0; x < r; ++x)
        {
            for(int y = 0; y < c; ++y)
            {
                m(id++, i) = matrix_set( i )(x, y);
            }
        }
    }


    return m;
}

Vector rowOfMatrixToVector(const Matrix& mx, int row)
{
    assert(row <= mx.rows());

    int cols_num = mx.cols();

    Vector v(cols_num);

    for(int i = 0; i < cols_num; ++i)
        v( i ) = mx(row, i);

    return v;
}

Vector columnOfMatrixToVector(const Matrix& mx, int column)
{
    assert(column <= mx.cols());

    int rows_num = mx.rows();

    Vector v(rows_num);

    for(int i = 0; i < rows_num; ++i)
        v( i ) = mx(i, column);

    return v;
}

ostream& operator<< (ostream& os, const Vector& lav)
{
    os << "(";

    for(auto el : lav.data)
        os << el << ", ";

    if(lav.data.empty())
        os << " )";
    else
        os << "bb b" << ")";

    return os;
}

Vector Vector::operator+(const Vector& lav) const
{
    size_t len = data.size();

    assert(len == lav.data.size());

    Vector addition;

    addition.data.resize(len, 0);

    for(size_t i = 0; i < len; ++i)
        addition(i) = at(i) + lav(i);

    return addition;
}

Vector& Vector::operator+=(const Vector& lav)
{
    return *this = *this + lav;
}

Vector Vector::operator-(const Vector& lav) const
{
    size_t len = data.size();

    assert(len == lav.data.size());

    Vector subtraction;

    subtraction.data.resize(data.size(), 0);

    for(size_t i = 0; i < len; ++i)
        subtraction(i) = at(i) - lav(i);

    return subtraction;
}

Vector& Vector::operator-=(const Vector& lav)
{
    return *this = *this - lav;
}

Fraction Vector::operator*(const Vector& lav) const // dot product
{
    size_t len = data.size();

    assert(len == lav.data.size());

    Fraction dot_prod;

    for(size_t i = 0; i < len; ++i)
        dot_prod += at(i) * lav(i);

    return dot_prod;
}

// vectorial product
Vector Vector::operator->*(const Vector& lav) const
{
    size_t len = data.size();

    assert( (len == lav.data.size()) && len == 3);

    return {at(1) * lav.at(2) - at(2) * lav.at(1),
            - (at(2) * lav.at(0) - at(0) * lav.at(2)),
            at(0) * lav.at(1) - at(1) * lav.at(0) };
}

Vector operator*(const Vector& lav, Fraction n)
{
    Vector mult;

    mult.data.resize(lav.data.size(), 0);

    int i = 0;
    for( auto el : lav.data)
        mult.at(i++) = el * n;

    return mult;
}

Vector operator*(Fraction n, const Vector& lav)
{
    Vector mult;

    mult.data.resize(lav.data.size(), 0);

    int i = 0;
    for( auto el : lav.data)
        mult.at(i++) = el * n;

    return mult;
}

double Vector::norm() const
{
    Fraction n;

    size_t len = dimension();

    for(size_t i = 0; i < len; ++i)
        n += pow_fract(at(i), 2);

    return sqrt(n.to_double());
}

Fraction Vector::norm_Power2() const
{
    Fraction n;

    size_t len = dimension();

    for(size_t i = 0; i < len; ++i)
        n += pow_fract(at(i), 2);

    return n;
}

bool is_orthogonal(std::initializer_list<Vector> vec_set)
{
    assert(vec_set.size() > 1);

    std::vector<Vector> vec(vec_set);

    size_t len = vec.size();

    for(size_t i = 0; i < len; ++i )
        assert(vec.at(i).dimension() == vec.at(0).dimension());

    for( size_t i = 0; i < len - 1; ++i)
        for( size_t j = i + 1; j < len; ++j)
            if (vec.at(i) * vec.at(j) == 0)
                return true;

    return false;
}

Vector proj(Vector u, Vector a)
{
    return Fraction(u*a, a.norm_Power2()) * a;
}

Vector proj_orthogonal(Vector u, Vector a)
{
    return u - proj(u, a);
}

bool is_linearly_dependent(std::initializer_list<Vector> vec_set)
{
    Matrix mx = vectorsToMatrix(vec_set).gaussJordanElimination();

    int rows_num = mx.rows();
    int cols_num = mx.cols();

    int num_non_zero_numbers = 0;
    for(int i = 0; i < rows_num; ++i)
    {
        for(int j = 0; j < cols_num; ++j)
        {
            if(mx(i, j) != 0)
                ++num_non_zero_numbers;
        }

        if(num_non_zero_numbers > 1)
            return true;

        num_non_zero_numbers = 0;
    }

    return false;
}

bool is_linearly_dependent(initializer_list<Matrix> matrices_set)
{
    assert(matrices_set.size() > 0);

    vector<Matrix> vecs(matrices_set);

    int len = vecs.size();
    for(int i = 0; i < len; ++i)
        assert(vecs(i).size() == vecs(0).size() && vecs(i).size() > 0);

    int r = vecs(0).rows();
    int c = vecs(0).cols();

    Matrix mx(r, c);

    vecs.push_back(mx);

    Matrix m = turnMatricesIntoLinearCombination(vecs);

    if( is_consistent(m))
        return false;
    else
        return true;
}

bool is_linearly_independent(std::initializer_list<Vector>vec_set)
{
    return ! is_linearly_dependent(vec_set);
}

bool is_linearly_independent(initializer_list<Matrix> matrices_set)
{
    return ! is_linearly_dependent(matrices_set);
}

bool is_linear_combination(std::initializer_list<Vector> vec_set, Vector vec)
{
    vector<Vector> vecs(vec_set);

    vecs.push_back(vec);

    Matrix mx = vectorsToMatrix(vecs);

    if( ! is_consistent(mx))
        return false;

    mx = mx.gaussJordanElimination();

    int rows_num = mx.rows();
    int cols_num = mx.cols();

    Vector results = columnOfMatrixToVector(mx, cols_num - 1);

    Vector combination(rows_num);

    for(int i = 0; i < rows_num; ++i)
    {
        for(int j = 0; j < cols_num - 1; ++j)
            combination(i) += results(j) * vecs.at(j)(i);
    }

    if(vec == combination)
        return true;
    else
        return false;
}

bool is_linear_combination(std::initializer_list<Matrix> matrices_set, Matrix mx)
{
    assert(matrices_set.size() > 0);

    vector<Matrix> vecs(matrices_set);
    vecs.push_back(mx);

    Matrix m = turnMatricesIntoLinearCombination(vecs);

    int cols_num = m.cols();

    vector<Vector> vec_lav(cols_num);

    for(int i = 0; i < cols_num; ++i)
        vec_lav(i) = columnOfMatrixToVector(m, i);

    if( ! is_consistent(m))
        return false;

    m = m.gaussJordanElimination();

    Vector results = columnOfMatrixToVector(m, cols_num - 1);

    Vector combination(m.rows());

    for(int i = 0; i < cols_num - 1; ++i)
        combination += results(i) * vec_lav.at(i);

    Vector lav = vec_lav(vec_lav.size() - 1);

    if(lav == combination)
        return true;
    else
        return false;
}

bool is_basis(std::initializer_list<Vector> vec_set)
{
    assert(vec_set.size() > 0);

    vector<Vector> vec(vec_set);

    int len = vec.size();
    for(int i = 0; i < len; ++i)
        assert(vec(i).dimension() == vec(0).dimension());

    if(vec.size() != vec(0).dimension())
        return false;

    return ! is_linearly_dependent(vec_set);
}

bool is_basis(std::initializer_list<Matrix> matrices_set)
{
    return ! is_linearly_dependent(matrices_set);
}

Vector change_basis(Vector vec, std::initializer_list<Vector> basis_from,
                    std::initializer_list<Vector> basis_to)
{
    assert(basis_to.size() == basis_from.size());
    assert(vec.dimension() == basis_from.size());

    Matrix from = vectorsToMatrix(basis_from);
    Matrix to = vectorsToMatrix(basis_to);

    Matrix transition_matrix = transitionMatrix(from, to);

    int vec_dimension = vec.dimension();

    Matrix vec_matrix(vec_dimension, 1);

    for(int i = 0; i < vec_dimension; ++i)
        vec_matrix(i,0) = vec(i);

    Matrix new_basis_vec_matrix = transition_matrix * vec_matrix;

    Vector vec_in_new_basis(vec_dimension);

    for(int i = 0; i < vec_dimension; ++i)
        vec_in_new_basis(i) = new_basis_vec_matrix(i,0);

    return vec_in_new_basis;
}

Vector change_basis(Vector vec_in_standard_basis, std::initializer_list<Vector> destination_basis)
{
    return coordinate_vector_relative_to_basis(destination_basis, vec_in_standard_basis);
}

bool spans_space(std::initializer_list<Vector> vec_set)
{
    return ! is_linearly_dependent(vec_set);
}

bool spans_space(std::initializer_list<Matrix> matrix_set)
{
    return ! is_linearly_dependent(matrix_set);
}

bool is_in_span(Vector vec, std::initializer_list<Vector> span)
{
    return is_linear_combination(span, vec);
}

Vector coordinate_vector_relative_to_basis(std::initializer_list<Vector> basis,
        Vector vec)
{
    assert(basis.size() == vec.dimension());

    vector<Vector> vecs(basis);

    vecs.push_back(vec);

    Matrix mx = vectorsToMatrix(vecs);

    mx = mx.gaussJordanElimination();

    int rows_num = mx.rows();
    int cols_num = mx.cols();

    if(! is_consistent(mx))
        throw runtime_error("the basis is linearly dependent");

    Vector coordinate_vector(rows_num);

    for(int i = 0; i < rows_num; ++i)
        coordinate_vector(i) = mx(i, cols_num - 1);

    return coordinate_vector;
}

Vector vector_with_coordinate_relative_to_basis(initializer_list<Vector> basis,
        Vector coordinate_vec)
{
    assert(basis.size() > 0);

    assert(coordinate_vec.dimension() == basis.size());

    vector<Vector> vecs(basis);

    int len = vecs.size();
    for(int i = 0; i < len; ++i)
        assert(vecs(i).dimension() == vecs(0).dimension());

    assert(coordinate_vec.dimension() == vecs(0).dimension());

    size_t basis_size = basis.size();
    size_t vec_size = vecs(0).dimension();

    Vector vec(vec_size);

    for(size_t i = 0; i < basis_size; ++i)
        for(size_t j = 0; j < vec_size; ++j)
            vec(i) += coordinate_vec(j) * vecs.at(j)(i);

    return vec;
}

std::vector<Vector> row_space_basis(Matrix mx)
{
    mx = mx.gaussJordanElimination();

    int rows_num = mx.rows();
    int cols_num = mx.cols();

    vector<Vector> space_basis;
    Vector lav(cols_num);

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            if(mx(i, j) != 0)
            {
                for(int j = 0; j < cols_num; ++j)
                    lav(j) = mx(i, j);

                space_basis.push_back(lav);

                break;
            }

    return space_basis;
}

vector<Vector> column_space_basis(Matrix mx)
{
    Matrix m = mx.gaussJordanElimination();

    int rows_num = mx.rows();
    int cols_num = mx.cols();

    vector<Vector> space_basis;

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
        {
            Vector temp(rows_num);

            if(m(i, j) != 0)
            {
                for(int k = 0; k < rows_num; ++k)
                    temp( k ) = mx(k, j);

                space_basis.push_back(temp);

                break;
            }
        }

    return space_basis;
}


vector<Vector> null_space(Matrix mx)
{
    Matrix m = mx.gaussJordanElimination();

    int rows_num = mx.rows();
    int cols_num = mx.cols();

    vector<int> pivot_cols;

    vector<Vector> free_variables(cols_num);

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            if(m(i, j) != 0)
            {
                // keeps all cols numbers so it is guaranteed that the column that contains a pivot won't
                // be used for the null space
                pivot_cols.push_back(j);

                break;
            }

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
        {
            if(m(i,j) != 0)
            {
                for(int k = 0; k < cols_num; ++k)
                {
                    // the j'th column is the one with pivot so it can not be used for the null space
                    // meaning that it has to be above or below

                    // if it is below it means that the k'th column might be one with free variable,
                    // it will be checked, if it is free it will be added zero because to get to the
                    // j'th column it had to get past only zeroes
                    if( k < j )
                    {
                        // starting from the second row, before immediately adding 0(zero), it will be checked
                        // whether the column is one that contains a pivot, in case it does the 0 won't be added
                        if(i > 0)
                        {
                            if(find(pivot_cols.cbegin(), pivot_cols.cend(), k) == pivot_cols.cend())
                                free_variables(j).push_back(0);
                        }
                        else
                            free_variables(j).push_back(0);
                    }
                    else if(k > j && find(pivot_cols.cbegin(), pivot_cols.cend(), k) == pivot_cols.cend())
                    {
                        free_variables(j).push_back( -m(i, k) );
                    }
                }
                break;
            }
        }

    int num_vectors = free_variables.size();
    int dimension;

    // get the dimension of the vector that will be of the null space
    for(int i = 0; i < num_vectors; ++i)
        if (free_variables(i).dimension() != 0)
        {
            dimension = free_variables(i).dimension();
            break;
        }


    // add the Identity Matrix to the rows in the new matrix which correspond to the 'free' columns
    // in the original matrix, making sure the number of rows equals the number of columns in the
    // original matrix (otherwise, we couldn't multiply the original matrix against our new matrix)
    int ind = 0;
    for(int i = 0; i < num_vectors; ++i)
    {
        if(free_variables(i).dimension() == 0)
        {
            for(int j = 0; j < dimension; ++j)
                if(j == ind)
                    free_variables(i).push_back(1);
                else
                    free_variables(i).push_back(0);

            ++ind;
        }
    }

    vector<Vector> space_basis(dimension, num_vectors);

    for(int i = 0; i < dimension; ++i)
        for(int j = 0; j < num_vectors; ++j)
            space_basis.at(i)( j ) = free_variables.at(j)(i);

    return space_basis;
}

std::size_t column_space_dim(Matrix mx)
{
    mx = mx.gaussJordanElimination();

    int rows_num = mx.rows();
    int cols_num = mx.cols();
    int dimension = 0;

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            if(mx(i, j) != 0)
            {
                ++dimension;

                break;
            }

    return dimension;
}

std::size_t row_space_dim(Matrix mx)
{
    mx = mx.gaussJordanElimination();

    int rows_num = mx.rows();
    int cols_num = mx.cols();
    int dimension = 0;

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            if(mx(i, j) != 0)
            {
                ++dimension;
                break;
            }

    return dimension;
}

std::size_t nullity(Matrix mx)
{
    Matrix m = mx.gaussJordanElimination();

    int rows_num = mx.rows();
    int cols_num = mx.cols();

    vector<int> pivot_cols;

    vector<Vector> free_variables(cols_num);

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            if(m(i, j) != 0)
            {
                pivot_cols.push_back(j);

                break;
            }

    int dimension = 0;

    for(int i = 0; i < rows_num; ++i)
        for(int j = 0; j < cols_num; ++j)
            if(m(i,j) != 0)
            {
                for(int k = 0; k < cols_num; ++k)
                {
                    if(k < j )
                    {
                        if(i > 0)
                        {
                            if(find(pivot_cols.cbegin(), pivot_cols.cend(), k) == pivot_cols.cend())
                                ++dimension;
                        }
                        else
                            ++dimension;
                    }
                    else if(k > j && find(pivot_cols.cbegin(), pivot_cols.cend(), k) == pivot_cols.cend())
                        ++dimension;
                }

                return dimension;
            }

    return 0;
}

}// L_Algebra namespace

main.cpp

#include <iostream>
#include <math.h>
//#include <boost/timer/timer.hpp>
#include "Matrix.h"
#include "LA_Vector.h"
#include <vector>
#include <boost/format.hpp>

using namespace L_Algebra;
using namespace std;

int main()
{

    vector<int> vec;
    vec.push_back(76);
    vec.push_back(76);
    vec.push_back(76);
    vec.push_back(76);
    vec.push_back(76);

    Vector vv(vec);

    int sd = 87, ds = 56;

    Fraction ffr = 10;


    Matrix b(3,4,3);
    Matrix c{5,5,3};

    Matrix a = {{-5, 5, -6, 1, 0}, {0, -5, 10, -3, 3}, {1, 11, 6, 1, 7}, {4, 5, -9, 9, -7}, {-5, 10, 0, -4, 4}};
    Matrix s = {{5, 5, -6, 1, 0}, {3, 4, 5, 7, 8}, {1, 11, 6, 1, 7}, {4, 5, -9, 9, -7}, {5, 10, 0, -4, 4}};
    Matrix s1 = {{5, 5, -6, 1, 0}, {3, 4, 5, 7, 8}, {1, 11, 6, 1, 7}, {4, 5, -9, 9, -7}, {5, 10, 0, -4, 4}};

    cout << a * 23;


    Matrix sw = {{-5}};

    Matrix d = {{1, 0, 2}, {2, 3, 7}};//, {-2, 2, 1, 7}, {-2, 3, 4, 1} };
    Matrix e = {{1, 1}, {0, 0} };
    Matrix g = {{0, 1}, {1, 0} };
    Matrix h = {{1, 0}, {0, 1} };
    Matrix i = {{1, 1}, {0, 1 } };

    // cout << turnMatricesIntoLinearCombination({e, g, h, i});

    try
    {
        //  cout << boost::format("%1% %3%") % 36 % 77 % 34;
    }
    catch (exception& e)
    {
        cout << e.what();
    }


    Matrix f = { {4, 0, 7, 6}, {1, 0, 7, 7}, {8, 0, 8, 8}};//, {-1, -4, -5, 0} };
    Matrix ff = { {4, 2, 7, 6, 5, 6}, {1, 7, 7, 7, 8, 0}, {8, 2, 8, 8, 9, 1}, {-1, -4, -5, 0, 1, 5} };

    Matrix mx1 = { {4, 1, 3, 1}, {3, 1, 3, 0}, {5, 1, 4, 1} };
    Matrix mx11 = { {1, 4, 8, 2}, {1, 4, 4, 9}, {1, 4, 4, 3}, {1, 4, 5, 5} };



   // cout << f << endl << endl;

   // vector<Vector> test = null_space(mx11);

    //cout << f.gaussJordanElimination();

//    for(auto e : test)
//        cout << e << endl;
//
//    cout << endl << nullity(f);


    b(0,2) = 4;
    b(1,2) = 5;
    b(1,3) = 2;
    b(2,0) = -8;
    b(2,3) = 9;
    b(0,0) = 1;
    b(0,1) = 2;



    //cout << mx11 << endl << endl;
    //vector<Vector> test3 = null_space(mx11);

//        for(auto e : test3)
//        cout << e << endl;


    //  cout << mx11.determinant();


    /*

     Vector lav1 = {1, 2, 1};
    Vector lav2 = {2, 9, 0};
    Vector lav = {3, 3, 4};

    Vector lav1 = {1, 5, 3};
    Vector lav2 = {-2, 6, 2};
    Vector lav = {3, -1, 1};

    Vector lav1 = {1, 2, -1};
    Vector lav2 = {6, 4, 2};
    Vector lav3 = {9, 2, 7};

    Vector lav1 = {3, 6, -9, -3};
    Vector lav2 = {6, -2, 5, 1};
    Vector lav3 = {-1, -2, 3, 1};
    Vector lav4 = {2, 3, 0, -2, 0};

    Vector lav3 = {3, 2, 1};
    */

    // cout << p.gaussJordanElimination();

    Matrix mx({ {3, 1, 1, 1}, {5, 2, 3, -2}});//,{-1, -2, 3, 1}});

    //  cout << mx.gaussJordanElimination();


    initializer_list<initializer_list<Fraction>> A = { {1, 3}, {1, -2} };
    initializer_list<Vector> B = { {3, 5}, {1, 2} };
    initializer_list<Vector> C = {{1, 0, 0, 0, }, {-2, 1, 0, 0, }, {5, 3, 0, 0}, {0, 0, 1, 0}, {3, 0, 0, 0} };
    //  Vector vec = {3, 2};

    Matrix gt(A);
    Matrix wz = { {0, 0, 0, 2, 9, 6}, {0, 0, 0, 4, 5, 8} };
    Matrix wzf = { {3, 2, 9, 2, 9, 6}, {6, 4, 5, 4, 5, 8} };
    Matrix z = { {1, 3, -2, 0, 2, 0}, {2, 6, -5, -2, 4, -3}, {0, 0, 5, 10, 0, 15}, {2, 6, 0, 8, 4, 18} };

//    cout << gt;

    Matrix dz = { {4, 1, 5, 1, 7, 8, 2}, {6, 3, 3, 5, 2, 3, 1}};//, {0, 0, 5, 10, 0, 15}, {2, 6, 0, 8, 4, 18} };

    Matrix fz = { {1, 3, 4, 4}, {2, 3, 5, 4}, {9, 1, 7, 2}};// {-1, -4, -5, 0} };
    Matrix tfz = { {1, 3, 4, 4, 1}, {2, 3, 5, 4, 5}, {9, 1, 7, 2, 3}};// {-1, -4, -5, 0} };

    Matrix khan = { {1, 1, 2, 3, 2}, {1, 1, 3, 1, 4} };
    Matrix kha = { {2, 0, 2}, {-1, 0, -1}, {-1, 0, -1} };

//    boost::timer::cpu_timer timer;
//    wz.gaussJordanElimination();
    //  timer.stop();


    //  cout << timer.format();

    Vector lav1 = {0, -2, 2};
    Vector lav2 = {1, 3, -1};
    Vector lav3 = {9, 0, 0};
    Vector lav4 = {4, 0, 2};
    Vector v = { 0, 0, 0};

    Matrix p = { {4, 0}, {-2, -2} };
    Matrix ph = { {1, -1}, {2, 3} };
    Matrix ph1 = { {0, 2}, {1, 4} };
    Matrix ph2 = { {-1, 5}, {7, 1} };
    Matrix ph21 = { {6, -8}, {-1, -8} };
    Matrix ph3 = { {6, 0}, {3, 8} };
    Matrix ph0 = { {0, 0}, {0, 0} };

    Fraction fr1(27, 17);
    Fraction fr2(43, 34);
    Fraction fr3(-29, 306);

    Matrix mcf(3, 3, {2, 3, 5, 6, 4, 5, 5, 8, 9});

    double db = 10.0 / 3;

    Fraction frt;



    // cout << frt;

    // cout << s << endl;


    try
    {
//        cout << s.main_diagonal() << endl;
//        cout << s.secondary_diagonal() << endl;

        //cout << coordinate_vector_relative_to_basis({ {0,1,0}, { {-4,5}, 0, {3,5}, }, { {3,5}, 0, {4,5} } }, {1,1,1});

        //cout << change_basis(vec, A, B);

        //cout << kha.gaussJordanElimination() << endl;

        //vector<Vector> v = null_space(kha);
        //  cout << coordinate_vector_relative_to_basis({ lav1, lav2,lav3}, lav4);

        // for(auto e : v)
        //    cout << e << endl;

        //  cout << endl << khan.rank();
    }
    catch(exception& e)
    {
        cout << e.what();
    }

//cout << lav2 * (lav ->* lav1);

}

What I am looking for is reviews on every possible aspect: C++ best practices (taking into account C++20), algorithms used, code simplicity/readability/organization, potential bugs, tips, tricks, warnings etc etc.

Worthy to note that I have tested every functionality as best as I could which I am quite sure it is not good enough.