Floating point approximately equal

Question

I implemented the following code for determining if two floating point numbers are equal is_float_equal(...). It handles the tests for my use cases well, but I was wondering if the community could provide feedback on the implementation and its correctness?

My objectives are to consider two floating pointer numbers (represented by value_t) equal if their difference is subnormal or their difference is less than machine epsilon scaled by the largest number.

// feb 16 417 PM
#include <iostream>
#include <iomanip>
#include <cassert>
#include <limits>
#include <cmath>
#include <fstream>

std::fstream llog("out.txt", std::ios::out);

// https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon
// https://stackoverflow.com/questions/17333/what-is-the-most-effective-way-for-float-and-double-comparison
// https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/

/**
 * @brief compares l and r floating point values for equality
 * @param l the left value
 * @param r the right value
 * @returns true if l is equal to r
 */
template <typename value_t>
bool is_float_equal(value_t l, value_t r)
{
    // https://en.cppreference.com/w/cpp/numeric/math/isnan
    if (std::isnan(l) || std::isnan(r))
    {
        return false;
    }
    // https://en.cppreference.com/w/cpp/numeric/math/isinf
    if (std::isinf(l) && std::isinf(r))
    {
        return l == r; // NOLINT
    }
    if (l < r)
    {
        value_t const t{l};
        l = r;
        r = t;
    }
    value_t const d{l - r};
    value_t constexpr e{std::numeric_limits<value_t>::epsilon()};
    value_t constexpr m{std::numeric_limits<value_t>::min()};
    return d <= m || d <= l * e;
}

void tests()
{
    float const inf{std::numeric_limits<float>::infinity()};
    assert(-inf == -inf && is_float_equal(-inf, -inf));
    assert(-inf != inf && !is_float_equal(-inf, inf));
    assert(inf != -inf && !is_float_equal(inf, -inf));
    assert(inf == inf && is_float_equal(inf, inf));

    float const zn{-0.0f};
    float const zp{0.0f};
    assert(zn == zn && is_float_equal(zn, zn));
    assert(zp == zn && is_float_equal(zp, zn));
    assert(zp == zp && is_float_equal(zp, zp));

    assert(std::nextafter(zn, zp) == zp);
}

int main()
{
    tests();
    using value_t = float;
    value_t constexpr e{std::numeric_limits<value_t>::epsilon()};
    value_t constexpr max{std::numeric_limits<value_t>::max()};
    value_t constexpr m{std::numeric_limits<value_t>::min()};
    value_t l{std::numeric_limits<value_t>::lowest()};
    value_t r{l};
    value_t d{0};
    unsigned long long ec{0};
    unsigned long long tc{0};
    std::cout << std::fixed << std::setprecision(50) << l << ' ' << r << ' ' << e << '\n';
    while (r < max)
    {
        // https://en.wikipedia.org/wiki/Floating-point_arithmetic
        // https://peps.python.org/pep-0485/
        float d{r - l};
        while (d <= m || std::abs((r - l) / r) <= e)
        {
            r = std::nextafter(r, max);
            d = r - l;
        }

        ++tc;
        double const p{static_cast<double>(ec) / static_cast<double>(tc)};
        std::cout << std::fixed << std::setprecision(50) << l << ' ' << r << ' ' << d << ' ' << p << '\n';
        if (!is_float_equal(l, l))
        {
            llog << std::endl
                 << "l not equal to l\n";
            is_float_equal(l, l);
            return EXIT_FAILURE;
        }
        if (is_float_equal(l, r))
        {
            ++ec;
            llog << std::endl
                 << ec << '\n'
                 << "l equal to r\n";
            is_float_equal(l, r);
            if (ec > 1000)
            {
                return EXIT_FAILURE;
            }
        }
        l = r;
    }
    return EXIT_SUCCESS;
}

Answer 1

When we multiply epsilon by the largest number (absolute value, as noted in chux's answer), we end up with a value that's about std::nextafter(l) - l. So we can replace all that work by simply testing whether the numbers are adjacent according to std::nextafter:

#include <cmath>
#include <limits>
#include <concepts>

template<std::floating_point value_t> [[nodiscard]]
constexpr bool is_float_equal(value_t l, value_t r)
{
    constexpr auto infinity = std::numeric_limits<value_t>::infinity();

    auto const min = std::nextafter(l, -infinity);
    auto const max = std::nextafter(l, infinity);
    return (min <= r && r <= max);
}

Simpler still (thanks to Ilmari Karonen, we can just take advantage of std::nextafter() doing what we need in terms of the "direction" argument and just write

template<std::floating_point value_t> [[nodiscard]]
constexpr bool is_float_equal(value_t l, value_t r)
{
    return r == std::nextafter(l, r);
}

• No need to std::swap() the arguments or to special-case NaNs and infinities.
• Works well close to ±0.
• No multiplications.
• Constexpr - which it should already have been.
• Full compliance with the test suite (until my timeout of 60 seconds expired, anyway).

---

Style comments:

There's a huge amount of output, which I had to disable to get a reasonable amount of the test suite to run. A good test should be silent when everything is correct, and only print the failures.

Failure messages should go to std::cerr, not std::cout.

I dislike the test program writing files on my system (doesn't matter that the file is empty - the fact it overwrites any existing file is bad enough).

Unused variable: value_t d{0};. Worse, it's shadowed by another d of the same type in smaller scope.

There's a couple of places we call is_float_equal() and don't use the result. Adding [[nodiscard]] identifies these.

Answer 2

Right goal?

"difference is less than machine epsilon scaled by the largest number."

l * e scales epsilon by the greatest number. (most positive).

I`d expect scaling by the largest value in magnitude (e.g. absolute value)

Consider two negative numbers like -1.0 and -1.000...001 that differ only in the last bit. d <= m || d <= l * e is never true as it is like 1e-17 < 1e-300 || 1e-17 <= some_negative.

Code needs ... || d <= max(abs(l), abs(r)) * e or the like.

Minor: Epsilon first

Consider reordering as d <= m is rarely true.

// return d <= m || d <= l * e;
return d <= l * e || d <= m;

Answer 3

The name of the function is wrong. The == operator returns whether two floating point numbers are equal. You return whether they are equal or very close together. Call your function "is_almost_equal" or "maybe_equal". Can be accompanied by "definitely_greater" (greater and not maybe_equal), "maybe_greater_equal" (greater or maybe_equal) etc.

I think you are much to strict. What you want to catch is numbers that should have been equal but are not due to rounding errors, and avoid false positives (independent values that are close together by coincidence). Rounding errors can be a lot larger. But being within 10^-15 by coincidence won't happen.

Swift has a function like this in its standard library, and it checks whether about half the bits match.