src/base/SkFloatingPoint.cpp

/*
 * Copyright 2023 Google LLC
 *
 * Use of this source code is governed by a BSD-style license that can be
 * found in the LICENSE file.
 */

#include "include/private/base/SkFloatingPoint.h"

#include <algorithm>
#include <cmath>
#include <limits>

// Return the positive magnitude of a double.
// * normalized - given 1.bbb...bbb x 2^e return 2^e.
// * subnormal - return 0.
// * nan & infinity - return infinity
static double magnitude(double a) {
    static constexpr int64_t extractMagnitude =
            0b0'11111111111'0000000000000000000000000000000000000000000000000000;
    int64_t bits;
    memcpy(&bits, &a, sizeof(bits));
    bits &= extractMagnitude;
    double out;
    memcpy(&out, &bits, sizeof(out));
    return out;
}

bool sk_doubles_nearly_equal_ulps(double a, double b, uint8_t maxUlpsDiff) {

    // The maximum magnitude to construct the ulp tolerance. The proper magnitude for
    // subnormal numbers is minMagnitude, which is 2^-1021, so if a and b are subnormal (having a
    // magnitude of 0) use minMagnitude. If a or b are infinity or nan, then maxMagnitude will be
    // +infinity. This means the tolerance will also be infinity, but the expression b - a below
    // will either be NaN or infinity, so a tolerance of infinity doesn't matter.
    static constexpr double minMagnitude = std::numeric_limits<double>::min();
    const double maxMagnitude = std::max(std::max(magnitude(a), minMagnitude), magnitude(b));

    // Given a magnitude, this is the factor that generates the ulp for that magnitude.
    // In numbers, 2 ^ (-precision + 1) = 2 ^ -52.
    static constexpr double ulpFactor = std::numeric_limits<double>::epsilon();

    // The tolerance in ULPs given the maxMagnitude. Because the return statement must use <
    // for comparison instead of <= to correctly handle infinities, bump maxUlpsDiff up to get
    // the full maxUlpsDiff range.
    const double tolerance = maxMagnitude * (ulpFactor * (maxUlpsDiff + 1));

    // The expression a == b is mainly for handling infinities, but it also catches the exact
    // equals.
    return a == b || std::abs(b - a) < tolerance;
}

bool sk_double_nearly_zero(double a) {
    return a == 0 || fabs(a) < std::numeric_limits<float>::epsilon();
}