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

They are compliant with the standard, but do not implement all of it. "Round to nearest, ties away from zero" is one of the IEEE-754 rounding modes. — Ben Voigt, Sep 23 '15 at 18:13
IEEE-754 notes that the schoolbook rounding mode exists, but is only intended for decimal floats. `An implementation of this standard shall provide roundTiesToEven and the three directed rounding attributes. A decimal format implementation of this standard shall provide roundTiesToAway as a userselectable rounding-direction attribute. The rounding attribute roundTiesToAway is not required for a binary format implementation.` Rounding of ties to even is the only portable (and sensible) mode. — 68ejxfcj5669, Sep 23 '15 at 18:17
Instead of at `round()` look at `rint()` used with a default rounding mode of "round to nearest or even". — njuffa, Sep 23 '15 at 18:32
You might be interested in CPython's solution to this, [here](https://github.com/python/cpython/blob/3.5/Objects/floatobject.c#L955-L958). — Mark Dickinson, Sep 23 '15 at 18:33
@njuffa Thanks a lot for the pointer to `rint` (extremely unintuitive name). I had to double check the C standard, but it appears to indeed do the right thing when `__STDC_IEC_559__` is set. The real mystery is why `round` does something different from `rint` and `FE_TONEAREST`. — 68ejxfcj5669, Sep 23 '15 at 19:09

Stephen Canon · Answer 1 · 2015-09-24T11:04:52.347

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.

I agree that (int)(x + 0.5f) is wrong for many reasons. However, MSVC in particular detects this idiom and replaces it with cvtss2si, which I would guess is why it is commonly used. — 68ejxfcj5669, Sep 24 '15 at 17:44
@68ejxfcj5669: I doubt that's the case, but I would love to be proved wrong. I don't have a windows machine to test on, and searching the web doesn't turn up any accounts of this transformation. Do you have an actual example of a program where MSVC did that? — Stephen Canon, Sep 24 '15 at 18:06
Any kind of transformation which changes numeric behavior can only happen under fast math. Which is not an option for anyone who wants well-defined behavior. — stgatilov, Jan 09 '23 at 08:24

chux - Reinstate Monica · Answer 2 · 2015-09-24T04:47:09.567

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

score 6 · Answer 3 · answered Sep 23 '15 at 18:13

6

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.

answered Sep 23 '15 at 18:13

gnasher729

51,477
5
75
98

Is it guaranteed that `(round(x) - x)` will return 0.5/-0.5 if the fractional component of x is exactly half, i.e. is this an exact operation? – 68ejxfcj5669 Sep 23 '15 at 18:19
3

Yes. The result of `round(x) - x` will be exactly representable (assuming an IEEE 754 format). – Mark Dickinson Sep 23 '15 at 18:21
3

An evil trick is to use `2.0 * round(0.5 * x)` in the case where `fabs(x - round(x)) == 0.5`, which saves the somewhat messy task of figuring out whether `round(x)` is even or odd. – Mark Dickinson Sep 23 '15 at 18:24

score 3 · Answer 4 · answered Sep 25 '15 at 22:44

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.

I checked this on all floats in range from INT_MIN to INT_MAX: it provides exactly the same results as `nearbyint` and `_mm_cvtps_epi32` with default rounding. I want to emulate _mm_cvtps_epi32 in non-SIMD code path, but nearbyint is slow as hell! Performance of my DXT1 compressor is 47 MP/s with this implementation, while only 18 MP/s with nearbyint (MSVC 2022). Thanks to C++ committee for providing yet another necessary but unusable math function =( — stgatilov, Jan 09 '23 at 08:19

score 3 · Answer 5 · answered Jul 15 '19 at 19:48

3

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

answered Jul 15 '19 at 19:48

David G

5,408
1
23
19

score 1 · Answer 6 · edited Jun 20 '20 at 09:12

1

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

Logic : error = 0.00001

number=2.5

temp = floor(2.5)%2 = 2%2 = 0

x = -1 + temp = -1

x*error + number = 2.40009

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

edited Jun 20 '20 at 09:12

Community

1
1

answered Dec 23 '17 at 08:07

Shivansh Jagga

1,541
1
15
24

2

If 2.500001 rounds to 2 instead of 3, then this isn't a standards-conforming round. – Mark Dickinson Dec 23 '17 at 08:50
The program is written for someone who wants this implemented in their own code. The error can be changed to any precision point needed in the code,i.e, it can even be the smallest float value possible. `error = sys.float_info.min` (for Python) , and a smiliar equivalent for C++ – Shivansh Jagga Dec 23 '17 at 08:55
1

Okay, but it doesn't answer the question: the OP is looking for an IEEE 754-compliant rounding operation. The code in this answer doesn't give that. And using `error = sys.float_info.min` isn't a solution. If you test your code, you'll find that even the current code doesn't give `4` from an input of `3.5`; it gives `3`. – Mark Dickinson Dec 23 '17 at 09:14

IEEE-754 compliant round-half-to-even

6 Answers6