// Reading and writing a matrix of dimension m*n
// m -- number of rows
// n -- number of columns
#include <stdio.h>
#include <math.h>

bool readMatrix(FILE* f, double **a, int* m, int* n); 
bool writeMatrix(FILE* f, const double *a, int m, int n);
int gaussMethod(double *a, int m, int n); // Return a matrix rank
double det(double *a, int m, int n);
bool solveLinearSystem(double *a, int n, double* x); // matrix n*(n+1)

const double EPS = 1e-8;

int main() {
    int m, n;
    double *a = 0;
    
    FILE* f = fopen("input.txt", "rt");
    if (f == NULL) {
        perror("Cannot open an input file");
        return (-1);
    }
    
    if (!readMatrix(f, &a, &m, &n)) {
        printf("Incorrect format of matrix.\n");
        delete[] a;
        fclose(f);
        return (-1);
    }
    fclose(f);
    
    printf("Matrix of size %d*%d:\n", m, n);
    writeMatrix(stdout, a, m, n);
    
    int rank = gaussMethod(a, m, n);
    
    printf("Row echelon form of matrix:\n");
    writeMatrix(stdout, a, m, n);

    printf("Rank of matrix: %d\n", rank);
    if (m == n) {
        printf("Determinant of matrix: %f\n", det(a, m, n));
    }
    
    // Writing the resulting matrix to the file "output.txt"
    FILE* g = fopen("output.txt", "wt");
    if (g == NULL) {
        perror("Cannot open an output file");
    } else {
        if (!writeMatrix(g, a, m, n))
            perror("Could not write a matrix to the output file");
        fclose(g);
    }

    delete[] a;
    
    return 0;
}

bool readMatrix(FILE* f, double **a, int* m, int* n) {
    // 2 3     
    // 11 12 13
    // 21 22 23
    int rows, cols;
    if (fscanf(f, "%d%d", &rows, &cols) < 2) {
        return false;
    }
    *m = rows; *n = cols;
    double *matrix = new double[rows*cols];
    *a = matrix;
    for (int i = 0; i < rows; ++i) {
        // Read i-th row of matrix
        for (int j = 0; j < cols; ++j) {
            if (fscanf(f, "%lf", matrix + i*cols + j) < 1) {
                // Read error
                return false;
            }
        }
    }
    return true;
}
    
bool writeMatrix(FILE* f, const double *a, int m, int n) {
    if (fprintf(f, "%d %d\n", m, n) <= 0) {
        return false;
    } 
    
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            if (fprintf(f, "%12.4f", a[i*n + j]) <= 0)
                return false;
        }
        fprintf(f, "\n");
    }    
    return true;
}

int gaussMethod(double *a, int m, int n) { // Return a matrix rank 
    int i = 0;
    int j = 0;
    while (i < m && j < n) {
        // 1. Find the maximal element in j-th column
        int k = i;
        double maxElem = fabs(a[k*n + j]);
        for (int l = i+1; l < m; ++l) {
            if (fabs(a[l*n + j]) > maxElem) {
                maxElem = fabs(a[l*n + j]);
                k = l;
            }
        }
        
        if (maxElem <= EPS) {
            // Zero colunm
            for (int l = i; l < m; ++l) {
                a[l*n + j] = 0.;
            }
            ++j;
            continue;
        }
        
        if (k != i) {
            // Swap rows i, k
            for (int l = j; l < n; ++l) {
                double tmp = a[i*n + l];
                a[i*n + l] = a[k*n + l];
                a[k*n + l] = (-tmp);
            }
        }            
            
        for (int l = i+1; l < m; ++l) {
            // Annihilate the element a_lj
            double r = a[l*n + j]/a[i*n + j];
            a[l*n + j] = 0.;    
            for (int s = j+1; s < n; ++s) {
                a[l*n + s] -= r*a[i*n + s];
            }
        }
        ++i; ++j;
    }
    return i;
}
            
double det(double *a, int m, int n) {
    double d = 1.;
    int k = n;
    if (m < n)
        k = m;
    for (int i = 0; i < k; ++i) {
        d *= a[i*n + i];
    }
    return d;
}