3 files changed, 693 insertions, 0 deletions
diff --git a/include/mm b/include/mm/mm.hpp
index 7e9f02b..7e9f02b 100644
--- a/include/mm
+++ b/include/mm/mm.hpp
diff --git a/include/mm/mmmatrix.hpp b/include/mm/mmmatrix.hpp
new file mode 100644
index 0000000..5285429
--- /dev/null
+++ b/include/mm/mmmatrix.hpp
@@ -0,0 +1,374 @@
+/* mmmatrix.hpp
+ * Part of Mathematical library built (ab)using Modern C++ 17 abstractions.
+ *
+ * This library is not intended to be _performant_, it does not contain
+ * hand written SMID / SSE / AVX optimizations. It is instead an example
+ * of highly inefficient (but abstract!) code, where matrices can contain any 
+ * data type.
+ *
+ * Naoki Pross <naopross@thearcway.org>
+ * 2018 ~ 2019
+ */
+#pragma once
+
+#include <iostream>
+#include <iomanip>
+#include <cstring>
+#include <cassert>
+#include <initializer_list>
+#include <array>
+
+namespace mm {
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    class basic_matrix;
+
+    /* specialization of basic_matrx for Cols = 1 */
+    template<typename T, std::size_t Rows>
+    class row_vec;
+
+    /* specialization of basic_matrx for Rows = 1 */
+    template<typename T, std::size_t Cols>
+    class col_vec;
+
+    /* shorter name for basic_matrix */
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    class matrix;
+
+    /* specialization of basic_matrix for Rows == Cols */
+    template<typename T, std::size_t N>
+    class square_matrix;
+
+    // template<typename T, std::size_t N>
+    // class diag_matrix;
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::basic_matrix {
+public:
+    using type = T;
+
+    template<typename U, std::size_t ORows, std::size_t OCols>
+    friend class mm::basic_matrix;
+
+    static constexpr std::size_t rows = Rows;
+    static constexpr std::size_t cols = Cols;
+
+    basic_matrix();
+
+    // from initializer_list
+    basic_matrix(std::initializer_list<std::initializer_list<T>> l);
+
+    // copyable and movable
+    basic_matrix(const basic_matrix<T, Rows, Cols>& other);
+    basic_matrix(basic_matrix<T, Rows, Cols>&& other);
+
+    // copy from another matrix
+    template<std::size_t ORows, std::size_t OCols>
+    basic_matrix(const basic_matrix<T, ORows, OCols>& other);
+
+    // access data
+    T& at(std::size_t row, std::size_t col);
+    const T& at(std::size_t row, std::size_t col) const;
+    // allows to access a matrix M at row j col k with M[j][k]
+    auto operator[](std::size_t index);
+
+    void swap_rows(std::size_t x, std::size_t y);
+    void swap_cols(std::size_t x, std::size_t y);
+
+    // mathematical operations
+    virtual basic_matrix<T, Cols, Rows> transposed() const;
+    inline basic_matrix<T, Cols, Rows> td() const { return transposed(); }
+
+    // bool is_invertible() const;
+    // basic_matrix<T, Rows, Cols> inverse() const;
+
+
+    /// downcast to square matrix
+    static inline constexpr bool is_square() { return (Rows == Cols); }
+    inline constexpr square_matrix<T, Rows> to_square() const {
+        static_assert(is_square());
+        return static_cast<square_matrix<T, Rows>>(*this);
+    }
+
+
+    /// downcast to row_vector
+    static inline constexpr bool is_row_vec() { return (Cols == 1); }
+    inline constexpr row_vec<T, Rows> to_row_vec() const {
+        static_assert(is_row_vec());
+        return static_cast<row_vec<T, Rows>>(*this);
+    }
+
+    /// downcast to col_vector
+    static inline constexpr bool is_col_vec() { return (Rows == 1); }
+    inline constexpr col_vec<T, Cols> to_col_vec() const {
+        static_assert(is_col_vec());
+        return static_cast<col_vec<T, Cols>>(*this);
+    }
+
+protected:
+    template<typename ConstIterator>
+    basic_matrix(ConstIterator begin, ConstIterator end);
+
+private:
+    std::array<T, Rows * Cols> data;
+};
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::basic_matrix<T, Rows, Cols>::basic_matrix() {
+    std::fill(data.begin(), data.end(), 0);
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::basic_matrix<T, Rows, Cols>::basic_matrix(
+    std::initializer_list<std::initializer_list<T>> l
+) {
+    assert(l.size() == Rows);
+    auto data_it = data.begin();
+
+    for (auto&& row : l) {
+        data_it = std::copy(row.begin(), row.end(), data_it);
+    }
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::basic_matrix<T, Rows, Cols>::basic_matrix(
+    const mm::basic_matrix<T, Rows, Cols>& other
+) : data(other.data) {}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::basic_matrix<T, Rows, Cols>::basic_matrix(
+    mm::basic_matrix<T, Rows, Cols>&& other
+) : data(std::forward<decltype(other.data)>(other.data)) {}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+template<std::size_t ORows, std::size_t OCols>
+mm::basic_matrix<T, Rows, Cols>::basic_matrix(
+    const mm::basic_matrix<T, ORows, OCols>& other
+) {
+    static_assert((ORows <= Rows),
+        "cannot copy a taller matrix into a smaller one"
+    );
+
+    static_assert((OCols <= Cols),
+        "cannot copy a larger matrix into a smaller one"
+    );
+
+    std::fill(data.begin(), data.end(), 0);
+    for (unsigned row = 0; row < Rows; row++)
+        for (unsigned col = 0; col < Cols; col++)
+            this->at(row, col) = other.at(row, col);
+}
+
+/* protected construtor */
+template<typename T, std::size_t Rows, std::size_t Cols>
+template<typename ConstIterator>
+mm::basic_matrix<T, Rows, Cols>::basic_matrix(
+    ConstIterator begin, ConstIterator end
+) {
+    assert(static_cast<unsigned>(std::distance(begin, end)) >= ((Rows * Cols)));
+    std::copy(begin, end, data.begin());
+}
+
+
+/* member functions */
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+T& mm::basic_matrix<T, Rows, Cols>::at(std::size_t row, std::size_t col) {
+    assert(row < Rows); // "out of row bound"
+    assert(col < Cols); // "out of column bound"
+
+    return data[row * Cols + col];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+const T& mm::basic_matrix<T, Rows, Cols>::at(std::size_t row, std::size_t col) const {
+    assert(row < Rows); // "out of row bound"
+    assert(col < Cols); // "out of column bound"
+
+    return data[row * Cols + col];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+auto mm::basic_matrix<T, Rows, Cols>::operator[](std::size_t index) {
+    if constexpr (is_row_vec() || is_col_vec()) {
+        return data.at(index);
+    } else {
+        return row_vec<T, Rows>(
+            data.cbegin() + (index * Cols),
+            data.cbegin() + ((index + 1) * Cols) + 1
+        );
+    }
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+void mm::basic_matrix<T, Rows, Cols>::swap_rows(std::size_t x, std::size_t y) {
+    if (x == y)
+        return;
+
+    for (unsigned col = 0; col < Cols; col++)
+        std::swap(this->at(x, col), this->at(y, col));
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+void mm::basic_matrix<T, Rows, Cols>::swap_cols(std::size_t x, std::size_t y) {
+    if (x == y)
+        return;
+
+    for (unsigned row = 0; row < rows; row++)
+        std::swap(this->at(row, x), this->at(row, y));
+}
+
+template<typename T, std::size_t M, std::size_t N>
+mm::basic_matrix<T, N, M> mm::basic_matrix<T, M, N>::transposed() const {
+    mm::basic_matrix<T, N, M> result;
+
+    for (unsigned row = 0; row < M; row++)
+        for (unsigned col = 0; col < N; col++)
+            result.at(col, row) = this->at(row, col);
+
+    return result;
+}
+
+
+/* operator overloading */
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::basic_matrix<T, Rows, Cols> operator+(
+    const mm::basic_matrix<T, Rows, Cols>& a,
+    const mm::basic_matrix<T, Rows, Cols>& b
+) {
+    mm::basic_matrix<T, Rows, Cols> result;
+
+    for (unsigned row = 0; row < Rows; row++)
+        for (unsigned col = 0; col < Cols; col++)
+            result.at(row, col) = a.at(row, col) + b.at(row, col);
+    
+    return result;
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::basic_matrix<T, Rows, Cols> operator*(
+    const mm::basic_matrix<T, Rows, Cols>& m,
+    const T& scalar
+) {
+    mm::basic_matrix<T, Rows, Cols> result;
+    for (unsigned row = 0; row < Rows; row++)
+        for (unsigned col = 0; col < Cols; col++)
+            result.at(row, col) = m.at(row, col) * scalar;
+
+    return result;
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::basic_matrix<T, Rows, Cols> operator*(
+    const T& scalar,
+    const mm::basic_matrix<T, Rows, Cols>& m
+) {
+    return m * scalar;
+}
+
+template<typename T, std::size_t M, std::size_t P1, std::size_t P2, std::size_t N>
+mm::basic_matrix<T, M, N> operator*(
+    const mm::basic_matrix<T, M, P1>& a,
+    const mm::basic_matrix<T, P2, N>& b
+) {
+    static_assert(P1 == P2, "invalid matrix multiplication");
+    mm::basic_matrix<T, M, N>  result;
+
+    // TODO: use a more efficient algorithm
+    for (unsigned row = 0; row < M; row++)
+        for (unsigned col = 0; col < N; col++)
+            for (unsigned k = 0; k < P1; k++)
+                result.at(row, col) = a.at(row, k) * b.at(k, col);
+
+    return result;
+}
+
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::basic_matrix<T, Rows, Cols> operator-(
+    const mm::basic_matrix<T, Rows, Cols>& a,
+    const mm::basic_matrix<T, Rows, Cols>& b
+) {
+    return a + (static_cast<T>(-1) * b);
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols, unsigned NumW = 3>
+std::ostream& operator<<(std::ostream& os, const mm::basic_matrix<T, Rows, Cols>& m) {
+    for (unsigned row = 0; row < Rows; row++) {
+        os << "[ ";
+        for (unsigned col = 0; col < (Cols -1); col++) {
+            os << std::setw(NumW) << m.at(row, col) << ", ";
+        }
+        os << std::setw(NumW) << m.at(row, (Cols -1)) << " ]\n";
+    }
+
+    return os;
+}
+
+
+
+
+/* row vector specialization */
+template<typename T, std::size_t Rows>
+class mm::row_vec : public mm::basic_matrix<T, Rows, 1> {
+public:
+    using mm::basic_matrix<T, Rows, 1>::basic_matrix;
+};
+
+/* column vector specialization */
+template<typename T, std::size_t Cols>
+class mm::col_vec : public mm::basic_matrix<T, 1, Cols> {
+public:
+    using mm::basic_matrix<T, 1, Cols>::basic_matrix;
+};
+
+/* general specialization (alias) */
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::matrix : public mm::basic_matrix<T, Rows, Cols> {
+public:
+    using mm::basic_matrix<T, Rows, Cols>::basic_matrix;
+};
+
+/* square matrix specializaiton */
+template<typename T, std::size_t N>
+class mm::square_matrix : public mm::basic_matrix<T, N, N> {
+public:
+    using mm::basic_matrix<T, N, N>::basic_matrix;
+
+    /// in place transpose
+    void transpose();  
+    inline void t() { transpose(); }
+
+    T trace();
+    inline T tr() { return trace(); }
+
+    /// in place inverse
+    void invert();
+
+    // get the identity of size N
+    static inline constexpr square_matrix<T, N> identity() {
+        square_matrix<T, N> i;
+        for (unsigned row = 0; row < N; row++)
+            for (unsigned col = 0; col < N; col++)
+                i.at(row, col) = (row == col) ? 1 : 0;
+
+        return i;
+    }
+};
+
+
+template<typename T, std::size_t N>
+void mm::square_matrix<T, N>::transpose() {
+    for (unsigned row = 0; row < N; row++)
+        for (unsigned col = 0; col < row; col++)
+            std::swap(this->at(row, col), this->at(col, row));
+}
+
+template<typename T, std::size_t N>
+T mm::square_matrix<T, N>::trace() {
+    T sum = 0;
+    for (unsigned i = 0; i < N; i++)
+        sum += this->at(i, i);
+
+    return sum;
+}
diff --git a/include/mm/mmvec.hpp b/include/mm/mmvec.hpp
new file mode 100644
index 0000000..1939388
--- /dev/null
+++ b/include/mm/mmvec.hpp
@@ -0,0 +1,319 @@
+/* mmvec.hpp
+ * Part of Mathematical library built (ab)using Modern C++ 17 abstractions.
+ *
+ * This library is not intended to be _performant_, it does not contain
+ * hand written SMID / SSE / AVX optimizations. It is instead an example
+ * of highly abstracted code, where Vectors can contain any data type.
+ *
+ * As a challenge, the vector data structure has been built on a container
+ * of static capacity. But if a dynamic base container is needed, the code
+ * should be easily modifiable to add further abstraction, by templating
+ * the container, and by consequence the allocator.
+ *
+ * Naoki Pross <naopross@thearcway.org>
+ * 2018 ~ 2019
+ */
+#pragma once
+
+#include <cassert>
+#include <cmath>
+
+#include <iostream>
+#include <array>
+#include <algorithm>
+#include <numeric>
+#include <complex>
+#include <initializer_list>
+
+namespace mm {
+    // generic implementation
+    template<typename T, std::size_t d>
+    struct basic_vec;
+
+    // usable specializations
+    template<typename T, std::size_t d>
+    struct vec;
+    template<typename T>
+    struct vec3;
+    template<typename T>
+    struct vec2;
+}
+
+template<typename T, std::size_t d>
+struct mm::basic_vec : public std::array<T, d> {
+    using type = T;
+    static constexpr std::size_t dimensions = d;
+
+    // TODO: template away these
+    static constexpr T null_element = static_cast<T>(0);
+    static constexpr T unit_element = static_cast<T>(1);
+    static constexpr T unit_additive_inverse_element = static_cast<T>(-1);
+
+    basic_vec();
+    basic_vec(const std::initializer_list<T> l);
+
+    // copyable to a vector of size n <= d
+    template<std::size_t n> basic_vec(const basic_vec<T, n>& other);
+    // movable to a vector of size n <= d
+    template<std::size_t n> basic_vec(basic_vec<T, n>&& other);
+
+    T length() const;
+
+    // copy operator=
+    template<std::size_t n>
+    basic_vec<T, d>& operator=(const mm::basic_vec<T, n>& other);
+
+    // move operator=
+    template<std::size_t n>
+    basic_vec<T, d>& operator=(mm::basic_vec<T, n>&& other);
+
+    template<std::size_t n>
+    basic_vec<T, d>& operator+=(const mm::basic_vec<T, n>& other);
+
+    template<std::size_t n>
+    basic_vec<T, d>& operator-=(const mm::basic_vec<T, n>& other);
+
+
+    basic_vec<T, d>& operator*=(const T& scalar);
+};
+
+
+// member functions for basic_vec
+
+template<typename T, std::size_t d>
+mm::basic_vec<T, d>::basic_vec() : std::array<T, d>() {
+    this->fill(basic_vec<T, d>::null_element);
+}
+
+template<typename T, std::size_t d>
+mm::basic_vec<T, d>::basic_vec(const std::initializer_list<T> l) {
+    // why can't this sh*t be a constexpr with static_assert???
+    assert(l.size() <= d);
+    std::copy(l.begin(), l.end(), this->begin());
+}
+
+template<typename T, std::size_t d>
+template<std::size_t n>
+mm::basic_vec<T, d>::basic_vec(const mm::basic_vec<T, n>& other) {
+    *this = other;
+}
+
+template<typename T, std::size_t d>
+template<std::size_t n>
+mm::basic_vec<T, d>::basic_vec(basic_vec<T, n>&& other) {
+    *this = other;
+}
+
+template<typename T, std::size_t d>
+T mm::basic_vec<T, d>::length() const {
+    return std::sqrt(std::accumulate(this->begin(), this->end(),
+        basic_vec<T, d>::null_element,
+        [](const T& init, const T& val) -> T {
+            return init + val * val;
+        }
+    ));
+}
+
+
+// memeber operator overloads for basic_vec
+
+template<typename T, std::size_t d>
+template<std::size_t n>
+mm::basic_vec<T, d>& mm::basic_vec<T, d>::operator=(const mm::basic_vec<T, n>& other) {
+    static_assert(
+        d >= n, "cannot copy higher dimensional vector into a smaller one"
+    );
+
+    std::copy(other.begin(), other.end(), this->begin());
+
+    return *this;
+}
+
+template<typename T, std::size_t d>
+template<std::size_t n>
+mm::basic_vec<T, d>& mm::basic_vec<T, d>::operator=(mm::basic_vec<T, n>&& other) {
+    static_assert(
+        d >= n, "cannot move a higher dimensional vector into a smaller one"
+    );
+
+    std::move(other.begin(), other.end(), this->begin());
+
+    return *this;
+}
+
+
+template<typename T, std::size_t d>
+template<std::size_t n>
+mm::basic_vec<T, d>& mm::basic_vec<T, d>::operator+=(const mm::basic_vec<T, n>& other) {
+    *this = *this + other;
+    return *this;
+}
+
+template<typename T, std::size_t d>
+template<std::size_t n>
+mm::basic_vec<T, d>& mm::basic_vec<T, d>::operator-=(const mm::basic_vec<T, n>& other) {
+    *this = *this - other;
+    return *this;
+}
+
+template<typename T, std::size_t d>
+mm::basic_vec<T, d>& mm::basic_vec<T, d>::operator*=(const T& scalar) {
+    *this = *this * scalar;
+    return *this;
+}
+
+
+// operator overloads for basic_vec
+
+template<typename T, std::size_t d>
+mm::basic_vec<T, d> operator+(const mm::basic_vec<T, d>& rhs, const mm::basic_vec<T, d>& lhs) {
+    mm::basic_vec<T, d> out;
+    
+    std::transform(rhs.begin(), rhs.end(), lhs.begin(), out.begin(),
+        [](const T& r, const T& l) -> T {
+            return r + l;
+        }
+    );
+
+    return out;
+}
+
+template<typename T, std::size_t d>
+mm::basic_vec<T, d> operator*(const mm::basic_vec<T, d>& rhs, const T& lhs) {
+    return lhs * rhs;
+}
+
+template<typename T, std::size_t d>
+mm::basic_vec<T, d> operator*(const T& rhs, const mm::basic_vec<T, d>& lhs) {
+    mm::basic_vec<T, d> out;
+
+    std::transform(lhs.begin(), lhs.end(), out.begin(), 
+        [rhs](const T& t) -> T {
+            return t * rhs;
+    });
+
+    return out;
+}
+
+template<typename T, std::size_t d>
+mm::basic_vec<T, d> operator-(const mm::basic_vec<T, d>& rhs, const mm::basic_vec<T, d>& lhs) {
+    return rhs + mm::basic_vec<T, d>::unit_additive_inverse_element * lhs;
+}
+
+template<typename T, std::size_t d>
+T operator*(const mm::basic_vec<T, d>& rhs, const mm::basic_vec<T, d>& lhs) {
+    return std::inner_product(rhs.begin(), rhs.end(), lhs.begin(), 0);
+}
+
+template<typename T, std::size_t d>
+std::ostream& operator<<(std::ostream& os, const mm::basic_vec<T, d>& v) {
+    os << "<";
+    std::for_each(v.begin(), v.end() -1, [&](const T& el) {
+        os << el << ", ";
+    });
+    os << v.back() << ">";
+
+    return os;
+}
+
+
+// actual vectors to use in your code
+
+template<typename T, std::size_t d>
+class mm::vec: public mm::basic_vec<T, d> {
+public:
+    using mm::basic_vec<T, d>::basic_vec;
+};
+
+
+// three dimensional specialization with a static cross product
+// TODO: specialize operator+ for spherical coordinates
+
+template<typename T>
+class mm::vec3 : public mm::basic_vec<T, 3> {
+public:
+    using mm::basic_vec<T, 3>::basic_vec;
+
+    T& x() { return this->at(0); }
+    T& y() { return this->at(1); }
+    T& z() { return this->at(2); }
+
+    const T& x() const { return this->at(0); }
+    const T& y() const { return this->at(1); }
+    const T& z() const { return this->at(2); }
+
+    T zenith() const;
+    T azimuth() const;
+    vec3<T> spherical() const;
+
+    static vec3<T> cross(const vec3<T>& rhs, const vec3<T>& lhs);
+};
+
+template<typename T>
+T mm::vec3<T>::zenith() const {
+    return std::acos(this->z() / this->length());
+}
+
+template<typename T>
+T mm::vec3<T>::azimuth() const {
+    return std::atan(this->y() / this->x());
+}
+
+template<typename T>
+mm::vec3<T> mm::vec3<T>::spherical() const {
+    return mm::vec3<T> {
+        this->length(),
+        this->zenith(),
+        this->azimuth(),
+    };
+}
+
+template<typename T>
+mm::vec3<T> mm::vec3<T>::cross(const vec3<T>& rhs, const vec3<T>& lhs) {
+    mm::vec3<T> res;
+
+    res.x() = (rhs.y() * lhs.z()) - (rhs.z() * lhs.y());
+    res.y() = (rhs.z() * lhs.x()) - (rhs.x() * lhs.z());
+    res.z() = (rhs.x() * lhs.y()) - (rhs.y() * lhs.x());
+
+    return res;
+}
+
+
+// two dimensional specialization with a polar conversion
+// TODO: specialize operator+ for polar coordinates
+
+template<typename T>
+class mm::vec2: public mm::basic_vec<T, 2> {
+public:
+    using mm::basic_vec<T, 2>::basic_vec;
+
+    T& x() { return this->at(0); }
+    T& y() { return this->at(1); }
+
+    const T& x() const { return this->at(0); }
+    const T& y() const { return this->at(1); }
+
+    T angle() const;
+    vec2<T> polar() const;
+
+    static vec3<T> cross(const vec2<T>& rhs, const vec2<T>& lhs);
+};
+
+template<typename T>
+T mm::vec2<T>::angle() const {
+    return std::atan(this->y() / this->x());
+}
+
+template<typename T>
+mm::vec2<T> mm::vec2<T>::polar() const {
+    return {
+        this->length(),
+        this->angle()
+    };
+}
+
+template<typename T>
+mm::vec3<T> mm::vec2<T>::cross(const mm::vec2<T>& rhs, const mm::vec2<T>& lhs) {
+    return mm::vec3<T>::cross(mm::vec3<T>(rhs), mm::vec3<T>(lhs));
+}