Implement Matlab's eps(x) function in C++

Question

I'm trying to implement Matlab's eps(x) function in C++

For example, in Matlab:

>> eps(587.3888)
ans = 1.1369e-13
>> eps(single(587.3888))
ans = 6.1035e-05

However, when I try to do this in C++ I am not able to get the correct single precision answer.

#include <limits>
#include <iostream>
#include <math.h>

#define DEBUG(x) do { std::cerr << x << std::endl; } while (0)
#define DEBUG2(x) do { std::cerr << #x << ": " << x << std::endl; } while (0)

int main() {

    float epsf = std::numeric_limits<float>::epsilon();
    DEBUG2(epsf);
    double epsd = std::numeric_limits<double>::epsilon();
    DEBUG2(epsd);

    float espxf = nextafter(float(587.3888), epsf) - float(587.3888);
    double espxd = nextafter(double(587.3888), epsd) - double(587.3888);
    DEBUG2(espxf);
    DEBUG2(espxd);

}

Running the program I get the following output:

$ ./a.out 
epsf: 1.19209e-07
epsd: 2.22045e-16
espxf: -1.13687e-13
espxd: -1.13687e-13

It seems as though for some reason even though the eps values for single and double precision are correct, the output using nextafter function only outputs the double precision value. My value for epsxf should be 6.1035e-05 as it is in Matlab.

Any thoughts?

Answer 1

Include <cmath> and call std::nextafter, and your code will work, provided you have a C++11 compiler.

Including <math.h> and calling ::nextafter invokes the C version of the function. The C implementation of nextafter obviously supports no overloads, so C provides a nextafterf for single-precision result, as well as nextafterl for quad-precision. (Simply calling double-precision nextafter with float fails because the argument gets converted to double.) If you don't have a C++11 compiler, you can fix your code by invoking ::nextafterf.

Answer 2

Use libraries. Matlab's eps function in other languages is called ULP, for unit in last place. According to Wikipedia article on ULP, the following function from the boost C++ library can be used to calculate the floating point distance between two doubles a and b:

boost::math::float_distance(a, b)

The documentation for float_distance is here.