From ca6fa4ad3efd94a138cb2bf0105902c1a378d980 Mon Sep 17 00:00:00 2001
From: Danny Robson <danny@nerdcruft.net>
Date: Wed, 21 Jan 2015 23:32:06 +1100
Subject: [PATCH] float: use ULP based comparison

---
 float.cpp | 56 +++++++++++++++++++++++++++++++++++++++++++++----------
 float.hpp |  2 ++
 2 files changed, 48 insertions(+), 10 deletions(-)

diff --git a/float.cpp b/float.cpp
index a2eccd84..500ecbb6 100644
--- a/float.cpp
+++ b/float.cpp
@@ -76,22 +76,58 @@ ieee_float<E, S>::operator==(floating_t _floating) const {
 #include <iostream>
 
 
+//template <unsigned int E, unsigned int S>
+//bool
+//ieee_float<E, S>::almost_equal (floating_t a,
+//                                floating_t b) { 
+//    // Static cast to avoid integer casting warnings when using uint16_t for half
+//    static const floating_t epsilon = static_cast<floating_t> (0.001);
+//    const floating_t diff = static_cast<floating_t> (std::fabs (a - b));
+//
+//    // * Use an exact equality first so that infinities are not indirectly compared. This would generate NaNs in the diff.
+//    // * Do not use gte or lte. This stops an infinite from making infinities on both sides of the inequality.
+//    return exactly_equal (a, b)           ||
+//           diff < epsilon * std::fabs (a) ||
+//           diff < epsilon * std::fabs (b);
+//}
+
+
 template <unsigned int E, unsigned int S>
 bool
 ieee_float<E, S>::almost_equal (floating_t a,
-                                floating_t b) { 
-    // Static cast to avoid integer casting warnings when using uint16_t for half
-    static const floating_t epsilon = static_cast<floating_t> (0.001);
-    const floating_t diff = static_cast<floating_t> (std::fabs (a - b));
-
-    // * Use an exact equality first so that infinities are not indirectly compared. This would generate NaNs in the diff.
-    // * Do not use gte or lte. This stops an infinite from making infinities on both sides of the inequality.
-    return exactly_equal (a, b)           ||
-           diff < epsilon * std::fabs (a) ||
-           diff < epsilon * std::fabs (b);
+                                floating_t b)
+{
+    return almost_equal (a, b, 10000);
 }
 
 
+template <unsigned int E, unsigned int S>
+bool
+ieee_float<E, S>::almost_equal (floating_t _a,
+                                floating_t _b,
+                                uint_t ulps)
+{
+    union {
+        floating_t f;
+        sint_t     s;
+    } a, b;
+
+    a.f = _a;
+    b.f = _b;
+
+    if (std::isnan (a.f) || std::isnan (b.f))
+        return false;
+
+    if (a.s == b.s)
+        return true;
+
+    uint_t diff = std::abs (a.s - b.s);
+    std::cerr << "diff: " << diff << '\n';
+    return diff <= ulps;
+}
+
+
+
 template class ieee_float< 5,  10>; // ieee_half
 template class ieee_float< 8,  23>; // ieee_single;
 template class ieee_float<11,  52>; // ieee_double;
diff --git a/float.hpp b/float.hpp
index e6167cd0..b3af11b6 100644
--- a/float.hpp
+++ b/float.hpp
@@ -14,6 +14,7 @@ class ieee_float {
 
         static const unsigned int BIAS = (1 << (EXPONENT - 1)) - 1;
 
+        typedef typename bits_type<TOTAL_BITS>::sint     sint_t;
         typedef typename bits_type<TOTAL_BITS>::uint     uint_t;
         typedef typename bits_type<TOTAL_BITS>::floating floating_t;
 
@@ -50,6 +51,7 @@ class ieee_float {
         bool operator==   (floating_t) const;
 
         static bool almost_equal (floating_t, floating_t);
+        static bool almost_equal (floating_t, floating_t, uint_t ulps);
 };