GMT_equal_double and GMT_equal_float fail for values close to zero
Consider this pseudocode example:
float someFloat = 67329.2348f; // arbitrarily chosen
float nextFloat = NearbyFloat(someFloat, 1); // exactly one ULP away from 'someFloat'
float diff = nextFloat - someFloat; // = .0078125000; diff is 1,006,632,960 ULPs away from zero!
Btw, FLT_EPSILON is 872,415,232 ULPs away from zero.
#1 Updated by Florian about 7 years ago
- Status changed from New to Resolved
- % Done changed from 0 to 100
- Affected version changed from all to 5.x-dev
Here is another typical catastrophic cancellation:
float degree = 360.0; float lessDegree = nearbyFloat(degree, -1); // = 359.9999695 float diff = degree - lessDegree; // = 0.0000305Now here comes the example why this is dangerous.
#define GMT_CONV_LIMIT 1.0e-8 #define GMT_IS_ZERO(x) (fabs (x) < GMT_CONV_LIMIT) // test if lessDegree ~= degree (very common gmt code): GMT_IS_ZERO (degree - lessDegree) // returns false even though lessDegree is only 1 ULP away from degree!This is probably not what we want. Though, I cheated in this example because I used float instead of double, this example shows that such code can lead to unexpected results. In gmt it worked so far because GMT_CONV_LIMIT is very very large in comparison to DBL_EPSILON. GMT_IS_ZERO---if used in such a way---would also fail for large doubles (10 million or so). Here are some new functions that could be used instead:
/* common_math.c */ * int32_abs Like abs but for int32_t * int64_abs Like abs but for int64_t * floatAlmostEqualUlps Compare 2 floats with ULPs (not safe for values near zero) * doubleAlmostEqualUlps Same as floatAlmostEqualUlps for double * floatAlmostEqualUlpsAndAbs Compare 2 floats with an absolute epsilon check (values near zero), then based on ULPs * doubleAlmostEqualUlpsAndAbs Same as floatAlmostEqualUlpsAndAbs for double * floatAlmostEqual Convenience macro for floatAlmostEqualUlps * doubleAlmostEqual Convenience macro for doubleAlmostEqualUlps * floatAlmostEqualZero Convenience macro for floatAlmostEqualUlpsAndAbs * doubleAlmostEqualZero Convenience macro for doubleAlmostEqualUlpsAndAbsThey are mostly inspired by these C++ functions and don't have the overflow bug #2. *AlmostEqualUlpsAndAbs and *AlmostEqualZero also can handle values close to zero, which GMT_equal_double could not. So instead of
it is safer to use
doubleAlmostEqual (lon, 360)
or if both parameters might be zero
doubleAlmostEqualZero (lon1, lon2)
This has been fixed in r9839.