#include <cmath>
#include "Matrix.h"

// Elementary transfors of matrix rows
// Swap two matrix rows and multiply one of them
// by -1 so that determinant be invariant.
void Matrix::swapRows(int i, int k) {
    for (int j = 0; j < n; ++j) {
        double tmp = (*this)[i][j];
        (*this)[i][j] = (*this)[k][j];
        (*this)[k][j] = (-tmp);
    }
}
        
// Add k-th row of matrix multiplied by the coeff. c
// to i-th row
void Matrix::addRows(int i, int k, double c) {
    for (int j = 0; j < n; ++j) {
        (*this)[i][j] += (*this)[k][j]*c;
    }
}

// Multiply the i-th row of matrix by the coeff. c
void Matrix::multiplyRow(int i, double c) {
    for (int j = 0; j < n; ++j) {
        (*this)[i][j] *= c;
    }
}

// Transform the matrix to the row echelon form
// using elementary transforms of matrix rows.
// Return the rank of matrix.
int Matrix::gaussElimination(double eps) {
    int i = 0;
    int j = 0;
    while (i < m && j < n) {
        // Find the maximal element in j-th column
        // starting from i-th row
        double maxVal = fabs((*this)[i][j]);
        int imax = i;
        for (int k = i+1; k < m; ++k) {
            if (fabs((*this)[k][j]) > maxVal) {
                maxVal = fabs((*this)[k][j]);
                imax = k;
            }
        }
        if (maxVal <= eps) {
            // Zero column
            ++j;
            continue;
        }
        
        if (i != imax) {
            swapRows(i, imax);
        }
        
        assert(fabs((*this)[i][j]) > eps);
        double r = (*this)[i][j]; // Resolving element
        for (int k = i+1; k < m; ++k) {
            double c = (-(*this)[k][j]/r);
            addRows(k, i, c);
        }
        ++i; 
        ++j;
    }       
    return i;   // Rank of matrix == i == 
                //     number of non-zero rows
}

double Matrix::det() const {
    if (m != n)
        throw std::invalid_argument(
            "determinant of non-square matrix"
        );
    Matrix b(*this);
    int r = b.gaussElimination();
    if (r < m)
        return 0.;
    double d = b[0][0];
    for (int i = 1; i < m; ++i) {
        d *= b[i][i];
    }
    return d;
}