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

int Matrix::gauss(double eps) {
    int i = 0;
    int j = 0;
    while (i < m && j < n) {
        // 1. Search for maximal element in j-th column
        //    starting from i-th row
        double amax = fabs((*this)[i][j]);
        int k = i;
        for (int l = i+1; l < m; ++l) {
            if (fabs((*this)[l][j]) > amax) {
                amax = fabs((*this)[l][j]);
                k = l;
            }
        }
        
        if (amax <= eps) {
            ++j; continue;
        }
        
        if (k != i) {
            // Swap rows i and k
            for (int s = j; s < n; ++s) {
                double tmp = (*this)[i][s];
                (*this)[i][s] = (*this)[k][s];
                (*this)[k][s] = (-tmp);
            }
        }
        
        // Annihilate j-th column starting from i+1-th row
        for (int k = i+1; k < m; ++k) {
            double coeff = (*this)[k][j]/(*this)[i][j];
            for (int s = j; s < n; ++s) {
                (*this)[k][s] -= coeff*(*this)[i][s];
            }
        }    
        ++i;
        ++j;        
    }
    return i;
}

int Matrix::rref(double eps) {  // Reduced row echelon form     
    int rank = gauss();
    int i = rank - 1;
    while (i >= 0) {
        // Find the first nonzero element in i-th row
        int j = 0;
        while (j < cols() && fabs((*this)[i][j]) <= eps) {
            ++j;
        }
        assert(j < cols());
        
        // Divide the i-th row by the first nonzero element a[i][j]
        double coeff = 1./(*this)[i][j]; 
        (*this)[i][j] = 1.;
        for (int k = j+1; k < cols(); ++k) {
            (*this)[i][k] *= coeff;
        }
        
        // Annihilate the j-th column above i-th row
        for (int l = 0; l < i; ++l) {
            double coeff = (*this)[l][j];
            (*this)[l][j] = 0.;
            for (int k = j+1; k < cols(); ++k) {
                (*this)[l][k] -= (*this)[i][k]*coeff;
            }
        }
        --i;    // To the next row in reverse order
    }
    return rank;
}

double Matrix::det() const {
    if (rows() != cols()) {
        throw std::invalid_argument(
            "determinant of non-square matrix"
        );
    }
    Matrix a = *this;
    int rank = a.gauss();
    if (rank < a.rows())
        return 0.;
    double d = 1.;
    for (int i = 0; i < a.rows(); ++i) {
        d *= a[i][i];
    }
    return d;
}


std::ostream& operator<<(std::ostream& s, const Matrix& a) {
    for (int i = 0; i < a.rows(); ++i) {
        for (int j = 0; j < a.cols(); ++j) {
            if (j > 0)
                s << ' ';   // Delimiter
            s << a[i][j];
        }
        s << std::endl;
    }
    return s;
}

std::istream& operator>>(std::istream& s, Matrix& a) {
    for (int i = 0; s.good() && i < a.rows(); ++i) {
        for (int j = 0; s.good() && j < a.cols(); ++j) {
            double x;
            if (s >> x) {
                a[i][j] = x;
            }            
        }
    }
    return s;
}