Math: better algorithm for comparing floating-point values.

Still not ideal (for values around zero), but can't do much better for
such general case.

src/Magnum/Math/Algorithms/Test/GramSchmidtTest.cpp

@@ -90,9 +90,9 @@ void GramSchmidtTest::orthonormalize() {
     CORRADE_COMPARE(Vector3::dot(orthonormalized[1], orthonormalized[2]), 0.0f);
 
     /* Just to be sure */
-    Matrix3x3 expected(Vector3( 0.303046f,  0.505076f,  0.808122f),
+    Matrix3x3 expected(Vector3( 0.3030458f,  0.5050763f,  0.8081220f),
-                       Vector3( 0.928316f, -0.348119f, -0.130544f),
+                       Vector3( 0.9283164f, -0.3481189f, -0.1305445f),
-                       Vector3(-0.215388f, -0.789754f,  0.574367f));
+                       Vector3(-0.2153877f, -0.7897540f,  0.5743665f));
     CORRADE_COMPARE(orthonormalized, expected);
 }
 

src/Magnum/Math/Test/Matrix4Test.cpp

@@ -263,10 +263,10 @@ void Matrix4Test::rotation() {
     CORRADE_COMPARE(Matrix4::rotation(Deg(-74.0f), {-1.0f, 2.0f, 2.0f}), Matrix4());
     CORRADE_COMPARE(o.str(), "Math::Matrix4::rotation(): axis must be normalized\n");
 
-    Matrix4 matrix({ 0.35612214f,  -0.80181062f, 0.47987163f, 0.0f},
+    Matrix4 matrix({ 0.35612202f, -0.80181062f, 0.47987163f, 0.0f},
-                   { 0.47987163f,   0.59757638f,  0.6423595f, 0.0f},
+                   { 0.47987163f,  0.59757626f,  0.6423596f, 0.0f},
-                   {-0.80181062f, 0.0015183985f, 0.59757638f, 0.0f},
+                   {-0.80181062f,  0.00151846f, 0.59757626f, 0.0f},
-                   {        0.0f,          0.0f,        0.0f, 1.0f});
+                   {        0.0f,         0.0f,        0.0f, 1.0f});
     CORRADE_COMPARE(Matrix4::rotation(Deg(-74.0f), Vector3(-1.0f, 2.0f, 2.0f).normalized()), matrix);
 }
 
@@ -388,9 +388,9 @@ void Matrix4Test::rotationNormalizedPart() {
 
 void Matrix4Test::rotationPart() {
     Matrix4 rotation = Matrix4::rotation(Deg(-74.0f), Vector3(-1.0f, 2.0f, 2.0f).normalized());
-    Matrix3x3 expectedRotationPart(Vector3( 0.35612214f,  -0.80181062f, 0.47987163f),
+    Matrix3x3 expectedRotationPart(Vector3( 0.35612206f, -0.80181074f, 0.47987169f),
-                                   Vector3( 0.47987163f,   0.59757638f,  0.6423595f),
+                                   Vector3( 0.47987163f,  0.59757626f, 0.64235962f),
-                                   Vector3(-0.80181062f, 0.0015183985f, 0.59757638f));
+                                   Vector3(-0.80181062f,  0.00151846f, 0.59757626f));
 
     /* For rotation and translation this is the same as rotationScaling() */
     Matrix4 rotationTranslation = rotation*Matrix4::translation({2.0f, 5.0f, -3.0f});

src/Magnum/Math/Test/QuaternionTest.cpp

@@ -347,7 +347,7 @@ void QuaternionTest::slerp() {
     CORRADE_COMPARE(o.str(), "Math::Quaternion::slerp(): quaternions must be normalized\n");
 
     Quaternion slerp = Quaternion::slerp(a, b, 0.35f);
-    CORRADE_COMPARE(slerp, Quaternion({0.119165f, 0.0491109f, 0.0491109f}, 0.990442f));
+    CORRADE_COMPARE(slerp, Quaternion({0.1191653f, 0.0491109f, 0.0491109f}, 0.9904423f));
 }
 
 void QuaternionTest::transformVector() {

src/Magnum/Math/Test/TypeTraitsTest.cpp

@@ -34,17 +34,29 @@ class TypeTraitsTest: public Corrade::TestSuite::Tester {
     public:
         TypeTraitsTest();
 
-        void equalsFloatingPoint();
         void equalsIntegral();
+        void equalsFloatingPoint0();
+        void equalsFloatingPoint1();
+        void equalsFloatingPointLarge();
+        void equalsFloatingPointInfinity();
+        void equalsFloatingPointNaN();
 
     private:
-        template<class T> void _equalsFloatingPoint();
+        template<class> void _equalsIntegral();
-        template<class T> void _equalsIntegral();
+        template<class> void _equalsFloatingPoint0();
+        template<class> void _equalsFloatingPoint1();
+        template<class> void _equalsFloatingPointLarge();
+        template<class> void _equalsFloatingPointInfinity();
+        template<class> void _equalsFloatingPointNaN();
 };
 
 TypeTraitsTest::TypeTraitsTest() {
     addTests({&TypeTraitsTest::equalsIntegral,
-              &TypeTraitsTest::equalsFloatingPoint});
+              &TypeTraitsTest::equalsFloatingPoint0,
+              &TypeTraitsTest::equalsFloatingPoint1,
+              &TypeTraitsTest::equalsFloatingPointLarge,
+              &TypeTraitsTest::equalsFloatingPointInfinity,
+              &TypeTraitsTest::equalsFloatingPointNaN});
 }
 
 void TypeTraitsTest::equalsIntegral() {
@@ -58,29 +70,68 @@ void TypeTraitsTest::equalsIntegral() {
     _equalsIntegral<Long>();
 }
 
-void TypeTraitsTest::equalsFloatingPoint() {
+template<class T> void TypeTraitsTest::_equalsIntegral() {
-    _equalsFloatingPoint<Float>();
+    CORRADE_VERIFY(!TypeTraits<T>::equals(1, 1+TypeTraits<T>::epsilon()));
+}
+
+void TypeTraitsTest::equalsFloatingPoint0() {
+    _equalsFloatingPoint0<Float>();
     #ifndef MAGNUM_TARGET_GLES
-    _equalsFloatingPoint<Double>();
+    _equalsFloatingPoint0<Double>();
     #endif
 }
 
-template<class T> void TypeTraitsTest::_equalsIntegral() {
+template<class T> void TypeTraitsTest::_equalsFloatingPoint0() {
-    CORRADE_VERIFY(!TypeTraits<T>::equals(1, 1+TypeTraits<T>::epsilon()));
+    CORRADE_VERIFY(TypeTraits<T>::equals(T(0)+TypeTraits<T>::epsilon()/T(2), T(0)));
+    CORRADE_VERIFY(!TypeTraits<T>::equals(T(0)+TypeTraits<T>::epsilon()*T(2), T(0)));
 }
 
-template<class T> void TypeTraitsTest::_equalsFloatingPoint() {
+void TypeTraitsTest::equalsFloatingPoint1() {
+    _equalsFloatingPoint1<Float>();
+    #ifndef MAGNUM_TARGET_GLES
+    _equalsFloatingPoint1<Double>();
+    #endif
+}
+
+template<class T> void TypeTraitsTest::_equalsFloatingPoint1() {
     CORRADE_VERIFY(TypeTraits<T>::equals(T(1)+TypeTraits<T>::epsilon()/T(2), T(1)));
-    CORRADE_VERIFY(!TypeTraits<T>::equals(T(1)+TypeTraits<T>::epsilon()*T(2), T(1)));
+    CORRADE_VERIFY(!TypeTraits<T>::equals(T(1)+TypeTraits<T>::epsilon()*T(3), T(1)));
+}
+
+void TypeTraitsTest::equalsFloatingPointLarge() {
+    _equalsFloatingPointLarge<Float>();
+    #ifndef MAGNUM_TARGET_GLES
+    _equalsFloatingPointLarge<Double>();
+    #endif
+}
+
+template<class T> void TypeTraitsTest::_equalsFloatingPointLarge() {
+    CORRADE_VERIFY(TypeTraits<T>::equals(T(25)+TypeTraits<T>::epsilon()*T(2), T(25)));
+    CORRADE_VERIFY(!TypeTraits<T>::equals(T(25)+TypeTraits<T>::epsilon()*T(75), T(25)));
+}
 
-    {
+void TypeTraitsTest::equalsFloatingPointInfinity() {
-        CORRADE_EXPECT_FAIL("Comparing to infinity is broken");
+    _equalsFloatingPointInfinity<Float>();
-        CORRADE_VERIFY(TypeTraits<T>::equals(std::numeric_limits<T>::infinity(),
+    #ifndef MAGNUM_TARGET_GLES
-                                                std::numeric_limits<T>::infinity()));
+    _equalsFloatingPointInfinity<Double>();
-    }
+    #endif
+}
+
+template<class T> void TypeTraitsTest::_equalsFloatingPointInfinity() {
+    CORRADE_VERIFY(TypeTraits<T>::equals(std::numeric_limits<T>::infinity(),
+                                         std::numeric_limits<T>::infinity()));
+}
+
+void TypeTraitsTest::equalsFloatingPointNaN() {
+    _equalsFloatingPointNaN<Float>();
+    #ifndef MAGNUM_TARGET_GLES
+    _equalsFloatingPointNaN<Double>();
+    #endif
+}
 
+template<class T> void TypeTraitsTest::_equalsFloatingPointNaN() {
     CORRADE_VERIFY(!TypeTraits<T>::equals(std::numeric_limits<T>::quiet_NaN(),
-                                              std::numeric_limits<T>::quiet_NaN()));
+                                          std::numeric_limits<T>::quiet_NaN()));
 }
 
 }}}

src/Magnum/Math/Test/UnitTest.cpp

@@ -86,8 +86,8 @@ void UnitTest::constructConversion() {
 }
 
 void UnitTest::compare() {
-    CORRADE_VERIFY(Sec(25.0f + TypeTraits<Float>::epsilon()/2) == Sec(25.0f));
+    CORRADE_VERIFY(Sec(25.0f + TypeTraits<Float>::epsilon()/2.0f) == Sec(25.0f));
-    CORRADE_VERIFY(Sec(25.0f + TypeTraits<Float>::epsilon()*2) != Sec(25.0f));
+    CORRADE_VERIFY(Sec(25.0f + TypeTraits<Float>::epsilon()*75.0f) != Sec(25.0f));
 
     constexpr bool c = Sec(3.0f) < Sec(3.0f);
     constexpr bool d = Sec(3.0f) <= Sec(3.0f);

src/Magnum/Math/TypeTraits.h

@@ -106,14 +106,6 @@ template<class T> struct TypeTraits: Implementation::TypeTraitsDefault<T> {
     #endif
 };
 
-/** @bug Infinity comparison! */
-
-/**
- * @todo Implement better fuzzy comparison algorithm, like at
- * http://floating-point-gui.de/errors/comparison/ or
- * http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
- */
-
 /* Integral scalar types */
 namespace Implementation {
     template<class T> struct TypeTraitsIntegral: TypeTraitsDefault<T> {
@@ -159,13 +151,31 @@ template<> struct TypeTraits<Long>: Implementation::TypeTraitsIntegral<Long> {
 
 /* Floating-point scalar types */
 namespace Implementation {
-    template<class T> struct TypeTraitsFloatingPoint {
-        TypeTraitsFloatingPoint() = delete;
 
-        static bool equals(T a, T b) {
+template<class T> struct TypeTraitsFloatingPoint {
-            return std::abs(a - b) < TypeTraits<T>::epsilon();
+    TypeTraitsFloatingPoint() = delete;
-        }
+
-    };
+    static bool equals(T a, T b);
+};
+
+/* Adapted from http://floating-point-gui.de/errors/comparison/ */
+template<class T> bool TypeTraitsFloatingPoint<T>::equals(const T a, const T b) {
+    /* Shortcut for binary equality (also infinites) */
+    if (a == b) return true;
+
+    const T absA = std::abs(a);
+    const T absB = std::abs(b);
+    const T difference = std::abs(a - b);
+
+    /* One of the numbers is zero or both are extremely close to it, relative
+       error is meaningless */
+    if (a == T{} || b == T{} || difference < TypeTraits<T>::epsilon())
+        return difference < TypeTraits<T>::epsilon();
+
+    /* Relative error */
+    return difference/(absA + absB) < TypeTraits<T>::epsilon();
+}
+
 }
 
 template<> struct TypeTraits<Float>: Implementation::TypeTraitsFloatingPoint<Float> {