
Date: Fri, 10 Jan 2020 18:30:23 +0100 From: Szabolcs Nagy <nsz@...t70.net> To: paul zimmermann <Paul.Zimmermann@...ia.fr> Cc: Szabolcs Nagy <Szabolcs.Nagy@....com>, nd@....com, jens.gustedt@...ia.fr, Vincent.Lefevre@...lyon.fr, musl@...ts.openwall.com Subject: Re: Re: musl mathematical functions * paul zimmermann <Paul.Zimmermann@...ia.fr> [20200110 17:01:43 +0100]: > Dear Szabolcs, > > thank you for your answer. > > I understand the issues of slowing down the code and/or breaking symmetry, > but in my opinion the ordering should be: > > accuracy >> symmetry >> speed > > where "x >> y" means that "x is more important than y". what do you think about directed rounding mode behaviour? i think libm functions are extremely rarely used with nonnearest rounding mode so i think NR accuracy >> DR accuracy >> NR symmetry >> NR speed >> DR symmetry >> DR speed where NR is nearest rounding and DR is directed rounding. and by accuracy i just mean correct behavirour with respect to exceptions and results (i.e. small ulp errors). > > Maybe you can find some tricks in the "Handbook of FloatingPoint Arithmetic"? > > Note that our mpcheck tool can also check for symmetry. > > Anyway, if you do some changes, I'll be happy to run mpcheck again and send > you the new results. thanks. > > Best regards, > Paul > > > From: Szabolcs Nagy <Szabolcs.Nagy@....com> > > 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> > > ThreadTopic: musl mathematical functions > > ThreadIndex: AQHVxieg5o3AZI5d3UWouuHlhctYy6fg5FoA > > Date: Wed, 8 Jan 2020 15:28:54 +0000 > > useragent: Mozilla/5.0 (X11; Linux aarch64; rv:60.0) Gecko/20100101 > > Thunderbird/60.9.0 > > nodisclaimer: True > > OriginalAuthenticationResults: spf=none (sender IP is ) > > smtp.mailfrom=Szabolcs.Nagy@....com; > > > > 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.