C++ Floating-Point Number Comparison

Introduction

Floating-point number comparison has always been a problem in computer applications. Because a real number cannot often be precisely represented by floating-point numbers and floating-pointer computations have rounding errors, the order of floating point calculation matters for bitwise equality. For example, a + b + c might not be bitwise equal to c + b + a. Because the floating-point errors could accumulate and propagate, solving a problem using different but correct mathematical recipes might result in different outcomes. Therefore, we would often need to compare floating-point numbers to determine if our certain algorithm implementation produces correct results.

In this blog post, I would like to discuss how to compare floating-point numbers in C++ empirically.

C++ Floating-Point Number Comparison

It is common to see that a absolute tolerance value and a relative tolerance value were used to compare floating-point numbers, such as the numpy.allclose function. However, the problem is that the relative error between two small close floating-point numbers might be larger than the relative error between two large close floating-point numbers. Similarly, the absolute error between two small close floating-point numbers might be smaller than the absolute error between two large close floating-point numbers.

Therefore, we could use different absolute tolerance value and relative tolerance value for floating numbers in different ranges.

test_float_diff.cpp

#include <cassert>
#include <cmath>
#include <functional>
#include <iostream>
#include <limits>
#include <type_traits>

// Commutative (symmetric)
template <typename T,
          std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
bool float_close(T a, T b, std::function<T(T, T)> const& rtol_func,
                 std::function<T(T, T)> const& atol_func)
{
    if (std::isinf(a) || std::isinf(b))
    {
        return false;
    }
    if (std::isnan(a) || std::isnan(b))
    {
        return false;
    }

    T const abs_a{std::abs(a)};
    T const abs_b{std::abs(b)};
    T const abs_diff{std::abs(a - b)};
    T const abs_sum{abs_a + abs_b};
    T const rtol{rtol_func(a, b)};
    T const atol{atol_func(a, b)};

    if (a == b)
    {
        return true;
    }
    else
    {
        return (abs_diff / abs_sum < rtol && abs_diff < atol);
    }
}

// This function can be customized to the user's needs.
template <typename T,
          std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
T dynamic_rtol_func(T a, T b)
{
    T const abs_a{std::abs(a)};
    T const abs_b{std::abs(b)};
    T const abs_sum{abs_a + abs_b};

    T rtol{0};

    if (a == 0 && b == 0)
    {
        rtol = std::numeric_limits<T>::infinity();
    }
    else if (abs_sum < 1e-8)
    {
        rtol = 5;
    }
    else if (abs_sum < 1e-5)
    {
        rtol = 1;
    }
    else if (abs_sum < 1e-2)
    {
        rtol = 1e-2;
    }
    else
    {
        rtol = 1e-4;
    }

    return rtol;
}

// This function can be customized to the user's needs.
template <typename T,
          std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
T dynamic_atol_func(T a, T b)
{
    T const abs_a{std::abs(a)};
    T const abs_b{std::abs(b)};
    T const abs_sum{abs_a + abs_b};

    T atol{0};

    if (a == 0 && b == 0)
    {
        atol = std::numeric_limits<T>::min();
    }
    else if (abs_sum < 1e-8)
    {
        atol = 1e-9;
    }
    else if (abs_sum < 1e-5)
    {
        atol = 1e-7;
    }
    else if (abs_sum < 1e-2)
    {
        atol = 1e-5;
    }
    else
    {
        atol = 1e-4;
    }

    return atol;
}

int main()
{
    float const a{std::numeric_limits<float>::min()};
    float const b{std::numeric_limits<float>::min() * 100};
    std::cout << "a: " << a << std::endl;
    std::cout << "b: " << b << std::endl;
    bool const is_close{float_close<float>(a, b, dynamic_rtol_func<float>,
                                           dynamic_atol_func<float>)};
    std::cout << "is_close: " << is_close << std::endl;
    assert(is_close);
}

$ g++ test_float_diff.cpp -o test_float_diff -std=c++14
$ ./test_float_diff
a: 1.17549e-38
b: 1.17549e-36
is_close: 1

References

• The Floating-Point Guide - Comparison