| 1 | /* |
|---|
| 2 | * Copyright 2023 Google LLC |
|---|
| 3 | * |
|---|
| 4 | * Use of this source code is governed by a BSD-style license that can be |
|---|
| 5 | * found in the LICENSE file. |
|---|
| 6 | */ |
|---|
| 7 | |
|---|
| 8 | #include "include/private/base/SkFloatingPoint.h" |
|---|
| 9 | |
|---|
| 10 | #include <algorithm> |
|---|
| 11 | #include <cmath> |
|---|
| 12 | #include <limits> |
|---|
| 13 | |
|---|
| 14 | // Return the positive magnitude of a double. |
|---|
| 15 | // * normalized - given 1.bbb...bbb x 2^e return 2^e. |
|---|
| 16 | // * subnormal - return 0. |
|---|
| 17 | // * nan & infinity - return infinity |
|---|
| 18 | static double magnitude(double a) { |
|---|
| 19 | static constexpr int64_t = |
|---|
| 20 | 0b0'11111111111'0000000000000000000000000000000000000000000000000000; |
|---|
| 21 | int64_t bits; |
|---|
| 22 | memcpy(dest: &bits, src: &a, n: sizeof(bits)); |
|---|
| 23 | bits &= extractMagnitude; |
|---|
| 24 | double out; |
|---|
| 25 | memcpy(dest: &out, src: &bits, n: sizeof(out)); |
|---|
| 26 | return out; |
|---|
| 27 | } |
|---|
| 28 | |
|---|
| 29 | bool sk_doubles_nearly_equal_ulps(double a, double b, uint8_t maxUlpsDiff) { |
|---|
| 30 | |
|---|
| 31 | // The maximum magnitude to construct the ulp tolerance. The proper magnitude for |
|---|
| 32 | // subnormal numbers is minMagnitude, which is 2^-1021, so if a and b are subnormal (having a |
|---|
| 33 | // magnitude of 0) use minMagnitude. If a or b are infinity or nan, then maxMagnitude will be |
|---|
| 34 | // +infinity. This means the tolerance will also be infinity, but the expression b - a below |
|---|
| 35 | // will either be NaN or infinity, so a tolerance of infinity doesn't matter. |
|---|
| 36 | static constexpr double minMagnitude = std::numeric_limits<double>::min(); |
|---|
| 37 | const double maxMagnitude = std::max(a: std::max(a: magnitude(a), b: minMagnitude), b: magnitude(a: b)); |
|---|
| 38 | |
|---|
| 39 | // Given a magnitude, this is the factor that generates the ulp for that magnitude. |
|---|
| 40 | // In numbers, 2 ^ (-precision + 1) = 2 ^ -52. |
|---|
| 41 | static constexpr double ulpFactor = std::numeric_limits<double>::epsilon(); |
|---|
| 42 | |
|---|
| 43 | // The tolerance in ULPs given the maxMagnitude. Because the return statement must use < |
|---|
| 44 | // for comparison instead of <= to correctly handle infinities, bump maxUlpsDiff up to get |
|---|
| 45 | // the full maxUlpsDiff range. |
|---|
| 46 | const double tolerance = maxMagnitude * (ulpFactor * (maxUlpsDiff + 1)); |
|---|
| 47 | |
|---|
| 48 | // The expression a == b is mainly for handling infinities, but it also catches the exact |
|---|
| 49 | // equals. |
|---|
| 50 | return a == b || std::abs(lcpp_x: b - a) < tolerance; |
|---|
| 51 | } |
|---|
| 52 | |
|---|
| 53 | bool sk_double_nearly_zero(double a) { |
|---|
| 54 | return a == 0 || fabs(x: a) < std::numeric_limits<float>::epsilon(); |
|---|
| 55 | } |
|---|
| 56 | |
|---|
