Re: Re: musl mathematical functions

[Szabolcs Nagy <nsz@port70.net>]

* paul zimmermann <Paul.Zimmermann@inria.fr> [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 <Szabolcs.Nagy@arm.com>
> > CC: nd <nd@arm.com>, "jens.gustedt@inria.fr" <jens.gustedt@inria.fr>,
> > 	"Vincent.Lefevre@ens-lyon.fr" <Vincent.Lefevre@ens-lyon.fr>,
> > 	"musl@lists.openwall.com" <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 <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.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
> > > 
> > > 
> > > 
> >