Tag: algorithm

permutations – What’s wrong with the following shuffle algorithm?

Suppose I have an array A of size n, with the initial state of:

A(0) == 1, A(1) = 2, ... A(n-1) = n

I know that one way to get a uniformly-distributed permutation of that array is to use Fisher-Yates Algorithm, but I am more interested in what’s flawed with the following naive approach:

for i = 0 to n-1:
   rand_i = random(0, n-1, UNIFORM_DIST)
   swap(A(i), A(rand_i))

When I tried to use that with n = 3, it seems like some permutations are more likely to be drawn than others.
Here’s the output of running the above algorithm 10,000,000 times, and then averaging on the results (shown in percentage of likelihood):

#> swap.out --iteartions 10000000
00: (1,2,3): 14.8237
01: (1,3,2): 18.5105
02: (2,1,3): 18.504
03: (2,3,1): 18.5216
04: (3,1,2): 14.7975
05: (3,2,1): 14.8427

To rule out the possibility that the Pseudo-Random-Generator is not uniform, I compared that with an algorithm that builds db of all permutations and then, using the same PRG – chooses one permutation:

#> choose_permutation.out --iterations 10000000
00: (1,2,3): 16.6847
01: (1,3,2): 16.6649
02: (2,1,3): 16.6731
03: (2,3,1): 16.6706
04: (3,1,2): 16.6516
05: (3,2,1): 16.655

This pattern is consistent; permutations 1,2,3 always come up with a higher likelihood than 0,4,5.
With n = 4:

#> swap.out --iteartions 10000000
00: (1,2,3,4): 3.90774
01: (1,2,4,3): 3.90958
02: (1,3,2,4): 3.91321
03: (1,3,4,2): 5.46167
04: (1,4,2,3): 4.29965
05: (1,4,3,2): 3.51932
06: (2,1,3,4): 3.89793
07: (2,1,4,3): 5.83975
08: (2,3,1,4): 5.45905
09: (2,3,4,1): 5.47537
10: (2,4,1,3): 4.30715
11: (2,4,3,1): 4.3
12: (3,1,2,4): 4.29691
13: (3,1,4,2): 4.3075
14: (3,2,1,4): 3.51411
15: (3,2,4,1): 4.29534
16: (3,4,1,2): 4.30452
17: (3,4,2,1): 3.90896
18: (4,1,2,3): 3.12621
19: (4,1,3,2): 3.50413
20: (4,2,1,3): 3.52526
21: (4,2,3,1): 3.12843
22: (4,3,1,2): 3.89564
23: (4,3,2,1): 3.90257

And comparing to randomly choosing from db of permutations:

#> choose_permutation.out --iterations 10000000  
00: (1,2,3,4): 4.16284
01: (1,2,4,3): 4.16721
02: (1,3,2,4): 4.15893
03: (1,3,4,2): 4.17306
04: (1,4,2,3): 4.15853
05: (1,4,3,2): 4.16169
06: (2,1,3,4): 4.16584
07: (2,1,4,3): 4.17245
08: (2,3,1,4): 4.17309
09: (2,3,4,1): 4.15519
10: (2,4,1,3): 4.17007
11: (2,4,3,1): 4.17163
12: (3,1,2,4): 4.16276
13: (3,1,4,2): 4.17367
14: (3,2,1,4): 4.17147
15: (3,2,4,1): 4.16955
16: (3,4,1,2): 4.16576
17: (3,4,2,1): 4.1659
18: (4,1,2,3): 4.16653
19: (4,1,3,2): 4.16749
20: (4,2,1,3): 4.17225
21: (4,2,3,1): 4.1668
22: (4,3,1,2): 4.16797
23: (4,3,2,1): 4.15932

What’s flawed with that method and why are those specific permutations come up more often?

Complexity in time and memory for graph search algorithm

I am working on an assignment where I have to write an algorithm to detect all vertices that lie in a cycle in a graph and then calculate its complexity. I have come up with an algorithm in pseudocode. My approach is to ‘delete’ all vertices that have only one connecting edge until no vertexes are left that have only one connecting edge. The remaining vertices, if any, are part of a cycle.

The graph is undirected.

/*
// vertices is an array containing every vertex
// e.g. vertices = (0,1,2,3,...)
*/
vertices()

/*
// edges is a two dimensional array defining each edge
// by the two vertices it connects
// e.g. edges = ((0,1),(2,3),(4,1),...)
*/
edges()

verticesInCycle = cycleDetection(vertices, edges, vertices)

verticesNotInCycle = vertices

foreach (verticesInCycle as vertexIC) {
    if (verticesNotInCycle.contains(vertexIC)) {
        verticesNotInCycle.remove(vertexIC);
    }
}


/*
// vertices() = vertices to test
// edges() = edges to test them against
// verticesLeft() = all untested vertices that are left
*/
function cycleDetection(vertices(), edges()), verticesLeft()) {
    recur = false
    nextVertices()

    foreach (vertices as vertex) {
        count = 0
        connectingEdge = null

        foreach (edges as edge) {
            if (vertex == edge(0) || vertex == edge(1)) {
                count++

                if (count > 1) {
                    break
                } else {
                    connectingEdge = edge
                }
            }
        }

        if (count == 1) {
            edges.remove(connectingEdge)
            verticesLeft.remove(vertex)

            if (connectingEdge(0) == vertex) {
                nextVertices.add(connectingEdge(1))
            } else {
                nextVertices.add(connectingEdge(0))
            }

            recur = true
        }
    }

    if (recur) {
        cycleDetection(nextVertices, edges, verticesLeft)
    } else {
        return verticesLeft
    }
}

Example Run

Example graph

10  6
  / 
  1   5
  |   |
  2   4-7
  | /
  8 3
   
    9

Initial call of cycleDetection

x   o
  / 
  o   o
  |   |
  o   o-x
  | /
  o o
   
    x

First recursive call

    o
   / 
  o   o
  |   |
  o   o
  | /
  x o

Second recursive call

    o
   / 
  o   o
  |   |
  o   o
    /
    o

Done!

Now I am calculating the complexity (time and memory)

I calculated that the complexity is approximately O(1.4(V*E)), where N is the number of vertices and E the number of edges.

This is based on the time (or iterations) because the maximum number of iterations (vertices and edges to check) is approximately 1.4(n*e), where n is the number of vertices and e the number of edges.

The worst case, which requires the most iterations, is a line of connected vertices with a cycle at one end.

o
|
| o-o-o-o-o-o-o-o-o-o
|/
o

In this case the complexity is the number of vertices n times the number of edges e multiplied by approximately 1.4. Initially all edges need to be iterated for every vertex but subsequently only the remaining edges need to be iterated for one vertex at a time. In the end about a fourth of n*e needs to be iterated through after the initial iteration of all edges for every vertex.

Memory usage

Memory usage is directly dependant on the amount of recursive calls, with the first call needing the most memory. Every call of cycleDetection uses v + e + vL memory, v + e + vL being the size of the three arrays passed as parameters.

The worst case for memory usage is a line of connected vertices with a cycle at one end (as above). In this case the amount of recursive calls required is v - 2, where v is the number of vertices given that v > 2. So memory complexity is approximately O(3) work memory.

A. I am unsure if my time complexity is formulated right? It doesn’t need to be spot on but I also saw that you can leave out the multiplier and write something like O(VxE**2).

B. I am pretty sure my memory complexity is wrong? I have read this and formulated my complexity accordingly based on the fact that I have three variables that I am passing to my function with every recursive call.

gr.group theory – Algorithm for Root System of Coxeter Group Generated by Permutations

Suppose we are given a group $G$ in terms of generators $t_1, …, t_n$ which are order 2 in $S_m$ (however we don’t assume anything other than that these elements generate $G$ and have order 2). Since $G$ is finite and generated by transpositions, it must have a root system. What is the best known algorithm for finding the root system?

computational geometry – An algorithm to split an area into multiple polygons based on other polygons intersection

I have a list of n polygons (A,B,C,D,E,...) which possibly intersect each other.

I need to find a new list of polygons (or multi-polygons) which looks like this:

(
Area only in A (A difference from the rest)
Area only in B
Area only in C
...

Area only in A and B (A intersection with B and difference from the rest)
Area only in A and C
Area only in B and C
...

Area only in A,B and C (A intersection with B and C and difference from the rest)
Area only in A,B and D 
Area only in A,B and E 
....

)```

The only way I can find is to iterate all the combinations and run a simple intersection/difference algorithms but the number rapidly increases.

algorithm – How to calculate best move in Match-3 game?

I am having a hard time identifying how to approach this problem.

If you have a match 3 game and the objects can only be swapped horizontally and vertically, and only swapped if it will result in a match of 3 or more jewels, which will be removed, with each removed jewel adding +1 to the score.

How would I find the best possible move? I’m confused at where to start and been searching for a while now and getting nowhere.

  • Board 6×6 or 8×8
  • Matches = 3 or more of same kind
  • write function to find best swap for maximum score.

Any help would be greatly appreciated !

EDIT: So I could iterate through the entire board with a nested for loop for the rows and columns, but im not sure what approach to take from here.

For example, say I start with top left which is “red” and I want to test if swapping this will result in a match. How would I go about searching and identifying all possibilities that would result in swapping this left equaling a match and score of 5.

**If I am being really unclear I apologies and will remove my question and possible reword.

Colour block example

Need Advise Machine Learning Algorithm for product Recommendation

enter image description here

I wanna give recommendation to buyer which product is should be bought by price, sold, shipping price, and location. If my location is C (sorry, it should be written “C” not “B”) which product is optimal for me from price, sold, shipping price, and location.

Since I’m new in machine learning, which one of machine learning algorithm is good to implement?

Thank you.

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.

time complexity – Modification to the Algorithm merge sort

I was trying to solve this problem but I’m struggling in the time complexity
Can u help me please:
Consider the following modification to the Algorithm mergesort. We apply
the algorithm on the input array A(1..n) and continue the recursive calls
until the size of a substance becomes relatively small, say m or less. At
this point, we switch to Algorithm insertion sort and apply it on the
small instance. So, the first test of the modified algorithm will look like
the following:
if high − low + 1 ≤ m then insertion sort(A(low..high)).
What is the largest value of m in terms of n such that the running time
of the modified algorithm will still be Θ(n log n)? You may assume for
simplicity that n is a power of 2

Optimization algorithm seems biased – Software Engineering Stack Exchange

I’ve written an algorithm designed to solve a minimization problem which can be demonstrated with the following simple example:

I have three observers (who’s locations are known) and one or more signal source(s) (who’s location(s) are unknown). Every time an observer “sees” a signal, the time at which the signal was seen is recorded. What I want to do is take these three sets of data (one from each observer), and produce a list of triplets, where each triplet has one timestamp contributed by each observer, such that the the difference between the maximum and minimum value is minimized, and each triplet is unique (meaning no two triplets share a timestamp).

Given three observers i,j and k, my algorithm accomplishes this by selecting one observer (this selection is arbitrary, let’s choose i), and for each timestamp recorded by that observer I find the closest timestamp from each of the other two observers and produce a triplet {t_i,t_j,t_k}. I record this triplet.

After iterating over all timestamps recorded by i, and producing a triplet for each, I sort this list of triplets by max(t_i,t_j,t_k)−min(t_i,t_j,t_k). I then iterate over this list, and with each iteration store t_i,t_j,t_k in a cache, and if t_i,t_j,t_k aren’t already in the cache, I add them to a second list. This second list ends up being a list of unique, minimized triplets.

The problem is, my algorithm seems to have some bias. When I reconstruct the location of each triplet, there is a clustering at the coordinates of the chosen observer i. This doesn’t necessarily indicate bias, but it is highly unlikely that such a clustering would occur. So, for now, I’m assuming (with a fair amount of confidence) that there’s an algorithm issue. (Note that I’m certain, beyond doubt, that there is no issue with the reconstruction algorithm.)

Is this a valid solution to this problem? Have I made a mistake? I feel like I must’ve done something really stupid here and haven’t noticed it. Is there a better way to solve this?