7

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?

Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
kyle
  • 197
  • 1
  • 2
  • 11
  • 1
    MATLAB's `eps` always gives positive results. The above code would give negative results if `x` is greater than `epsf`. [Here](http://coliru.stacked-crooked.com/a/68546c8c401c0610)'s the fixed code: `double eps(float x) { float xp = std::abs(x); double x1 = std::nextafter(xp, xp + 1.0f); return x1 - xp; }` – legends2k Dec 28 '15 at 12:02

2 Answers2

8

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.

user4815162342
  • 141,790
  • 18
  • 296
  • 355
3

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.

horchler
  • 18,384
  • 4
  • 37
  • 73
  • yes, thank you I was aware of the boost implementation, however for our class we are not able to use the boost libs. – kyle Nov 23 '13 at 23:20
  • I really wonder how in the world this answer rates a down vote. Not only did I explain the concept of ULP, but I suggested a robust alternative solution that works for both single and double precision. The OP had not indicated the his wish not to use boost. – horchler Nov 24 '13 at 16:36
  • 1
    There's a +1 for you to balance out the -1. Thanks for the background info on ULP. Helpful to understand as I try to implement this in Go. – Nathaniel Jones Nov 17 '21 at 19:40