Feature suggestion: Add support for testing equality of floats with precision specified in ulps #864
Description
As discussed in #863 , there are cases where one want to check approximate equality of floats in ulps instead of relative error. This is useful to check that results are rounded according to expectations.
I will use floats to denote all floating point types.
I have implemented a draft and tested it with doctest.
The draft defines two functions ulpeq to check equality between floats within specified ulps, and ulpadd which can add a positive number of ulps to a float. The implementation is tested both with random numbers and edge cases.
The random test case have two configurable parameters: maximum number of ulps MAX_ULPS to test the implementation for, and the number of ITERATIONS .
#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
#include <doctest/doctest.h>
#include <cassert>
#include <limits>
#include <random>
#define MAX_ULPS 100
#define ITERATIONS 1'000'000
template <class T>
std::enable_if_t<not std::numeric_limits<T>::is_integer, bool>
ulpeq(T x, T y, int n)
{
assert(n >= 0);
if (x == y) {
return true;
}
if (x > y) {
return ulpeq(y, x, n);
}
if (n == 0) {
return x == y;
}
return ulpeq(std::nextafter(x, y), y, n - 1);
}
template <typename T>
std::enable_if_t<not std::numeric_limits<T>::is_integer, T>
ulpadd(T x, int n)
{
assert(n >= 0);
if (n <= 0) {
return x;
}
return ulpadd(std::nextafter(x, std::numeric_limits<T>::infinity()), n - 1);
}
TEST_CASE_TEMPLATE("ulpeq-random", T, float, double, long double)
{
std::mt19937 g(0);
std::uniform_int_distribution<int> ulp_distrib(0, MAX_ULPS);
std::uniform_real_distribution<T> value_distrib(std::numeric_limits<T>::min(), std::numeric_limits<T>::max());
for (int i = 0; i < ITERATIONS; ++i)
{
const int n = ulp_distrib(g);
const T x = value_distrib(g);
const T y = ulpadd(x, n);
REQUIRE(ulpeq(x, y, n));
REQUIRE(ulpeq(y, x, n));
}
}
TEST_CASE_TEMPLATE("ulpeq-edge-cases", T, float, double, long double)
{
// Extreme values values
REQUIRE(ulpeq(std::numeric_limits<T>::max(), std::numeric_limits<T>::max(), 0));
REQUIRE(ulpeq(+std::numeric_limits<T>::min(), +std::numeric_limits<T>::min(), 0));
REQUIRE(ulpeq(T{0}, T{0}, 0));
REQUIRE(ulpeq(-std::numeric_limits<T>::min(), -std::numeric_limits<T>::min(), 0));
REQUIRE(ulpeq(std::numeric_limits<T>::lowest(), std::numeric_limits<T>::lowest(), 0));
// Infinities
REQUIRE(ulpeq(+std::numeric_limits<T>::infinity(), +std::numeric_limits<T>::infinity(), 0));
REQUIRE(ulpeq(-std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity(), 0));
REQUIRE(ulpeq(-std::numeric_limits<T>::infinity(), std::numeric_limits<T>::lowest(), 1));
REQUIRE(ulpeq(std::numeric_limits<T>::lowest(), -std::numeric_limits<T>::infinity(), 1));
REQUIRE(ulpeq(std::numeric_limits<T>::infinity(), std::numeric_limits<T>::max(), 1));
REQUIRE(ulpeq(std::numeric_limits<T>::max(), std::numeric_limits<T>::infinity(), 1));
// Subnormals
for (int i = 1; i <= MAX_ULPS; ++i) {
REQUIRE(ulpeq(T{0}, +ulpadd(T{0}, i), i));
REQUIRE(ulpeq(T{0}, -ulpadd(T{0}, i), i));
}
}
Interestingly, in the docs for std::numeric_limits<T>::epsilon() here https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon there is even an example for checking equality within ulps. That example implementation is problematic since it performs floating point computations to check the floating points... I am pretty sure that implementation is either doomed to fail for certain floating point values, both in terms of false positives and false negatives.
Given the ulpeq function above another class like doctest::Approx could be defined for checking approximate equality within ulps instead of relative error. I don't think adding a method .ulp(int n) as suggested in the referenced thread is a good idea, better add a completely new class which has it owns semantics and guarantees?
I intend the draft here to be 100% correct, would be great if this can be examined by several people to find more edge cases.Happy to discuss this and receive input from as many as possible before proceeding to final implementation (if the community agree that this is a good idea to implement).