Date: Sun, 19 Mar 2017 15:39:53 +0100 From: Szabolcs Nagy <nsz@...t70.net> To: musl@...ts.openwall.com Subject: Re: [PATCH] fix underflow exception in fma and fmal * Rich Felker <dalias@...c.org> [2017-03-19 10:12:49 -0400]: > On Sun, Mar 19, 2017 at 04:36:14AM +0100, Szabolcs Nagy wrote: > > another corner case in the freebsd fma code where signaling underflow > > may be missed for an inexact subnormal result. > > (fmaf and x86 fma are not affected) > > --- > > src/math/fma.c | 7 +++++++ > > src/math/fmal.c | 8 ++++++++ > > 2 files changed, 15 insertions(+) > > > > diff --git a/src/math/fma.c b/src/math/fma.c > > index 741ccd75..c69918d1 100644 > > --- a/src/math/fma.c > > +++ b/src/math/fma.c > > @@ -279,6 +279,13 @@ static inline double add_and_denormalize(double a, double b, int scale) > > uhi.i += 1 - (((uhi.i ^ ulo.i) >> 62) & 2); > > sum.hi = uhi.f; > > } > > +#ifdef FE_UNDERFLOW > > + /* > > + * Raise underflow manually because scalbn won't do it if all > > + * lost bits are 0: fma(-0x1p-1000, 0x1.000001p-74, 0x1p-1022) > > + */ > > + feraiseexcept(FE_UNDERFLOW); > > +#endif > > Can you explain why it should happen if all lost bits are zero > (usually that's an exact case). I imagine it's something specific to > fma or its implementation but it's not obvious to me. > this case is for nearest rounding mode when the result is in the subnormal range, at this point the result is represented as hi,lo,scale but the final returned value is computed as scalbn(hi,scale) (the last bits of hi are adjusted if required for correct rounding), however scalbn fails to raise underflow if lo!=0 and all lost bits of hi are 0. the example is such a case: 0x1p-1022 - 0x1.000001p-1074 then hi=1-eps,lo=-0x1p-76,scale=-1022 or maybe with shifted scale and exponents, but in the end only one bit is lost from hi which is zero, alternatively i could do scalbn(lo,scale) too to raise underflow.
Powered by blists - more mailing lists
Powered by Openwall GNU/*/Linux - Powered by OpenVZ