[<prev] [next>] [<thread-prev] [thread-next>] [day] [month] [year] [list]
Date: Wed, 8 Jan 2020 10:46:28 -0500
From: Rich Felker <dalias@...c.org>
To: musl@...ts.openwall.com
Subject: Re: Re: musl mathematical functions
On Wed, Jan 08, 2020 at 03:28:54PM +0000, Szabolcs Nagy wrote:
> On 08/01/2020 13:29, paul zimmermann wrote:
> > Dear Szabolcs,
> >
> > my colleague Jens Gustedt told me that you lead the development of mathematical
> > functions in musl.
> >
> > I just tried our mpcheck tool (https://gforge.inria.fr/projects/mpcheck) which
> > checks the accuracy of mathematical functions, by comparing them to MPFR (which
> > is supposed to give correct rounding).
>
> thanks!
>
> CCing the musl list as it should be discussed there.
>
> >
> > For the GNU libc here is what I get for example for double precision
> > (with 10000 random inputs per function):
> >
> > zimmerma@...ate:~/svn/mpcheck$ ./mpcheck-double --seed=588493
> > GCC: 9.2.1 20200104
> > GNU libc version: 2.29
> > GNU libc release: stable
> > MPFR 3.1.6
> > ...
> > Max. errors : 3.59 (nearest), 5.80 (directed) [seed=588493]
> > Incorrect roundings: 483084 (basic 0)
> > Wrong side of directed rounding: 245029
> > Wrong monotonicity: 31701
> > Wrong errno: 992 (suppressed 992)
> > Wrong inexact: 11 (suppressed 11)
> > Wrong underflow: 42 (suppressed 42)
> >
> > This means (among other things) that the maximal error found on those random
> > inputs is 3.59 ulps for rounding to nearest, and 5.80 ulps for directed
> > rounding.
> >
> > With musl (revision 70d8060) I get:
> >
> > zimmerma@...ate:~/svn/mpcheck$ ./mpcheck-double --seed=588493
> > GCC: 9.2.1 20200104
> > MPFR 3.1.6
> > ...
> > Max. errors : 5.30 (nearest), 1.44e19 (directed) [seed=588493]
> > Incorrect roundings: 407422 (basic 0)
> > Wrong side of directed rounding: 130496
> > Wrong errno: 131411 (suppressed 10901)
> > Wrong inexact: 125 (suppressed 125)
> > Wrong overflow: 16 (suppressed 0)
> > Wrong underflow: 178 (suppressed 108)
> >
> > We get a slightly larger maximal error for rounding to nearest (5.30 instead
> > of 3.59 for the GNU libc) but a huge maximal error for directed rounding.
> >
> > The 1.44e19 error is obtained for the "sin" function, with input
> > x=4.2725660088821189e2 and rounding upwards.
>
> yes, this is a known issue (the math tests i use with
> musl finds this, but it's suppressed for now
> https://repo.or.cz/w/libc-test.git
> https://github.com/ARM-software/optimized-routines
> )
>
> these issues come from fdlibm via freebsd which
> does not support non-nearest rounding in the trig
> arg reduction code (and possibly in other places).
> http://git.musl-libc.org/cgit/musl/tree/src/math/__rem_pio2.c#n120
> (note the comment: assume round-to-nearest)
>
> i haven't fixed this because i don't have a good
> solution: the key broken part is something like
>
> y = round(x/p)
> z -= y*p
> /* i.e. z = x mod p, and z in [-p/2,p/2] */
> return poly(z)
>
> the problem is that the fast and portable way to
> do round relies on the current rounding mode and
> z can end up in the range [-p,p] with directed
> rounding, but the poly approx only works on
> [-p/2,p/2].
>
> some targets have single instruction round that's
> independent of the rounding mode, but most targets
> don't.
>
> changing fenv is slower than just calling round or
> rint, and doing an external call is already too
> expensive.
>
> one can do tricks such that rounding mode has
> less effect on arg reduction, e.g. add
>
> if (z > p/2 || z < -p/2) /* do something */
>
> or if branches are too expensive, instead of
>
> [...]
>
> now i lean towards fixing it in a way that's
> least expensive in the nearest-rounding case
> (at least for |x| < 100, beyond that performance
> does not matter much) and only care about
> symmetry in nearest rounding mode, this should
> be doable by adding a few ifs in the critical
> path that never trigger with nearest rounding.
>
> but other ideas are welcome.
I'm in favor of a fix that's an entirely predictable branch on default
rounding mode; on most modern cpus it should have zero measurable
performance cost. If you're concerned it might matter, we can look at
some measurements.
Rich
Powered by blists - more mailing lists
Confused about mailing lists and their use?
Read about mailing lists on Wikipedia
and check out these
guidelines on proper formatting of your messages.
