Linear Algebra Library in C++;

Question

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 INVERTIBLE\n");

    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 << "\b\b \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.

Answer 1

Some things to consider:

1. Fundamental types do not have move constructors, so num(std::move(_num)) is just the equivalent of num(_num)

2. If you're not doing template code, move definitions out of header files. This can cause naming conflicts if multiple files include Fraction.h

3. Having a ++ and -- operator for a Fraction doesn't make sense. What does it mean to increment a fraction. It seems that you've chosen to do frac + 1 but that if I wanted (num+1)/den

4. You can write num = num / n; as num /= n; which behaves like += or -=

5. For - operator, you call your intermediate variable sub, but in + and * operator you call them addition and multiplication. Keep it consistent. Also in / operator, you call the result sub when I think you meant division.

6. Your Matrix only takes a std::initializer_list<>. What is somebody wants to pass a std::vector<>? It seems that they would be out of luck

7. Use a for each loop instead of iterators in your Matrix constructor:

   for (const auto& row: values)
   {
       assert(row.length() != 0);
   }
   

There's probably some other things, but that was what I was able to find

Answer 2

How long is a long long?

It depends on the CPU architecture and the operating system how long a long long really is. It might help to be more specific, and specify that a Fraction is a fraction of 64-bit integers, and then use int64_t. Also, instead of writing long long, consider creating a type alias:

using Integer = long long;

And use that everywhere. That makes changing the type of the integers used very easy.

Unnecessary calls to std::move

There is no need to use std::move when copying an integer into another integer, it just clutters the code. Just write:

Fraction(Integer _num = 0, Integer _den = 1) : num{_num}, den{_den}

Avoid names starting with an underscore

There are certain rules for using underscores in identifiers. While the above use is actually OK, I recommend you don't start any name with an underscore, as that is an easier rule to remember. You also don't need the underscores in the above function definition, you can write:

Fraction(Integer num = 0, Integer den = 1) : num{num}, den{den}

Handling zero denominators

Your code assert()s that the denominator is not zero. Be aware that in release builds, assert() macros might be disabled. If you want to ensure you always report an error if the denominator is zero, consider throwing a std::domain_error.

However, consider that the following is perfectly fine code when dealing with floating point numbers:

float foo = 1.0 / 0.0;

The value of foo is well-defined in this case: it is positive infinity. You might want to support the denominator being zero. Just be aware of this inside simplify(), and just don't do anything if den == 0.

Reduce the amount of overloads you need to write

You have a lot of code duplication that can be reduced. Take for example Fraction's operator+: you have three variants:

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

You only need to write one variant:

friend Fraction operator+(const Fraction& lhs, const Fraction& rhs);

Since a Fraction can be implicitly constructed from a single long long, the above statement will handle any combination of long long and Fraction arguments.

Cast to the right type

The function to_double() lies and returns a long double instead. Note that double is not the same as long double, on x86 and x86_64 a long double is 80 bits instead of 64 bits, and there are even architectures where long double is 128 bits.

The implementation of the function to_float() casts numerator to double. Why not cast it to float instead?

Remove unused code

There's a lot of unused code. Some of it is commented out, but for example in Fraction::operator-(const Fraction &), there are two variables nu and de that are not used at all (and if they would have been they would have the wrong type).

Avoid unnecessary parentheses:

return (Fraction(n) + fr);

Can be written as:

return Fraction(n) + fr;

Optimize pow_fract()

There are more optimal ways to implement integer power functions, see this StackOverflow question.

Consider allowing zero size vectors and matrices

The constructors of Matrix and Vector all assert() that the constructor object has a non-zero size. But is it really necessary to limit that case? Most functions work perfectly well with zero-sized vectors and matrices, and you avoid the overhead of the check every time you construct an object. You only need this check in the rare cases where a function would cause a crash or undefined behaviour if the size is zero.

Use range-for and STL algorithms where applicable

I see a lot of old-style for-loops where you could have used a range-for or even an STL algorithm. For example, Matrix::Matrix() can be rewritten as:

Matrix::Matrix(std::initializer_list<std::initializer_list<Fraction>> values )
{
    rows_num = values.size();
    cols_num = 0;

    for(auto &row: values) {
        cols_num = row.size();
        break;
    }

    data.reserve(rows_num * cols_num);

    for(auto &row: values)
    {
        assert(row.size() == cols_num);
        std::copy(row.begin(), row.end(), std::back_inserter(data));
    }
}

As another example, Matrix::operator+(const Matrix &) can be written as:

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

    Matrix result(rows_num, cols_num);

    std::transform(data.begin(), data.end(), mx.data.begin(), result.data.begin(), std::plus);

    return result;
}

Note that the result matrix is initialized unnecessarily; consider adding a (possibly private) constructor that allows creating a Matrix of a given size with data not being initialized.

Use a std::vector in Vector

Why does the Vector class store its data in a std::deque? You don't need the functionality of a deque (like \$\mathcal{O}(1)\$ insertion and removal at both ends), but now you pay the price in performance and storage overhead.

Reduce the number of constructors of Vector

You have overloaded the constructor of Vector to handle std::vectors and std::deques of ints and Fractions as input. But what if I want to pass it a std::array<unsigned int>? Surely you can see that you cannot support everything this way unless you write hundreds of overloads, and even then you will miss some cases. If you really want to handle arbitrary contains being passed to the constructor, do what the STL does in its container classes, and write a template that takes a pair of iterators, like so:

template<class InputIt>
Vector(InputIt first, InputIt last): data(first, last) {}

That's all there is to it. Now you can do something like:

std::list<unsigned long> foo{1, 2, 3, 4, 5};
Vector vec(foo.begin(), foo.end());

Answer 3

I think I would change the Fraction::simplify, personally. I would maybe let it be a static class method that took a Fraction object and returned it's reduced representation without modifying the original object.

There are times when it might be useful to calculate a proportion without clobbering the original object, such as when calculating the Binomial proportion confidence interval, for example.