#include <iostream>
#include <iomanip>
using namespace std;

#define SIZE 5
#define fabs(x) (((x) > 0) ? (x) : -(x))
#define fnonzero(x) (fabs(x) > 0.0001f)

void SwapRows(double rc[SIZE][SIZE + 1], int r1, int r2)
{
    for (int i = 0; i < SIZE + 1; i++)
    {
        float t = rc[r1][i];
        rc[r1][i] = rc[r2][i];
        rc[r2][i] = t;
    }
}

void Print(double rc[SIZE][SIZE])
{
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE; j++)
        {
            if (fnonzero(rc[i][j]))
            {
                cout << rc[i][j] << '\t';
            }
            else
            {
                cout << 0 << '\t';
            }
        }
        cout << endl;
    }
    cout << endl;
}

void Print(double rc[SIZE][SIZE + 1])
{
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = 0; j < SIZE + 1; j++)
        {
            if (fnonzero(rc[i][j]))
            {
                cout << rc[i][j] << '\t';
            }
            else
            {
                cout << 0 << '\t';
            }
        }
        cout << endl;
    }
    cout << endl;
}

void AddRow(double rc[SIZE][SIZE + 1], int SrcR, double Scalar, int DestR)
{
    //cout << 'R' << (DestR + 1) << " += " << Scalar << 'R' << (SrcR + 1) << endl;
    for (int i = 0; i < SIZE + 1; i++)
    {
        rc[DestR][i] += Scalar * rc[SrcR][i];
    }
}

void AddRow(double rc[SIZE][SIZE], int SrcR, double Scalar, int DestR)
{
    for (int i = 0; i < SIZE; i++)
    {
        rc[DestR][i] += Scalar * rc[SrcR][i];
    }
}

void ScaleRow(double rc[SIZE][SIZE + 1], int SrcR, double Scalar)
{
    if (Scalar != 1)
    {
        //cout << 'R' << (SrcR + 1) << " *= " << Scalar << endl;
    }
    for (int i = 0; i < SIZE + 1; i++)
    {
        rc[SrcR][i] = rc[SrcR][i] * Scalar;
    }
}

void ScaleRow(double rc[SIZE][SIZE], int SrcR, double Scalar)
{
    for (int i = 0; i < SIZE; i++)
    {
        rc[SrcR][i] = rc[SrcR][i] * Scalar;
    }
}

int main()
{
    double rc[SIZE][SIZE + 1];
    
    double solutions[SIZE];
    cout << "Building matrix to find solutions" << endl;
    for (int i = 0; i < SIZE; i++)
    {
        solutions[i] = (rand() % 20) - 10;
        cout << solutions[i] << ' ';
    }
    cout << endl;
    for (int i = 0; i < SIZE; i++)
    {
        double sum = 0;
        for (int j = 0; j < SIZE; j++)
        {
            rc[i][j] = (rand() % 20) - 10;
            sum += rc[i][j] * solutions[j];
        }
        rc[i][SIZE] = sum;
    }
    
    Print(rc);
    
    //AddRow(rc, 0, -(rc[1][0] / rc[0][0]), 1);
    //AddRow(rc, 0, -(rc[2][0] / rc[0][0]), 2);
    //AddRow(rc, 1, -(rc[2][1] / rc[1][1]), 2);
    
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = i + 1; j < SIZE; j++)
        {
            AddRow(rc, i, -(rc[j][i] / rc[i][i]), j);
        }
    }
    
    //ScaleRow(rc, 0, 1.0f / rc[0][0]);
    //ScaleRow(rc, 1, 1.0f / rc[1][1]);
    //ScaleRow(rc, 2, 1.0f / rc[2][2]);
    
    for (int i = 0; i < SIZE; i++)
    {
        if (fnonzero(rc[i][i]))
        {
            ScaleRow(rc, i, 1.0f / rc[i][i]);
        }
    }
    
    //AddRow(rc, 1, -rc[0][1], 0);
    //AddRow(rc, 2, -rc[0][2], 0);
    //AddRow(rc, 2, -rc[1][2], 1);
    for (int i = 0; i < SIZE; i++)
    {
        for (int j = i + 1; j < SIZE; j++)
        {
            AddRow(rc, j, -rc[i][j], i);
        }
    }
    
    Print(rc);
    return 0;
}