
Date: Wed, 8 Jan 2020 15:28:54 +0000 From: Szabolcs Nagy <Szabolcs.Nagy@....com> To: paul zimmermann <Paul.Zimmermann@...ia.fr> CC: nd <nd@....com>, "jens.gustedt@...ia.fr" <jens.gustedt@...ia.fr>, "Vincent.Lefevre@...lyon.fr" <Vincent.Lefevre@...lyon.fr>, "musl@...ts.openwall.com" <musl@...ts.openwall.com> Subject: Re: musl mathematical functions 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$ ./mpcheckdouble 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$ ./mpcheckdouble 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/libctest.git https://github.com/ARMsoftware/optimizedroutines ) these issues come from fdlibm via freebsd which does not support nonnearest rounding in the trig arg reduction code (and possibly in other places). http://git.musllibc.org/cgit/musl/tree/src/math/__rem_pio2.c#n120 (note the comment: assume roundtonearest) 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 Shift = 0x1.8p52 y = x/p + Shift  Shift implement round as e.g. Shift = 0x1800000000.8p0 t = x/p + Shift tbits = representation_as_uint64(t) y = (double)(int32_t)(tbits >> 16) then z is in [p/2  p/2^16, p/2 + p/2^16] in all rounding modes and the polynomial can be made to work on that interval. the downside is that these tricks make the code slower and more importantly all such tricks break symmetry: x and x can have different arg reduction result. now i lean towards fixing it in a way that's least expensive in the nearestrounding 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. thanks. > > Indeed with the following program: > > #include <stdio.h> > #include <stdlib.h> > #include <math.h> > #include <fenv.h> > > int > main (int argc, char *argv[]) > { > double x = atof (argv[1]), y; > fesetround (FE_UPWARD); > y = sin (x); > printf ("sin(%.16e) = %.16e\n", x, y); > } > > I get with the GNU libc: > > $ ./a.out 4.2725660088821189e2 > sin(4.2725660088821190e+02) = 1.1766512962000004e14 > > and with musl: > > $ ./a.out 4.2725660088821189e2 > sin(4.2725660088821190e+02) = 2.2563645396544984e11 > > which is indeed very far from the correctly rounded result. > > Best regards, > Paul Zimmermann > > >
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.