5

Consider a "normal" real number TREAL x in C++ (not subnormal and not NaN/Infinite) (TREAL = float, double, long double)
Is the following the good solution to find the previous and next x from a floating-point point of view ?

TREAL xprev = (((TREAL)(1.)) - std::numeric_limits<TREAL>::epsilon()) * x;
TREAL xnext = (((TREAL)(1.)) + std::numeric_limits<TREAL>::epsilon()) * x;

Thank you very much.

Eitan T
  • 32,660
  • 14
  • 72
  • 109
Vincent
  • 57,703
  • 61
  • 205
  • 388
  • 2
    You notice that x prev next != x? – Daniel Jun 22 '12 at 15:56
  • Are you saying that you don't want (x+1) but rather the value if you increment the mantissa? – Skizz Jun 22 '12 at 15:59
  • Yes, I dont want x+1 but x +/- epsilon (it is to check boundaries taking in account possible precision problems) – Vincent Jun 22 '12 at 16:15
  • Same question there http://stackoverflow.com/questions/1336767/what-is-the-next-normalised-floating-point-number-afterbefore-a-normalised-flo/11063887#11063887 – aka.nice Jun 22 '12 at 17:04

4 Answers4

13

C99 and C++11 have nextafter, nextafterl and nextafterf functions in <math.h> and <cmath>. Implementing them with basic arithmetic and epsilon would be tedious as you'd need to take rounding into account. Working on the binary representation is probably easier, but I wonder about the effect of the sign and magnitude representation and the existence of -0.0 (see Fred's answer for what is needed).

Community
  • 1
  • 1
AProgrammer
  • 51,233
  • 8
  • 91
  • 143
5

Getting the next floating point number is a lot easier on the binary level:

float next(float f)
{
    unsigned x;
    memcpy(&x, &f, 4);
    ++x;
    memcpy(&f, &x, 4);
    return f;
}

Of course this will only work for systems where floating point numbers are stored "in ascending order", which happens to be the case for IEEE754.

Negative numbers will go towards negative infinity. Want them to go to zero instead? Use this:

float next(float f)
{
    int x;
    memcpy(&x, &f, 4);
    x += x >> 31 | 1;   // this will add 1 for positive f and -1 for negative f
    memcpy(&f, &x, 4);
    return f;
}
fredoverflow
  • 256,549
  • 94
  • 388
  • 662
  • 1
    doesn't work for -0 (it gives a NaN, nextafterf gives 1.4012984643e-45). – AProgrammer Jun 22 '12 at 18:42
  • @AProgrammer I had a small bug in there, shifting by 30 instead of 31. Does it work now? – fredoverflow Jun 22 '12 at 19:07
  • same problem. -0 is represented as 0x800000000 and next() should returns 0x00000001, I don't see a way to avoid a test for it. (then we will look at the qNaN and sNaN ;)) – AProgrammer Jun 23 '12 at 10:37
  • Brilliant solution! I translated it into C#. Thanks to a hack, it's not only safe code, but it's actually much more efficient than the C code, believe it or not. – IamIC Jun 03 '13 at 11:20
2

No, the ratio between "consecutive" floating point values is not uniform; this approach may miss some out, or leave you stuck at a point where xnext == x.

To move from one value to the next-largest value, you'd have to:

  • extract the mantissa and exponent;
  • increment the mantissa;
  • if that overflows, reset it and increment the exponent;
  • reconstruct the value from the exponent and mantissa.

The details are quite fiddly, and will probably require some knowledge of the floating point representation.

However, assuming a representation similar to IEEE, you could achieve this by reinterpreting the bit pattern as a large enough integer, and incrementing that integer. That will increment the mantissa, with any overflow going into the exponent, just as we want.

Mike Seymour
  • 249,747
  • 28
  • 448
  • 644
0

The following nextFloat function gives the correct result for all floats of at least minV (that is, floats for which the intermediate values in nextFloat stay out of the denormal range). I tested it for all floats starting at minV, up to FLT_MAX and the result was always equal to the result of nextFloatRef.

float nextFloatRef(float v)
{
    uint32_t vBits = reinterpret_cast<uint32_t&>(v);
    vBits++;
    return reinterpret_cast<float&>(vBits);
}

float nextFloat(float v)
{
    return v + v * nextFloatRef(FLT_EPSILON / 2);
}

float minV = FLT_MIN / (FLT_EPSILON / 2);

nextFloatRef(FLT_EPSILON / 2) is a constant so can be precomputed.

Lieven
  • 101
  • 2