3 files changed, 522 insertions, 0 deletions
diff --git a/include/mm/experiments/iterators.hpp b/include/mm/experiments/iterators.hpp
index 4afd58f..fb5f568 100644
--- a/include/mm/experiments/iterators.hpp
+++ b/include/mm/experiments/iterators.hpp
@@ -619,4 +619,184 @@ T& mm::const_diag_iterator<T, N>::operator*() const
         M.data[index * Cols + (index - position)];
 }
 
+/*
+ * SECOND IMPLEMENTATION
+ *
+ */
+
+// TODO, short term solution
+#define MM_ROW_ITER 0
+#define MM_COL_ITER 1
+#define MM_DIAG_ITER 2
+
+template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
+class mm::vector_iterator
+{
+    std::size_t index; // variable index
+
+    Grid& M;
+
+    const int position; // fixed index, negative too for diagonal iterator
+
+public:
+
+    vector_iterator(Grid& M, int position, std::size_t index = 0);
+
+    operator T&()
+    {
+        return *(*this);
+    }
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid> operator++()
+    {
+        vector_iterator<T, Rows, Cols, IterType, Grid> it = *this;
+        ++index;
+        return it;
+    }
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid> operator--()
+    {
+        vector_iterator<T, Rows, Cols, IterType, Grid> it = *this;
+        --index;
+        return it;
+    }
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid>& operator++(int)
+    {
+        ++index;
+        return *this;
+    }
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid>& operator--(int)
+    {
+        --index;
+        return *this;
+    }
+
+    bool operator==(const mm::vector_iterator<T, Rows, Cols, IterType, Grid>& other) const
+    {
+        return index == other.index;
+    }
+
+    bool operator!=(const mm::vector_iterator<T, Rows, Cols, IterType, Grid>& other) const
+    {
+        return index != other.index;
+    }
+
+    bool ok() const
+    {
+        if constexpr(IterType == MM_ROW_ITER)
+            return index < Cols;
+        else 
+            return index < Rows;
+    }
+
+    T& operator*();
+    T& operator[](std::size_t);
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid> begin()
+    {
+        return mm::vector_iterator<T, Rows, Cols, IterType, Grid>(M, position, 0);
+    }
+
+    mm::vector_iterator<T, Rows, Cols, IterType, Grid> end()
+    {
+        if constexpr(IterType == MM_ROW_ITER)
+            return mm::vector_iterator<T, Rows, Cols, IterType, Grid>(M, position, Cols);
+        else 
+            return mm::vector_iterator<T, Rows, Cols, IterType, Grid>(M, position, Rows);
+    }
+
+    /*
+     * Scalar product 
+     */
+
+    template<std::size_t P>
+    T operator*(const mm::vector_iterator<T, Rows, P, IterType, Grid>& v)
+    {
+        T out(0);
+
+        for (unsigned k(0); k < Rows; ++k)
+            out += (*this)[k] * v[k];
+
+        return out;
+    }
+
+    template<std::size_t P>
+    T operator*(const mm::vector_iterator<T, P, Cols, IterType, Grid>& v)
+    {
+        T out(0);
+
+        for (unsigned k(0); k < Cols; ++k)
+            out += (*this)[k] * v[k];
+
+        return out;
+    }
+};
+
+
+/* Row Iterators */
+
+namespace mm {
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    using row_iterator = vector_iterator<T, Rows, Cols, MM_ROW_ITER, mm::basic_matrix<T, Rows, Cols>>;
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    using col_iterator = vector_iterator<T, Rows, Cols, MM_COL_ITER, mm::basic_matrix<T, Rows, Cols>>;
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    using const_row_iterator = vector_iterator<typename std::add_const<T>::type, Rows, Cols, MM_ROW_ITER, typename std::add_const<mm::basic_matrix<T, Rows, Cols>>::type>;
+
+    template<typename T, std::size_t Rows, std::size_t Cols>
+    using const_col_iterator = vector_iterator<typename std::add_const<T>::type, Rows, Cols, MM_COL_ITER, typename std::add_const<mm::basic_matrix<T, Rows, Cols>>::type>;
+
+    template<typename T, std::size_t N>
+    using diag_iterator = vector_iterator<T, N, N, MM_DIAG_ITER, mm::basic_matrix<T, N, N>>;
+
+    template<typename T, std::size_t N>
+    using const_diag_iterator = vector_iterator<typename std::add_const<T>::type, N, N, MM_DIAG_ITER, typename std::add_const<mm::basic_matrix<T, N, N>>::type>;
+}
+
+
+/* Iterators implementation */
+
+template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
+mm::vector_iterator<T, Rows, Cols, IterType, Grid>::vector_iterator(Grid& _M, int pos, std::size_t i)
+    : index(i), M(_M), position(pos)
+{
+    if constexpr (IterType == MM_ROW_ITER) {
+        assert(pos < Rows);
+    } else if constexpr (IterType == MM_COL_ITER) {
+        assert(pos < Cols);
+    } else if constexpr (IterType == MM_DIAG_ITER) {
+        assert(abs(pos) < Rows);
+    }
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
+T& mm::vector_iterator<T, Rows, Cols, IterType, Grid>::operator*()
+{
+    if constexpr (IterType == MM_ROW_ITER)
+        return M.data[position * Cols + index];
+    else if constexpr (IterType == MM_COL_ITER)
+        return M.data[index * Cols + position];
+    else if constexpr (IterType == MM_DIAG_ITER)
+        return (position > 0) ?
+            M.data[(index + position) * Cols + index] :
+            M.data[index * Cols + (index - position)];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols, int IterType, class Grid>
+T& mm::vector_iterator<T, Rows, Cols, IterType, Grid>::operator[](std::size_t i)
+{
+    if constexpr (IterType == MM_ROW_ITER)
+        return M.data[position * Cols + i];
+    else if constexpr (IterType == MM_COL_ITER)
+        return M.data[i * Cols + position];
+    else if constexpr (IterType == MM_DIAG_ITER)
+        return (position > 0) ?
+            M.data[(i + position) * Cols + i] :
+            M.data[i * Cols + (i - position)];
+}
 
diff --git a/include/mm/experiments/old.hpp b/include/mm/experiments/old.hpp
new file mode 100644
index 0000000..939f335
--- /dev/null
+++ b/include/mm/experiments/old.hpp
@@ -0,0 +1,58 @@
+#pragma once
+
+
+/* 
+
+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)) {}*/
+
+
+/* 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 mm::row_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+
+        //return row_vec<T, Rows>(
+         //   data.cbegin() + (index * Cols),
+          //  data.cbegin() + ((index + 1) * Cols) + 1
+        );
+    }
+}*/
+
+/*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;
+}*/
diff --git a/include/mm/experiments/tranpositions.hpp b/include/mm/experiments/tranpositions.hpp
new file mode 100644
index 0000000..f124a5c
--- /dev/null
+++ b/include/mm/experiments/tranpositions.hpp
@@ -0,0 +1,284 @@
+#pragma once
+
+template<typename T, std::size_t Rows, std::size_t Cols, bool Regular>
+class mm::access : virtual public mm::basic_matrix<T, Rows, Cols>
+{
+public:
+    template<typename U, std::size_t ORows, std::size_t OCols>
+    friend class mm::access;
+
+    using mm::basic_matrix<T, Rows, Cols>::basic_matrix;
+
+    access(mm::basic_matrix<T, Rows, Cols>&& m) 
+        : basic_matrix<T, Rows, Cols>(m) {}
+
+    virtual T& at(std::size_t row, std::size_t col) = 0;
+    virtual const T& at(std::size_t row, std::size_t col) const = 0;
+
+    static constexpr std::size_t rows = Regular ? Rows : Cols;
+    static constexpr std::size_t cols = Regular ? Cols : Rows;
+};
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::row_access : public mm::access<T, Rows, Cols, true>
+{
+public:
+    template<typename U, std::size_t ORows, std::size_t OCols>
+    friend class mm::row_access;
+
+    using mm::access<T, Rows, Cols, true>::access;
+
+    row_access(mm::basic_matrix<T, Rows, Cols>&& m) 
+        : access<T, Rows, Cols, true>(m) {}
+
+    virtual T& at(std::size_t row, std::size_t col) override;
+    virtual const T& at(std::size_t row, std::size_t col) const override;
+
+    row_iterator<T, Rows, Cols> operator[](std::size_t index);
+    const_row_iterator<T, Rows, Cols> operator[](std::size_t index) const;
+};
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::col_access : public mm::access<T, Rows, Cols, false>
+{
+public:
+    template<typename U, std::size_t ORows, std::size_t OCols>
+    friend class mm::access;
+
+    using mm::access<T, Rows, Cols, false>::access;
+
+    col_access(mm::basic_matrix<T, Rows, Cols>&& m) 
+        : access<T, Rows, Cols, false>(m) {}
+
+    virtual T& at(std::size_t row, std::size_t col) override;
+    virtual const T& at(std::size_t row, std::size_t col) const override;
+
+    col_iterator<T, Rows, Cols> operator[](std::size_t index);
+    const_col_iterator<T, Rows, Cols> operator[](std::size_t index) const;
+};
+
+/*
+ * Square interface
+ */
+
+template<typename T, std::size_t N>
+class mm::square_interface : virtual public mm::basic_matrix<T, N, N>
+{
+public:
+
+    template<typename U, std::size_t ON>
+    friend class mm::square_interface;
+
+    using mm::basic_matrix<T, N, N>::basic_matrix;
+
+    T trace();
+    inline T tr() { return trace(); }
+
+    mm::diag_iterator<T, N> diagonal_iter(int index = 0)
+    {
+        return mm::diag_iterator<T, N>(*this, index);
+    }
+
+    mm::diag_iterator<T, N> diagonal_iter(int index = 0) const
+    {
+        return mm::const_diag_iterator<T, N>(*this, index);
+    }
+
+    // TODO, determinant
+
+    /// in place inverse
+    // TODO, det != 0
+    // TODO, use gauss jordan for invertible ones
+    //void invert();, TODO, section algorithm
+};
+
+/* 
+ * derivated classes 
+ */
+
+// simple format matrix
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::matrix : public mm::row_access<T, Rows, Cols>
+{
+public:
+    using mm::row_access<T, Rows, Cols>::row_access;
+
+    matrix(mm::t_matrix<T, Rows, Cols>&& m)
+        : row_access<T, Rows, Cols>(std::move<mm::basic_matrix<T, Rows, Cols>>(m))
+    {
+                
+    }
+
+    operator mm::t_matrix<T, Rows, Cols>()
+    {
+        return mm::t_matrix<T, Rows, Cols>(*this); 
+    }
+
+    mm::t_matrix<T, Rows, Cols>& t()
+    {
+        return static_cast<mm::t_matrix<T, Rows, Cols>>(*this);
+    }
+
+    const mm::t_matrix<T, Rows, Cols>& t() const
+    {
+        return static_cast<mm::t_matrix<T, Rows, Cols>>(*this);
+    }
+};
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::t_matrix : public mm::col_access<T, Rows, Cols>
+{
+public:
+    using mm::col_access<T, Rows, Cols>::col_access;
+
+    t_matrix(mm::matrix<T, Rows, Cols>&& m)
+        : col_access<T, Rows, Cols>(std::move<mm::basic_matrix<T, Rows, Cols>>(m))
+    {
+                
+    }
+
+    operator mm::matrix<T, Rows, Cols>()
+    {
+        return mm::matrix<T, Rows, Cols>(*this); 
+    }
+
+    mm::matrix<T, Rows, Cols>& t()
+    {
+        return static_cast<mm::matrix<T, Rows, Cols>>(*this);
+    }
+
+    const mm::matrix<T, Rows, Cols>& t() const
+    {
+        return static_cast<mm::matrix<T, Rows, Cols>>(*this);
+    }
+};
+
+// transposed matrix
+/*template<typename T, std::size_t Rows, std::size_t Cols>
+class mm::t_matrix : public mm::basic_matrix<T, Rows, Cols>, virtual public mm::access<T, Rows, Cols, false>
+{
+public:
+    using mm::basic_matrix<T, Rows, Cols>::basic_matrix;
+};*/
+
+/* row vector specialization */
+template<typename T, std::size_t Rows>
+class mm::row_vec : public mm::matrix<T, Rows, 1> {
+public:
+    using mm::matrix<T, Rows, 1>::matrix;
+
+    // TODO, begin, end
+};
+
+/* column vector specialization */
+template<typename T, std::size_t Cols>
+class mm::col_vec : public mm::matrix<T, 1, Cols> {
+public:
+    using mm::matrix<T, 1, Cols>::matrix;
+
+    // TODO, begin, end
+};
+
+/* square matrix specialization */
+template<typename T, std::size_t N>
+class mm::square_matrix : public mm::matrix<T, N, N>, public mm::square_interface<T, N> {
+public:
+    using mm::matrix<T, N, N>::matrix;
+
+    // get the identity of size N
+    static inline constexpr square_matrix<T, N> identity() {
+        mm::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>
+class mm::t_square_matrix : public mm::t_matrix<T, N, N>, public mm::square_interface<T, N> {
+public:
+    using mm::t_matrix<T, N, N>::t_matrix;
+};
+
+/* 
+ * K-diagonal square matrix format 
+ * K is bounded between ]-N, N[
+ */
+
+/*template<typename T, std::size_t N, std::size_t K>
+class mm::diagonal_matrix : public mm::square_matrix<T, N>
+{
+public:
+    using mm::square_matrix<T, N>::square_matrix;
+
+    // TODO, redefine at, operator[]
+    // TODO, matrix multiplication
+};*/
+
+/*
+ * Accessors implementation
+ */
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+T& mm::row_access<T, Rows, Cols>::at(std::size_t row, std::size_t col)
+{
+    return mm::basic_matrix<T, Rows, Cols>::data[row * Cols + col];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+const T& mm::row_access<T, Rows, Cols>::at(std::size_t row, std::size_t col) const
+{
+    return mm::basic_matrix<T, Rows, Cols>::data[row * Cols + col];
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+T& mm::col_access<T, Rows, Cols>::at(std::size_t row, std::size_t col)
+{
+    return mm::basic_matrix<T, Rows, Cols>::data[col * Cols + row]; // transpose
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+const T& mm::col_access<T, Rows, Cols>::at(std::size_t row, std::size_t col) const
+{
+    return mm::basic_matrix<T, Rows, Cols>::data[col * Cols + row]; // transpose
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::row_iterator<T, Rows, Cols> mm::row_access<T, Rows, Cols>::operator[](std::size_t index)
+{
+    return mm::row_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::const_row_iterator<T, Rows, Cols> mm::row_access<T, Rows, Cols>::operator[](std::size_t index) const
+{
+    return mm::const_row_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::col_iterator<T, Rows, Cols> mm::col_access<T, Rows, Cols>::operator[](std::size_t index)
+{
+    return mm::col_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+}
+
+template<typename T, std::size_t Rows, std::size_t Cols>
+mm::const_col_iterator<T, Rows, Cols> mm::col_access<T, Rows, Cols>::operator[](std::size_t index) const
+{
+    return mm::const_col_iterator<T, Rows, Cols>(*this, static_cast<int>(index));
+}
+
+/* Square interface implementation */
+
+template<typename T, std::size_t N>
+T mm::square_interface<T, N>::trace() 
+{
+    T sum = 0;
+    for (const auto& x : diagonal_iter())
+        sum += x;
+
+    return sum;
+}
+
+