Re: [PATCH] fix underflow exception in fma and fmal

[Szabolcs Nagy <nsz@port70.net>]

* Rich Felker <dalias@libc.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.