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

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);