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

I can never remember which of the popular sites is more untrustworthy:
  http://stackoverflow.com/questions/6520052/whats-wrong-with-cplusplus-com
so I confirmed it here:

http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1256.pdf
|
| The round functions round their argument to the nearest integer
| value in floating-point format, rounding halfway cases away from
| zero, regardless of the current rounding direction.

That means we don't want cephes, readily accessible e.g. here:
  https://github.com/deepmind/torch-cephes/blob/master/cephes/cmath/round.c
| Round double to nearest or even integer valued double
                  ^^^^^^^^^^^^^^^
(unless we specifically want that rounding direction even though
the C committee chose a different one).

As this article
  https://lists.freebsd.org/pipermail/freebsd-standards/2004-January/000359.html
points out, conformance is guaranteed if the round() implementation is
copied from the standard [F.9.6.6/2]...as long as #pragma FENV_ACCESS
and the library functions used in that reference implementation are
themselves correctly implemented; but I wouldn't assume those
assumptions are reliable for MinGW-w64.

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

That sounds like our best idea, and I'll point out that the newlib
implementation
  
https://sourceware.org/git/gitweb.cgi?p=newlib-cygwin.git;a=blob;f=newlib/libm/common/sf_round.c
(free source from sun) looks very similar, while what BSD does:
  
https://android.googlesource.com/platform/bionic/+/6861c6f/libm/upstream-freebsd/lib/msun/src/s_roundl.c
  
https://github.com/openbsd-unofficial/openbsd-src/blob/master/lib/libm/src/s_roundl.c
seems to be less carefully thought out. (For example, they don't
have a special case for floating-point values whose exponent is
large enough to shift all mantissa bits well to the left of the
binary point, in which case they're integers and adding 0.5 does
nothing.)

BTW, there's some interesting discussion here:

http://blog.frama-c.com/index.php?post/2013/05/02/nearbyintf1
http://blog.frama-c.com/index.php?post/2013/05/03/nearbyintf2
http://blog.frama-c.com/index.php?post/2013/05/04/nearbyintf3