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

using namespace std;

bool readMatrix(istream& s, int& m, int& n, double* &a);
bool readMatrix(string path, int& m, int& n, double* &a);

void writeMatrix(ostream& s, int m, int n, const double* a);
void writeMatrix(string path, int m, int n, const double* a);

// Transform a matrix a with m rows and n columns
// to the row echelon form. Return the rank of matrix.
int gauss(
    int m, int n, double *a,
    double eps = 1e-12
); 

// Reduced Row Echelon Form
int rref(
    int m, int n, double *a, 
    double eps = 1e-12
);

bool inverseMatrix(
    int n, const double* a,
    double* b, double eps = 1e-12
);    

int main() {
    int m, n;
    double *a;
    
    cout << fixed << setprecision(4);
    if (!readMatrix("input.txt", m, n, a)) {
        cerr << "Could not read the matrix" << endl;
        return -1;
    }
    cout << "Initial matrix:" << endl;
    writeMatrix(cout, m, n, a);
    
    double* a0 = new double[m*n];
    for (int i = 0; i < m*n; ++i) {
        a0[i] = a[i];
    }    
    
    cout << "Row echelon form of matrix:" << endl;
    int r = gauss(m, n, a);
    writeMatrix(cout, m, n, a);
    
    // Determinant
    if (m == n) {
        cout << "Determinant = ";
        if (r < m) {
            cout << 0 << endl;
        } else {
            double d = 1.;
            for (int i = 0; i < m; ++i) {
                d *= a[i*n + i];
            }
            cout << d << endl;
        }
    }   

    // Save the initial matrix
    for (int i = 0; i < m*n; ++i) {
        a[i] = a0[i];
    }    

    cout << "Reduced row echelon form:" << endl;
    r = rref(m, n, a0);
    writeMatrix(cout, m, n, a0);

    if (m == n && r == m) {
        double* b = new double[m*n];
        if (inverseMatrix(m, a, b)) {
            cout << "Inverse matrix:" << endl;
            writeMatrix(cout, m, m, b);
        }
        delete[] b;
    }
    delete[] a;
    delete[] a0;

    return 0;
}

bool readMatrix(istream& s, int& m, int& n, double* &a) {
    s >> m >> n;
    if (!s || m <= 0 || n <= 0)
        return false;
    a = new double[m*n];
    for (int i = 0; i < m*n; ++i) {
        if (!(s >> a[i])) {
            delete[] a;
            return false;
        }
    }
    return true;
}

bool readMatrix(string path, int& m, int& n, double* &a) {
    ifstream s(path);
    if (!s.is_open())
        return false;
    return readMatrix(s, m, n, a);
}

void writeMatrix(ostream& s, int m, int n, const double* a) {
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            if (j > 0)
                s << " ";
            s << a[i*n + j];
        }
        s << endl;
    }
}

int gauss(int m, int n, double *a, double eps) {
    int i = 0;
    int j = 0;
    while (i < m && j < n) {
        // 1. Find nonzero element in j-th columns
        //    starting with i-th row
        int imax = i;
        double amax = fabs(a[i*n + j]);
        for (int k = i+1; k < m; ++k) {
            if (fabs(a[k*n + j]) > amax) {
                amax = fabs(a[k*n + j]);
                imax = k;
            }
        }
        if (amax <= eps) {
            ++j;
            continue;
        }
        
        assert(fabs(a[imax*n + j]) > eps);
        if (imax != i) {
            // Swap rows i and imax
            for (int l = j; l < n; ++l) {
                double tmp = a[i*n + l];
                a[i*n + l] = a[imax*n + l];
                a[imax*n + l] = (-tmp);
            }
        }        
            
        assert(fabs(a[i*n + j]) > eps);
        double r = a[i*n + j];
        
        // Annihilate j-th column starting with i-th row
        for (int k = i+1; k < m; ++k) {
            double c = a[k*n + j]/r;
            for (int l = j; l < n; ++l) {
                a[k*n + l] -= a[i*n + l]*c;
            }
        }
        ++j; ++i;
    }
    return i;
}

// Reduced Row Echelon Form
int rref(int m, int n, double *a, double eps) {        
    int i = 0;
    int j = 0;
    while (i < m && j <= n) {
        // Find maximal element in j-th
        // column starting from i-th row
        double amax = fabs(a[i*n + j]);
        int imax = i;
        for (int k = i+1; k < m; ++k) {
            double acurr = fabs(a[k*n + j]);
            if (acurr > amax) {
                amax = acurr;
                imax = k;
            }
        }
        
        if (amax <= eps) {
            ++j; continue;
        }
        
        if (imax != i) {
            // Swap line i and imax
            for (int l = j; l < n; ++l) {
                double tmp = a[i*n + l];
                a[i*n + l] = a[imax*n + l];
                a[imax*n + l] = tmp;
            }
        }
        
        double r = a[i*n + j];
        for (int l = j; l < n; ++l) {
            a[i*n + l] /= r;
        }
        
        // Annihilate j-th column, except
        // i-th element in it
        for (int k = 0; k < m; ++k) {
            if (k == i)
                continue;
            r = a[k*n + j];
            for (int l = j; l < n; ++l) {
                a[k*n + l] -=
                    a[i*n + l]*r;
            }
        }
        ++i; ++j;
    }
    return i;
}            

bool inverseMatrix(
    int m, const double* a,
    double* b, double eps
) {
    int n = 2*m;
    double* aExt = new double[m*n];
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < m; ++j) {
            aExt[i*n + j] = a[i*m + j];
        }    
        for (int j = m; j < n; ++j) {
            aExt[i*n + j] = 0.;
        }    
        aExt[i*n + (m + i)] = 1.;
    }
    
    cout << "Extended matrix:" << endl;
    writeMatrix(cout, m, n, aExt);
    
    int r = rref(m, n, aExt, eps);
    if (r < m) {
        cerr << "Singular matrix, no inverse"
             << endl;
        return false;
    }    
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < m; ++j) {
            b[i*m + j] =
                aExt[i*n + (m + j)];
        }
    }        
    return true;
}