#ifndef MATRIX_H
#define MATRIX_H

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

class Matrix {
private:
    int m;  // Number of rows
    int n;  // Number of columns
    std::vector<double> elems;
public:
    Matrix(int num_rows = 1, int num_cols = 1):
        m(num_rows),
        n(num_cols),
        elems(m*n)
    {}
    
    Matrix(const Matrix& a):    // Copy constructor
        m(a.m),
        n(a.n),
        elems(a.elems)
    {}
    
    Matrix(Matrix&& a):     // Move constructor
        m(a.m),
        n(a.n),
        elems(std::move(a.elems))
    {
        a.m = 0; a.n = 0;
    }
    
    Matrix& operator=(const Matrix& a) {
        m = a.m;
        n = a.n;
        elems = a.elems;
        return *this;
    }
    
    Matrix& operator=(Matrix&& a) {
        if (this != &a) {
            m = a.m;
            n = a.n;
            a.m = 0;
            a.n = 0;
            elems = std::move(a.elems);
        }
        return *this;
    }
    
    int rows() const { 
        return m;
    }
    int columns() const {
        return n;
    }
    int cols() const {
        return columns();
    }

    const double* operator[](int i) const {
        return elems.data() + i*n;
    }
    
    // a[i][j] = y;
    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 index out of bounds");
        return elems[i*n + j];
    }

    // a.at(i, j) = y;
    double& at(int i, int j) {
        if (
            i < 0 || i >= m ||
            j < 0 || j >= n
        )
            throw std::out_of_range("Matrix index out of bounds");
        return elems[i*n + j];
    }

    Matrix operator+(const Matrix& b) const;
    Matrix operator-(const Matrix& b) const;
    Matrix operator*(const Matrix& b) const;
    Matrix operator/(const Matrix& b) const;

    // Gauss elimination (transform the matrix to 
    // the row echelon form)
    int gauss(double eps = 1e-12); // Return the rank of matrix
    int rref(double eps = 1e-12);  // Reduced row echelon form 

    double determinant(double eps = 1e-12) const;    
    Matrix inverse() const;
    void solveLinearSystem(
        const double* freeTerms,
        double* solution
    ) const;
};

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

#endif