#ifndef MATRIX_H
#define MATRIX_H

#include <iostream>
#include <vector>
#include <cmath>
#include <stdexcept>

class Matrix {
private:
    int m;  // Number of rows
    int n;  // Number of columns
    std::vector<double> elems;
public:
    Matrix(int r = 1, int c = 1):
        m(r),
        n(c),
        elems(m*n)
    {}
    
    /* Not necessary!
    Copy constructor, move constructor, 
    operator=, destructor
    */
    
    // Copy-constructor
    Matrix(const Matrix& a):
        m(a.m),
        n(a.n),
        elems(a.elems)
    {}
    
    // Move-constructor
    Matrix(Matrix&& a):
        m(a.m),
        n(a.n),
        elems(std::move(a.elems))
    {
        a.m = 0; a.n = 0;
    }
    
    // Copy-operator=
    Matrix& operator=(const Matrix& a) {
        m = a.m;
        n = a.n;
        elems = a.elems;
        return *this;
    }
    
    // Move-operator=
    Matrix& operator=(Matrix&& a) {
        m = a.m;
        n = a.n;
        elems = std::move(a.elems);
        a.m = 0;
        a.n = 0;
        return *this;
    }
    
    ~Matrix() {}
    
    int rows() const { return m; }
    int cols() const { return n; }
    
    void resize(int r, int c) {
        elems.resize(r*c);
        m = r; n = c;
        fill(0.);
    }
    
    void fill(double x = 0.) {
        for (int i = 0; i < m*n; ++i)
            elems[i] = x;
    }

    const double* operator[](int i) const {
        return elems.data() + i*n;    
    }

    double* operator[](int i) {
        return elems.data() + i*n;    
    }
    
    double at(int i, int j) const {
        if (i < 0 || i >= m ||
            j < 0 || j >= n
        )
            throw std::out_of_range(
                "matrix indices out of range"
            );
        return elems[i*n + j];
    }

    double& at(int i, int j) {
        if (i < 0 || i >= m ||
            j < 0 || j >= n
        )
            throw std::out_of_range(
                "matrix indices out of range"
            );
        return elems[i*n + j];
    }

    int gauss(double eps = 1e-10);  // Returns rank of matrix    
    int rref(double eps = 1e-10);   // Reduced row echelon form
    double det() const;
};

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

#endif