1 files changed, 473 insertions, 63 deletions

diff --git a/meowpp/math/Matrix.h b/meowpp/math/Matrix.h
index 0a52d8c..b4b4853 100644
--- a/meowpp/math/Matrix.h
+++ b/meowpp/math/Matrix.h
@@ -1,74 +1,484 @@
 #ifndef   math_Matrix_H__
 #define   math_Matrix_H__
 
+#include "utility.h"
 
-#include <cstdlib>
+#include "../Self.h"
 
 #include <vector>
+#include <algorithm>
+
+#include <cstdlib>
 
-namespace meow{
-  template<class Entry>
-  class Matrix{
-    private:
-      size_t             _rows;
-      size_t             _cols;
-      std::vector<Entry> _entries;
-      //
-      size_t index(size_t __r, size_t __c) const;
-    public:
-      Matrix();
-      Matrix(Matrix const& __m);
-      Matrix(size_t __rows, size_t __cols, Entry const& __entry);
-      ~Matrix();
-      
-      Matrix& copy(Matrix const& __m);
-      
-      bool valid() const;
-      
-      size_t rows() const;
-      size_t cols() const;
-      size_t size(size_t __rows, size_t __cols);
-      
-      Entry const& entry(size_t __i, size_t __j) const;
-      Entry const& entry(size_t __i, size_t __j, Entry const& __entry);
-      void entries(size_t __rFirst, size_t __rLast,
-                   size_t __cFirst, size_t __cLast,
-                   Entry const& __entry);
-      
-      Matrix subMatrix(size_t __rFirst, size_t __rLast,
-                       size_t __cFirst, size_t __cLast) const;
-      Matrix row(size_t __row) const;
-      Matrix col(size_t __col) const;
-      
-      Matrix positive() const;
-      Matrix negative() const;
-      
-      Matrix add(Matrix const& __m) const;
-      Matrix sub(Matrix const& __m) const;
-      Matrix mul(Matrix const& __m) const;
-      
-      Matrix mul(Entry const& __s) const;
-      Matrix div(Entry const& __s) const;
-      
-      Matrix  identity  () const;
-      Matrix& identitied();
-      
-      Matrix  transpose () const;
-      Matrix& transposed();
-      
-      Matrix& operator=(Matrix const& __m);
-      Entry const& operator()(size_t __i, size_t __j) const;
-      Entry      & operator()(size_t __i, size_t __j);
-      Matrix operator+() const;
-      Matrix operator-() const;
-      Matrix operator+(Matrix const& __m) const;
-      Matrix operator-(Matrix const& __m) const;
-      Matrix operator*(Matrix const& __m) const;
-      Matrix operator*(Entry const& __s) const;
-      Matrix operator/(Entry const& __s) const;
+namespace meow {
+/*!
+ * @brief \b matrix
+ *
+ * @author cat_leopard
+ */
+template<class Entry>
+class Matrix {
+private:
+  struct Myself {
+    size_t                rows_;
+    size_t                cols_;
+    std::vector<Entry> entries_;
+    Myself(): rows_(0), cols_(0), entries_(0) {
+    }
+    ~Myself() {
+    }
+    size_t index(size_t r, size_t c) const {
+      return r * cols_ + c;
+    }
+    Myself& copyFrom(Myself const& m) {
+      rows_    = m.   rows_;
+      cols_    = m.   cols_;
+      entries_ = m.entries_;
+      return *this;
+    }
   };
-}
+  Self<Myself> const self;
+public:
+  /*!
+   * @brief constructor
+   *
+   * Create an empty matrix with size \b 0x0.
+   * In other world, create an \b invalid matrix
+   */
+  Matrix(): self(true) { }
+
+  /*!
+   * @brief constructor
+   *
+   * Copy data from another one
+   *
+   * @param [in] m another matrix
+   */
+  Matrix(Matrix const& m): self(false) { self().copyFrom(m.self); }
+
+  /*!
+   * @brief constructor
+   *
+   * Create an \a r x \a c matrix with all entry be \a e
+   *
+   * @param [in] r number of rows
+   * @param [in] c number of columns
+   * @param [in] e inital entry
+   */
+  Matrix(size_t r, size_t c, Entry const& e): self(true) { reset(r, c, e); }
+
+  //! @brief destructor
+  ~Matrix() { }
+
+  /*!
+   * @brief copy
+   *
+   * Copy data from another matrix
+   *
+   * @param [in] m matrix
+   * @return *this
+   */
+  Matrix& copyFrom(Matrix const& m) {
+    self().copyFrom(m.self);
+    return *this;
+  }
+
+  /*!
+   * @brief reference
+   *
+   * Reference itself to another matrix
+   *
+   * @param [in] m matrix
+   * @return *this
+   */
+  Matrix& referenceFrom(Matrix const& m) {
+    self().referenceFrom(m.self);
+    return *this;
+  }
+
+  //! @brief reset the size of the matrix to \a r x \a c with entry all be \a e
+  void reset(size_t r, size_t c, Entry const& e) {
+    self()->rows_ = r;
+    self()->cols_ = c;
+    self()->entries_.clear();
+    self()->entries_.resize(r * c, e);
+  }
+
+  //! @brief Return whether it is a \b valid matrix
+  bool valid() const {
+    return (rows() > 0 && cols() > 0);
+  }
+
+  //! @brief Return number of rows
+  size_t rows() const {
+    return self->rows_;
+  }
+
+  //! @brief Return number of cols
+  size_t cols() const {
+    return self->cols_;
+  }
+
+  //! @brief Return number of rows times number of cols
+  size_t size() const {
+    return rows() * cols();
+  }
+
+  /*!
+   * @brief resize the matrix such that number of rows become \a r.
+   *
+   * New created entry will be \a e
+   *
+   * @param [in] r new number of rows
+   * @param [in] e inital entry
+   * @return new number of rows
+   */
+  size_t rows(size_t r, Entry const& e) {
+    if (r != rows()) {
+      self()->entries_.resize(r * cols(), e);
+      self()->rows_ = r;
+    }
+    return rows();
+  }
+
+  /*!
+   * @brief resize the matrix such that number of cols become \a c
+   *
+   * New created entry will be \a e
+   *
+   * @param [in] c new number of columns
+   * @param [in] e inital entry
+   * @return new number of columns
+   */
+  size_t cols(size_t c, Entry const& e) {
+    if (c != cols()) {
+      Self<Myself> const old(false);
+      old().copyFrom(self);
+      self()->entries_.resize(rows() * c);
+      self()->cols_ = c;
+      for (size_t i = 0, I = rows(); i < I; i++) {
+        size_t j, J1 = std::min(old->cols_, cols()), J2 = cols();
+        for (j =  0; j < J1; j++)
+          self()->entries_[self->index(i, j)] = old->entries_[old->index(i, j)];
+        for (j = J1; j < J2; j++)
+          self()->entries_[self->index(i, j)] = e;
+      }
+    }
+    return cols();
+  }
+
+  /*!
+   * @brief resize
+   * 
+   * Resize to \a r x \a c, with new created entry be \a e
+   *
+   * @param [in] r number of rows
+   * @param [in] c number of rows
+   * @param [in] e inital entry
+   * @return \a r * \a c
+   */
+  size_t size(size_t r, size_t c, Entry const& e) {
+    cols(c, e);
+    rows(r, e);
+    return rows() * cols();
+  }
+
+  //! @brief Access the entry at \a r x \a c
+  Entry entry(size_t r, size_t c) const {
+    return self->entries_[self->index(r, c)];
+  }
+
+  //! @brief Change the entry at \a r x \a c
+  Entry entry(size_t r, size_t c, Entry const& e) {
+    self()->entries_[self->index(r, c)] = e;
+    return entry(r, c);
+  }
 
-#include "Matrix.hpp"
+  /*!
+   * @brief Change the entries from \a rFirst x \a cFirst to \a rLast x \a cLast
+   * 
+   * @param [in] rFirst 
+   * @param [in] rLast 
+   * @param [in] cFirst
+   * @param [in] cLast 
+   * @param [in] e value
+   * @return void
+   */
+  void entries(ssize_t rFirst, ssize_t rLast,
+               ssize_t cFirst, ssize_t cLast,
+               Entry const& e) {
+    for (ssize_t r = rFirst; r <= rLast; r++) {
+      for (ssize_t c = cFirst; c <=cFirst; c++) {
+        entry(r, c, e);
+      }
+    }
+  }
+
+  /*!
+   * @brief Return a \a rLast-rFirst+1 x \a cLast-cFirst+1 matrix
+   *
+   * With value be the entries from \a rFirst x \a cFirst to \a rLast x \a cLast
+   * 
+   * @param [in] rFirst
+   * @param [in] rLast
+   * @param [in] cFirst
+   * @param [in] cLast
+   * @return a matrix
+   */
+  Matrix subMatrix(size_t rFirst, size_t rLast,
+                   size_t cFirst, size_t cLast) const {
+    if (rFirst > rLast || cFirst > cLast) return Matrix();
+    if (rFirst == 0 || cFirst == 0) {
+      Matrix ret(*this);
+      ret.size(rLast + 1, cLast + 1, Entry(0));
+      return ret;
+    }
+    Matrix ret(rLast - rFirst + 1, cLast - cFirst + 1, entry(rFirst, cFirst));
+    for (size_t r = rFirst; r <= rLast; r++)
+      for (size_t c = cFirst; c <= cLast; c++)
+        ret.entry(r - rFirst, c - cFirst, entry(r, c));
+    return ret;
+  }
+
+  //! @brief Return the \a r -th row
+  Matrix row(size_t r) const {
+    return subMatrix(r, r, 0, cols() - 1);
+  }
+
+  //! @brief Return the \a c -th column
+  Matrix col(size_t c) const {
+    return subMatrix(0, rows() - 1, c, c);
+  }
+
+  //! @brief return +\a (*this)
+  Matrix positive() const {
+    return *this;
+  }
+
+  //! @brief return -\a (*this)
+  Matrix negative() const {
+    Matrix ret(*this);
+    for (size_t r = 0, R = rows(); r < R; r++)
+      for (size_t c = 0, C = cols(); c < C; c++)
+        ret.entry(r, c, -ret.entry(r, c));
+    return ret;
+  }
+
+  /*! @brief return \a (*this) + \a m.
+   *
+   * If the size not match, it will return an invalid matrix
+   */
+  Matrix add(Matrix const& m) const {
+    if (rows() != m.rows() || cols() != m.cols()) return Matrix();
+    Matrix ret(*this);
+    for (size_t r = 0, R = rows(); r < R; r++)
+      for (size_t c = 0, C = cols(); c < C; c++)
+        ret.entry(r, c, ret.entry(r, c) + m.entry(r, c));
+    return ret;
+  }
+
+  /*! @brief return \a (*this) - \a m.
+   *
+   * If the size not match, it will return an invalid matrix
+   */
+  Matrix sub(Matrix const& m) const {
+    if (rows() != m.rows() || cols() != m.cols()) return Matrix();
+    Matrix ret(*this);
+    for (size_t r = 0, R = rows(); r < R; r++)
+      for (size_t c = 0, C = cols(); c < C; c++)
+        ret.entry(r, c, ret.entry(r, c) - m.entry(r, c));
+    return ret;
+  }
+
+  /*! @brief return \a (*this) times \a m.
+   *
+   * If the size not match, it will return an invalid matrix
+   */
+  Matrix mul(Matrix const& m) const {
+    if (cols() != m.rows()) return Matrix();
+    Matrix ret(rows(), m.cols(), Entry(0));
+    for (size_t r = 0, R = rows(); r < R; r++)
+      for (size_t c = 0, C = m.cols(); c < C; c++)
+        for (size_t k = 0, K = cols(); k < K; k++)
+          ret.entry(r, c, ret.entry(r, c) + entry(r, k) * m.entry(k, c));
+    return ret;
+  }
+
+  //! @brief return \a (*this) times \a s. \a s is a scalar
+  Matrix mul(Entry const& s) const {
+    Matrix ret(*this);
+    for (size_t r = 0, R = rows(); r < R; r++)
+      for (size_t c = 0, C = cols(); c < C; c++)
+        ret.entry(r, c, ret.entry(r, c) * s);
+    return ret;
+  }
+
+  //! @brief return \a (*this) / \a s. \a s is a scalar
+  Matrix div(Entry const& s) const {
+    Matrix ret(*this);
+    for (size_t r = 0, R = rows(); r < R; r++)
+      for (size_t c = 0, C = cols(); c < C; c++)
+        ret.entry(r, c, ret.entry(r, c) / s);
+    return ret;
+  }
+
+  //! @brief Return a identity matrix with size equal to itself
+  Matrix identity() const {
+    Matrix ret(*this);
+    ret.identitied();
+    return ret;
+  }
+
+  /*!
+   * @brief Let itself be an identity matrix
+   *
+   * Our definition of Identity matrix is 1 for entry(i, i) and 0 otherwise.
+   */
+  Matrix& identitied() {
+    for (size_t r = 0, R = rows(); r < R; r++)
+      for (size_t c = 0, C = cols(); c < C; c++)
+        entry(r, c, (r == c ? Entry(1) : Entry(0)));
+    return *this;
+  }
+
+  /*!
+   * @brief Return a matrix which is an inverse matrix of \a (*this)
+   *
+   * If inverse matrix doesn't exist, it will return a invalid matrix
+   */
+  Matrix inverse() const {
+    if (rows() != cols() || rows() == 0) return Matrix<Entry>();
+    Matrix tmp(rows(), cols() * 2, Entry(0));
+    for (size_t r = 0, R = rows(); r < R; r++) {
+      for (size_t c = 0, C = cols(); c < C; c++) {
+        tmp.entry(r, c, entry(r, c));
+        tmp.entry(r, c + cols(), (r == c ? Entry(1) : Entry(0)));
+      }
+    }
+    tmp.triangulared();
+    for (ssize_t r = rows() - 1; r >= 0; r--) {
+      if (tmp(r, r) == Entry(0)) return Matrix<Entry>();
+      for (ssize_t r2 = r - 1; r2 >= 0; r2--) {
+        Entry rat(-tmp.entry(r2, r) / tmp.entry(r, r));
+        for (size_t c = r, C = tmp.cols(); c < C; c++) {
+          tmp.entry(r2, c, tmp.entry(r2, c) + rat * tmp(r, c));
+        }
+      }
+      Entry rat(tmp.entry(r, r));
+      for (size_t c = cols(), C = tmp.cols(); c < C; c++) {
+        tmp.entry(r, c - cols(), tmp.entry(r, c) / rat);
+      }
+    }
+    tmp.size(cols(), rows(), Entry(0));
+    return tmp;
+  }
+  
+  //! @brief let itself become itself's inverse matrix
+  Matrix& inversed() {
+    copyFrom(inverse());
+    return *this;
+  }
+
+  //! @brief return itself's transpose matrix
+  Matrix  transpose () const {
+    Matrix ret(cols(), rows(), Entry(0));
+    for (size_t r = 0, R = cols(); r < R; r++)
+      for (size_t c = 0, C = rows(); c < C; c++)
+        ret.entry(r, c, entry(c, r));
+    return ret;
+  }
+  
+  //! @brief Let itself become itself's transpose matrix
+  Matrix& transposed() {
+    copyFrom(transpose());
+    return *this;
+  }
+
+  //! @brief return a matrix which is the triangular form of \a (*this)
+  Matrix triangular() const {
+    Matrix<Entry> ret(*this);
+    ret.triangulared();
+    return ret;
+  }
+  
+  //! @brief triangluar itself
+  Matrix& triangulared() {
+    for (size_t r = 0, c = 0, R = rows(), C = cols(); r < R && c < C; r++) {
+      ssize_t maxR;
+      for ( ; c < C; c++) {
+        maxR = -1;
+        for (size_t r2 = r; r2 < R; r2++)
+          if (maxR == -1 || tAbs(entry(r2, c)) > tAbs(entry(maxR, c)))
+            maxR = r2;
+        if (entry(maxR, c) != Entry(0)) break;
+      }
+      if (c >= C) break;
+      if (maxR != (ssize_t)r) {
+        for (size_t c2 = c; c2 < C; c2++)
+          std::swap(self()->entries_[self->index(   r, c2)],
+                    self()->entries_[self->index(maxR, c2)]);
+      }
+      for (size_t r2 = r + 1; r2 < R; r2++) {
+        Entry rati = -entry(r2, c) / entry(r, c);
+        entry(r2, c, Entry(0));
+        for (size_t c2 = c + 1; c2 < C; c2++)
+          entry(r2, c2, entry(r2, c2) + entry(r, c2) * rati);
+      }
+    }
+    return *this;
+  }
+
+  //! @brief same as \a copyFrom
+  Matrix& operator=(Matrix const& m) {
+    return copyFrom(m);
+  }
+  
+  //! @brief same as \a entry(r,c)
+  Entry operator()(size_t r, size_t c) const {
+    return entry(r, c);
+  }
+  
+  //! @brief same as \a entry(r,c,e)
+  Entry operator()(size_t r, size_t c, Entry const& e) {
+    return entry(r, c, e);
+  }
+  
+  //! @brief same as \a positive()
+  Matrix operator+() const {
+    return positive();
+  }
+  
+  //! @brief same as \a negative()
+  Matrix operator-() const {
+    return negative();
+  }
+  
+  //! @brief same as \a add(m)
+  Matrix operator+(Matrix const& m) const {
+    return add(m);
+  }
+  
+  //! @brief same as \a sub(m)
+  Matrix operator-(Matrix const& m) const {
+    return sub(m);
+  }
+  
+  //! @brief same as \a mul(m)
+  Matrix operator*(Matrix const& m) const {
+    return mul(m);
+  }
+  
+  //! @brief same as \a mul(m)
+  Matrix operator*(Entry const& s) const {
+    return mul(s);
+  }
+  
+  //! @brief same as \a div(s)
+  Matrix operator/(Entry const& s) const {
+    return div(s);
+  }
+};
+
+}
 
 #endif // math_Matrix_H__