| 394 | |
| 395 | template <typename FloatType, typename IntType> |
| 396 | static inline bool GDALFloatAlmostEquals(FloatType A, FloatType B, |
| 397 | unsigned maxUlps) |
| 398 | { |
| 399 | static_assert(sizeof(FloatType) == sizeof(IntType)); |
| 400 | // This function will allow maxUlps-1 floats between A and B. |
| 401 | |
| 402 | // Make sure maxUlps is non-negative and small enough that the default NAN |
| 403 | // won't compare as equal to anything. |
| 404 | if (maxUlps >= 4 * 1024 * 1024) |
| 405 | { |
| 406 | CPLError(CE_Failure, CPLE_IllegalArg, "Invalid maxUlps"); |
| 407 | return false; |
| 408 | } |
| 409 | |
| 410 | const auto MapToInteger = [](FloatType x) |
| 411 | { |
| 412 | IntType i = 0; |
| 413 | memcpy(&i, &x, sizeof(i)); |
| 414 | |
| 415 | constexpr int NBITS = 8 * static_cast<int>(sizeof(i)) - 1; |
| 416 | constexpr IntType SHIFT = (static_cast<IntType>(1) << NBITS); |
| 417 | |
| 418 | // Make i lexicographically ordered with negative values |
| 419 | // remapped to the [0, SHIFT[ range and positive values in the |
| 420 | // [SHIFT, UINT_MAX[ range |
| 421 | if ((i >> NBITS) != 0) |
| 422 | i = SHIFT - (i & ~SHIFT); |
| 423 | else |
| 424 | i += SHIFT; |
| 425 | return i; |
| 426 | }; |
| 427 | |
| 428 | const auto aInt = MapToInteger(A); |
| 429 | const auto bInt = MapToInteger(B); |
| 430 | |
| 431 | return ((aInt > bInt) ? aInt - bInt : bInt - aInt) <= maxUlps; |
| 432 | } |
| 433 | |
| 434 | bool GDALFloatAlmostEquals(float A, float B, unsigned maxUlps) |
| 435 | { |