Is there a standard way to get the nth `nextafter` floating-point value in C++

Question

C++ has std::nextafter(), which returns the next representable value after a given floating-point value f. In my case, I'd like to allow for n bits of slop in the lower mantissa bits, so 3 bits of slop would require getting the 8th next value after some given value f. I could call nextafter() eight times, but is there a better way to handle this?

For most values, you could get by with casting the FP value to uint_64, adding the tolerance (1<<3 for 3 bits of slop), and then casting back to double, thanks to the layout of IEEE 754. However, that relies on IEEE 754 floating point (a good assumption, but not rock solid either).

(For background, I want to use this to hoist ray-surface intersection points, which are occasionally located inside the surface due to FP imprecision. Those familiar with robust floating-point will understand why epsilon is a terrible solution.)

Answer 1

As far as I know, there is no standard function for this. Boost has you covered though: See boost::math::float_advance. If you're using this for comparing two floats, you probably want boost::math::float_distance instead.

Answer 2

A naive approach could be to multiply by 8 the distance between a value and the next representable float, instead of calling 8 times std::nextafter

double advance_float(double x, int d)
{
    double step = std::copysign((std::nextafter(x, x + d) - x) * d, d);
    return x + step;
}

Here there are some tests, but it's up to you to determine if this is suitable for your use case.

Edit

As noted by Steve Hollash, x may be so big that x + d == d. Daniel Jour suggested to take advantage of frexp (and ldexp), but in the following attempt, I'll use a different approach to determine the direction.

double advance_float(double x, int d)
{
    const double to = std::copysign(std::numeric_limits<double>::infinity(), d);
    const double next = std::nextafter(x, to);
    return x + std::copysign(d * (next - x), d);
}

Note that it assumes that std::numeric_limits<double>::has_infinity == true, otherwise ::lowest() and ::max() must be used.

Those are some results

         x    d             previous                     x                   next
------------------------------------------------------------------------------------------
           1  1     0x1.fffffffffffffp-1                   0x1p+0     0x1.0000000000001p+0
           1  8     0x1.ffffffffffff8p-1                   0x1p+0     0x1.0000000000008p+0
     3.14159  8      0x1.921fb54442d1p+1     0x1.921fb54442d18p+1      0x1.921fb54442d2p+1
      100.01  8     0x1.900a3d70a3d69p+6     0x1.900a3d70a3d71p+6     0x1.900a3d70a3d79p+6
     -100.01  8    -0x1.900a3d70a3d79p+6    -0x1.900a3d70a3d71p+6    -0x1.900a3d70a3d69p+6
       1e+67  8   0x1.7bd29d1c87a11p+222   0x1.7bd29d1c87a19p+222   0x1.7bd29d1c87a21p+222
       1e-59  8    0x1.011c2eaabe7dp-196   0x1.011c2eaabe7d8p-196    0x1.011c2eaabe7ep-196
           0  8 -0x0.0000000000008p-1022                   0x0p+0  0x0.0000000000008p-1022
4.94066e-324  8 -0x0.0000000000007p-1022  0x0.0000000000001p-1022  0x0.0000000000009p-1022