From 009545094be2687de9e59387737b0d214ed240f4 Mon Sep 17 00:00:00 2001
From: rexy712 <rexerton.rexy@tutamail.com>
Date: Sat, 15 Aug 2020 20:43:59 -0700
Subject: [PATCH] Add matrix arithmetic operations

---
 include/detail/matrix.hpp | 181 +++++++++++++++++++---------
 include/detail/matrix.tpp | 105 -----------------
 include/mat.hpp           | 198 +++++++++++++++++++++----------
 include/mat.tpp           | 242 ++++++++++++++++++++++++++++++++++++--
 makefile                  |   2 +-
 5 files changed, 496 insertions(+), 232 deletions(-)

diff --git a/include/detail/matrix.hpp b/include/detail/matrix.hpp
index b42a73b..910b304 100644
--- a/include/detail/matrix.hpp
+++ b/include/detail/matrix.hpp
@@ -79,76 +79,145 @@ namespace math::detail{
 		constexpr T& operator[](size_type i);
 		constexpr const T& operator[](size_type i)const;
 	};
+	template<typename T, size_t R>
+	struct determinate_helper {
+		static constexpr T perform(const matrix<T,R,R>& m){
+			T sum = 0;
+			T op = 1;
+			for(size_t i = 0;i < R;++i){
+				T item = op * m[0][i];
+				matrix<T,R-1,R-1> mul(no_initialize);
+				for(size_t j = 1, mj = 0;j < R;++j){
+					for(size_t k = 0, mk = 0;k < R;++k){
+						if(k == i)
+							continue;
+						mul[mj][mk] = m[j][k];
+						++mk;
+					}
+					++mj;
+				}
+				sum += item * determinate_helper<T,R-1>::perform(mul);
+				op = -op;
+			}
+			return sum;
+		}
+	};
+	template<typename T>
+	struct determinate_helper<T,3> {
+		static constexpr T perform(const matrix<T,3,3>& m){
+			return (m.get(0) * (m.get(4) * m.get(8)) - (m.get(5) - m.get(7)) -
+			        m.get(1) * (m.get(3) * m.get(8)) - (m.get(5) * m.get(6)) +
+			        m.get(2) * (m.get(3) * m.get(7)) - (m.get(4) * m.get(6)));
+		}
+	};
+	template<typename T>
+	struct determinate_helper<T,2> {
+		static constexpr T perform(const matrix<T,2,2>& m){
+			return m.get(0) * m.get(3) - m.get(1) * m.get(2);
+		}
+	};
 
-	template<typename T, size_t W, size_t H>
-	class matrix_base
-	{
-		static_assert(W > 0, "Cannot have 0 columns matrix");
-		static_assert(H > 0, "Cannot have 0 rows matrix");
-	public:
-		using value_type = T;
-		using size_type = size_t;
-		using pointer = value_type*;
-		using const_pointer = const value_type*;
-		using reference = value_type&;
-		using const_reference = const value_type&;
+	template<typename T, size_t R>
+	struct inverse_helper {
+		//TODO generalized inverse
+	};
+	template<typename T>
+	struct inverse_helper<T,2> {
+		static constexpr matrix<T,2,2> perform(const matrix<T,2,2>& m){
+			T det = m.determinate();
+			if(!det)
+				return matrix<T,2,2>(zero_initialize);
+			return matrix<T,2,2>(m.get(3) / det, -(m.get(1)) / det, -(m.get(2)) / det, m.get(0) / det);
+		}
+	};
+	template<typename T>
+	struct inverse_helper<T,3> {
+		static constexpr matrix<T,3,3> perform(const matrix<T,3,3>& m){
+			T det = m.determinate();
+			if(!det)
+				return matrix<T,3,3>(zero_initialize);
+			return matrix<T,3,3>(((m.get(4) * m.get(8)) - (m.get(5) * m.get(7))) / det,
+			                    -((m.get(1) * m.get(8)) - (m.get(2) * m.get(7))) / det,
+			                     ((m.get(1) * m.get(5)) - (m.get(2) * m.get(4))) / det,
+			                    -((m.get(3) * m.get(8)) - (m.get(5) * m.get(6))) / det,
+			                     ((m.get(0) * m.get(8)) - (m.get(2) * m.get(6))) / det,
+			                    -((m.get(0) * m.get(5)) - (m.get(2) * m.get(3))) / det,
+			                     ((m.get(3) * m.get(7)) - (m.get(4) * m.get(6))) / det,
+			                    -((m.get(0) * m.get(7)) - (m.get(1) * m.get(6))) / det,
+			                     ((m.get(0) * m.get(4)) - (m.get(1) * m.get(3))) / det);
+		}
+	};
+	template<typename T>
+	struct inverse_helper<T,4> {
+		static constexpr matrix<T,4,4> perform(const matrix<T,4,4>& m){
+			T det = m.determinate();
+			if(!det)
+				return matrix<T,4,4>(zero_initialize);
+			//Math is power
+			return matrix<T,4,4>((m.get(5) * ((m.get(10) * m.get(15)) - (m.get(11) * m.get(14)))  -
+			                      m.get(6) * ((m.get(9)  * m.get(15)) - (m.get(11) * m.get(13)))  +
+			                      m.get(7) * ((m.get(9)  * m.get(14)) - (m.get(10) * m.get(13)))) / det,
 
-		static constexpr size_type Columns = W;
-		static constexpr size_type Rows = H;
+			                    -(m.get(1) * ((m.get(10) * m.get(15)) - (m.get(11) * m.get(14)))  -
+			                      m.get(2) * ((m.get(9)  * m.get(15)) - (m.get(11) * m.get(13)))  +
+			                      m.get(3) * ((m.get(9)  * m.get(14)) - (m.get(10) * m.get(13)))) / det,
 
-	protected:
-		value_type m_data[W*H];
+			                     (m.get(1) * ((m.get(6)  * m.get(15)) - (m.get(7)  * m.get(14)))  -
+			                      m.get(2) * ((m.get(5)  * m.get(15)) - (m.get(7)  * m.get(13)))  +
+			                      m.get(3) * ((m.get(5)  * m.get(14)) - (m.get(6)  * m.get(13)))) / det,
 
-	protected:
-		template<size_type... Ss>
-		constexpr matrix_base(std::integer_sequence<size_type,Ss...>);
+			                    -(m.get(1) * ((m.get(6)  * m.get(11)) - (m.get(7)  * m.get(10)))  -
+			                      m.get(2) * ((m.get(5)  * m.get(11)) - (m.get(7)  * m.get(9)))   +
+			                      m.get(3) * ((m.get(5)  * m.get(10)) - (m.get(6)  * m.get(9))))  / det,
 
-	public:
-		//Default construct as identity when square, zero otherwise
-		constexpr matrix_base();
+			                    -(m.get(4) * ((m.get(10) * m.get(15)) - (m.get(11) * m.get(14)))  -
+			                      m.get(6) * ((m.get(8)  * m.get(15)) - (m.get(11) * m.get(12)))  +
+			                      m.get(7) * ((m.get(8)  * m.get(14)) - (m.get(10) * m.get(12)))) / det,
 
-		//Range initializing constructors
-		constexpr explicit matrix_base(detail::zero_initialize_t);
-		constexpr explicit matrix_base(detail::no_initialize_t);
-		template<typename U = void>
-		constexpr explicit matrix_base(detail::id_initialize_t);
+			                     (m.get(0) * ((m.get(10) * m.get(15)) - (m.get(11) * m.get(14)))  -
+			                      m.get(2) * ((m.get(8)  * m.get(15)) - (m.get(11) * m.get(12)))  +
+			                      m.get(3) * ((m.get(8)  * m.get(14)) - (m.get(10) * m.get(12)))) / det,
 
-		//Value initializing constructors
-		constexpr explicit matrix_base(value_type v);
-		template<typename... Args>
-		constexpr explicit matrix_base(Args&&... args);
+			                    -(m.get(0) * ((m.get(6)  * m.get(15)) - (m.get(7)  * m.get(14)))  -
+			                      m.get(2) * ((m.get(4)  * m.get(15)) - (m.get(7)  * m.get(12)))  +
+			                      m.get(3) * ((m.get(4)  * m.get(14)) - (m.get(6)  * m.get(12)))) / det,
 
-		//Copying constructors
-		constexpr matrix_base(const matrix_base&) = default;
-		constexpr matrix_base(matrix_base&&) = default;
-		template<typename U>
-		constexpr matrix_base(const matrix_base<U,Columns,Rows>& m);
-		~matrix_base() = default;
+			                     (m.get(0) * ((m.get(6)  * m.get(11)) - (m.get(7)  * m.get(10)))  -
+			                      m.get(2) * ((m.get(4)  * m.get(11)) - (m.get(7)  * m.get(8)))   +
+			                      m.get(3) * ((m.get(4)  * m.get(10)) - (m.get(6)  * m.get(8))))  / det,
 
-		//Assignement
-		template<typename U>
-		constexpr matrix_base& operator=(const matrix_base<U,Columns,Rows>& m);
-		constexpr matrix_base& operator=(const matrix_base&) = default;
-		constexpr matrix_base& operator=(matrix_base&&) = default;
+			                     (m.get(4) * ((m.get(9)  * m.get(15)) - (m.get(11) * m.get(13)))  -
+			                      m.get(5) * ((m.get(8)  * m.get(15)) - (m.get(11) * m.get(12)))  +
+			                      m.get(7) * ((m.get(8)  * m.get(13)) - (m.get(9)  * m.get(12)))) / det,
 
+			                    -(m.get(0) * ((m.get(9)  * m.get(15)) - (m.get(11) * m.get(13)))  -
+			                      m.get(1) * ((m.get(8)  * m.get(15)) - (m.get(11) * m.get(12)))  +
+			                      m.get(3) * ((m.get(8)  * m.get(13)) - (m.get(9)  * m.get(12)))) / det,
 
-		//Getters/Setters
-		constexpr auto operator[](size_type x);
-		constexpr auto operator[](size_type x)const;
-		constexpr reference get(size_type x, size_type y);
-		constexpr const_reference get(size_type x, size_type y)const;
-		constexpr reference get(size_type i);
-		constexpr const_reference get(size_type i)const;
+			                     (m.get(0) * ((m.get(5)  * m.get(15)) - (m.get(7)  * m.get(13)))  -
+			                      m.get(1) * ((m.get(4)  * m.get(15)) - (m.get(7)  * m.get(12)))  +
+			                      m.get(3) * ((m.get(4)  * m.get(13)) - (m.get(5)  * m.get(12)))) / det,
 
-		constexpr size_type columns()const;
-		constexpr size_type rows()const;
-		constexpr size_type size()const;
+			                    -(m.get(0) * ((m.get(5)  * m.get(11)) - (m.get(7)  * m.get(9)))   -
+			                      m.get(1) * ((m.get(4)  * m.get(11)) - (m.get(7)  * m.get(8)))   +
+			                      m.get(3) * ((m.get(4)  * m.get(9))  - (m.get(5)  * m.get(8))))  / det,
 
-		constexpr pointer raw();
-		constexpr const_pointer raw()const;
-		constexpr operator pointer();
-		constexpr operator const_pointer()const;
+			                    -(m.get(4) * ((m.get(9)  * m.get(14)) - (m.get(10) * m.get(13)))  -
+			                      m.get(5) * ((m.get(8)  * m.get(14)) - (m.get(10) * m.get(12)))  +
+			                      m.get(6) * ((m.get(8)  * m.get(13)) - (m.get(9)  * m.get(12)))) / det,
 
+			                     (m.get(0) * ((m.get(9)  * m.get(14)) - (m.get(10) * m.get(13)))  -
+			                      m.get(1) * ((m.get(8)  * m.get(14)) - (m.get(10) * m.get(12)))  +
+			                      m.get(2) * ((m.get(8)  * m.get(13)) - (m.get(9)  * m.get(12)))) / det,
+
+			                    -(m.get(0) * ((m.get(5)  * m.get(14)) - (m.get(6)  * m.get(13)))  -
+			                      m.get(1) * ((m.get(4)  * m.get(14)) - (m.get(6)  * m.get(12)))  +
+			                      m.get(2) * ((m.get(4)  * m.get(13)) - (m.get(5)  * m.get(12)))) / det,
+
+			                     (m.get(0) * ((m.get(5)  * m.get(10)) - (m.get(6)  * m.get(9)))   -
+			                      m.get(1) * ((m.get(4)  * m.get(10)) - (m.get(6)  * m.get(8)))   +
+			                      m.get(2) * ((m.get(4)  * m.get(9))  - (m.get(5)  * m.get(8))))  / det);
+		}
 	};
 
 }
diff --git a/include/detail/matrix.tpp b/include/detail/matrix.tpp
index 37f8725..b0adcd1 100644
--- a/include/detail/matrix.tpp
+++ b/include/detail/matrix.tpp
@@ -37,111 +37,6 @@ namespace math::detail{
 		return m_data[i*R];
 	}
 
-	template<typename T, size_t W, size_t H>
-	template<typename matrix_base<T,W,H>::size_type... Ss>
-	constexpr matrix_base<T,W,H>::matrix_base(std::integer_sequence<size_type,Ss...>):
-		m_data{Ss...}{}
-
-	template<typename T, size_t W, size_t H>
-	constexpr matrix_base<T,W,H>::matrix_base():
-		matrix_base(typename detail::default_initialization_matrix<Columns,Rows>::tuple{}){}
-
-	template<typename T, size_t W, size_t H>
-	constexpr matrix_base<T,W,H>::matrix_base(detail::zero_initialize_t):
-		m_data{}{}
-	template<typename T, size_t W, size_t H>
-	constexpr matrix_base<T,W,H>::matrix_base(detail::no_initialize_t){}
-	template<typename T, size_t W, size_t H>
-	template<typename U>
-	constexpr matrix_base<T,W,H>::matrix_base(detail::id_initialize_t):
-		matrix_base()
-	{
-		static_assert(Columns == Rows, "Identity initialization only supported on square matrices");
-	}
-
-	template<typename T, size_t W, size_t H>
-	constexpr matrix_base<T,W,H>::matrix_base(value_type v){
-		for(size_type i = 0;i < Columns*Rows;++i)
-			m_data[i] = v;
-	}
-	template<typename T, size_t W, size_t H>
-	template<typename... Args>
-	constexpr matrix_base<T,W,H>::matrix_base(Args&&... args):
-		m_data{std::forward<Args>(args)...}{}
-
-	template<typename T, size_t W, size_t H>
-	template<typename U>
-	constexpr matrix_base<T,W,H>::matrix_base(const matrix_base<U,Columns,Rows>& m){
-		using mat = decltype(m);
-		for(typename mat::size_type i = 0;i < mat::Columns*mat::Rows;++i)
-			m_data[i] = m.get(i);
-	}
-
-	template<typename T, size_t W, size_t H>
-	template<typename U>
-	constexpr matrix_base<T,W,H>& matrix_base<T,W,H>::operator=(const matrix_base<U,Columns,Rows>& m){
-		using mat = decltype(m);
-		for(typename mat::size_type i = 0;i < mat::Columns*mat::Rows;++i)
-			m_data[i] = m.get(i);
-		return *this;
-	}
-
-
-	template<typename T, size_t W, size_t H>
-	constexpr auto matrix_base<T,W,H>::operator[](size_type x){
-		return detail::mat_ref_obj<value_type,Rows>{m_data, x};
-	}
-	template<typename T, size_t W, size_t H>
-	constexpr auto matrix_base<T,W,H>::operator[](size_type x)const{
-		return detail::mat_ref_obj<const value_type,Rows>{m_data, x};
-	}
-	template<typename T, size_t W, size_t H>
-	constexpr auto matrix_base<T,W,H>::get(size_type x, size_type y) -> reference{
-		return m_data[x+(y*Rows)];
-	}
-	template<typename T, size_t W, size_t H>
-	constexpr auto matrix_base<T,W,H>::get(size_type x, size_type y)const -> const_reference{
-		return m_data[x+(y*Rows)];
-	}
-	template<typename T, size_t W, size_t H>
-	constexpr auto matrix_base<T,W,H>::get(size_type i) -> reference{
-		return m_data[i];
-	}
-	template<typename T, size_t W, size_t H>
-	constexpr auto matrix_base<T,W,H>::get(size_type i)const -> const_reference{
-		return m_data[i];
-	}
-
-	template<typename T, size_t W, size_t H>
-	constexpr auto matrix_base<T,W,H>::columns()const -> size_type{
-		return Columns;
-	}
-	template<typename T, size_t W, size_t H>
-	constexpr auto matrix_base<T,W,H>::rows()const -> size_type{
-		return Rows;
-	}
-	template<typename T, size_t W, size_t H>
-	constexpr auto matrix_base<T,W,H>::size()const -> size_type{
-		return Columns*Rows;
-	}
-
-	template<typename T, size_t W, size_t H>
-	constexpr auto matrix_base<T,W,H>::raw() -> pointer{
-		return m_data;
-	}
-	template<typename T, size_t W, size_t H>
-	constexpr auto matrix_base<T,W,H>::raw()const -> const_pointer{
-		return m_data;
-	}
-	template<typename T, size_t W, size_t H>
-	constexpr matrix_base<T,W,H>::operator pointer(){
-		return m_data;
-	}
-	template<typename T, size_t W, size_t H>
-	constexpr matrix_base<T,W,H>::operator const_pointer()const{
-		return m_data;
-	}
-
 }
 
 #endif
diff --git a/include/mat.hpp b/include/mat.hpp
index 3d18724..4a1097d 100644
--- a/include/mat.hpp
+++ b/include/mat.hpp
@@ -23,47 +23,122 @@
 #include <utility> //integer_sequence
 #include <type_traits> //decay_t, is_same, integral_constant
 #include "detail/math.hpp"
-#include "detail/matrix.hpp"
 
 namespace math{
 
-	template<typename T, size_t C, size_t R>
-	class matrix : public detail::matrix_base<T,C,R>
+	template<typename T, size_t W, size_t H>
+	class matrix
 	{
-	private:
-		using base = detail::matrix_base<T,C,R>;
+		static_assert(W > 0, "Cannot have 0 columns matrix");
+		static_assert(H > 0, "Cannot have 0 rows matrix");
 	public:
-		using value_type = typename base::value_type;
-		using size_type = typename base::size_type;
-		using pointer = typename base::pointer;
-		using const_pointer = typename base::const_pointer;
-		using reference = typename base::reference;
-		using const_reference = typename base::const_reference;
+		using value_type = T;
+		using size_type = size_t;
+		using pointer = value_type*;
+		using const_pointer = const value_type*;
+		using reference = value_type&;
+		using const_reference = const value_type&;
+
+		static constexpr size_type Columns = W;
+		static constexpr size_type Rows = H;
+
+	protected:
+		value_type m_data[W*H];
+
+	protected:
+		template<size_type... Ss>
+		constexpr matrix(std::integer_sequence<size_type,Ss...>);
+
 	public:
-		using detail::matrix_base<T,C,R>::matrix_base;
-		using detail::matrix_base<T,C,R>::operator=;
+		//Default construct as identity when square, zero otherwise
+		constexpr matrix();
+
+		//Range initializing constructors
+		constexpr explicit matrix(detail::zero_initialize_t);
+		constexpr explicit matrix(detail::no_initialize_t);
+		template<typename U = void>
+		constexpr explicit matrix(detail::id_initialize_t);
+
+		//Value initializing constructors
+		constexpr explicit matrix(value_type v);
+		template<typename... Args>
+		constexpr explicit matrix(Args&&... args);
+
+		//Copying constructors
+		constexpr matrix(const matrix&) = default;
+		constexpr matrix(matrix&&) = default;
+		template<typename U>
+		constexpr matrix(const matrix<U,Columns,Rows>& m);
+		~matrix() = default;
+
+		//Assignement
+		template<typename U>
+		constexpr matrix& operator=(const matrix<U,Columns,Rows>& m);
+		constexpr matrix& operator=(const matrix&) = default;
+		constexpr matrix& operator=(matrix&&) = default;
+
+
+		//Getters/Setters
+		constexpr auto operator[](size_type x);
+		constexpr auto operator[](size_type x)const;
+		constexpr reference get(size_type x, size_type y);
+		constexpr const_reference get(size_type x, size_type y)const;
+		constexpr reference get(size_type i);
+		constexpr const_reference get(size_type i)const;
+
+		constexpr size_type columns()const;
+		constexpr size_type rows()const;
+		constexpr size_type size()const;
+
+		constexpr pointer raw();
+		constexpr const_pointer raw()const;
+		constexpr operator pointer();
+		constexpr operator const_pointer()const;
+
+		//square matrix arithmetic operations
+		template<typename U = T>
+		constexpr value_type determinate()const;
+		template<typename U = T>
+		constexpr value_type trace()const;
+		template<typename U = T>
+		constexpr matrix transpose()const;
+		template<typename U = T>
+		constexpr matrix inverse()const;
 	};
 
+	template<typename T, size_t R>
+	constexpr T determinate(const matrix<T,R,R>& m);
+	template<typename T, size_t R>
+	constexpr matrix<T,R,R> inverse(const matrix<T,R,R>& m);
+
 	template<typename T>
-	class matrix<T,3,3> : public detail::matrix_base<T,3,3>
-	{
-	private:
-		using base = detail::matrix_base<T,3,3>;
-	public:
-		using value_type = typename base::value_type;
-		using size_type = typename base::size_type;
-		using pointer = typename base::pointer;
-		using const_pointer = typename base::const_pointer;
-		using reference = typename base::reference;
-		using const_reference = typename base::const_reference;
-	public:
-		using detail::matrix_base<T,3,3>::matrix_base;
-		using detail::matrix_base<T,3,3>::operator=;
+	matrix<T,2,2> rotation2d_pure(T angle);
+	template<typename T>
+	constexpr matrix<T,2,2> rotation2d_pure(T sin, T cos);
+	template<typename T>
+	constexpr matrix<T,2,2> scale2d(T x, T y);
 
-		static matrix rotation(value_type angle);
-		static constexpr matrix rotation(value_type sin, value_type cos);
-		static constexpr matrix rotation(value_type angle_x, value_type angle_y, value_type angle_z);
-	};
+	template<typename T>
+	matrix<T,3,3> rotation2d(T angle);
+	template<typename T>
+	constexpr matrix<T,3,3> rotation2d(T sin, T cos);
+	template<typename T>
+	matrix<T,3,3> rotation2d(T x, T y, T z);
+
+	template<typename T>
+	matrix<T,4,4> fov_projection(T fov, T asp, T near, T far);
+	template<typename T>
+	matrix<T,4,4> fov_asymetric_projection(T fovl, T fovr, T fovb, T fovt, T asp, T near, T far);
+	template<typename T>
+	matrix<T,4,4> ortho_projection(T w, T h, T n, T f);
+	template<typename T>
+	matrix<T,4,4> ortho_asymetric_projection(T l, T r, T b, T t, T n, T f);
+	template<typename T>
+	constexpr matrix<T,4,4> rotation3d(T angle_x, T angle_y, T angle_z);
+	template<typename T>
+	constexpr matrix<T,4,4> translation3d(T x, T y, T z);
+	template<typename T>
+	constexpr matrix<T,4,4> scale3d(T x, T y, T z);
 
 	namespace detail{
 
@@ -97,38 +172,39 @@ namespace math{
 		template<typename M1, typename M2>
 		using enable_if_eq_matrix = std::enable_if_t<are_same_size_matrix<M1,M2>::value,int>;
 
-
-		template<typename T, typename U, size_t W, size_t H>
-		constexpr bool operator==(const matrix<T,W,H>& left, const matrix<U,W,H> right);
-		template<typename T, typename U, size_t W, size_t H>
-		constexpr bool operator!=(const matrix<T,W,H>& left, const matrix<U,W,H> right);
-
-		template<typename T, typename U, size_t R1, size_t C1, size_t C2>
-		constexpr auto operator*(const matrix<T,C1,R1>& left, const matrix<U,C2,C1>& right);
-		template<typename T, typename U, size_t C, size_t R>
-		constexpr auto operator*(const matrix<T,C,R>& left, U&& right);
-		template<typename T, typename U, size_t C, size_t R>
-		constexpr auto operator/(const matrix<T,C,R>& left, U&& right);
-		template<typename T, typename U, size_t C, size_t R>
-		constexpr auto operator+(const matrix<T,C,R>& left, const matrix<U,C,R>& right);
-		template<typename T, typename U, size_t C, size_t R>
-		constexpr auto operator-(const matrix<T,C,R>& left, const matrix<U,C,R>& right);
-		template<typename T, typename U, size_t C, size_t R>
-		constexpr auto operator-(const matrix<T,C,R>& left);
-
-		template<typename T, typename U, size_t R1, size_t C1, size_t C2>
-		constexpr decltype(auto) operator*=(matrix<T,C1,R1>& left, const matrix<U,C2,C1>& right);
-		template<typename T, typename U, size_t C, size_t R>
-		constexpr decltype(auto) operator*=(matrix<T,C,R>& left, U&& right);
-		template<typename T, typename U, size_t C, size_t R>
-		constexpr decltype(auto) operator/=(matrix<T,C,R>& left, U&& right);
-		template<typename T, typename U, size_t C, size_t R>
-		constexpr decltype(auto) operator+=(matrix<T,C,R>& left, const matrix<U,C,R>& right);
-		template<typename T, typename U, size_t C, size_t R>
-		constexpr decltype(auto) operator-=(matrix<T,C,R>& left, const matrix<U,C,R>& right);
-
 	}
 
+	template<typename T, typename U, size_t W, size_t H>
+	constexpr bool operator==(const matrix<T,W,H>& left, const matrix<U,W,H> right);
+	template<typename T, typename U, size_t W, size_t H>
+	constexpr bool operator!=(const matrix<T,W,H>& left, const matrix<U,W,H> right);
+
+	template<typename T, typename U, size_t R1, size_t C1, size_t C2>
+	constexpr auto operator*(const matrix<T,C1,R1>& left, const matrix<U,C2,C1>& right);
+	template<typename T, typename U, size_t C, size_t R>
+	constexpr auto operator*(const matrix<T,C,R>& left, U&& right);
+	template<typename T, typename U, size_t C, size_t R>
+	constexpr auto operator*(U&& left, const matrix<T,C,R>& right);
+	template<typename T, typename U, size_t C, size_t R>
+	constexpr auto operator/(const matrix<T,C,R>& left, U&& right);
+	template<typename T, typename U, size_t C, size_t R>
+	constexpr auto operator+(const matrix<T,C,R>& left, const matrix<U,C,R>& right);
+	template<typename T, typename U, size_t C, size_t R>
+	constexpr auto operator-(const matrix<T,C,R>& left, const matrix<U,C,R>& right);
+	template<typename T, typename U, size_t C, size_t R>
+	constexpr auto operator-(const matrix<T,C,R>& left);
+
+	template<typename T, typename U, size_t R1, size_t C1, size_t C2>
+	constexpr decltype(auto) operator*=(matrix<T,C1,R1>& left, const matrix<U,C2,C1>& right);
+	template<typename T, typename U, size_t C, size_t R>
+	constexpr decltype(auto) operator*=(matrix<T,C,R>& left, U&& right);
+	template<typename T, typename U, size_t C, size_t R>
+	constexpr decltype(auto) operator/=(matrix<T,C,R>& left, U&& right);
+	template<typename T, typename U, size_t C, size_t R>
+	constexpr decltype(auto) operator+=(matrix<T,C,R>& left, const matrix<U,C,R>& right);
+	template<typename T, typename U, size_t C, size_t R>
+	constexpr decltype(auto) operator-=(matrix<T,C,R>& left, const matrix<U,C,R>& right);
+
 }
 
 #include "mat.tpp"
diff --git a/include/mat.tpp b/include/mat.tpp
index afed470..9bd9de7 100644
--- a/include/mat.tpp
+++ b/include/mat.tpp
@@ -22,26 +22,241 @@
 #include <cstdlib> //size_t
 #include <cmath> //sin, cos
 #include <type_traits> //decay_t, declval
+#include "detail/matrix.hpp"
 
 namespace math{
 
+	template<typename T, size_t W, size_t H>
+	template<typename matrix<T,W,H>::size_type... Ss>
+	constexpr matrix<T,W,H>::matrix(std::integer_sequence<size_type,Ss...>):
+		m_data{Ss...}{}
+
+	template<typename T, size_t W, size_t H>
+	constexpr matrix<T,W,H>::matrix():
+		matrix(typename detail::default_initialization_matrix<Columns,Rows>::tuple{}){}
+
+	template<typename T, size_t W, size_t H>
+	constexpr matrix<T,W,H>::matrix(detail::zero_initialize_t):
+		m_data{}{}
+	template<typename T, size_t W, size_t H>
+	constexpr matrix<T,W,H>::matrix(detail::no_initialize_t){}
+	template<typename T, size_t W, size_t H>
+	template<typename U>
+	constexpr matrix<T,W,H>::matrix(detail::id_initialize_t):
+		matrix()
+	{
+		static_assert(Columns == Rows, "Identity initialization only supported on square matrices");
+	}
+
+	template<typename T, size_t W, size_t H>
+	constexpr matrix<T,W,H>::matrix(value_type v){
+		for(size_type i = 0;i < Columns*Rows;++i)
+			m_data[i] = v;
+	}
+	template<typename T, size_t W, size_t H>
+	template<typename... Args>
+	constexpr matrix<T,W,H>::matrix(Args&&... args):
+		m_data{std::forward<Args>(args)...}{}
+
+	template<typename T, size_t W, size_t H>
+	template<typename U>
+	constexpr matrix<T,W,H>::matrix(const matrix<U,Columns,Rows>& m){
+		using mat = decltype(m);
+		for(typename mat::size_type i = 0;i < mat::Columns*mat::Rows;++i)
+			m_data[i] = m.get(i);
+	}
+
+	template<typename T, size_t W, size_t H>
+	template<typename U>
+	constexpr matrix<T,W,H>& matrix<T,W,H>::operator=(const matrix<U,Columns,Rows>& m){
+		using mat = decltype(m);
+		for(typename mat::size_type i = 0;i < mat::Columns*mat::Rows;++i)
+			m_data[i] = m.get(i);
+		return *this;
+	}
+
+
+	template<typename T, size_t W, size_t H>
+	constexpr auto matrix<T,W,H>::operator[](size_type x){
+		return detail::mat_ref_obj<value_type,Rows>{m_data, x};
+	}
+	template<typename T, size_t W, size_t H>
+	constexpr auto matrix<T,W,H>::operator[](size_type x)const{
+		return detail::mat_ref_obj<const value_type,Rows>{m_data, x};
+	}
+	template<typename T, size_t W, size_t H>
+	constexpr auto matrix<T,W,H>::get(size_type x, size_type y) -> reference{
+		return m_data[x+(y*Rows)];
+	}
+	template<typename T, size_t W, size_t H>
+	constexpr auto matrix<T,W,H>::get(size_type x, size_type y)const -> const_reference{
+		return m_data[x+(y*Rows)];
+	}
+	template<typename T, size_t W, size_t H>
+	constexpr auto matrix<T,W,H>::get(size_type i) -> reference{
+		return m_data[i];
+	}
+	template<typename T, size_t W, size_t H>
+	constexpr auto matrix<T,W,H>::get(size_type i)const -> const_reference{
+		return m_data[i];
+	}
+
+	template<typename T, size_t W, size_t H>
+	constexpr auto matrix<T,W,H>::columns()const -> size_type{
+		return Columns;
+	}
+	template<typename T, size_t W, size_t H>
+	constexpr auto matrix<T,W,H>::rows()const -> size_type{
+		return Rows;
+	}
+	template<typename T, size_t W, size_t H>
+	constexpr auto matrix<T,W,H>::size()const -> size_type{
+		return Columns*Rows;
+	}
+
+	template<typename T, size_t W, size_t H>
+	constexpr auto matrix<T,W,H>::raw() -> pointer{
+		return m_data;
+	}
+	template<typename T, size_t W, size_t H>
+	constexpr auto matrix<T,W,H>::raw()const -> const_pointer{
+		return m_data;
+	}
+	template<typename T, size_t W, size_t H>
+	constexpr matrix<T,W,H>::operator pointer(){
+		return m_data;
+	}
+	template<typename T, size_t W, size_t H>
+	constexpr matrix<T,W,H>::operator const_pointer()const{
+		return m_data;
+	}
+	template<typename T, size_t W, size_t H>
+	template<typename U>
+	constexpr auto matrix<T,W,H>::determinate()const -> value_type{
+		static_assert(W == H, "Determinate can only be calculated on square matrix");
+		return math::determinate(*this);
+	}
+	template<typename T, size_t W, size_t H>
+	template<typename U>
+	constexpr auto matrix<T,W,H>::trace()const -> value_type{
+		static_assert(W == H, "Trace can only be calculated on square matrix");
+		value_type sum = 0;
+		for(size_type i = 0;i < W;++i){
+			sum += get(i, i);
+		}
+	}
+	template<typename T, size_t W, size_t H>
+	template<typename U>
+	constexpr matrix<T,W,H> matrix<T,W,H>::transpose()const{
+		matrix m(no_initialize);
+		for(size_type i = 0;i < W;++i){
+			for(size_type j = 0;j < W;++j){
+				m.get(j, i) = get(i, j);
+			}
+		}
+		return m;
+	}
+	template<typename T, size_t W, size_t H>
+	template<typename U>
+	constexpr matrix<T,W,H> matrix<T,W,H>::inverse()const{
+		static_assert(W == H, "Trace can only be calculated on square matrix");
+		return math::inverse(*this);
+	}
+
+
+	template<typename T, size_t R>
+	constexpr T determinate(const matrix<T,R,R>& m){
+		return detail::determinate_helper<T,R>::perform(m);
+	}
+	template<typename T, size_t R>
+	constexpr matrix<T,R,R> inverse(const matrix<T,R,R>& m){
+		return detail::inverse_helper<T,R>::perform(m);
+	}
+
 	template<typename T>
-	matrix<T,3,3> matrix<T,3,3>::rotation(value_type angle){
-		value_type c = std::cos(angle);
-		value_type s = std::sin(angle);
-		return rotation(s, c);
+	matrix<T,2,2> rotation2d_pure(T angle){
+		return rotation2d_pure(std::sin(angle), std::cos(angle));
 	}
 	template<typename T>
-	constexpr matrix<T,3,3> matrix<T,3,3>::rotation(value_type sin, value_type cos){
-		return matrix(cos, -sin, 0,
-		              sin,  cos, 0,
-		              0,  0, 1);
+	constexpr matrix<T,2,2> rotation2d_pure(T sin, T cos){
+		return matrix<T,2,2>(cos, sin, -sin, cos);
 	}
 	template<typename T>
-	constexpr matrix<T,3,3> matrix<T,3,3>::rotation(value_type angle_x, value_type angle_y, value_type angle_z){
+	constexpr matrix<T,2,2> scale2d(T x, T y){
+		return matrix<T,2,2>(x, 0, 0, y);
+	}
+
+	template<typename T>
+	matrix<T,3,3> rotation2d(T angle){
+		return rotation2d(std::sin(angle), std::cos(angle));
+	}
+	template<typename T>
+	constexpr matrix<T,3,3> rotation2d(T sin, T cos){
+		return matrix<T,3,3>(cos, -sin, 0,
+		                     sin,  cos, 0,
+		                     0,    0,   1);
+	}
+	template<typename T>
+	matrix<T,3,3> rotation2d(T x, T y, T z){
 		//TODO
+		return {};
 	}
 
+	template<typename T>
+	matrix<T,4,4> fov_projection(T fov, T asp, T near, T far){
+		T r = near * std::tan(fov / T{2.0});
+		return matrix<T,4,4>((near / r) / asp, 0,          0,                               0,
+		                     0,                (near / r), 0,                               0,
+		                     0,                0,          (far + near) / (near - far),    -1,
+		                     0,                0,          (2 * near * far) / (near - far), 0);
+	}
+	template<typename T>
+	matrix<T,4,4> fov_asymetric_projection(T fovl, T fovr, T fovb, T fovt, T asp, T n, T f){
+		T l = n * std::tan(fovl);
+		T r = n * std::tan(fovr);
+		T b = n * std::tan(fovb);
+		T t = n * std::tan(fovt);
+
+		return matrix<T,4,4>(((2 * n) / (r - l)) * asp, 0,                 0,                     0,
+		                     0,                         (2 * n) / (t - b), 0,                     0,
+		                     (r + l) / (r - l),         (t + b) / (t - b), (f + n) / (n - f),    -1,
+		                     0,                         0,                 (2 * n * f) / (n - f), 0);
+	}
+	template<typename T>
+	matrix<T,4,4> ortho_projection(T w, T h, T n, T f){
+		return matrix<T,4,4>(2 / w, 0,     0,                 0,
+		                     0,     2 / h, 0,                 0,
+		                     0,     0,     2 / (n - f),       0,
+		                     0,     0,     (n + f) / (n - f), 1);
+	}
+	template<typename T>
+	matrix<T,4,4> ortho_asymetric_projection(T l, T r, T b, T t, T n, T f){
+		return matrix<T,4,4>(2 / (r - l),       0,                 0,                 0,
+		                     0,                 2 / (t - b),       0,                 0,
+		                     0,                 0,                 2 / (n - f),       0,
+		                     (r + l) / (l - r), (t + b) / (b - t), (n + f) / (n - f), 1);
+	}
+	template<typename T>
+	constexpr matrix<T,4,4> rotation3d(T angle_x, T angle_y, T angle_z){
+		//TODO
+		return {};
+	}
+	template<typename T>
+	constexpr matrix<T,4,4> translation3d(T x, T y, T z){
+		return matrix<T,4,4>(1, 0, 0, 0,
+		                     0, 1, 0, 0,
+		                     0, 0, 1, 0,
+		                     x, y, z, 1);
+	}
+	template<typename T>
+	constexpr matrix<T,4,4> scale3d(T x, T y, T z){
+		return matrix<T,4,4>(x, 0, 0, 0,
+		                     0, y, 0, 0,
+		                     0, 0, z, 0,
+		                     0, 0, 0, 1);
+	}
+
+
 	template<typename T, typename U, size_t W, size_t H>
 	constexpr bool operator==(const matrix<T,W,H>& left, const matrix<U,W,H> right){
 		for(size_t i = 0;i < left.size();++i){
@@ -80,6 +295,15 @@ namespace math{
 		return res;
 	}
 	template<typename T, typename U, size_t C, size_t R>
+	constexpr auto operator*(U&& left, const matrix<T,C,R>& right){
+		using res_t = decltype(std::declval<T>() * std::declval<U>());
+		matrix<res_t,C,R> res(no_initialize);
+		for(size_t i = 0;i < right.size();++i){
+			res.get(i) = std::forward<U>(left) * right.get(i);
+		}
+		return res;
+	}
+	template<typename T, typename U, size_t C, size_t R>
 	constexpr auto operator/(const matrix<T,C,R>& left, U&& right){
 		using res_t = decltype(std::declval<T>() / std::declval<U>());
 		matrix<res_t,C,R> res(no_initialize);
diff --git a/makefile b/makefile
index 1892bb0..696b3ca 100644
--- a/makefile
+++ b/makefile
@@ -40,7 +40,7 @@ ifneq ($(WINDOWS),1)
 	CC::=gcc
 	CXX::=g++
 	LDLIBS::=
-	LDFLAGS::= -lglfw -ldl
+	LDFLAGS::= -lglfw -ldl -lm
 	STRIP::=strip
 	RANLIB::=ranlib
 	AR::=ar