[<prev] [next>] [<thread-prev] [thread-next>] [day] [month] [year] [list]
Message-ID: <20200110173023.GS23985@port70.net>
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> [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@....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>
> > 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@....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$ ./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@...ate:~/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
> > >
> > >
> > >
> >
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.
