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

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 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