#ifndef MATRIX_H
#define MATRIX_H 1

#include <cmath>

class MatrixException {
private:
    const char* errorReason;
    
public:
    MatrixException(const char* reason = ""):
        errorReason(reason)
    {}
    
    const char* what() const {
        return errorReason;
    }
};    

// Rectangular matrix of type double
class Matrix {
private:
    int m;      // Number if rows
    int n;      // Number of columns
    double *a;  // Array of matrix elements
public:
    Matrix(int numRows = 1, int numCols = 1):
        m(numRows),
        n(numCols),
        a(new double[m*n])
    {
        erase();
    }
    
    void erase(double v = 0.) {
        for (int i = 0; i < m*n; ++i) {
            a[i] = v;
        }
    }
    
    // Copy-constructor
    Matrix(const Matrix& matr):
        m(matr.m),
        n(matr.n),
        a(new double[m*n])
    {
        for (int i = 0; i < m*n; ++i) {
            a[i] = matr.a[i];
        }
    }
    
    // Destructor: release memory
    ~Matrix() { delete[] a; }
    
    // copy-operator
    Matrix& operator=(const Matrix& matr) {
        if (m*n < matr.m*matr.n) {
            // Not enougth memory
            delete[] a;
            a = new double[matr.m*matr.n];
        }
        m = matr.m;
        n = matr.n;
        for (int i = 0; i < m*n; ++i) {
            a[i] = matr.a[i];
        }
        return *this;
    }
           
    void resize(int numRows, int numCols) {
        if (m*n < numRows*numCols) {
            // Not enougth memory
            delete[] a;
            a = new double[numRows*numCols];
        }
        m = numRows;
        n = numCols;
        erase();
    }
    
    int rows() const { return m; }
    int cols() const { return n; }
    
    double at(int i, int j) const { // x = abc.at[3, 4];
        if (
            i < 0 || i >= m ||
            j < 0 || j >= n
        ) {
            throw MatrixException("indices out of bounds");
        }
        return a[i*n + j];
    }   

    double& at(int i, int j) {  // abc.at(1, 2) = 1.5;
        if (
            i < 0 || i >= m ||
            j < 0 || j >= n
        ) {
            throw MatrixException("indices out of bounds");
        }
        return a[i*n + j];
    }
    
    // x = abc[i][j];
    // abc[l][k] = 17.3;
    const double* operator[](int i) const {
        return a + i*n;
    }

    double* operator[](int i) {
        return a + i*n;
    }
    
    Matrix operator+(const Matrix& matr) const;
    
    Matrix operator*(const Matrix& mart) const;

    // Row echelon form
    int gauss(double eps = 1e-8); // Return rank of matrix
    
    double det() const;
    
    void inverse(Matrix& inverted) const;
};

std::ostream& operator<<(std::ostream& s, const Matrix& matrix);
std::istream& operator>>(std::istream& s, Matrix& matrix);       

#endif /* MATRIX_H */