[<prev] [next>] [<thread-prev] [day] [month] [year] [list]
Message-ID: <20200118201435.GH23985@port70.net>
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> [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;
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.
