math: fp_compare: Adapt fp_compare_epsilon API to actual use

UHD had an issue where the design of fp_compare_epsilon and its usage
differed. In fact, the *only* usage of fp_compare_epsilon outside of
unit tests was to do a fuzzy frequency comparison, and it always took
a form like this:

```cpp
// The argument EPSILON may be implied, i.e., using the default
if (fp_compare_epsilon<double>(test_freq, EPSILON) < boundary_freq) {
    // ...
}
```

However, the API of fp_compare_epsilon was such that it would apply
DOUBLE_PRECISION_EPSILON to part of the frequency comparison, thus
rendering the argument EPSILON obsolete. When the default EPSILON was
used, this was OK, but only when the floating point type of
fp_compare_epsilon<> was `double`, and not `float`.

As an example, consider the following:
```
if (fp_compare_epsilon<double>(1e9 + x, LITTLE_EPSILON) == 1e9) {
    // ....
}
double BIG_EPSILON = x * 10;
if (fp_compare_epsilon<double>(1e9 + x, BIG_EPSILON) == 1e9) {
    // ....
}
```

If you expect the second comparison to pass even if the first failed,
then you are not alone. However, that's not what UHD would do. Because
of the aforementioned behaviour, it would use DOUBLE_PRECISION_EPSILON
for the right hand comparison, which would fail again.

Instead of fixing the instances of fp_compare_epsilon throughout UHD,
this patch changes the comparison algorithm from "very close with
tolerance epsilon" to "close enough with tolerance epsilon". This
requires only one side to be close to the other, using its own epsilon,
so the aforementioned example would always pass on the second check.

However, this exposed a second bug in fp_compare_epsilon. For
greater-/less-than comparisons, it would use epsilon like a delta value,
i.e., it would check if

    a + epsilon < b - epsilon

That means that if a < b, but (b-a) < 2*epsilon, this check would return
"false", i.e., it would report that a >= b, which is incorrect. These
operators are now changed such that they first check equality of a and
b using the algorithm described in the code, and then compare the values
of a and b (ignoring epsilon) directly. A unit test for this case was
added.

2 files changed, 40 insertions, 32 deletions

diff --git a/host/include/uhd/utils/fp_compare_epsilon.ipp b/host/include/uhd/utils/fp_compare_epsilon.ipp
index 99e467df4..a00cb0696 100644
--- a/host/include/uhd/utils/fp_compare_epsilon.ipp
+++ b/host/include/uhd/utils/fp_compare_epsilon.ipp
@@ -9,6 +9,7 @@
 #include <cmath>
 #include <typeinfo>
 
+
 #pragma once
 
 
@@ -49,15 +50,32 @@ namespace uhd { namespace math { namespace fp_compare {
         _epsilon = copy._epsilon;
     }
 
-    template<typename float_t> UHD_INLINE
-    bool operator==(fp_compare_epsilon<float_t> lhs, fp_compare_epsilon<float_t> rhs) {
-
-        bool lhs_compare = ((std::abs(lhs._value - rhs._value) / std::abs(lhs._value))
+    template <typename float_t>
+    UHD_INLINE bool operator==(
+        fp_compare_epsilon<float_t> lhs, fp_compare_epsilon<float_t> rhs)
+    {
+        // If raw bit values are equal, then they're equal. This also catches
+        // the case where both are 0.0!
+        if (lhs._value == rhs._value) {
+            return true;
+        }
+
+        // If one of them is within epsilon of zero, but the other is not, then
+        // they're also not equal.
+        if ((std::abs(lhs._value) <= lhs._epsilon && std::abs(rhs._value) > rhs._epsilon)
+            || (std::abs(lhs._value) > lhs._epsilon
+                && std::abs(rhs._value) <= rhs._epsilon)) {
+            return false;
+        }
+
+        // In all other cases, we use the "close enough with tolerance epsilon"
+        // algorithm as described in math.hpp.
+        const bool lhs_compare = ((std::abs(lhs._value - rhs._value) / std::abs(lhs._value))
                     <= lhs._epsilon);
-        bool rhs_compare = ((std::abs(lhs._value - rhs._value) / std::abs(rhs._value))
+        const bool rhs_compare = ((std::abs(lhs._value - rhs._value) / std::abs(rhs._value))
                     <= rhs._epsilon);
 
-        return (lhs_compare && rhs_compare);
+        return (lhs_compare || rhs_compare);
     }
 
     template<typename float_t> UHD_INLINE
@@ -67,7 +85,7 @@ namespace uhd { namespace math { namespace fp_compare {
 
     template<typename float_t> UHD_INLINE
     bool operator<(fp_compare_epsilon<float_t> lhs, fp_compare_epsilon<float_t> rhs) {
-        return (lhs._value + lhs._epsilon) <  (rhs._value - rhs._epsilon);
+        return (lhs != rhs) && (lhs._value < rhs._value);
     }
 
     template<typename float_t> UHD_INLINE
@@ -77,7 +95,7 @@ namespace uhd { namespace math { namespace fp_compare {
 
     template<typename float_t> UHD_INLINE
     bool operator>(fp_compare_epsilon<float_t> lhs, fp_compare_epsilon<float_t> rhs) {
-        return (lhs._value - lhs._epsilon) > (rhs._value + rhs._epsilon);
+        return (lhs != rhs) && (lhs._value > rhs._value);
     }
 
     template<typename float_t> UHD_INLINE
@@ -85,15 +103,10 @@ namespace uhd { namespace math { namespace fp_compare {
         return !(lhs < rhs);
     }
 
-    template<typename float_t> UHD_INLINE
-    bool operator==(fp_compare_epsilon<float_t> lhs, double rhs) {
-
-        bool lhs_compare = ((std::abs(lhs._value - rhs) / std::abs(lhs._value))
-                    <= lhs._epsilon);
-        bool rhs_compare = ((std::abs(lhs._value - rhs) / std::abs(rhs))
-                    <= DOUBLE_PRECISION_EPSILON);
-
-        return (lhs_compare && rhs_compare);
+    template <typename float_t>
+    UHD_INLINE bool operator==(fp_compare_epsilon<float_t> lhs, double rhs)
+    {
+        return lhs == fp_compare_epsilon<float_t>(static_cast<float_t>(rhs));
     }
 
     template<typename float_t> UHD_INLINE
@@ -104,7 +117,7 @@ namespace uhd { namespace math { namespace fp_compare {
     template<typename float_t> UHD_INLINE
     bool operator<(fp_compare_epsilon<float_t> lhs, double rhs) {
 
-        return (lhs._value + lhs._epsilon) <  (rhs - DOUBLE_PRECISION_EPSILON);
+        return (lhs != rhs) && (lhs._value < rhs);
     }
 
     template<typename float_t> UHD_INLINE
@@ -115,7 +128,7 @@ namespace uhd { namespace math { namespace fp_compare {
     template<typename float_t> UHD_INLINE
     bool operator>(fp_compare_epsilon<float_t> lhs, double rhs) {
 
-        return (lhs._value - lhs._epsilon) > (rhs + DOUBLE_PRECISION_EPSILON);
+        return (lhs != rhs) && (lhs._value > rhs);
     }
 
     template<typename float_t> UHD_INLINE
@@ -123,15 +136,10 @@ namespace uhd { namespace math { namespace fp_compare {
         return !(lhs < rhs);
     }
 
-    template<typename float_t> UHD_INLINE
-    bool operator==(double lhs, fp_compare_epsilon<float_t> rhs) {
-
-        bool lhs_compare = ((std::abs(lhs - rhs._value) / std::abs(lhs))
-                    <= DOUBLE_PRECISION_EPSILON);
-        bool rhs_compare = ((std::abs(lhs - rhs._value) / std::abs(rhs._value))
-                    <= rhs._epsilon);
-
-        return (lhs_compare && rhs_compare);
+    template <typename float_t>
+    UHD_INLINE bool operator==(double lhs, fp_compare_epsilon<float_t> rhs)
+    {
+        return fp_compare_epsilon<float_t>(static_cast<float_t>(lhs)) == rhs;
     }
 
     template<typename float_t> UHD_INLINE
@@ -142,7 +150,7 @@ namespace uhd { namespace math { namespace fp_compare {
     template<typename float_t> UHD_INLINE
     bool operator<(double lhs, fp_compare_epsilon<float_t> rhs) {
 
-        return (lhs + DOUBLE_PRECISION_EPSILON) <  (rhs._value - rhs._epsilon);
+        return (lhs != rhs) && (lhs < rhs._value);
     }
 
     template<typename float_t> UHD_INLINE
@@ -153,7 +161,7 @@ namespace uhd { namespace math { namespace fp_compare {
     template<typename float_t> UHD_INLINE
     bool operator>(double lhs, fp_compare_epsilon<float_t> rhs) {
 
-        return (lhs - DOUBLE_PRECISION_EPSILON) > (rhs._value + rhs._epsilon);
+        return (lhs != rhs) && (lhs > rhs._value);
     }
 
     template<typename float_t> UHD_INLINE
diff --git a/host/include/uhd/utils/math.hpp b/host/include/uhd/utils/math.hpp
index 6c8fceae9..6ee46e98a 100644
--- a/host/include/uhd/utils/math.hpp
+++ b/host/include/uhd/utils/math.hpp
@@ -85,10 +85,10 @@ public:
  * There are obviously a lot of strategies for defining floating point
  * equality, and in the end it all comes down to the domain at hand. UHD's
  * floating-point-with-epsilon comparison algorithm is based on the method
- * presented in Knuth's "The Art of Computer Science" called "very close
+ * presented in Knuth's "The Art of Computer Science" called "close enough
  * with tolerance epsilon".
  *
- *      [(|u - v| / |u|) <= e] && [(|u - v| / |v|) <= e]
+ *      [(|u - v| / |u|) <= e] || [(|u - v| / |v|) <= e]
  *
  * UHD's modification to this algorithm is using the denominator's epsilon
  * value (since each float_t object has its own epsilon) for each