Re: musl mathematical functions

Submitted by Szabolcs Nagy on Jan. 18, 2020, 8:14 p.m.

Details

Message ID 20200118201435.GH23985@port70.net
State New
Series "Re: musl mathematical functions"
Headers show

Commit Message

Szabolcs Nagy Jan. 18, 2020, 8:14 p.m.
* 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)

Patch hide | download patch | download mbox

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;