CODEWITHSUNDEEP
floating-point
pgast
pgast

Reputation: 1677

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

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 ?

Upvotes: 6

Views: 2302

Answers (5)

Stephen Canon
Stephen Canon

Reputation: 106207

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:

double x = 1.0;
uint64_t rep;
assert(sizeof x == sizeof rep);
memcpy(&rep, &x, sizeof x);
rep += 1;
memcpy(&x, &rep, sizeof x);

This adds one to the significand of x by operating on the bitwise-representation of the floating-point value; 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.

Upvotes: 6

aka.nice
aka.nice

Reputation: 9382

1.0 - epsilon is not the predecessor of 1.0, so the negative counter part does not work at all...
The predecessor of 1.0 is 1.0-epsilon/2.0

Upvotes: 0

pgast
pgast

Reputation: 1677

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.

Upvotes: 1

Jitse Niesen
Jitse Niesen

Reputation: 4532

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.

Upvotes: 3

Robert Kern
Robert Kern

Reputation: 13430

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.

Upvotes: 8

Related Questions