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.