From 3e496cbda1d809d2ca15f69cfa231424bf47179f Mon Sep 17 00:00:00 2001 From: Martin Braun Date: Wed, 12 Jan 2022 14:38:51 +0100 Subject: 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(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(1e9 + x, LITTLE_EPSILON) == 1e9) { // .... } double BIG_EPSILON = x * 10; if (fp_compare_epsilon(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. --- host/include/uhd/utils/fp_compare_epsilon.ipp | 68 +++++++++++++++------------ host/include/uhd/utils/math.hpp | 4 +- 2 files changed, 40 insertions(+), 32 deletions(-) (limited to 'host/include') 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 #include + #pragma once @@ -49,15 +50,32 @@ namespace uhd { namespace math { namespace fp_compare { _epsilon = copy._epsilon; } - template UHD_INLINE - bool operator==(fp_compare_epsilon lhs, fp_compare_epsilon rhs) { - - bool lhs_compare = ((std::abs(lhs._value - rhs._value) / std::abs(lhs._value)) + template + UHD_INLINE bool operator==( + fp_compare_epsilon lhs, fp_compare_epsilon 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 UHD_INLINE @@ -67,7 +85,7 @@ namespace uhd { namespace math { namespace fp_compare { template UHD_INLINE bool operator<(fp_compare_epsilon lhs, fp_compare_epsilon rhs) { - return (lhs._value + lhs._epsilon) < (rhs._value - rhs._epsilon); + return (lhs != rhs) && (lhs._value < rhs._value); } template UHD_INLINE @@ -77,7 +95,7 @@ namespace uhd { namespace math { namespace fp_compare { template UHD_INLINE bool operator>(fp_compare_epsilon lhs, fp_compare_epsilon rhs) { - return (lhs._value - lhs._epsilon) > (rhs._value + rhs._epsilon); + return (lhs != rhs) && (lhs._value > rhs._value); } template UHD_INLINE @@ -85,15 +103,10 @@ namespace uhd { namespace math { namespace fp_compare { return !(lhs < rhs); } - template UHD_INLINE - bool operator==(fp_compare_epsilon 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 + UHD_INLINE bool operator==(fp_compare_epsilon lhs, double rhs) + { + return lhs == fp_compare_epsilon(static_cast(rhs)); } template UHD_INLINE @@ -104,7 +117,7 @@ namespace uhd { namespace math { namespace fp_compare { template UHD_INLINE bool operator<(fp_compare_epsilon lhs, double rhs) { - return (lhs._value + lhs._epsilon) < (rhs - DOUBLE_PRECISION_EPSILON); + return (lhs != rhs) && (lhs._value < rhs); } template UHD_INLINE @@ -115,7 +128,7 @@ namespace uhd { namespace math { namespace fp_compare { template UHD_INLINE bool operator>(fp_compare_epsilon lhs, double rhs) { - return (lhs._value - lhs._epsilon) > (rhs + DOUBLE_PRECISION_EPSILON); + return (lhs != rhs) && (lhs._value > rhs); } template UHD_INLINE @@ -123,15 +136,10 @@ namespace uhd { namespace math { namespace fp_compare { return !(lhs < rhs); } - template UHD_INLINE - bool operator==(double lhs, fp_compare_epsilon 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 + UHD_INLINE bool operator==(double lhs, fp_compare_epsilon rhs) + { + return fp_compare_epsilon(static_cast(lhs)) == rhs; } template UHD_INLINE @@ -142,7 +150,7 @@ namespace uhd { namespace math { namespace fp_compare { template UHD_INLINE bool operator<(double lhs, fp_compare_epsilon rhs) { - return (lhs + DOUBLE_PRECISION_EPSILON) < (rhs._value - rhs._epsilon); + return (lhs != rhs) && (lhs < rhs._value); } template UHD_INLINE @@ -153,7 +161,7 @@ namespace uhd { namespace math { namespace fp_compare { template UHD_INLINE bool operator>(double lhs, fp_compare_epsilon rhs) { - return (lhs - DOUBLE_PRECISION_EPSILON) > (rhs._value + rhs._epsilon); + return (lhs != rhs) && (lhs > rhs._value); } template 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 -- cgit v1.2.3