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

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

template<class T>
matrix<T>::matrix(int r, int c) : rows(r), cols(c), contents(new T[r*c]) {
    for (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 (int i = 0; i < m.rows*m.cols; i++)
        contents[i] = m.contents[i];
}

template<class T>
matrix<T>::matrix(int r, 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) {
    //Assume dimensions match, because we're using C++
    T* out = new T[rows*cols];
    for (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) {
    //Assume dimensions match, because we're using C++
    //We're using row-major order
    matrix<T> out(rows, m.cols);
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < m.cols; j++) {
            T sum = 0;
            for (int x = 0; x < cols; x++)
                sum += (*this)[i][x] * m.contents[x*m.cols+j];
                //todo: figure out why m[x][j] fails to compile here
            out[i][j] = sum;
        }
    }
    return out;
}