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

[Greg Chicares]

On 2016-12-16 01:39, Vadim Zeitlin wrote:
> On Thu, 15 Dec 2016 16:02:54 +0000 Greg Chicares <address@hidden> wrote:
[...]
> 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
[...snip example...]
>  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.

Thanks, that's interesting to know.

>  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()
[...]
>  Looking at its source code:
> 
> https://sourceforge.net/p/mingw-w64/mingw-w64/ci/master/tree/mingw-w64-crt/math/roundl.c

A one-character change to that URL gives their round() and roundf()
implementations, which are essentially the same. I think there's a good
chance that their round() would fail for some argument; maybe roundf()
would also.

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

We don't actually call roundl() anywhere in lmi, except to unit-test the
RTL implementation. I don't think we call round() either in production
today, though we do in 'currency.hpp'. Now that we see the MinGW-w64
implementation, maybe we should replace that with lmi's 'round_to'.

Alternatively, we could use mingw.org's implementation:
  
https://sourceforge.net/u/cstrauss/mingwrt/ci/master/tree/mingwex/math/round_internal.h
but I would suspect that glibc's is better tested.

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

All the more reason to avoid division by zero:
  https://sourceforge.net/p/mingw-w64/bugs/270/

> I don't know if the code can
> be just copied from glibc (is MinGW-w64 licence compatible?)

AFAICT, MinGW-w64 disavows copyright, and glibc is LGPL, so both are
compatible with lmi (and we're already using Drepper's md5 code), but I
don't know whether MinGW-w64 could use LGPL code; mingw.org refused to:
  http://www.mingw.org/wiki/c99
| Unfortunately using glibc is not an option because of licence issues

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

I think that's a good idea. Maybe I'll do that over Christmas, just
because I enjoy this sort of thing.