From 90a249b5ffb06cbb8303f2d6b4f039f8eb567046 Mon Sep 17 00:00:00 2001
From: Christophe Riccio <mail@g-truc.net>
Date: Sat, 11 Jan 2014 16:44:15 +0100
Subject: [PATCH] Optimized matrix inverse and division code (#149)

---
 glm/detail/func_matrix.inl     | 113 ++-----------------------------
 glm/detail/type_mat.hpp        |  22 +-----
 glm/detail/type_mat2x2.inl     |  23 ++++++-
 glm/detail/type_mat3x3.inl     |  32 ++++++++-
 glm/detail/type_mat4x4.inl     | 118 ++++++++++++++++++---------------
 glm/gtc/matrix_transform.inl   |   5 +-
 readme.txt                     |   1 +
 test/core/core_func_matrix.cpp |  72 ++++++++++++++++++--
 test/core/core_type_mat4x4.cpp |  76 ++++++++++++++++++---
 9 files changed, 255 insertions(+), 207 deletions(-)

diff --git a/glm/detail/func_matrix.inl b/glm/detail/func_matrix.inl
index 761006a8..c6a2c025 100644
--- a/glm/detail/func_matrix.inl
+++ b/glm/detail/func_matrix.inl
@@ -416,7 +416,7 @@ namespace detail
 	{
 		static T call(detail::tmat3x3<T, P> const & m)
 		{
-			return 
+			return
 				+ m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2])
 				- m[1][0] * (m[0][1] * m[2][2] - m[2][1] * m[0][2])
 				+ m[2][0] * (m[0][1] * m[1][2] - m[1][1] * m[0][2]);
@@ -441,114 +441,9 @@ namespace detail
 				+ (m[1][0] * SubFactor01 - m[1][1] * SubFactor03 + m[1][3] * SubFactor05),
 				- (m[1][0] * SubFactor02 - m[1][1] * SubFactor04 + m[1][2] * SubFactor05));
 
-			return m[0][0] * DetCof[0]
-					+ m[0][1] * DetCof[1]
-					+ m[0][2] * DetCof[2]
-					+ m[0][3] * DetCof[3];
-		}
-	};
-
-	template <template <class, precision> class matType, typename T, precision P>
-	struct compute_inverse{};
-
-	template <typename T, precision P>
-	struct compute_inverse<detail::tmat2x2, T, P>
-	{
-		static detail::tmat2x2<T, P> call(detail::tmat2x2<T, P> const & m)
-		{
-			T Determinant = determinant(m);
-
-			detail::tmat2x2<T, P> Inverse(
-				+ m[1][1] / Determinant,
-				- m[0][1] / Determinant,
-				- m[1][0] / Determinant,
-				+ m[0][0] / Determinant);
-
-			return Inverse;
-		}
-	};
-
-	template <typename T, precision P>
-	struct compute_inverse<detail::tmat3x3, T, P>
-	{
-		static detail::tmat3x3<T, P> call(detail::tmat3x3<T, P> const & m)
-		{
-			T Determinant = determinant(m);
-
-			detail::tmat3x3<T, P> Inverse(detail::tmat3x3<T, P>::_null);
-			Inverse[0][0] = + (m[1][1] * m[2][2] - m[2][1] * m[1][2]);
-			Inverse[1][0] = - (m[1][0] * m[2][2] - m[2][0] * m[1][2]);
-			Inverse[2][0] = + (m[1][0] * m[2][1] - m[2][0] * m[1][1]);
-			Inverse[0][1] = - (m[0][1] * m[2][2] - m[2][1] * m[0][2]);
-			Inverse[1][1] = + (m[0][0] * m[2][2] - m[2][0] * m[0][2]);
-			Inverse[2][1] = - (m[0][0] * m[2][1] - m[2][0] * m[0][1]);
-			Inverse[0][2] = + (m[0][1] * m[1][2] - m[1][1] * m[0][2]);
-			Inverse[1][2] = - (m[0][0] * m[1][2] - m[1][0] * m[0][2]);
-			Inverse[2][2] = + (m[0][0] * m[1][1] - m[1][0] * m[0][1]);
-			Inverse /= Determinant;
-
-			return Inverse;
-		}
-	};
-
-	template <typename T, precision P>
-	struct compute_inverse<detail::tmat4x4, T, P>
-	{
-		static detail::tmat4x4<T, P> call(detail::tmat4x4<T, P> const & m)
-		{
-			T Coef00 = m[2][2] * m[3][3] - m[3][2] * m[2][3];
-			T Coef02 = m[1][2] * m[3][3] - m[3][2] * m[1][3];
-			T Coef03 = m[1][2] * m[2][3] - m[2][2] * m[1][3];
-
-			T Coef04 = m[2][1] * m[3][3] - m[3][1] * m[2][3];
-			T Coef06 = m[1][1] * m[3][3] - m[3][1] * m[1][3];
-			T Coef07 = m[1][1] * m[2][3] - m[2][1] * m[1][3];
-
-			T Coef08 = m[2][1] * m[3][2] - m[3][1] * m[2][2];
-			T Coef10 = m[1][1] * m[3][2] - m[3][1] * m[1][2];
-			T Coef11 = m[1][1] * m[2][2] - m[2][1] * m[1][2];
-
-			T Coef12 = m[2][0] * m[3][3] - m[3][0] * m[2][3];
-			T Coef14 = m[1][0] * m[3][3] - m[3][0] * m[1][3];
-			T Coef15 = m[1][0] * m[2][3] - m[2][0] * m[1][3];
-
-			T Coef16 = m[2][0] * m[3][2] - m[3][0] * m[2][2];
-			T Coef18 = m[1][0] * m[3][2] - m[3][0] * m[1][2];
-			T Coef19 = m[1][0] * m[2][2] - m[2][0] * m[1][2];
-
-			T Coef20 = m[2][0] * m[3][1] - m[3][0] * m[2][1];
-			T Coef22 = m[1][0] * m[3][1] - m[3][0] * m[1][1];
-			T Coef23 = m[1][0] * m[2][1] - m[2][0] * m[1][1];
-
-			detail::tvec4<T, P> const SignA(+1, -1, +1, -1);
-			detail::tvec4<T, P> const SignB(-1, +1, -1, +1);
-
-			detail::tvec4<T, P> Fac0(Coef00, Coef00, Coef02, Coef03);
-			detail::tvec4<T, P> Fac1(Coef04, Coef04, Coef06, Coef07);
-			detail::tvec4<T, P> Fac2(Coef08, Coef08, Coef10, Coef11);
-			detail::tvec4<T, P> Fac3(Coef12, Coef12, Coef14, Coef15);
-			detail::tvec4<T, P> Fac4(Coef16, Coef16, Coef18, Coef19);
-			detail::tvec4<T, P> Fac5(Coef20, Coef20, Coef22, Coef23);
-
-			detail::tvec4<T, P> Vec0(m[1][0], m[0][0], m[0][0], m[0][0]);
-			detail::tvec4<T, P> Vec1(m[1][1], m[0][1], m[0][1], m[0][1]);
-			detail::tvec4<T, P> Vec2(m[1][2], m[0][2], m[0][2], m[0][2]);
-			detail::tvec4<T, P> Vec3(m[1][3], m[0][3], m[0][3], m[0][3]);
-
-			detail::tvec4<T, P> Inv0 = SignA * (Vec1 * Fac0 - Vec2 * Fac1 + Vec3 * Fac2);
-			detail::tvec4<T, P> Inv1 = SignB * (Vec0 * Fac0 - Vec2 * Fac3 + Vec3 * Fac4);
-			detail::tvec4<T, P> Inv2 = SignA * (Vec0 * Fac1 - Vec1 * Fac3 + Vec3 * Fac5);
-			detail::tvec4<T, P> Inv3 = SignB * (Vec0 * Fac2 - Vec1 * Fac4 + Vec2 * Fac5);
-
-			detail::tmat4x4<T, P> Inverse(Inv0, Inv1, Inv2, Inv3);
-
-			detail::tvec4<T, P> Row0(Inverse[0][0], Inverse[1][0], Inverse[2][0], Inverse[3][0]);
-
-			T Determinant = dot(m[0], Row0);
-
-			Inverse /= Determinant;
-
-			return Inverse;
+			return
+				m[0][0] * DetCof[0] + m[0][1] * DetCof[1] +
+				m[0][2] * DetCof[2] + m[0][3] * DetCof[3];
 		}
 	};
 }//namespace detail
diff --git a/glm/detail/type_mat.hpp b/glm/detail/type_mat.hpp
index 3d213bfb..25095382 100644
--- a/glm/detail/type_mat.hpp
+++ b/glm/detail/type_mat.hpp
@@ -44,26 +44,8 @@ namespace detail
 	template <typename T, precision P> struct tmat4x3;
 	template <typename T, precision P> struct tmat4x4;
 	
-	template <typename T>
-	struct is_matrix
-	{
-		enum is_matrix_enum
-		{
-			_YES = 0,
-			_NO = 1
-		};
-	};
-	
-	#define GLM_DETAIL_IS_MATRIX(T)	\
-		template <>					\
-		struct is_matrix			\
-		{							\
-			enum is_matrix_enum		\
-			{						\
-				_YES = 1,			\
-				_NO = 0				\
-			};						\
-		}
+	template <template <class, precision> class matType, typename T, precision P>
+	struct compute_inverse{};
 }//namespace detail
 
 	/// @addtogroup core_precision
diff --git a/glm/detail/type_mat2x2.inl b/glm/detail/type_mat2x2.inl
index 72e3e1ec..fc7414d8 100644
--- a/glm/detail/type_mat2x2.inl
+++ b/glm/detail/type_mat2x2.inl
@@ -358,7 +358,7 @@ namespace detail
 	template <typename U>
 	GLM_FUNC_QUALIFIER tmat2x2<T, P>& tmat2x2<T, P>::operator/= (tmat2x2<U, P> const & m)
 	{
-		return (*this = *this * compute_inverse_mat2(m));
+		return (*this = *this * detail::compute_inverse<detail::tmat2x2, T, P>::call(m));
 	}
 
 	template <typename T, precision P>
@@ -393,6 +393,23 @@ namespace detail
 		return Result;
 	}
 
+	template <typename T, precision P>
+	struct compute_inverse<detail::tmat2x2, T, P>
+	{
+		static detail::tmat2x2<T, P> call(detail::tmat2x2<T, P> const & m)
+		{
+			T Determinant = determinant(m);
+
+			detail::tmat2x2<T, P> Inverse(
+				+ m[1][1] / Determinant,
+				- m[0][1] / Determinant,
+				- m[1][0] / Determinant,
+				+ m[0][0] / Determinant);
+
+			return Inverse;
+		}
+	};
+
 	template <typename T, precision P>
 	GLM_FUNC_QUALIFIER tmat2x2<T, P> compute_inverse_mat2(tmat2x2<T, P> const & m)
 	{
@@ -607,7 +624,7 @@ namespace detail
 		typename tmat2x2<T, P>::row_type & v
 	)
 	{
-		return detail::compute_inverse_mat2(m) * v;
+		return detail::compute_inverse<detail::tmat2x2, T, P>::call(m) * v;
 	}
 
 	template <typename T, precision P>
@@ -617,7 +634,7 @@ namespace detail
 		tmat2x2<T, P> const & m
 	)
 	{
-		return v * detail::compute_inverse_mat2(m);
+		return v * detail::compute_inverse<detail::tmat2x2, T, P>::call(m);
 	}
 
 	template <typename T, precision P>
diff --git a/glm/detail/type_mat3x3.inl b/glm/detail/type_mat3x3.inl
index 9c473e12..81a15874 100644
--- a/glm/detail/type_mat3x3.inl
+++ b/glm/detail/type_mat3x3.inl
@@ -393,7 +393,7 @@ namespace detail
 	template <typename U>
 	GLM_FUNC_QUALIFIER tmat3x3<T, P> & tmat3x3<T, P>::operator/= (tmat3x3<U, P> const & m)
 	{
-		return (*this = *this * detail::compute_inverse_mat3(m));
+		return (*this = *this * detail::compute_inverse<detail::tmat3x3, T, P>::call(m));
 	}
 
 	template <typename T, precision P>
@@ -430,6 +430,32 @@ namespace detail
 		return Result;
 	}
 
+	template <typename T, precision P>
+	struct compute_inverse<detail::tmat3x3, T, P>
+	{
+		static detail::tmat3x3<T, P> call(detail::tmat3x3<T, P> const & m)
+		{
+			T Determinant =
+				+ m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2])
+				- m[1][0] * (m[0][1] * m[2][2] - m[2][1] * m[0][2])
+				+ m[2][0] * (m[0][1] * m[1][2] - m[1][1] * m[0][2]);
+
+			detail::tmat3x3<T, P> Inverse(detail::tmat3x3<T, P>::_null);
+			Inverse[0][0] = + (m[1][1] * m[2][2] - m[2][1] * m[1][2]);
+			Inverse[1][0] = - (m[1][0] * m[2][2] - m[2][0] * m[1][2]);
+			Inverse[2][0] = + (m[1][0] * m[2][1] - m[2][0] * m[1][1]);
+			Inverse[0][1] = - (m[0][1] * m[2][2] - m[2][1] * m[0][2]);
+			Inverse[1][1] = + (m[0][0] * m[2][2] - m[2][0] * m[0][2]);
+			Inverse[2][1] = - (m[0][0] * m[2][1] - m[2][0] * m[0][1]);
+			Inverse[0][2] = + (m[0][1] * m[1][2] - m[1][1] * m[0][2]);
+			Inverse[1][2] = - (m[0][0] * m[1][2] - m[1][0] * m[0][2]);
+			Inverse[2][2] = + (m[0][0] * m[1][1] - m[1][0] * m[0][1]);
+			Inverse /= Determinant;
+
+			return Inverse;
+		}
+	};
+
 	template <typename T, precision P>
 	GLM_FUNC_QUALIFIER tmat3x3<T, P> compute_inverse_mat3(tmat3x3<T, P> const & m)
 	{
@@ -720,7 +746,7 @@ namespace detail
 		typename tmat3x3<T, P>::row_type const & v
 	)
 	{
-		return detail::compute_inverse_mat3(m) * v;
+		return detail::compute_inverse<detail::tmat3x3, T, P>::call(m) * v;
 	}
 
 	template <typename T, precision P>
@@ -730,7 +756,7 @@ namespace detail
 		tmat3x3<T, P> const & m
 	)
 	{
-		return v * detail::compute_inverse_mat3(m);
+		return v * detail::compute_inverse<detail::tmat3x3, T, P>::call(m);
 	}
 
 	template <typename T, precision P>
diff --git a/glm/detail/type_mat4x4.inl b/glm/detail/type_mat4x4.inl
index 3de0d7af..db14796d 100644
--- a/glm/detail/type_mat4x4.inl
+++ b/glm/detail/type_mat4x4.inl
@@ -80,10 +80,10 @@ namespace detail
 		tmat4x4<T, P> const & m
 	)
 	{
-		this->value[0] = m.value[0];
-		this->value[1] = m.value[1];
-		this->value[2] = m.value[2];
-		this->value[3] = m.value[3];
+		this->value[0] = m[0];
+		this->value[1] = m[1];
+		this->value[2] = m[2];
+		this->value[3] = m[3];
 	}
 
 	template <typename T, precision P>
@@ -93,10 +93,10 @@ namespace detail
 		tmat4x4<T, Q> const & m
 	)
 	{
-		this->value[0] = m.value[0];
-		this->value[1] = m.value[1];
-		this->value[2] = m.value[2];
-		this->value[3] = m.value[3];
+		this->value[0] = m[0];
+		this->value[1] = m[1];
+		this->value[2] = m[2];
+		this->value[3] = m[3];
 	}
 
 	template <typename T, precision P>
@@ -461,7 +461,7 @@ namespace detail
 	template <typename U>
 	GLM_FUNC_QUALIFIER tmat4x4<T, P> & tmat4x4<T, P>::operator/= (tmat4x4<U, P> const & m)
 	{
-		return (*this = *this * detail::compute_inverse_mat4(m));
+		return (*this = *this * detail::compute_inverse<detail::tmat4x4, T, P>::call(m));
 	}
 
 	template <typename T, precision P>
@@ -501,52 +501,62 @@ namespace detail
 	}
 
 	template <typename T, precision P>
-	GLM_FUNC_QUALIFIER tmat4x4<T, P> compute_inverse_mat4(tmat4x4<T, P> const & m)
+	struct compute_inverse<detail::tmat4x4, T, P>
 	{
-		// Calculate all mat2 determinants
-		T const SubFactor00 = m[2][2] * m[3][3] - m[3][2] * m[2][3];
-		T const SubFactor01 = m[2][1] * m[3][3] - m[3][1] * m[2][3];
-		T const SubFactor02 = m[2][1] * m[3][2] - m[3][1] * m[2][2];
-		T const SubFactor03 = m[2][0] * m[3][3] - m[3][0] * m[2][3];
-		T const SubFactor04 = m[2][0] * m[3][2] - m[3][0] * m[2][2];
-		T const SubFactor05 = m[2][0] * m[3][1] - m[3][0] * m[2][1];
-		T const SubFactor06 = m[1][2] * m[3][3] - m[3][2] * m[1][3];
-		T const SubFactor07 = m[1][1] * m[3][3] - m[3][1] * m[1][3];
-		T const SubFactor08 = m[1][1] * m[3][2] - m[3][1] * m[1][2];
-		T const SubFactor09 = m[1][0] * m[3][3] - m[3][0] * m[1][3];
-		T const SubFactor10 = m[1][0] * m[3][2] - m[3][0] * m[1][2];
-		T const SubFactor11 = m[1][1] * m[3][3] - m[3][1] * m[1][3];
-		T const SubFactor12 = m[1][0] * m[3][1] - m[3][0] * m[1][1];
-		T const SubFactor13 = m[1][2] * m[2][3] - m[2][2] * m[1][3];
-		T const SubFactor14 = m[1][1] * m[2][3] - m[2][1] * m[1][3];
-		T const SubFactor15 = m[1][1] * m[2][2] - m[2][1] * m[1][2];
-		T const SubFactor16 = m[1][0] * m[2][3] - m[2][0] * m[1][3];
-		T const SubFactor17 = m[1][0] * m[2][2] - m[2][0] * m[1][2];
-		T const SubFactor18 = m[1][0] * m[2][1] - m[2][0] * m[1][1];
+		static detail::tmat4x4<T, lowp> call(detail::tmat4x4<T, P> const & m)
+		{
+			T Coef00 = m[2][2] * m[3][3] - m[3][2] * m[2][3];
+			T Coef02 = m[1][2] * m[3][3] - m[3][2] * m[1][3];
+			T Coef03 = m[1][2] * m[2][3] - m[2][2] * m[1][3];
 
-		tmat4x4<T, P> Inverse(
-			+ m[1][1] * SubFactor00 - m[1][2] * SubFactor01 + m[1][3] * SubFactor02,
-			- m[1][0] * SubFactor00 + m[1][2] * SubFactor03 - m[1][3] * SubFactor04,
-			+ m[1][0] * SubFactor01 - m[1][1] * SubFactor03 + m[1][3] * SubFactor05,
-			- m[1][0] * SubFactor02 + m[1][1] * SubFactor04 - m[1][2] * SubFactor05,
-			- m[0][1] * SubFactor00 + m[0][2] * SubFactor01 - m[0][3] * SubFactor02,
-			+ m[0][0] * SubFactor00 - m[0][2] * SubFactor03 + m[0][3] * SubFactor04,
-			- m[0][0] * SubFactor01 + m[0][1] * SubFactor03 - m[0][3] * SubFactor05,
-			+ m[0][0] * SubFactor02 - m[0][1] * SubFactor04 + m[0][2] * SubFactor05,
-			+ m[0][1] * SubFactor06 - m[0][2] * SubFactor07 + m[0][3] * SubFactor08,
-			- m[0][0] * SubFactor06 + m[0][2] * SubFactor09 - m[0][3] * SubFactor10,
-			+ m[0][0] * SubFactor11 - m[0][1] * SubFactor09 + m[0][3] * SubFactor12,
-			- m[0][0] * SubFactor08 + m[0][1] * SubFactor10 - m[0][2] * SubFactor12,
-			- m[0][1] * SubFactor13 + m[0][2] * SubFactor14 - m[0][3] * SubFactor15,
-			+ m[0][0] * SubFactor13 - m[0][2] * SubFactor16 + m[0][3] * SubFactor17,
-			- m[0][0] * SubFactor14 + m[0][1] * SubFactor16 - m[0][3] * SubFactor18,
-			+ m[0][0] * SubFactor15 - m[0][1] * SubFactor17 + m[0][2] * SubFactor18);
-		
-		T Determinant = static_cast<T>(+ m[0][0] * Inverse[0][0] + m[0][1] * Inverse[1][0] + m[0][2] * Inverse[2][0] + m[0][3] * Inverse[3][0]);
+			T Coef04 = m[2][1] * m[3][3] - m[3][1] * m[2][3];
+			T Coef06 = m[1][1] * m[3][3] - m[3][1] * m[1][3];
+			T Coef07 = m[1][1] * m[2][3] - m[2][1] * m[1][3];
 
-		Inverse /= Determinant;
-		return Inverse;
-	}
+			T Coef08 = m[2][1] * m[3][2] - m[3][1] * m[2][2];
+			T Coef10 = m[1][1] * m[3][2] - m[3][1] * m[1][2];
+			T Coef11 = m[1][1] * m[2][2] - m[2][1] * m[1][2];
+
+			T Coef12 = m[2][0] * m[3][3] - m[3][0] * m[2][3];
+			T Coef14 = m[1][0] * m[3][3] - m[3][0] * m[1][3];
+			T Coef15 = m[1][0] * m[2][3] - m[2][0] * m[1][3];
+
+			T Coef16 = m[2][0] * m[3][2] - m[3][0] * m[2][2];
+			T Coef18 = m[1][0] * m[3][2] - m[3][0] * m[1][2];
+			T Coef19 = m[1][0] * m[2][2] - m[2][0] * m[1][2];
+
+			T Coef20 = m[2][0] * m[3][1] - m[3][0] * m[2][1];
+			T Coef22 = m[1][0] * m[3][1] - m[3][0] * m[1][1];
+			T Coef23 = m[1][0] * m[2][1] - m[2][0] * m[1][1];
+
+			detail::tvec4<T, P> Fac0(Coef00, Coef00, Coef02, Coef03);
+			detail::tvec4<T, P> Fac1(Coef04, Coef04, Coef06, Coef07);
+			detail::tvec4<T, P> Fac2(Coef08, Coef08, Coef10, Coef11);
+			detail::tvec4<T, P> Fac3(Coef12, Coef12, Coef14, Coef15);
+			detail::tvec4<T, P> Fac4(Coef16, Coef16, Coef18, Coef19);
+			detail::tvec4<T, P> Fac5(Coef20, Coef20, Coef22, Coef23);
+
+			detail::tvec4<T, P> Vec0(m[1][0], m[0][0], m[0][0], m[0][0]);
+			detail::tvec4<T, P> Vec1(m[1][1], m[0][1], m[0][1], m[0][1]);
+			detail::tvec4<T, P> Vec2(m[1][2], m[0][2], m[0][2], m[0][2]);
+			detail::tvec4<T, P> Vec3(m[1][3], m[0][3], m[0][3], m[0][3]);
+
+			detail::tvec4<T, P> Inv0(Vec1 * Fac0 - Vec2 * Fac1 + Vec3 * Fac2);
+			detail::tvec4<T, P> Inv1(Vec0 * Fac0 - Vec2 * Fac3 + Vec3 * Fac4);
+			detail::tvec4<T, P> Inv2(Vec0 * Fac1 - Vec1 * Fac3 + Vec3 * Fac5);
+			detail::tvec4<T, P> Inv3(Vec0 * Fac2 - Vec1 * Fac4 + Vec2 * Fac5);
+
+			detail::tvec4<T, P> SignA(+1, -1, +1, -1);
+			detail::tvec4<T, P> SignB(-1, +1, -1, +1);
+			detail::tmat4x4<T, P> Inverse(Inv0 * SignA, Inv1 * SignB, Inv2 * SignA, Inv3 * SignB);
+
+			detail::tvec4<T, P> Row0(Inverse[0][0], Inverse[1][0], Inverse[2][0], Inverse[3][0]);
+
+			T OneOverDeterminant = static_cast<T>(1) / dot(m[0], Row0);
+
+			return Inverse * OneOverDeterminant;
+		}
+	};
 
 	// Binary operators
 	template <typename T, precision P>
@@ -823,7 +833,7 @@ namespace detail
 		typename tmat4x4<T, P>::row_type const & v
 	)
 	{
-		return detail::compute_inverse_mat4(m) * v;
+		return detail::compute_inverse<detail::tmat4x4, T, P>::call(m) * v;
 	}
 
 	template <typename T, precision P>
@@ -833,7 +843,7 @@ namespace detail
 		tmat4x4<T, P> const & m
 	)
 	{
-		return v * detail::compute_inverse_mat4(m);
+		return v * detail::compute_inverse<detail::tmat4x4, T, P>::call(m);
 	}
 
 	template <typename T, precision P>
diff --git a/glm/gtc/matrix_transform.inl b/glm/gtc/matrix_transform.inl
index 9bec7091..78c9ce85 100644
--- a/glm/gtc/matrix_transform.inl
+++ b/glm/gtc/matrix_transform.inl
@@ -82,9 +82,8 @@ namespace glm
 		T c = cos(a);
 		T s = sin(a);
 
-		detail::tvec3<T, P> axis = normalize(v);
-
-		detail::tvec3<T, P> temp = (T(1) - c) * axis;
+		detail::tvec3<T, P> axis(normalize(v));
+		detail::tvec3<T, P> temp((T(1) - c) * axis);
 
 		detail::tmat4x4<T, P> Rotate(detail::tmat4x4<T, P>::_null);
 		Rotate[0][0] = c + temp[0] * axis[0];
diff --git a/readme.txt b/readme.txt
index b515bf21..d73c18bb 100644
--- a/readme.txt
+++ b/readme.txt
@@ -45,6 +45,7 @@ GLM 0.9.5.1: 2014-XX-XX
 - Fixed error 'inverse' is not a member of 'glm' from glm::unProject (#146)
 - Fixed mismatch between some declarations and definitions
 - Fixed inverse link error when using namespace glm; (#147)
+- Optimized matrix inverse and division code (#149)
 
 ================================================================================
 GLM 0.9.5.0: 2013-12-25
diff --git a/test/core/core_func_matrix.cpp b/test/core/core_func_matrix.cpp
index ebe33d5f..1f117956 100644
--- a/test/core/core_func_matrix.cpp
+++ b/test/core/core_func_matrix.cpp
@@ -8,6 +8,11 @@
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 
 #include <glm/matrix.hpp>
+#include <glm/gtc/matrix_transform.hpp>
+#include <glm/gtc/ulp.hpp>
+#include <vector>
+#include <ctime>
+#include <cstdio>
 
 using namespace glm;
 
@@ -175,18 +180,71 @@ int test_inverse()
 	glm::mat2x2 I2x2 = A2x2 * B2x2;
 	Failed += I2x2 == glm::mat2x2(1) ? 0 : 1;
 
+
+
 	return Failed;
 }
 
+std::size_t const Count(10000000);
+
+template <typename VEC3, typename MAT4>
+int test_inverse_perf(std::size_t Instance, char const * Message)
+{
+	std::vector<MAT4> TestInputs;
+	TestInputs.resize(Count);
+	std::vector<MAT4> TestOutputs;
+	TestOutputs.resize(TestInputs.size());
+
+	VEC3 Axis(glm::normalize(VEC3(1.0f, 2.0f, 3.0f)));
+
+	for(std::size_t i = 0; i < TestInputs.size(); ++i)
+	{
+		typename MAT4::value_type f = static_cast<typename MAT4::value_type>(i + Instance) * typename MAT4::value_type(0.1) + typename MAT4::value_type(0.1);
+		TestInputs[i] = glm::rotate(glm::translate(MAT4(1), Axis * f), f, Axis);
+		//TestInputs[i] = glm::translate(MAT4(1), Axis * f);
+	}
+
+	std::clock_t StartTime = std::clock();
+
+	for(std::size_t i = 0; i < TestInputs.size(); ++i)
+		TestOutputs[i] = glm::inverse(TestInputs[i]);
+
+	std::clock_t EndTime = std::clock();
+
+	for(std::size_t i = 0; i < TestInputs.size(); ++i)
+		TestOutputs[i] = TestOutputs[i] * TestInputs[i];
+
+	typename MAT4::value_type Diff(0);
+	for(std::size_t Entry = 0; Entry < TestOutputs.size(); ++Entry)
+	{
+		MAT4 i(1.0);
+		MAT4 m(TestOutputs[Entry]);
+		for(glm::length_t y = 0; y < m.length(); ++y)
+		for(glm::length_t x = 0; x < m[y].length(); ++x)
+			Diff = glm::max(m[y][x], i[y][x]);
+	}
+
+	//glm::uint Ulp = 0;
+	//Ulp = glm::max(glm::float_distance(*Dst, *Src), Ulp);
+
+	printf("inverse<%s>(%f): %d\n", Message, Diff, EndTime - StartTime);
+
+	return 0;
+};
 
 int main()
 {
-	int Failed = 0;
-	Failed += test_matrixCompMult();
-	Failed += test_outerProduct();
-	Failed += test_transpose();
-	Failed += test_determinant();
-	Failed += test_inverse();
-	return Failed;
+	int Error(0);
+	Error += test_matrixCompMult();
+	Error += test_outerProduct();
+	Error += test_transpose();
+	Error += test_determinant();
+	Error += test_inverse();
+	for(std::size_t i = 0; i < 1; ++i)
+	{
+		Error += test_inverse_perf<glm::vec3, glm::mat4>(i, "mat4");
+		Error += test_inverse_perf<glm::dvec3, glm::dmat4>(i, "dmat4");
+	}
+	return Error;
 }
 
diff --git a/test/core/core_type_mat4x4.cpp b/test/core/core_type_mat4x4.cpp
index e4a32638..6c6abc2f 100644
--- a/test/core/core_type_mat4x4.cpp
+++ b/test/core/core_type_mat4x4.cpp
@@ -16,19 +16,19 @@
 void print(glm::dmat4 const & Mat0)
 {
 	printf("mat4(\n");
-	printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[0][0], Mat0[0][1], Mat0[0][2], Mat0[0][3]);
-	printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[1][0], Mat0[1][1], Mat0[1][2], Mat0[1][3]);
-	printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[2][0], Mat0[2][1], Mat0[2][2], Mat0[2][3]);
-	printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f))\n\n", Mat0[3][0], Mat0[3][1], Mat0[3][2], Mat0[3][3]);
+	printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[0][0], Mat0[0][1], Mat0[0][2], Mat0[0][3]);
+	printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[1][0], Mat0[1][1], Mat0[1][2], Mat0[1][3]);
+	printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[2][0], Mat0[2][1], Mat0[2][2], Mat0[2][3]);
+	printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f))\n\n", Mat0[3][0], Mat0[3][1], Mat0[3][2], Mat0[3][3]);
 }
 
 void print(glm::mat4 const & Mat0)
 {
 	printf("mat4(\n");
-	printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[0][0], Mat0[0][1], Mat0[0][2], Mat0[0][3]);
-	printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[1][0], Mat0[1][1], Mat0[1][2], Mat0[1][3]);
-	printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[2][0], Mat0[2][1], Mat0[2][2], Mat0[2][3]);
-	printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f))\n\n", Mat0[3][0], Mat0[3][1], Mat0[3][2], Mat0[3][3]);
+	printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[0][0], Mat0[0][1], Mat0[0][2], Mat0[0][3]);
+	printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[1][0], Mat0[1][1], Mat0[1][2], Mat0[1][3]);
+	printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[2][0], Mat0[2][1], Mat0[2][2], Mat0[2][3]);
+	printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f))\n\n", Mat0[3][0], Mat0[3][1], Mat0[3][2], Mat0[3][3]);
 }
 
 int test_inverse_mat4x4()
@@ -107,6 +107,66 @@ int test_inverse()
 		Error += glm::all(glm::epsilonEqual(Identity[3], glm::vec4(0.0f, 0.0f, 0.0f, 1.0f), glm::vec4(0.01f))) ? 0 : 1;
 	}
 
+	{
+		glm::highp_mat4 const Matrix(
+			glm::highp_vec4(0.6f, 0.2f, 0.3f, 0.4f), 
+			glm::highp_vec4(0.2f, 0.7f, 0.5f, 0.3f), 
+			glm::highp_vec4(0.3f, 0.5f, 0.7f, 0.2f), 
+			glm::highp_vec4(0.4f, 0.3f, 0.2f, 0.6f));
+		glm::highp_mat4 const Inverse = glm::inverse(Matrix);
+		glm::highp_mat4 const Identity = Matrix * Inverse;
+
+		printf("highp_mat4 inverse\n");
+		print(Matrix);
+		print(Inverse);
+		print(Identity);
+
+		Error += glm::all(glm::epsilonEqual(Identity[0], glm::highp_vec4(1.0f, 0.0f, 0.0f, 0.0f), glm::highp_vec4(0.01f))) ? 0 : 1;
+		Error += glm::all(glm::epsilonEqual(Identity[1], glm::highp_vec4(0.0f, 1.0f, 0.0f, 0.0f), glm::highp_vec4(0.01f))) ? 0 : 1;
+		Error += glm::all(glm::epsilonEqual(Identity[2], glm::highp_vec4(0.0f, 0.0f, 1.0f, 0.0f), glm::highp_vec4(0.01f))) ? 0 : 1;
+		Error += glm::all(glm::epsilonEqual(Identity[3], glm::highp_vec4(0.0f, 0.0f, 0.0f, 1.0f), glm::highp_vec4(0.01f))) ? 0 : 1;
+	}
+
+	{
+		glm::mediump_mat4 const Matrix(
+			glm::mediump_vec4(0.6f, 0.2f, 0.3f, 0.4f), 
+			glm::mediump_vec4(0.2f, 0.7f, 0.5f, 0.3f), 
+			glm::mediump_vec4(0.3f, 0.5f, 0.7f, 0.2f), 
+			glm::mediump_vec4(0.4f, 0.3f, 0.2f, 0.6f));
+		glm::mediump_mat4 const Inverse = glm::inverse(Matrix);
+		glm::mediump_mat4 const Identity = Matrix * Inverse;
+
+		printf("mediump_mat4 inverse\n");
+		print(Matrix);
+		print(Inverse);
+		print(Identity);
+
+		Error += glm::all(glm::epsilonEqual(Identity[0], glm::mediump_vec4(1.0f, 0.0f, 0.0f, 0.0f), glm::mediump_vec4(0.01f))) ? 0 : 1;
+		Error += glm::all(glm::epsilonEqual(Identity[1], glm::mediump_vec4(0.0f, 1.0f, 0.0f, 0.0f), glm::mediump_vec4(0.01f))) ? 0 : 1;
+		Error += glm::all(glm::epsilonEqual(Identity[2], glm::mediump_vec4(0.0f, 0.0f, 1.0f, 0.0f), glm::mediump_vec4(0.01f))) ? 0 : 1;
+		Error += glm::all(glm::epsilonEqual(Identity[3], glm::mediump_vec4(0.0f, 0.0f, 0.0f, 1.0f), glm::mediump_vec4(0.01f))) ? 0 : 1;
+	}
+
+	{
+		glm::lowp_mat4 const Matrix(
+			glm::lowp_vec4(0.6f, 0.2f, 0.3f, 0.4f), 
+			glm::lowp_vec4(0.2f, 0.7f, 0.5f, 0.3f), 
+			glm::lowp_vec4(0.3f, 0.5f, 0.7f, 0.2f), 
+			glm::lowp_vec4(0.4f, 0.3f, 0.2f, 0.6f));
+		glm::lowp_mat4 const Inverse = glm::inverse(Matrix);
+		glm::lowp_mat4 const Identity = Matrix * Inverse;
+
+		printf("lowp_mat4 inverse\n");
+		print(Matrix);
+		print(Inverse);
+		print(Identity);
+
+		Error += glm::all(glm::epsilonEqual(Identity[0], glm::lowp_vec4(1.0f, 0.0f, 0.0f, 0.0f), glm::lowp_vec4(0.01f))) ? 0 : 1;
+		Error += glm::all(glm::epsilonEqual(Identity[1], glm::lowp_vec4(0.0f, 1.0f, 0.0f, 0.0f), glm::lowp_vec4(0.01f))) ? 0 : 1;
+		Error += glm::all(glm::epsilonEqual(Identity[2], glm::lowp_vec4(0.0f, 0.0f, 1.0f, 0.0f), glm::lowp_vec4(0.01f))) ? 0 : 1;
+		Error += glm::all(glm::epsilonEqual(Identity[3], glm::lowp_vec4(0.0f, 0.0f, 0.0f, 1.0f), glm::lowp_vec4(0.01f))) ? 0 : 1;
+	}
+
 	{
 		glm::mat4 const Matrix(
 			glm::vec4(0.6f, 0.2f, 0.3f, 0.4f),