#ifndef __UNIT3_EXAMPLE_MATRIX_H__
#define __UNIT3_EXAMPLE_MATRIX_H__

#include <cstdlib>

class Matrix
{
public:
    Matrix(size_t rows, size_t cols);
    ~Matrix();

    void resize(size_t rows, size_t cols);

    double get_elem(size_t r, size_t c);
    void set_elem(size_t r, size_t c, double val);

    size_t rows();
    size_t cols();

private:
    double *_data;
    size_t _rows, _cols;
};

// computes C = A * B;
bool Matrix_multiply(Matrix *C, Matrix *A, Matrix *B);

#endif // ! defined __UNIT3_EXAMPLE_MATRIX_H__

#ifndef __UNIT3_EXAMPLE_VECTOR_H__
#define __UNIT3_EXAMPLE_VECTOR_H__

#include <cstdlib>

class Vector
{
public:
    Vector(size_t len);
    ~Vector();

    void resize(size_t len);

    double norm();
    double norm(double p = 2);

    double get_elem(size_t n);
    void set_elem(size_t n, double val);

    size_t length();

private:
    double *_data;
    size_t _len;
};

#endif // ! defined __UNIT3_EXAMPLE_VECTOR_H__

#include "matrix.h"

#include <cmath>

Matrix::Matrix(size_t rows, size_t cols)
{
    _rows = rows;
    _cols = cols;

    if (rows * cols == 0) {
        _data = NULL;
    } else {
        _data = new double[rows * cols];
    }
}

Matrix::~Matrix()
{
    delete[] _data;
}

double
Matrix::get_elem(size_t r, size_t c)
{
    if (r >= _rows || c >= _cols)
        return NAN;

    return _data[c + r * _cols];
}

void
Matrix::set_elem(size_t r, size_t c, double val)
{
    if (r >= _rows || c >= _cols)
        return;

    _data[c + r * _cols] = val;
}

void
Matrix::resize(size_t rows, size_t cols)
{
    _rows = rows;
    _cols = cols;

    delete[] _data;

    if (rows * cols == 0) {
        _data = NULL;
    } else {
        _data = new double[rows * cols];
    }
}

size_t
Matrix::rows()
{
    return _rows;
}

size_t
Matrix::cols()
{
    return _cols;
}

bool
Matrix_multiply(Matrix *C, Matrix *A, Matrix *B)
{
    if (A->cols() != B->rows())
        return false;

    C->resize(A->rows(), B->cols());

    for (int i = 0; i < A->rows(); ++i) {
        for (int j = 0; j < B->cols(); ++j) {
            double c = 0;
            for (int k = 0; k < A->cols(); ++k) {
                c += A->get_elem(i, k) * B->get_elem(k, j);
            }
            C->set_elem(i, j, c);
        }
    }

    return true;
}

#include "vector.h"

#include <cmath>

using namespace std;

Vector::Vector(size_t len)
{
    _len = len;

    if (len == 0) {
        _data = NULL;
    } else {
        _data = new double[len];
    }
}

Vector::~Vector()
{
    delete[] _data;
}

double
Vector::norm()
{
    norm(2);
}

double
Vector::norm(double p)
{
    if (p < 1) {
        return NAN;
    } else if (p == 1) {
        double sum = 0;
        for (int i = 0; i < _len; ++i) {
            sum += abs(_data[i]);
        }
        return sum;
    } else if (p == 2) {
        double sum = 0;
        for (int i = 0; i < _len; ++i) {
            sum += _data[i] * _data[i];
        }
        return sqrt(sum);
    } else if (! isfinite(p)) {
        double max = 0;
        for (int i = 0; i < _len; ++i) {
            if (abs(_data[i]) > max)
                max = abs(_data[i]);
        }
        return max;
    } else {
        double sum = 0;
        for (int i = 0; i < _len; ++i) {
            sum += pow(_data[i], p);
        }
        return pow(sum, 1/p);
    }
}

double
Vector::get_elem(size_t n)
{
    return _data[n];
}

void
Vector::set_elem(size_t n, double val)
{
    _data[n] = val;
}

size_t
Vector::length()
{
    return _len;
}