#573 roundl() doesn't work correctly for numbers close to halfway point

[Vadim Zeitlin]

roundl() depends on the current rounding mode and rounds incorrectly by default: this function is not supposed to depend on it in the first place (see http://en.cppreference.com/w/cpp/numeric/math/round) and it should round the numbers smaller than, even if close to, 0.5 down and not up as it does: the result of roundl(nextafterl(0.5L, 0.0L)) is 1 in almost all rounding modes and is only the correct 0 when using FE_UPWARD, ironically enough.

Here is a small test program that can be used to reproduce this bug:

#include <cfenv>
#include <cmath>
#include <stdio.h>

void round_with_mode(int mode, char const* name, long double x)
{
    fesetround(mode);
    printf("Using %-13s rounding mode roundl(%.20Lf) = %.20Lf\n", name, x, roundl(x));
}

int main()
{
    const auto x = nextafterl(0.5L, 0.0L);

#define CALL_ROUND_WITH_MODE(mode, x) round_with_mode(mode, #mode, x)

    CALL_ROUND_WITH_MODE(FE_DOWNWARD, x);
    CALL_ROUND_WITH_MODE(FE_TONEAREST, x);
    CALL_ROUND_WITH_MODE(FE_TOWARDZERO, x);
    CALL_ROUND_WITH_MODE(FE_UPWARD, x);

#undef CALL_ROUND_WITH_MODE
}

It gives the same, incorrect, results with both 4.9.1 that we usually use and the latest 6.2.0:

Using FE_DOWNWARD   rounding mode roundl(0.49999999999999999997) = 1.00000000000000000000
Using FE_TONEAREST  rounding mode roundl(0.49999999999999999997) = 1.00000000000000000000
Using FE_TOWARDZERO rounding mode roundl(0.49999999999999999997) = 1.00000000000000000000
Using FE_UPWARD     rounding mode roundl(0.49999999999999999997) = 0.00000000000000000000

Notice that the problem only exists with long double, everything works correctly with just plain double. The same program also shows the expected results when using the same compiler versions under Linux.

Looking at this function implementation it seems to be quite naive, especially compared to the implementation in glibc. I don't know how does this work, but wouldn't it be possible to reuse glibc code?

If not, there is the trivial solution of always setting the rounding more to FE_UPWARD at the beginning of this function and restoring it at the end, but this is clearly not going to be optimal.

Please let me know if I can do anything to help with fixing this, TIA!

[Doug Semler]

The runtime can't "just use" the glibc implementation because of licensing. In the versions (6.3 and 7.1) I have tested, the problem exists at point nextafter (0.5, 0) for all versions of float, double, and long double.

Anyway, the reason for the bug is the (res - x) computation...non constant floating point operations take the rounding mode into account. So you end up with internal calculation of 1.0 - 0.4999999999997 in the FPU, which depending on the rounding mode can get you exactly 0.50000000000000 (round down, nearest, and towards zero). This is also why other "trivial" solutions such as trunc(x + 0.5) doesn't work for all positive values, etc. It's the MATH operation that screws things up with its rounding mode if the result happens to round differently in the least significant bit of the converted value.

Isn't floating point fun?

Sandboxing without looking at the glibc, it seems to me that the mingw runtime needs to avoid doing math operations for all the floating point functions that do rounding, and look at the various bits of the floating point value to get the sign, exponent, mantissa, etc and mask/bit-bang things to be correct across all values. Otherwise, floating point values one epsilon away from midpoints may violate IEEE and std rounding.