template <class T>
class matrix {
  public:
    unsigned int rows, cols;
    T* contents;
    matrix(unsigned int, unsigned int);
    matrix(const matrix<T>&);
    matrix(unsigned int, unsigned int, T*);
    matrix<T> operator*(const matrix<T>&) const;
	matrix<T> operator+(const matrix<T>&) const;
    ~matrix();
    T* operator[](unsigned int) const;
	const matrix<T>& operator=(const matrix<T>&);
};

template<class T>
const matrix<T>& matrix<T>::operator=(const matrix<T>& m) {
    rows = m.rows;
    cols = m.cols;
    delete [] contents;
    contents = new T[rows*cols];
    for (unsigned int i = 0; i<rows*cols; i++)
        contents[i] = m.contents[i];
    return m;
}

template<class T>
T* matrix<T>::operator[](unsigned int i) const {
    return &contents[i*cols];
}

template<class T>
matrix<T>::matrix(unsigned int r, unsigned int c) : rows(r), cols(c), contents(new T[r*c]) {
    for (unsigned int i = 0; i < r*c; i++)
        contents[i] = 0;
}

template<class T>
matrix<T>::matrix(const matrix<T>& m)
  : rows(m.rows), cols(m.cols), contents(new T[m.rows*m.cols]) {
    for (unsigned int i = 0; i < m.rows*m.cols; i++)
        contents[i] = m.contents[i];
}

template<class T>
matrix<T>::matrix(unsigned int r, unsigned int c, T* t) : rows(r), cols(c), contents(t) { }

template<class T>
matrix<T>::~matrix() {
    delete [] contents;
}

template<class T>
matrix<T> matrix<T>::operator+(const matrix<T>& m) const {
    //Assume dimensions match, because we're using C++
    T* out = new T[rows*cols];
    for (unsigned int i = 0; i < rows*cols; i++)
        out[i] = contents[i] + m.contents[i];
    return matrix(rows, cols, out);
}

template<class T>
matrix<T> matrix<T>::operator*(const matrix<T>& m) const {
    //Assume dimensions match, because we're using C++
    //We're using row-major order
    matrix<T> out(rows, m.cols);
    for (unsigned int i = 0; i < rows; i++) {
        for (unsigned int j = 0; j < m.cols; j++) {
            T sum = 0;
            for (unsigned int x = 0; x < cols; x++)
                sum += (*this)[i][x] * m[x][j];
            out[i][j] = sum;
        }
    }
    return out;
}