From mboxrd@z Thu Jan 1 00:00:00 1970 X-Msuck: nntp://news.gmane.org/gmane.linux.lib.musl.general/15117 Path: news.gmane.org!.POSTED.blaine.gmane.org!not-for-mail From: Szabolcs Nagy Newsgroups: gmane.linux.lib.musl.general Subject: Re: Re: musl mathematical functions Date: Fri, 10 Jan 2020 18:30:23 +0100 Message-ID: <20200110173023.GS23985@port70.net> References: <5cf6c0cf-6988-daac-8b74-43bb0e2c625f@arm.com> Reply-To: musl@lists.openwall.com Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Injection-Info: blaine.gmane.org; posting-host="blaine.gmane.org:195.159.176.226"; logging-data="36492"; mail-complaints-to="usenet@blaine.gmane.org" User-Agent: Mutt/1.10.1 (2018-07-13) Cc: Szabolcs Nagy , nd@arm.com, jens.gustedt@inria.fr, Vincent.Lefevre@ens-lyon.fr, musl@lists.openwall.com To: paul zimmermann Original-X-From: musl-return-15133-gllmg-musl=m.gmane-mx.org@lists.openwall.com Fri Jan 10 18:35:36 2020 Return-path: Envelope-to: gllmg-musl@m.gmane-mx.org Original-Received: from mother.openwall.net ([195.42.179.200]) by blaine.gmane.org with smtp (Exim 4.89) (envelope-from ) id 1ipy7M-000liv-4u for gllmg-musl@m.gmane-mx.org; Fri, 10 Jan 2020 18:30:40 +0100 Original-Received: (qmail 11647 invoked by uid 550); 10 Jan 2020 17:30:35 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: List-ID: Original-Received: (qmail 11626 invoked from network); 10 Jan 2020 17:30:35 -0000 Mail-Followup-To: paul zimmermann , Szabolcs Nagy , nd@arm.com, jens.gustedt@inria.fr, Vincent.Lefevre@ens-lyon.fr, musl@lists.openwall.com Content-Disposition: inline In-Reply-To: Xref: news.gmane.org gmane.linux.lib.musl.general:15117 Archived-At: * paul zimmermann [2020-01-10 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 non-nearest 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 Floating-Point 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 > > CC: nd , "jens.gustedt@inria.fr" , > > "Vincent.Lefevre@ens-lyon.fr" , > > "musl@lists.openwall.com" > > Thread-Topic: musl mathematical functions > > Thread-Index: AQHVxieg5o3AZI5d3UWouuHlhctYy6fg5FoA > > Date: Wed, 8 Jan 2020 15:28:54 +0000 > > user-agent: Mozilla/5.0 (X11; Linux aarch64; rv:60.0) Gecko/20100101 > > Thunderbird/60.9.0 > > nodisclaimer: True > > Original-Authentication-Results: spf=none (sender IP is ) > > smtp.mailfrom=Szabolcs.Nagy@arm.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@tomate:~/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@tomate:~/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 > > > > 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 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. > > > > thanks. > > > > > > > > Indeed with the following program: > > > > > > #include > > > #include > > > #include > > > #include > > > > > > 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.1766512962000004e-14 > > > > > > and with musl: > > > > > > $ ./a.out 4.2725660088821189e2 > > > sin(4.2725660088821190e+02) = -2.2563645396544984e-11 > > > > > > which is indeed very far from the correctly rounded result. > > > > > > Best regards, > > > Paul Zimmermann > > > > > > > > > > >