
Date: Sat, 18 Jan 2020 21:14:35 +0100 From: Szabolcs Nagy <nsz@...t70.net> To: paul zimmermann <Paul.Zimmermann@...ia.fr> Cc: 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> [20200110 19:35:08 +0100]: > > i think libm functions are extremely rarely used with > > nonnearest 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 186614 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.921fb54442d18p1, invpio2 = 6.36619772367581382433e01, /* 0x3FE45F30, 0x6DC9C883 */ pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ pio2_1t = 6.07710050650619224932e11, /* 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;
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.