Matrix.h

/*
 * Gamedev Framework (gf)
 * Copyright (C) 2016-2017 Julien Bernard
 *
 * This software is provided 'as-is', without any express or implied
 * warranty.  In no event will the authors be held liable for any damages
 * arising from the use of this software.
 *
 * Permission is granted to anyone to use this software for any purpose,
 * including commercial applications, and to alter it and redistribute it
 * freely, subject to the following restrictions:
 *
 * 1. The origin of this software must not be misrepresented; you must not
 *    claim that you wrote the original software. If you use this software
 *    in a product, an acknowledgment in the product documentation would be
 *    appreciated but is not required.
 * 2. Altered source versions must be plainly marked as such, and must not be
 *    misrepresented as being the original software.
 * 3. This notice may not be removed or altered from any source distribution.
 */
#ifndef GF_MATRIX_H
#define GF_MATRIX_H

#include <cstddef>
#include <algorithm>

#include "Portability.h"
#include "Vector.h"

namespace gf {
#ifndef DOXYGEN_SHOULD_SKIP_THIS
inline namespace v1 {
#endif

  template<typename T, std::size_t ROWS, std::size_t COLS>
  struct Matrix {
    #ifndef DOXYGEN_SHOULD_SKIP_THIS
    static_assert(ROWS > 0, "ROWS must be strictly positive");
    static_assert(COLS > 0, "COLS must be strictly positive");
    #endif

    using value_type = T;

    static constexpr std::size_t Rows = ROWS;

    static constexpr std::size_t Cols = COLS;

    Matrix() = default;

    explicit Matrix(T val) {
      std::fill_n(data, ROWS * COLS, val);
    }

    union {
      T grid[ROWS][COLS]; 
      T data[ROWS * COLS]; 
    };
  };

  template<typename T>
  struct Matrix<T, 2, 2> {
    using value_type = T;

    static constexpr std::size_t Rows = 2;

    static constexpr std::size_t Cols = 2;

    Matrix() = default;

    constexpr Matrix(T xx, T xy, T yx, T yy)
    : data{xx, xy, yx, yy}
    {

    }

    union {
      T grid[2][2]; 
      T data[4]; 
      struct {
        T xx, xy;
        T yx, yy;
      };
    };
  };

  template<typename T>
  struct Matrix<T, 3, 3> {
    using value_type = T;

    static constexpr std::size_t Rows = 3;

    static constexpr std::size_t Cols = 3;

    Matrix() = default;

    constexpr Matrix(T xx, T xy, T xz, T yx, T yy, T yz, T zx, T zy, T zz)
    : data{xx, xy, xz, yx, yy, yz, zx, zy, zz}
    {

    }

    union {
      T grid[3][3]; 
      T data[9]; 
      struct {
        T xx, xy, xz;
        T yx, yy, yz;
        T zx, zy, zz;
      };
    };
  };

  template<typename T>
  struct Matrix<T, 4, 4> {
    using value_type = T;

    static constexpr std::size_t Rows = 4;

    static constexpr std::size_t Cols = 4;

    constexpr Matrix(T xx, T xy, T xz, T xw, T yx, T yy, T yz, T yw, T zx, T zy, T zz, T zw, T wx, T wy, T wz, T ww)
    : data{xx, xy, xz, xw, yx, yy, yz, yw, zx, zy, zz, zw, wx, wy, wz, ww}
    {

    }

    union {
      T grid[4][4]; 
      T data[16]; 
      struct {
        T xx, xy, xz, xw;
        T yx, yy, yz, yw;
        T zx, zy, zz, zw;
        T wx, wy, wz, ww;
      };
    };
  };

  using Matrix2f = Matrix<float, 2, 2>;

  using Matrix3f = Matrix<float, 3, 3>;

  using Matrix4f = Matrix<float, 4, 4>;

// MSVC does not like extern template
#ifndef _MSC_VER
  extern template struct Matrix<float, 2, 2>;
  extern template struct Matrix<float, 3, 3>;
  extern template struct Matrix<float, 4, 4>;
#endif

  /*
   * equality operations
   */

  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  bool operator==(const Matrix<T, ROWS, COLS>& lhs, const Matrix<T, ROWS, COLS>& rhs) {
    for (std::size_t i = 0; i < ROWS; ++i) {
      for (std::size_t j = 0; j < COLS; ++j) {
        if (lhs.grid[i][j] != rhs.grid[i][j]) {
          return false;
        }
      }
    }

    return true;
  }

  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  bool operator!=(const Matrix<T, ROWS, COLS>& lhs, const Matrix<T, ROWS, COLS>& rhs) {
    return !(lhs == rhs);
  }

  /*
   * unary operators
   */

  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  Matrix<T, ROWS, COLS> operator-(const Matrix<T, ROWS, COLS>& val) {
    Matrix<T, ROWS, COLS> out;

    for (std::size_t i = 0; i < ROWS; ++i) {
      for (std::size_t j = 0; j < COLS; ++j) {
        out.grid[i][j] = - val.grid[i][j];
      }
    }

    return out;
  }

  /*
   * arithmetic component-wise operations
   */

  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  Matrix<T, ROWS, COLS> operator+(const Matrix<T, ROWS, COLS>& lhs, const Matrix<T, ROWS, COLS>& rhs) {
    Matrix<T, ROWS, COLS> out;

    for (std::size_t i = 0; i < ROWS; ++i) {
      for (std::size_t j = 0; j < COLS; ++j) {
        out.grid[i][j] = lhs.grid[i][j] + rhs.grid[i][j];
      }
    }

    return out;
  }

  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  Matrix<T, ROWS, COLS>& operator+=(Matrix<T, ROWS, COLS>& lhs, const Matrix<T, ROWS, COLS>& rhs) {
    for (std::size_t i = 0; i < ROWS; ++i) {
      for (std::size_t j = 0; j < COLS; ++j) {
        lhs.grid[i][j] += rhs.grid[i][j];
      }
    }

    return lhs;
  }


  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  Matrix<T, ROWS, COLS> operator-(const Matrix<T, ROWS, COLS>& lhs, const Matrix<T, ROWS, COLS>& rhs) {
    Matrix<T, ROWS, COLS> out;

    for (std::size_t i = 0; i < ROWS; ++i) {
      for (std::size_t j = 0; j < COLS; ++j) {
        out.grid[i][j] = lhs.grid[i][j] - rhs.grid[i][j];
      }
    }

    return out;
  }

  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  Matrix<T, ROWS, COLS>& operator-=(Matrix<T, ROWS, COLS>& lhs, const Matrix<T, ROWS, COLS>& rhs) {
    for (std::size_t i = 0; i < ROWS; ++i) {
      for (std::size_t j = 0; j < COLS; ++j) {
        lhs.grid[i][j] -= rhs.grid[i][j];
      }
    }

    return lhs;
  }

  /*
   * multiplication and division
   */

  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  Matrix<T, ROWS, COLS> operator*(const Matrix<T, ROWS, COLS>& lhs, T rhs) {
    Matrix<T, ROWS, COLS> out;

    for (std::size_t i = 0; i < ROWS; ++i) {
      for (std::size_t j = 0; j < COLS; ++j) {
        out.grid[i][j] = lhs.grid[i][j] * rhs;
      }
    }

    return out;
  }

  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  Matrix<T, ROWS, COLS>& operator*=(Matrix<T, ROWS, COLS>& lhs, T rhs) {
    for (std::size_t i = 0; i < ROWS; ++i) {
      for (std::size_t j = 0; j < COLS; ++j) {
        lhs.grid[i][j] *= rhs;
      }
    }

    return lhs;
  }

  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  Matrix<T, ROWS, COLS> operator*(T lhs, const Matrix<T, ROWS, COLS>& rhs) {
    Matrix<T, ROWS, COLS> out;

    for (std::size_t i = 0; i < ROWS; ++i) {
      for (std::size_t j = 0; j < COLS; ++j) {
        out.grid[i][j] = lhs * rhs.grid[i][j];
      }
    }

    return out;
  }

  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  Matrix<T, ROWS, COLS> operator/(const Matrix<T, ROWS, COLS>& lhs, T rhs) {
    Matrix<T, ROWS, COLS> out;

    for (std::size_t i = 0; i < ROWS; ++i) {
      for (std::size_t j = 0; j < COLS; ++j) {
        out.grid[i][j] = lhs.grid[i][j] / rhs;
      }
    }

    return out;
  }

  template<typename T, std::size_t ROWS, std::size_t COLS>
  inline
  Matrix<T, ROWS, COLS>& operator/=(Matrix<T, ROWS, COLS>& lhs, T rhs) {
    for (std::size_t i = 0; i < ROWS; ++i) {
      for (std::size_t j = 0; j < COLS; ++j) {
        lhs.grid[i][j] /= rhs;
      }
    }

    return lhs;
  }

  template<typename T, std::size_t S1, std::size_t S2>
  inline
  Vector<T, S1> operator*(const Matrix<T, S1, S2>& lhs, const Vector<T, S2>& rhs) {
    Vector<T, S1> out;

    for (std::size_t i = 0; i < S1; ++i) {
      T val{0};

      for (std::size_t j = 0; j < S2; ++j) {
        val += lhs.grid[i][j] * rhs.data[j];
      }

      out.data[i] = val;
    }

    return out;
  }

  template<typename T, std::size_t S1, std::size_t S2>
  inline
  Vector<T, S2> operator*(const Vector<T, S1>& lhs, const Matrix<T, S1, S2>& rhs) {
    Vector<T, S2> out;

    for (std::size_t j = 0; j < S2; ++j) {
      T val{0};

      for (std::size_t i = 0; i < S1; ++i) {
        val += lhs.data[i] * rhs.grid[i][j];
      }

      out.data[j] = val;
    }

    return out;
  }


  template<typename T, std::size_t S1, std::size_t S2, std::size_t S3>
  inline
  Matrix<T, S1, S3> operator*(const Matrix<T, S1, S2>& lhs, const Matrix<T, S2, S3>& rhs) {
    Matrix<T, S1, S3> out;

    for (std::size_t i = 0; i < S1; ++i) {
      for (std::size_t j = 0; j < S3; ++j) {
        T val{0};

        for (std::size_t k = 0; k < S2; ++k) {
          val += lhs.grid[i][k] * rhs.grid[k][j];
        }

        out.grid[i][j] = val;
      }
    }

    return out;
  }

  template<typename T, std::size_t N>
  inline
  Matrix<T, N, N>& operator*=(Matrix<T, N, N>& lhs, const Matrix<T, N, N>& rhs) {
    lhs = lhs * rhs;
    return lhs;
  }

  /*
   * usual operations
   */

  template<typename MatrixType>
  inline
  MatrixType identity() {
    static_assert(MatrixType::Rows == MatrixType::Cols, "identity() is only defined for square matrices.");

    MatrixType out;

    typedef typename MatrixType::value_type value_type;

    for (std::size_t i = 0; i < MatrixType::Rows; ++i) {
      for (std::size_t j = 0; j < MatrixType::Cols; ++j) {
        out.grid[i][j] = (i == j) ? value_type{1} : value_type{0};
      }
    }

    return out;
  }

  template<typename T, std::size_t S1, std::size_t S2>
  inline
  Matrix<T, S2, S1> transpose(const Matrix<T, S1, S2>& mat) {
    Matrix<T, S2, S1> out;

    for (std::size_t i = 0; i < S1; ++i) {
      for (std::size_t j = 0; j < S2; ++j) {
        out.grid[j][i] = mat.grid[i][j];
      }
    }

    return out;
  }

  // https://en.wikipedia.org/wiki/Invertible_matrix

  template<typename T>
  inline
  Matrix<T, 2, 2> invert(const Matrix<T, 2, 2>& mat) {
    Matrix<T, 2, 2> out;

    out.xx = mat.yy;
    out.xy = - mat.xy;
    out.yx = - mat.yx;
    out.yy = mat.xx;

    T det = mat.xx * mat.yy - mat.yx * mat.xy;
    out /= det;
    return out;
  }

  template<typename T>
  inline
  Matrix<T, 3, 3> invert(const Matrix<T, 3, 3>& mat) {
    Matrix<T, 3, 3> out;

    out.xx = mat.yy * mat.zz - mat.zy * mat.yz;
    out.xy = - (mat.xy * mat.zz - mat.zy * mat.xz);
    out.xz = mat.xy * mat.yz - mat.yy * mat.xz;
    out.yx = - (mat.yx * mat.zz - mat.zx * mat.yz);
    out.yy = mat.xx * mat.zz - mat.zx * mat.xz;
    out.yz = - (mat.xx * mat.yz - mat.yx * mat.xz);
    out.zx = mat.yx * mat.zy - mat.zx * mat.yy;
    out.zy = - (mat.xx * mat.zy - mat.zx * mat.xy);
    out.zz = mat.xx * mat.yy - mat.yx * mat.xy;

    T det = mat.xx * out.xx + mat.xy * out.yx + mat.xz * out.zx;
    out /= det;
    return out;
  }

#ifndef DOXYGEN_SHOULD_SKIP_THIS
}
#endif
}

#endif // GAME_MATRIX_H