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