*Re: musl mathematical functions[not found] <mwlfqit6tx.fsf@tomate.loria.fr>@ 2020-01-08 15:28 ` Szabolcs Nagy2020-01-08 15:46 ` Rich Felker 2020-01-10 16:01 ` paul zimmermann 0 siblings, 2 replies; 6+ messages in thread From: Szabolcs Nagy @ 2020-01-08 15:28 UTC (permalink / raw) To: paul zimmermann;+Cc:nd, jens.gustedt, Vincent.Lefevre, musl 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 > > > ^ permalink raw reply [flat|nested] 6+ messages in thread

*Re: Re: musl mathematical functions2020-01-08 15:28 ` musl mathematical functions Szabolcs Nagy@ 2020-01-08 15:46 ` Rich Felker2020-01-10 16:01 ` paul zimmermann 1 sibling, 0 replies; 6+ messages in thread From: Rich Felker @ 2020-01-08 15:46 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 6+ messages in thread

*2020-01-08 15:28 ` musl mathematical functions Szabolcs Nagy 2020-01-08 15:46 ` Rich FelkerRe: musl mathematical functions@ 2020-01-10 16:01 ` paul zimmermann2020-01-10 17:30 ` Szabolcs Nagy 1 sibling, 1 reply; 6+ messages in thread From: paul zimmermann @ 2020-01-10 16:01 UTC (permalink / raw) To: Szabolcs Nagy;+Cc:nd, jens.gustedt, Vincent.Lefevre, musl 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". 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. 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 > > > > > > > ^ permalink raw reply [flat|nested] 6+ messages in thread

*Re: Re: musl mathematical functions2020-01-10 16:01 ` paul zimmermann@ 2020-01-10 17:30 ` Szabolcs Nagy2020-01-10 18:35 ` paul zimmermann 0 siblings, 1 reply; 6+ messages in thread From: Szabolcs Nagy @ 2020-01-10 17:30 UTC (permalink / raw) To: paul zimmermann;+Cc:Szabolcs Nagy, nd, jens.gustedt, Vincent.Lefevre, musl * 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 > > > > > > > > > > > ^ permalink raw reply [flat|nested] 6+ messages in thread

*Re: Re: musl mathematical functions2020-01-10 17:30 ` Szabolcs Nagy@ 2020-01-10 18:35 ` paul zimmermann2020-01-18 20:14 ` [musl] " Szabolcs Nagy 0 siblings, 1 reply; 6+ messages in thread From: paul zimmermann @ 2020-01-10 18:35 UTC (permalink / raw) To: Szabolcs Nagy;+Cc:Szabolcs.Nagy, nd, jens.gustedt, Vincent.Lefevre, musl Dear Szabolcs, > 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. yes this makes sense. > and by accuracy i just mean correct behavirour with respect > to exceptions and results (i.e. small ulp errors). note that if directed rounding is used to implement interval arithmetic, it is very important to have the return value on the right side with respect to the exact value (at the cost of a few ulps of accuracy). Paul ^ permalink raw reply [flat|nested] 6+ messages in thread

*Re: [musl] Re: musl mathematical functions2020-01-10 18:35 ` paul zimmermann@ 2020-01-18 20:14 ` " Szabolcs Nagy0 siblings, 0 replies; 6+ messages in thread From: Szabolcs Nagy @ 2020-01-18 20:14 UTC (permalink / raw) To: paul zimmermann;+Cc:Szabolcs.Nagy, nd, jens.gustedt, Vincent.Lefevre, musl * paul zimmermann <Paul.Zimmermann@inria.fr> [2020-01-10 19:35:08 +0100]: > > 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. > > yes this makes sense. > > > and by accuracy i just mean correct behavirour with respect > > to exceptions and results (i.e. small ulp errors). > > note that if directed rounding is used to implement interval > arithmetic, it is very important to have the return value on > the right side with respect to the exact value (at the cost > of a few ulps of accuracy). getting on the right side would regress the performance for all users for something theoretical (existing math libraries don't support it). at least i think it would require accessing fenv (changing the rounding mode) or other expensive operation in the hot code path. (expensive rounding mode change is one of the reasons interval arithmetics is not practical, lack of compiler support for fenv access is another, the instruction set and language semantics should be designed differently for it to be practical.) i'd recommend using an arbitrary precision library or a correctly rounded math library (e.g. tr 18661-4 reserves separate cr prefixed symbols for that), not libc functions for interval arithmetic. large ulp errors can usually be fixed without significant penalty for nearest rounding. e.g. i'm considering a fix for trig arg reduction with two additional branches (i think this can be simplified with more code changes and the cost can be eliminated on targets with rounding mode independent rounding instructions) diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c index d403f81c..80fd72c8 100644 --- a/src/math/__rem_pio2.c +++ b/src/math/__rem_pio2.c @@ -36,6 +36,7 @@ */ static const double toint = 1.5/EPS, +pio4 = 0x1.921fb54442d18p-1, invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ @@ -122,6 +123,17 @@ medium: n = (int32_t)fn; r = x - fn*pio2_1; w = fn*pio2_1t; /* 1st round, good to 85 bits */ + if (predict_false(r - w < -pio4)) { + n--; + fn--; + r = x - fn*pio2_1; + w = fn*pio2_1t; + } else if (predict_false(r - w > pio4)) { + n++; + fn++; + r = x - fn*pio2_1; + w = fn*pio2_1t; + } y[0] = r - w; u.f = y[0]; ey = u.i>>52 & 0x7ff; ^ permalink raw reply [flat|nested] 6+ messages in thread

end of thread, back to indexThread overview:6+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- [not found] <mwlfqit6tx.fsf@tomate.loria.fr> 2020-01-08 15:28 ` musl mathematical functions Szabolcs Nagy 2020-01-08 15:46 ` Rich Felker 2020-01-10 16:01 ` paul zimmermann 2020-01-10 17:30 ` Szabolcs Nagy 2020-01-10 18:35 ` paul zimmermann 2020-01-18 20:14 ` [musl] " Szabolcs Nagy

mailing list of musl libc Archives are clonable: git clone --mirror http://inbox.vuxu.org/musl Example config snippet for mirrors Newsgroup available over NNTP: nntp://inbox.vuxu.org/vuxu.archive.musl AGPL code for this site: git clone https://public-inbox.org/public-inbox.git