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

glm/ext/quaternion_exponential.inl

@@ -46,17 +46,34 @@ 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);
+		if(abs(x.w / magnitude) > static_cast<T>(.5))
+		{
+			//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
 
-		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));
+			//Prevent a division by 0 error later on
+			T VectorMagnitude = x.x * x.x + x.y * x.y + x.z * x.z;
+			if (VectorMagnitude == 0) {
+				//Equivalent to raising a real number to a power
+				return qua<T, Q>(pow(x.w, y), 0, 0, 0);
+			}
 
-		return qua<T, Q>(cos(NewAngle) * magnitude * Mag, x.x * Div * Mag, x.y * Div * Mag, x.z * Div * Mag);
+			T SinAngle = sqrt(VectorMagnitude) / magnitude;
+			T Angle = asin(SinAngle);
+			T NewAngle = Angle * y;
+			T Div = sin(NewAngle) / SinAngle;
+			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);
+		}
+		else
+		{
+			//Scalar component is small, shouldn't cause loss of precision
+			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);
+		}
 	}
 
 	template<typename T, qualifier Q>