#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <iomanip>
#include <cmath>
#include <cassert>

using namespace std;

// Gaussian Elimination: Transform a Matrix to Row Echelon Form
//     using elementary transformations of matrix rows.
//     Return the rank of matrix (i.e. the number of nonzero rows
//     in Row Echelon Form)
// Parameters:
//     double *a  -- array of matrix elements in row major order
//     int rows   -- number of matrix rows
//     int cols   -- number of matrix columns
//     double eps -- precision of algorithm (default eps = 1e-8)
int gauss(double *a, int rows, int cols, double eps = 1e-8);

// Reverse pass after the gaussian elimination.
//     Transform the Row Echelon Form of matrix to the
//     Reduced Row Echelon Form
//     using elementary transformations of matrix rows.
//     Example:
//         Initial Matrix    Row Echelon Form
//          1  2  3  4       13.0000 14.0000 15.0000 16.0000
//          5  6  7  8  -->   0.0000 -0.9231 -1.8462 -2.7692
//          9 10 11 12        0.0000  0.0000  0.0000  0.0000
//         13 14 15 16        0.0000  0.0000  0.0000  0.0000
//
//                            Reduced Row Echelon Form
//                            1.0000 0.0000 -1.0000 -2.0000
//                      -->   0.0000 1.0000  2.0000  3.0000
//                            0.0000 0.0000  0.0000  0.0000
//                            0.0000 0.0000  0.0000  0.0000
void reversePass(
    double *a, int rows, int cols, int rank, 
    double eps = 1e-8
);
bool readMatrix(const char *path, double* &a, int& rows, int& cols);
bool writeMatrix(const char *path, const double *a, int rows, int cols);
bool printMatrix(ostream& f, const double *a, int rows, int cols);

int main() {
    double *a;
    int rows, cols;
    if (!readMatrix("input.txt", a, rows, cols)) {
        cout << "Cannot read a matrix" << endl;
        return (-1);
    }
    cout << "Initial matrix:" << endl;
    printMatrix(cout, a, rows, cols);
    
    double *a0 = new double[rows*cols];
    for (int i = 0; i < rows*cols; ++i)
        a0[i] = a[i];
    
    int rank = gauss(a, rows, cols);
    cout << "Rank of matrix = " << rank << endl;
    if (rows == cols) {
        double d = 1.;
        for (int i = 0; i < rows; ++i) {
            d *= a[i*cols + i];
        }
        cout <<"Determinant = " << d << endl;
    }

    cout << "Row Echelon Form of matrix:" << endl;
    printMatrix(cout, a, rows, cols);
    writeMatrix("output.txt", a, rows, cols);    
    
    reversePass(a, rows, cols, rank);
    cout << "After Reverse Pass (Reduced REF):" << endl;
    printMatrix(cout, a, rows, cols);
    writeMatrix("rref.txt", a, rows, cols);
    
    return 0;
}

int gauss(double *a, int rows, int cols, double eps) {
    int i = 0;
    int j = 0;
    while (i < rows && j < cols) {
        double maxElem = fabs(a[i*cols + j]);
        int imax = i;
        for (int k = i+1; k < rows; ++k) {
            // Search for the maximal element in j-th column
            if (fabs(a[k*cols + j]) > maxElem) {
                maxElem = fabs(a[k*cols + j]);
                imax = k;
            } 
        }    
        if (maxElem <= eps) {
            ++j;
            continue;
        }    
        if (i != imax) {
            // Swap i-th and imax-th rows
            for (int l = j; l < cols; ++l) {
                double tmp = a[i*cols + l];
                a[i*cols + l] = a[imax*cols + l];
                a[imax*cols + l] = (-tmp);
            }
        }        
        double r = a[i*cols + j];
        for (int k = i+1; k < rows; ++k) {
            double coeff = a[k*cols + j]/r;
            a[k*cols + j] = 0.;
            for (int l = j+1; l < cols; ++l) {
                a[k*cols + l] -= coeff*a[i*cols + l];
            }
        }
        ++i; ++j;    
    }
    return i;    
}

bool readMatrix(const char *path, double* &a, int& rows, int& cols) {
    ifstream f(path);
    if (!f.is_open())
        return false;

    if (!(f >> rows >> cols))
        return false;
    
    a = new double[rows*cols];
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            if (!(f >> a[i*cols + j]))
                return false;
        }    
    }
    return true;
}

bool writeMatrix(const char *path, const double *a, int rows, int cols) {
    ofstream f(path);
    if (!f.is_open())
        return false;
    f << rows << ' ' << cols << endl;
    bool res = printMatrix(f, a, rows, cols);
    f.close();
    return res;
}

bool printMatrix(ostream& f, const double *a, int rows, int cols) {
    // Define the maximal number of positions 
    // needed for printing a matrix element
    int lmax = 0;
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            ostringstream str;
            str << setprecision(4) << fixed;
            str << a[i*cols + j];
            int l = int(str.str().length());
            if (l > lmax)
                lmax = l;
        }
    }
    
    // Print the matrix
    f << setprecision(4) << fixed;
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            ostringstream str;
            str << setprecision(4) << fixed;
            str << a[i*cols + j];
            int l = int(str.str().length());
            
            if (j > 0)
                f << ' ';   // Delimiter
            
            // Padding at the left side
            for (int k = 0; k < lmax - l; ++k)
                f << ' ';
            f << str.str(); // Print the matrix element 
        }
        f << endl;
    }
    return true;
}

void reversePass(
    double *a, int rows, int cols, int rank, 
    double eps
) {
    assert(rank <= rows);    
    for (int i = rank-1; i >= 0; --i) {
        // Search for the first nonzero element in i-th row
        int j = 0;
        while (j < cols && fabs(a[i*cols + j]) <= eps)
            ++j;
        assert(j < cols && fabs(a[i*cols + j]) > eps);
        int j0 = j; // The colunm of the first non-zero element in i-th row
        double r = a[i*cols + j]; // The first non-zero element in i-th row
        
        // Divide i-th row by r
        a[i*cols + j] = 1.; ++j;
        while (j < cols) {
            a[i*cols + j] /= r;
            ++j;
        } 
        
        // Annihilate j-th column above i-th row
        for (int k = 0; k < i; ++k) {
            double coeff = a[k*cols + j0];
            a[k*cols + j0] = 0.;
            for (int s = j0+1; s < cols; ++s)
                a[k*cols + s] -= a[i*cols + s]*coeff;
        }    
    }
}