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;
}