diff --git a/glm/gtx/matrix_interpolation.inl b/glm/gtx/matrix_interpolation.inl
index 10972974..f4ba3a6f 100644
--- a/glm/gtx/matrix_interpolation.inl
+++ b/glm/gtx/matrix_interpolation.inl
@@ -2,17 +2,29 @@
 
 #include "../ext/scalar_constants.hpp"
 
+#include <limits>
+
 namespace glm
 {
 	template<typename T, qualifier Q>
 	GLM_FUNC_QUALIFIER void axisAngle(mat<4, 4, T, Q> const& m, vec<3, T, Q>& axis, T& angle)
 	{
-		T epsilon = static_cast<T>(0.01);
-		T epsilon2 = static_cast<T>(0.1);
+		T const epsilon =
+		    std::numeric_limits<T>::epsilon() * static_cast<T>(1e2);
 
-		if((abs(m[1][0] - m[0][1]) < epsilon) && (abs(m[2][0] - m[0][2]) < epsilon) && (abs(m[2][1] - m[1][2]) < epsilon))
+        bool const nearSymmetrical =
+            abs(m[1][0] - m[0][1]) < epsilon &&
+            abs(m[2][0] - m[0][2]) < epsilon &&
+            abs(m[2][1] - m[1][2]) < epsilon;
+
+		if(nearSymmetrical)
 		{
-			if ((abs(m[1][0] + m[0][1]) < epsilon2) && (abs(m[2][0] + m[0][2]) < epsilon2) && (abs(m[2][1] + m[1][2]) < epsilon2) && (abs(m[0][0] + m[1][1] + m[2][2] - static_cast<T>(3.0)) < epsilon2))
+            bool const nearIdentity =
+                abs(m[1][0] + m[0][1]) < epsilon &&
+                abs(m[2][0] + m[0][2]) < epsilon &&
+                abs(m[2][1] + m[1][2]) < epsilon &&
+                abs(m[0][0] + m[1][1] + m[2][2] - T(3.0)) < epsilon;
+			if (nearIdentity)
 			{
 				angle = static_cast<T>(0.0);
 				axis = vec<3, T, Q>(