[PATCH] fix underflow exception in fma and fmal

[PATCH] fix underflow exception in fma and fmal

[Szabolcs Nagy]

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
 	}
 	return scalbn(sum.hi, scale);
 }
diff --git a/src/math/fmal.c b/src/math/fmal.c
index 4506aac6..feb19ec1 100644
--- a/src/math/fmal.c
+++ b/src/math/fmal.c
@@ -121,6 +121,14 @@ static inline long double add_and_denormalize(long double a, long double b, int
 		bits_lost = -u.i.se - scale + 1;
 		if ((bits_lost != 1) ^ LASTBIT(u))
 			sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
+#ifdef FE_UNDERFLOW
+		/*
+		 * Raise underflow manually because scalbnl won't do it if all
+		 * lost bits are 0, e.g.:
+		 * fmal(-0x1p-10000L, 0x1.0000000000001p-6445L, 0x1p-16382L)
+		 */
+		feraiseexcept(FE_UNDERFLOW);
+#endif
 	}
 	return scalbnl(sum.hi, scale);
 }
-- 
2.11.0

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

[Rich Felker]

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.

Rich

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

[Szabolcs Nagy]

* 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.

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

[Rich Felker]

On Sun, Mar 19, 2017 at 03:39:53PM +0100, Szabolcs Nagy wrote:
> * 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.

That makes sense. I tend to prefer the scalbn(lo,scale) approach if
there aren't good reasons (performance?) against it, simply because
it's more self-documenting and less special-cased. But whichever you
like is fine. BTW we should probably check that scalbn raises inexact
in all cases it should; I'm not sure what it (especially asm versions)
does in cases where the scale is smaller than the min exponent.

Rich

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

[Szabolcs Nagy]

* Rich Felker <dalias@libc.org> [2017-03-19 10:53:40 -0400]:
> On Sun, Mar 19, 2017 at 03:39:53PM +0100, Szabolcs Nagy wrote:
> > 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.
> 
> That makes sense. I tend to prefer the scalbn(lo,scale) approach if
> there aren't good reasons (performance?) against it, simply because
> it's more self-documenting and less special-cased. But whichever you
> like is fine. BTW we should probably check that scalbn raises inexact
> in all cases it should; I'm not sure what it (especially asm versions)
> does in cases where the scale is smaller than the min exponent.

ok i can do it with scalbn.

generic scalbn is correct and i386 should be correct too:
it does a mul at the end which should raise the flags
(i added tests for this now).

i think we should fix the fma underflow before- vs
after-rounding issue, the current code is only correct
for after-rounding archs (x86, mips, sh) even with the
fix, but doing it correctly is tricky.