Re: [lmi] MinGW-w64 roundl() anomaly

[Vadim Zeitlin]

On Thu, 15 Dec 2016 16:02:54 +0000 Greg Chicares <address@hidden> wrote:

GC> 'round_test.cpp' doesn't pretend to be a comprehensive validation suite
GC> for the round() family, but it always succeeded with mingw.org's RTL
GC> supplement (libmingwex), yet it fails with MinGW-w64. AFAICS, the only
GC> problem it detects is
GC>   roundl(0.499999999999999999973)
GC> which I believe should always return 0.0L regardless of the current
GC> hardware rounding mode. It seems to do that iff the current mode is
GC> upward rounding--not if it's downward, toward zero, or to nearest.

 Yes, I can confirm all of the above and I can also easily reproduce the
problem outside of lmi with

---------------------------------- >8 --------------------------------------
#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
}
---------------------------------- >8 --------------------------------------

 BTW, initially I used "%La" instead of "%.20Lf" as format specifier and
was very surprised by the output, until I found the existing (and still
open) bug report at https://sourceforge.net/p/mingw-w64/bugs/459/ which
confirms that "%a" outputs "1.0" completely incorrectly when using this
compiler CRT. So just a word of warning: it's better to avoid "%a" with it,
even if I usually find it very useful when debugging floating point
problems.

 I also checked that the problem still happened with the latest MinGW-w64
version (6.2.0) and that it only happened with "long double" and roundl()
(whether called explicitly like this or as std::round() via overload
resolution), but this doesn't really help, of course...

GC> Because this function is supposed to round to nearest always,

 Yes, http://en.cppreference.com/w/cpp/numeric/math/round confirms it:

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

GC> I'd expect its implementation to start by setting that mode, so an
GC> error like this seems strange; but perhaps it's easy to fix.

 Looking at its source code:

https://sourceforge.net/p/mingw-w64/mingw-w64/ci/master/tree/mingw-w64-crt/math/roundl.c

it doesn't do anything with the rounding mode, but it's pretty naïve, e.g.
compared with the implementation in glibc (I _think_ this is the one which
is used for i386, although it's hard to be sure with the number of aliases
and indirections in glibc sources):

https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/ldbl-96/s_roundl.c

 FWIW we could probably use the version above as a workaround for now.


GC> Vadim, would you be interested in verifying my analysis and, if it's
GC> correct, following up with the MinGW-w64 people to get this fixed?

 I've opened a bug report at https://sourceforge.net/p/mingw-w64/bugs/573/
but I'm not sure if it's going to get fixed soon because, looking at the
history of changes for the files in math subdirectory, there doesn't seem
to be much active maintenance happening there. I don't know if the code can
be just copied from glibc (is MinGW-w64 licence compatible?) and it risks
being nontrivial redoing it on my own otherwise... Let's see if there are
going to be any replies (the bug tracker doesn't seem to very active
neither...) but if not, I think using the code from glibc in lmi itself
could be a good enough solution for now.

 What do you think?
VZ