Merge pull request #946 from qsantos/patch-1

Fix loss of precision on small angles in qua's pow #946

glm/ext/quaternion_exponential.inl

@@ -1,3 +1,5 @@
+#include "scalar_constants.hpp"
+
 namespace glm
 {
 	template<typename T, qualifier Q>
@@ -46,16 +48,30 @@ namespace glm
 		//To deal with non-unit quaternions
 		T magnitude = sqrt(x.x * x.x + x.y * x.y + x.z * x.z + x.w *x.w);
 
-		//Equivalent to raising a real number to a power
-		//Needed to prevent a division by 0 error later on
-		if(abs(x.w / magnitude) > static_cast<T>(1) - epsilon<T>() && abs(x.w / magnitude) < static_cast<T>(1) + epsilon<T>())
-			return qua<T, Q>(pow(x.w, y), 0, 0, 0);
+		T Angle;
+		if(abs(x.w / magnitude) > cos_one_over_two<T>())
+		{
+			//Scalar component is close to 1; using it to recover angle would lose precision
+			//Instead, we use the non-scalar components since sin() is accurate around 0
+
+			//Prevent a division by 0 error later on
+			T VectorMagnitude = x.x * x.x + x.y * x.y + x.z * x.z;
+			if (glm::abs(VectorMagnitude - static_cast<T>(0)) < glm::epsilon<T>()) {
+				//Equivalent to raising a real number to a power
+				return qua<T, Q>(pow(x.w, y), 0, 0, 0);
+			}
+
+			Angle = asin(sqrt(VectorMagnitude) / magnitude);
+		}
+		else
+		{
+			//Scalar component is small, shouldn't cause loss of precision
+			Angle = acos(x.w / magnitude);
+		}
 
-		T Angle = acos(x.w / magnitude);
 		T NewAngle = Angle * y;
 		T Div = sin(NewAngle) / sin(Angle);
 		T Mag = pow(magnitude, y - static_cast<T>(1));
-
 		return qua<T, Q>(cos(NewAngle) * magnitude * Mag, x.x * Div * Mag, x.y * Div * Mag, x.z * Div * Mag);
 	}
 

glm/ext/quaternion_trigonometric.inl

@@ -1,8 +1,15 @@
+#include "scalar_constants.hpp"
+
 namespace glm
 {
 	template<typename T, qualifier Q>
 	GLM_FUNC_QUALIFIER T angle(qua<T, Q> const& x)
 	{
+		if (abs(x.w) > cos_one_over_two<T>())
+		{
+			return asin(sqrt(x.x * x.x + x.y * x.y + x.z * x.z)) * static_cast<T>(2);
+		}
+
 		return acos(x.w) * static_cast<T>(2);
 	}
 

glm/ext/scalar_constants.hpp

@@ -30,6 +30,10 @@ namespace glm
 	template<typename genType>
 	GLM_FUNC_DECL GLM_CONSTEXPR genType pi();
 
+	/// Return the value of cos(1 / 2) for floating point types.
+	template<typename genType>
+	GLM_FUNC_DECL GLM_CONSTEXPR genType cos_one_over_two();
+
 	/// @}
 } //namespace glm
 

glm/ext/scalar_constants.inl

@@ -15,4 +15,10 @@ namespace glm
 		GLM_STATIC_ASSERT(std::numeric_limits<genType>::is_iec559, "'pi' only accepts floating-point inputs");
 		return static_cast<genType>(3.14159265358979323846264338327950288);
 	}
+
+	template<typename genType>
+	GLM_FUNC_QUALIFIER GLM_CONSTEXPR genType cos_one_over_two()
+	{
+		return genType(0.877582561890372716130286068203503191);
+	}
 } //namespace glm

glm/gtc/constants.inl

@@ -163,4 +163,5 @@ namespace glm
 	{
 		return genType(1.61803398874989484820458683436563811);
 	}
+
 } //namespace glm