IEEE-754 compliant round-half-to-even

Question

The C standard library provides the round, lround, and llround family of functions in C99. However, these functions are not IEEE-754 compliant, because they do not implement the "banker's rounding" of half-to-even as mandated by IEEE. Half-to-even rounding requires the result to be rounded to the nearest even value if the fractional component is exactly 0.5. The C99 standard instead mandates half-away-from-zero as noted on cppreference.com

1-3) Computes the nearest integer value to arg (in floating-point format), rounding halfway cases away from zero, regardless of the current rounding mode.

The usual ad-hoc way to implement rounding in C is the expression (int)(x + 0.5f) which, despite being incorrect in strict IEEE-754 math, is usually translated by compilers into the correct cvtss2si instruction. However, this is certainly not a portable assumption.

How can I implement a function that will round any floating point value with half-to-even semantics? If possible, the function should only rely upon language and standard library semantic, so that it can operate on non-IEEE floating point types. If this is not possible, an answer defined in terms of IEEE-754 bit representations is also acceptable. Please characterize any constants in terms of <limits.h> or <limits>.

Answer 1

Round the number x, and if the difference between x and round (x) is exactly +0.5 or -0.5, and round (x) is odd, then round (x) was rounded in the wrong direction, so you subtract the difference from x.

Answer 2

Use remainder(double x, 1.0) from the C standard library. This is independent of the current rounding mode.

The remainder functions compute the remainder x REM y required by IEC 60559

remainder() is useful here as is meets OP's ties to even requirement.

---

double round_to_nearest_ties_to_even(double x) {
  x -= remainder(x, 1.0);
  return x;
}

Test code

void rtest(double x) {
  double round_half_to_even = round_to_nearest_ties_to_even(x);
  printf("x:%25.17le   z:%25.17le \n", x, round_half_to_even);
}

void rtest3(double x) {
  rtest(nextafter(x, -1.0/0.0));
  rtest(x);
  rtest(nextafter(x, +1.0/0.0));
}

int main(void) {
  rtest3(-DBL_MAX);
  rtest3(-2.0);
  rtest3(-1.5);
  rtest3(-1.0);
  rtest3(-0.5);
  rtest(nextafter(-0.0, -DBL_MAX));
  rtest(-0.0);
  rtest(0.0);
  rtest(nextafter(0.0, +DBL_MAX));
  rtest3(0.5);
  rtest3(1.0);
  rtest3(1.5);
  rtest3(2.0);
  rtest3(DBL_MAX);
  rtest3(0.0/0.0);
  return 0;
}

Output

x:                     -inf   z:                     -inf 
x:-1.79769313486231571e+308   z:-1.79769313486231571e+308 
x:-1.79769313486231551e+308   z:-1.79769313486231551e+308 
x: -2.00000000000000044e+00   z: -2.00000000000000000e+00 
x: -2.00000000000000000e+00   z: -2.00000000000000000e+00 
x: -1.99999999999999978e+00   z: -2.00000000000000000e+00 
x: -1.50000000000000022e+00   z: -2.00000000000000000e+00 
x: -1.50000000000000000e+00   z: -2.00000000000000000e+00 tie to even
x: -1.49999999999999978e+00   z: -1.00000000000000000e+00 
x: -1.00000000000000022e+00   z: -1.00000000000000000e+00 
x: -1.00000000000000000e+00   z: -1.00000000000000000e+00 
x: -9.99999999999999889e-01   z: -1.00000000000000000e+00 
x: -5.00000000000000111e-01   z: -1.00000000000000000e+00 
x: -5.00000000000000000e-01   z:  0.00000000000000000e+00 tie to even 
x: -4.99999999999999944e-01   z:  0.00000000000000000e+00 
x:-4.94065645841246544e-324   z:  0.00000000000000000e+00 
x: -0.00000000000000000e+00   z:  0.00000000000000000e+00 
x:  0.00000000000000000e+00   z:  0.00000000000000000e+00 
x: 4.94065645841246544e-324   z:  0.00000000000000000e+00 
x:  4.99999999999999944e-01   z:  0.00000000000000000e+00 
x:  5.00000000000000000e-01   z:  0.00000000000000000e+00 tie to even 
x:  5.00000000000000111e-01   z:  1.00000000000000000e+00 
x:  9.99999999999999889e-01   z:  1.00000000000000000e+00 
x:  1.00000000000000000e+00   z:  1.00000000000000000e+00 
x:  1.00000000000000022e+00   z:  1.00000000000000000e+00 
x:  1.49999999999999978e+00   z:  1.00000000000000000e+00 
x:  1.50000000000000000e+00   z:  2.00000000000000000e+00 tie to even 
x:  1.50000000000000022e+00   z:  2.00000000000000000e+00 
x:  1.99999999999999978e+00   z:  2.00000000000000000e+00 
x:  2.00000000000000000e+00   z:  2.00000000000000000e+00 
x:  2.00000000000000044e+00   z:  2.00000000000000000e+00 
x: 1.79769313486231551e+308   z: 1.79769313486231551e+308 
x: 1.79769313486231571e+308   z: 1.79769313486231571e+308 
x:                      inf   z:                      inf 
x:                      nan   z:                      nan 
x:                      nan   z:                      nan 
x:                      nan   z:                      nan 

Answer 3

The C standard library provides the round, lround, and llround family of functions in C99. However, these functions are not IEEE-754 compliant, because they do not implement the "banker's rounding" of half-to-even as mandated by IEEE...

It doesn't make sense to talk about whether or not an individual function is "IEEE-754 compliant". IEEE-754 compliance requires that a set of data types operations with defined semantics be available. It does not require that those types or operations have specific names, nor does it require that only those operations be available. An implementation can provide whatever additional functions it wants and still be compliant. If an implementation wants to provide round-to-odd, round-random, round-away-from-zero, and trap-if-inexact, it can do so.

What IEEE-754 actually requires for rounding is that the following six operations are provided:

convertToIntegerTiesToEven(x)

convertToIntegerTowardZero(x)

convertToIntegerTowardPositive(x)

convertToIntegerTowardNegative(x)

convertToIntegerTiesToAway(x)

convertToIntegerExact(x)

In C and C++, the last five of these operations are bound to the trunc, ceil, floor, round, and rint functions, respectively. C11 and C++14 do not have a binding for the first, but future revisions will use roundeven. As you can see, round actually is one of the required operations.

However, roundeven is not available in current implementations, which brings us to the next part of your question:

The usual ad-hoc way to implement rounding in C is the expression (int)(x + 0.5f) which, despite being incorrect in strict IEEE-754 math, is usually translated by compilers into the correct cvtss2si instruction. However, this is certainly not a portable assumption.

The problems with that expression extend well-beyond "strict IEEE-754 math". It's totally incorrect for negative x, gives the wrong answer for nextDown(0.5), and turns all odd integers in the 2**23 binade into even integers. Any compiler that translates it into cvtss2si is horribly, horribly broken. If you have an example of that happening, I would love to see it.

How can I implement a function that will round any floating point value with half-to-even semantics?

As njuffa noted in a comment, you can ensure that the default rounding mode is set and use rint (or lrint, as it sounds like you actually want an integer result), or you can implement your own rounding function by calling round and then fixing up halfway cases like gnasher729 suggests. Once the n1778 bindings for C are adopted, you'll be able to use the roundeven or fromfp functions to do this operation without needing to control the rounding mode.

Answer 4

The float data type can represent all whole numbers, but no fractions, within the range 8388608.0f to 16777216.0f. Any float numbers which are larger than 8388607.5f are whole numbers, and no rounding will be necessary. Adding 8388608.0f to any non-negative float which is smaller than that will yield a whole number which will be rounded according to the current rounding mode (typically round-half-to-even). Subtracting 8388608.0f will then yield a properly-rounded version of the original (assuming it was in a suitable range).

Thus, it should be possible to do something like:

float round(float f)
{
  if (!(f > -8388608.0f && f < 8388608.0f)) // Return true for NaN
    return f;
  else if (f > 0)
    return (float)(f+8388608.0f)-8388608.0f;
  else
    return (float)(f-8388608.0f)+8388608.0f;
}

and take advantage of the natural rounding behavior of addition, without having to use any other "round-to-integer" facility.

Answer 5

The following is a simple implementation of round-half to even program which follows the IEEE standard of rounding.

Logic : error = 0.00001

1. number=2.5
2. temp = floor(2.5)%2 = 2%2 = 0
3. x = -1 + temp = -1
4. x*error + number = 2.40009
5. round(2.40009) = 2

Note: The error here is of 0.00001, i.e, if 2.500001 occurs, then it will round off to 2 instead of 3

Python 2.7 implementation:

temp = (number)     
rounded_number = int( round(-1+ temp%2)*0.00001 + temp )

C++ implementation: (Use math.h for floor function)

float temp = (number)     
int rounded_number = (int)( (-1+ temp%2)*0.00001 + temp + 0.5)

The output this would give would be as follows acc. to standards :

(3.5) -> 4

(2.5) -> 2

---

Edit 1 : As pointed out by @Mark Dickinson in the comments. The error can be modified as required in your code to standardize it. For python, to turn it into the smallest float value possible , you may do the following.

import sys
error = sys.float_info.min

Answer 6

Since C++11, there has a been a function in the STL that does half-even rounding. If the floating-point rounding mode is set to FE_TONEAREST (the default), then std::nearbyint will do the trick.

C++ Reference for std::nearbyint