Fix ldbl-128ibm fma (Inf, Inf, finite) (bug 23272)

Message ID 20180613220136.11438-1-tuliom@linux.ibm.com
State New
Headers show
Series
  • Fix ldbl-128ibm fma (Inf, Inf, finite) (bug 23272)
Related show

Commit Message

Tulio Magno Quites Machado Filho June 13, 2018, 10:01 p.m.
This solution is mainly based on Joseph's fix at commit
ca121b117f2c9c97a4c121334481a96c94fef3a0.

It has extra differences because:
 - part of the non-finite arguments were already being treated;
 - when x and y are +-Inf and z if finite, an overflow can be
   generated.

Tested on powerpc[|64|64le].

2018-06-13  Tulio Magno Quites Machado Filho  <tuliom@linux.ibm.com>

	[BZ 23272]
	* sysdeps/ieee754/ldbl-128ibm/s_fmal.c (__fmal): Handle all
	cases of non-finite arguments.

Signed-off-by: Tulio Magno Quites Machado Filho <tuliom@linux.ibm.com>
---
 sysdeps/ieee754/ldbl-128ibm/s_fmal.c | 31 ++++++++++++-------------------
 1 file changed, 12 insertions(+), 19 deletions(-)

Comments

Joseph Myers June 13, 2018, 11:07 p.m. | #1
On Wed, 13 Jun 2018, Tulio Magno Quites Machado Filho wrote:

> This solution is mainly based on Joseph's fix at commit
> ca121b117f2c9c97a4c121334481a96c94fef3a0.
> 
> It has extra differences because:
>  - part of the non-finite arguments were already being treated;
>  - when x and y are +-Inf and z if finite, an overflow can be
>    generated.

I'm not clear what the observed issues you are fixing are.  Is this patch 
to avoid running into libgcc generating spurious overflows for the case 
you mention, or is it intended to fix wrong results?

> +    if (!isfinite (x) && !isfinite (y) && isfinite(z))
> +      /* Compute the result as x * y to avoid an overflow.  */
> +      return x * y;

Why wouldn't this be needed if just one of x and y is not finite?  If x * 
y + z generates a spurious overflow in some cases where x * y is infinite, 
I'd expect that to apply equally to the case where the infinity comes from 
(Inf * finite), not just from (Inf * Inf).
Florian Weimer Sept. 21, 2018, 11:59 a.m. | #2
* Joseph Myers:

> On Wed, 13 Jun 2018, Tulio Magno Quites Machado Filho wrote:
>
>> This solution is mainly based on Joseph's fix at commit
>> ca121b117f2c9c97a4c121334481a96c94fef3a0.
>> 
>> It has extra differences because:
>>  - part of the non-finite arguments were already being treated;
>>  - when x and y are +-Inf and z if finite, an overflow can be
>>    generated.
>
> I'm not clear what the observed issues you are fixing are.  Is this patch 
> to avoid running into libgcc generating spurious overflows for the case 
> you mention, or is it intended to fix wrong results?
>
>> +    if (!isfinite (x) && !isfinite (y) && isfinite(z))
>> +      /* Compute the result as x * y to avoid an overflow.  */
>> +      return x * y;
>
> Why wouldn't this be needed if just one of x and y is not finite?  If x * 
> y + z generates a spurious overflow in some cases where x * y is infinite, 
> I'd expect that to apply equally to the case where the infinity comes from 
> (Inf * finite), not just from (Inf * Inf).

We still see a failures (see bug 23272).  How can we move this forward?

Thanks,
Florian
Joseph Myers Sept. 21, 2018, 12:46 p.m. | #3
On Fri, 21 Sep 2018, Florian Weimer wrote:

> * Joseph Myers:
> 
> > On Wed, 13 Jun 2018, Tulio Magno Quites Machado Filho wrote:
> >
> >> This solution is mainly based on Joseph's fix at commit
> >> ca121b117f2c9c97a4c121334481a96c94fef3a0.
> >> 
> >> It has extra differences because:
> >>  - part of the non-finite arguments were already being treated;
> >>  - when x and y are +-Inf and z if finite, an overflow can be
> >>    generated.
> >
> > I'm not clear what the observed issues you are fixing are.  Is this patch 
> > to avoid running into libgcc generating spurious overflows for the case 
> > you mention, or is it intended to fix wrong results?
> >
> >> +    if (!isfinite (x) && !isfinite (y) && isfinite(z))
> >> +      /* Compute the result as x * y to avoid an overflow.  */
> >> +      return x * y;
> >
> > Why wouldn't this be needed if just one of x and y is not finite?  If x * 
> > y + z generates a spurious overflow in some cases where x * y is infinite, 
> > I'd expect that to apply equally to the case where the infinity comes from 
> > (Inf * finite), not just from (Inf * Inf).
> 
> We still see a failures (see bug 23272).  How can we move this forward?

Do you mean bug 23584?  That answers the first question about the symptoms 
(both spurious overflow exceptions and bad results, apparently coming from 
__gcc_qadd mishandling certain cases involving infinities), but doesn't 
answer the second question about what happens when the infinite x * y 
comes from (Inf * finite) (I'd guess the same problem, in which case the 
proposed patch needs fixing further, but it looks like tests for that case 
also need adding to libm-test-fma.inc).
Florian Weimer Sept. 21, 2018, 12:51 p.m. | #4
* Joseph Myers:

> On Fri, 21 Sep 2018, Florian Weimer wrote:
>
>> * Joseph Myers:
>> 
>> > On Wed, 13 Jun 2018, Tulio Magno Quites Machado Filho wrote:
>> >
>> >> This solution is mainly based on Joseph's fix at commit
>> >> ca121b117f2c9c97a4c121334481a96c94fef3a0.
>> >> 
>> >> It has extra differences because:
>> >>  - part of the non-finite arguments were already being treated;
>> >>  - when x and y are +-Inf and z if finite, an overflow can be
>> >>    generated.
>> >
>> > I'm not clear what the observed issues you are fixing are.  Is this patch 
>> > to avoid running into libgcc generating spurious overflows for the case 
>> > you mention, or is it intended to fix wrong results?
>> >
>> >> +    if (!isfinite (x) && !isfinite (y) && isfinite(z))
>> >> +      /* Compute the result as x * y to avoid an overflow.  */
>> >> +      return x * y;
>> >
>> > Why wouldn't this be needed if just one of x and y is not finite?  If x * 
>> > y + z generates a spurious overflow in some cases where x * y is infinite, 
>> > I'd expect that to apply equally to the case where the infinity comes from 
>> > (Inf * finite), not just from (Inf * Inf).
>> 
>> We still see a failures (see bug 23272).  How can we move this forward?
>
> Do you mean bug 23584?

Yes, sorry I used the wrong bug number.

> That answers the first question about the symptoms (both spurious
> overflow exceptions and bad results, apparently coming from __gcc_qadd
> mishandling certain cases involving infinities), but doesn't answer
> the second question about what happens when the infinite x * y comes
> from (Inf * finite) (I'd guess the same problem, in which case the
> proposed patch needs fixing further, but it looks like tests for that
> case also need adding to libm-test-fma.inc).

Thanks.  I have to admit this is a bit over my head.

Tulio, can you take care of this, or should someone else pick this up?

Thanks,
Florian
Tulio Magno Quites Machado Filho Sept. 21, 2018, 5:27 p.m. | #5
Florian Weimer <fweimer@redhat.com> writes:

> * Joseph Myers:
>> That answers the first question about the symptoms (both spurious
>> overflow exceptions and bad results, apparently coming from __gcc_qadd
>> mishandling certain cases involving infinities), but doesn't answer
>> the second question about what happens when the infinite x * y comes
>> from (Inf * finite) (I'd guess the same problem, in which case the
>> proposed patch needs fixing further, but it looks like tests for that
>> case also need adding to libm-test-fma.inc).
>
> Thanks.  I have to admit this is a bit over my head.
>
> Tulio, can you take care of this, or should someone else pick this up?

Yes, it's in my task list.
I still have to confirm Joseph's points and will get back to you.
I still need more time, though.

Patch

diff --git a/sysdeps/ieee754/ldbl-128ibm/s_fmal.c b/sysdeps/ieee754/ldbl-128ibm/s_fmal.c
index e72a3e4d59..afba1ca727 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_fmal.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_fmal.c
@@ -127,37 +127,30 @@  long double
 __fmal (long double x, long double y, long double z)
 {
   double xhi, xlo, yhi, ylo, zhi, zlo;
-  int64_t hx, hy, hz;
-  int xexp, yexp, zexp;
   double scale_val;
   int scale_exp;
   ldbl_unpack (x, &xhi, &xlo);
-  EXTRACT_WORDS64 (hx, xhi);
-  xexp = (hx & 0x7ff0000000000000LL) >> 52;
   ldbl_unpack (y, &yhi, &ylo);
-  EXTRACT_WORDS64 (hy, yhi);
-  yexp = (hy & 0x7ff0000000000000LL) >> 52;
   ldbl_unpack (z, &zhi, &zlo);
-  EXTRACT_WORDS64 (hz, zhi);
-  zexp = (hz & 0x7ff0000000000000LL) >> 52;
 
-  /* If z is Inf or NaN, but x and y are finite, avoid any exceptions
-     from computing x * y.  */
-  if (zexp == 0x7ff && xexp != 0x7ff && yexp != 0x7ff)
+  if (__glibc_unlikely (!isfinite (x) || !isfinite (y) || x == 0 || y == 0))
+    if (!isfinite (x) && !isfinite (y) && isfinite(z))
+      /* Compute the result as x * y to avoid an overflow.  */
+      return x * y;
+    else
+      /* If x or y is Inf, NaN or zero compute as x * y + z.  */
+      return (x * y) + z;
+  else if (__glibc_unlikely (!isfinite (z)))
+    /* If z is Inf, but x and y are finite, the result should be z
+       rather than NaN.  */
     return (z + x) + y;
 
-  /* If z is zero and x are y are nonzero, compute the result as x * y
+  /* If z is zero and x and y are nonzero, compute the result as x * y
      to avoid the wrong sign of a zero result if x * y underflows to
      0.  */
-  if (z == 0 && x != 0 && y != 0)
+  if (z == 0)
     return x * y;
 
-  /* If x or y or z is Inf/NaN, or if x * y is zero, compute as x * y
-     + z.  */
-  if (xexp == 0x7ff || yexp == 0x7ff || zexp == 0x7ff
-      || x == 0 || y == 0)
-    return (x * y) + z;
-
   {
     SET_RESTORE_ROUND (FE_TONEAREST);