What is the next normalised floating point number after(before) a normalised floating point number f?

Question

Given a normalized floating point number f what is the next normalized floating point number after/before f.

With bit twiddling, extracting mantissa and exponent I have:

next_normalized(double&){
      if mantissa is not all ones
          maximally denormalize while maintaining equality 
          add 1 to mantissa
          normalize
      else 
          check overflow
          set mantissa to 1  
          add (mantissa size in bits) to exponent.
      endif
 }

But rather than do that can it be done with floating point operations?

As

std::numeric_limits<double>::epsilon()

is only an error difference in a "neighborhood" of 1. - e.g.:

normalized(d+=std::numeric_limits<double>::epsilon()) = d for d large

it seems more an error ratio than an error difference, thus my naive intuition is

(1.+std::numeric_limits<double>::epsilon())*f //should be the next.

And

(1.-std::numeric_limits<double>::epsilon())*f //should be the previous.

In particular I have 3 questions has anyone done any of the following (for IEEE754):

1)done the error analysis on this issue?

2)proved (or can prove) that for any normalized double d

    (1.+std::numeric_limits<double>::epsilon())*d != d ?

3)proved that for any normalized double number d no double f exists such that

    d < f < (1.+std::numeric_limits<double>::epsilon())*d ?

Answer 1

I’m not sure what you mean by “normalized double number”, but getting the next representable double number is done with the nextafter() function in most C standard math libraries.

Answer 2

The statement under 3) is false. If d is slightly smaller than 2, then there is exactly 1 floating-point number between d and (1+eps) * d. Here is a program to show it:

#include <limits>
#include <iostream>

int main(int, char**)
{
  using namespace std;
  double d = 1.875;
  cout.precision(18);
  cout << "d = " << d << "\n";
  double d2 = (1.+numeric_limits<double>::epsilon())*d;
  cout << "d2 = " << d2 << "\n";
  double f = d + (d2-d)/2;
  cout << "f = " << f << "\n";
}

The reason is that (1+eps) * 1.875 equals 1.875 + 1.875 * eps, which is rounded to 1.875 + 2 * eps. However, the difference between consecutive floating-point numbers between 1 and 2 is eps, so there is one floating-point number between 1.875 and 1.875 + 2 * eps, namely 1.875 + eps.

The statement under 2) is true, I think. And Robert Kern probably answered your real question.

Answer 3

As Robert Kern noted, you want the C nextafter( ) function, or the IEEE754 nextUp( ) and nextDown( ) functions, though those two are not widely implemented just yet.

If you want to avoid nextafter for some reason, you can do:

double next = x + scalbn(1.0, ilogb(x) - 52);

This adds 2^(exponent of x - 52) to x, which is exactly one unit in the last place (ULP).

If you don't have the usual cmath functions available, and don't mind type-system abuse that the compiler may object to, you can do:

union { double f; uint64_t u; } xbits = { .f = x };
xbits.u += UINT64_C(1);
double next = xbits.f;

(or, I suppose, use a reinterpret cast instead of a union in C++. I'm not a C++ guru.) This adds one to the mantissa of x directly; if the next value is in the next binade, this will carry into the exponent, returning the correct value. If you want it to work for negative values, you'll need to tweak that.

Answer 4

As stated below it turns out after a tiny bit of investigation that for positive floats in intel IEEE754 format of size n-bits that are < +infinity treating the concatenated exponent and significand as a n-1 bit unsigned integer adding one gets the next higher and (subtracting one get next lower)

And vice versa if negative. In particular one can interpret the n-1 bit integer as representing the absolute magnitude independent of the sign. And thus when negative one must subtract one to get the next floating point number closer to zero after the negative floating point number f.

Answer 5

Is this a duplicate of StackOverflow: next higher/lower IEEE double precision number ?