From mboxrd@z Thu Jan 1 00:00:00 1970 X-Msuck: nntp://news.gmane.org/gmane.linux.lib.musl.general/15104 Path: news.gmane.org!.POSTED.blaine.gmane.org!not-for-mail From: Rich Felker Newsgroups: gmane.linux.lib.musl.general Subject: Re: Re: musl mathematical functions Date: Wed, 8 Jan 2020 10:46:28 -0500 Message-ID: <20200108154628.GS30412@brightrain.aerifal.cx> 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="252668"; mail-complaints-to="usenet@blaine.gmane.org" User-Agent: Mutt/1.5.21 (2010-09-15) To: musl@lists.openwall.com Original-X-From: musl-return-15120-gllmg-musl=m.gmane.org@lists.openwall.com Wed Jan 08 16:47:43 2020 Return-path: Envelope-to: gllmg-musl@m.gmane.org Original-Received: from mother.openwall.net ([195.42.179.200]) by blaine.gmane.org with smtp (Exim 4.89) (envelope-from ) id 1ipDXf-000Ofi-Fg for gllmg-musl@m.gmane.org; Wed, 08 Jan 2020 16:46:43 +0100 Original-Received: (qmail 27747 invoked by uid 550); 8 Jan 2020 15:46:41 -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 27729 invoked from network); 8 Jan 2020 15:46:40 -0000 Content-Disposition: inline In-Reply-To: <5cf6c0cf-6988-daac-8b74-43bb0e2c625f@arm.com> Original-Sender: Rich Felker Xref: news.gmane.org gmane.linux.lib.musl.general:15104 Archived-At: 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@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 > > [...] > > 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