Is there a standard function for computing units in the last place in C++?

Question

For a number x in the binary floating point number system, units in the last place is defined as ulp(x) = 2^{E(x)} * 2^{1 - p} = 2^{E(x)}*epsilon, where the exponent is the exponent of the base 2 used to represent x, p is the precision of the number format (53 for the double precision format), and epsilon = 2^{1-p} is the machine tolerance.

Looking at std::numeric_limits I expected to find something like std::numeric_limits<double>::ulp(x), but this seems to be unavailable. Or something like std::numeric_limits<double>::exponent(x) that I can multiply with std::numeric_limits<double>::epsilon() but it is not there. The ulp function seems to be available in boost.

I can code something myself, but I would much rather use standardized code, if there is one. Did I miss something in the documentation?

Edit: to make it more clear, I am talking about the calculation of the exponent of a floating point number x, not getting min_ or max_ exponent range for a floating point number format (type).

Answer 1

The standard library indeed does not have an ulp function.

2^{E(x)}*epsilon

Can be calculated with:

int exp;
std::frexp(std::fabs(x), &exp);
double result = std::ldexp(std::numeric_limits<double>::epsilon(), exp - 1);

Caveat: This is only correct for normal floating point values. For denormal values there might not even exist a floating point value that can represent ULP.

---

I expected to find something like std::numeric_limits<double>::ulp(x), but this seems to be unavailable. Or something like std::numeric_limits<double>::exponent(x)

numeric_limits contains only constants (or functions which return constants). It doesn't have any functions that perform calculations and thus no functions which take input values. <cmath> header has functions for calculations.

---

As a sidenote to readers in general: If your intention is to find the next/previous representable value, use nextafter family of functions instead of adding/subtracting ULP.

Answer 2

I found an alternative, it doesn't fit into comments, so I'm posting an answer.

I think it can also be done by using std::nextafter and std::nexttoward:

template<typename FT> 
std::enable_if_t<std::is_floating_point_v<FT>, FT> 
ulp(FT x)
{
    if (x > 0)
        return std::nexttoward(x, std::numeric_limits<FT>::infinity()) - x;
    else 
        return x - std::nexttoward(x, -std::numeric_limits<FT>::infinity());
}

This is valid if nexttoward results in a floating point number that is next toward INF if the number is positive, and toward -INF if the number is negative. Because this number is a floating point number, the distance from x is going to be ulp(x) and it doesn't have a rounding error.

Here is a test program for this idea, it seems to work without using power functions and computing exponents manually (it's probably doing something like that under the hood). I tested it by checking:

// Check ulp by adding 0.5ulp to x and let the rounding to nearest round to x. 
assert((x + (0.5 *ulpx) == x));

Here is the MWE for it:

#include <cmath>
#include <cassert>
#include <iostream>
#include <iomanip>

using namespace std;

template<typename FT> 
std::enable_if_t<std::is_floating_point_v<FT>, FT> 
ulp(FT x)
{
    if (x > 0)
        return std::nexttoward(x, numeric_limits<FT>::infinity()) - x;
    else 
        return x - std::nexttoward(x, -numeric_limits<FT>::infinity());
}

void test_ulp(double x)
{
    double ulpx = ulp(x);

    std::cout << "ulpx = " << std::setprecision(25) << ulpx << std::endl;

    double xNext = x + ulpx; 

    // Check ulp by adding 0.5ulp to x and let the rounding to nearest round to x. 
    if (x > 0)
        assert((x + (0.5 *ulpx) == x));
    else 
        assert((x - (0.5 *ulpx) == x));

    std::cout << std::setprecision(25) <<  x << " + " << ulpx 
        << " == " << xNext << std::endl; 
}

int main()
{
    // Denormalized test. 
    test_ulp(0.); 

    // Epsilon test. 
    test_ulp(1.); 

    // Random number test. :) 
    test_ulp(42.); 

    // Epsilon test. 
    test_ulp(-1.); 
}

Edit: tried to find what nexttoward does, looked into gcc STL source code and only found this:

#ifndef __CORRECT_ISO_CPP11_MATH_H_PROTO_FP
   constexpr float
   nexttoward(float __x, long double __y)
   { return __builtin_nexttowardf(__x, __y); }

   constexpr long double
   nexttoward(long double __x, long double __y)
   { return __builtin_nexttowardl(__x, __y); }
 #endif

 #ifndef __CORRECT_ISO_CPP11_MATH_H_PROTO_INT
   template<typename _Tp>
     constexpr typename __gnu_cxx::__enable_if<__is_integer<_Tp>::__value, 
                                               double>::__type
     nexttoward(_Tp __x, long double __y)
     { return __builtin_nexttoward(__x, __y); }
 #endif

Answer 3

Technically, solution could be to extract the exponent bits, and subtract 1074 from it. This gives the exponent of ulp. Then scale the sign by this exponent if you want a signed ulp. Something like

int exp = extract_exponent_bits(x):
Int sign = 1 - 2* extract_sign_bit(x);
return ldexp((double)sign,exp-1074);

Bit extraction left as an exercize, pointer aliasing like in sun fdlibm is not an option nowadays, so just memcpy to an uint64, shift >> 52, etc...

Generalization to non ieee might be more involved, but are you going to need it?