Fix y0 and y1 exception handling for zero input [BZ #21134]

The Bessel functions of the second type (Yn) should raise the "divide
by zero" exception when input is zero (both positive and negative).
Current code gives the right output, but fails to set the exception.
This error is exposed for float, double, and long double when linking
with -lieee.  Without this flag, the error is not exposed, because the
wrappers for these functions, which use __kernel_standard
functionality, set the exception as expected.

Tested for powerpc64le.

	[BZ #21134]
	* sysdeps/ieee754/dbl-64/e_j0.c (__ieee754_y0): Raise the
	"divide by zero" exception when the input is zero.
	* sysdeps/ieee754/dbl-64/e_j1.c (__ieee754_y1): Likewise.
	* sysdeps/ieee754/flt-32/e_j0f.c (__ieee754_y0f): Likewise.
	* sysdeps/ieee754/flt-32/e_j1f.c (__ieee754_y1f): Likewise.
	* sysdeps/ieee754/ldbl-128/e_j0l.c (__ieee754_y0l): Likewise.
	* sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_y1l): Likewise.

diff --git a/ChangeLog b/ChangeLog
index b040928a3339eedaba2d37c4a9b9466de513deba..62961a2940a50166713ab115ca6d7ac8f77e01df 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,14 @@
+2017-02-15  Gabriel F. T. Gomes  <gftg@linux.vnet.ibm.com>
+
+       [BZ #21134]
+       * sysdeps/ieee754/dbl-64/e_j0.c (__ieee754_y0): Raise the
+       "divide by zero" exception when the input is zero.
+       * sysdeps/ieee754/dbl-64/e_j1.c (__ieee754_y1): Likewise.
+       * sysdeps/ieee754/flt-32/e_j0f.c (__ieee754_y0f): Likewise.
+       * sysdeps/ieee754/flt-32/e_j1f.c (__ieee754_y1f): Likewise.
+       * sysdeps/ieee754/ldbl-128/e_j0l.c (__ieee754_y0l): Likewise.
+       * sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_y1l): Likewise.
+
 2017-02-15  Joseph Myers  <joseph@codesourcery.com>
 
        * sysdeps/x86_64/fpu/test-double-vlen2.c: Move most contents to,

diff --git a/sysdeps/ieee754/dbl-64/e_j0.c b/sysdeps/ieee754/dbl-64/e_j0.c
index 9f25aa855e8683ed43f58db31e9606dfcff8b56b..4b440cf0d045b94778ba889b9309c46d2c5e7a4a 100644 (file)
--- a/sysdeps/ieee754/dbl-64/e_j0.c
+++ b/sysdeps/ieee754/dbl-64/e_j0.c
@@ -169,7 +169,7 @@ __ieee754_y0 (double x)
   if (ix >= 0x7ff00000)
     return one / (x + x * x);
   if ((ix | lx) == 0)
-    return -HUGE_VAL + x;                  /* -inf and overflow exception.  */
+    return -1 / zero; /* -inf and divide by zero exception.  */
   if (hx < 0)
     return zero / (zero * x);
   if (ix >= 0x40000000)         /* |x| >= 2.0 */

diff --git a/sysdeps/ieee754/dbl-64/e_j1.c b/sysdeps/ieee754/dbl-64/e_j1.c
index 4827fbf3d39bec50b396edfb8a5d152f4e3da25f..eb446fd102d91f9ab98f9ed979569e7ede7e876e 100644 (file)
--- a/sysdeps/ieee754/dbl-64/e_j1.c
+++ b/sysdeps/ieee754/dbl-64/e_j1.c
@@ -174,7 +174,7 @@ __ieee754_y1 (double x)
   if (__glibc_unlikely (ix >= 0x7ff00000))
     return one / (x + x * x);
   if (__glibc_unlikely ((ix | lx) == 0))
-    return -HUGE_VAL + x;
+    return -1 / zero; /* -inf and divide by zero exception.  */
   /* -inf and overflow exception.  */;
   if (__glibc_unlikely (hx < 0))
     return zero / (zero * x);

diff --git a/sysdeps/ieee754/flt-32/e_j0f.c b/sysdeps/ieee754/flt-32/e_j0f.c
index bd0b80fdb0b338812dae56a40b280400702be484..b783dd069dded429fa41ce8c4a93df6544797ab6 100644 (file)
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -105,7 +105,7 @@ __ieee754_y0f(float x)
        ix = 0x7fffffff&hx;
     /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf.  */
        if(ix>=0x7f800000) return  one/(x+x*x);
-       if(ix==0) return -HUGE_VALF+x;  /* -inf and overflow exception.  */
+       if(ix==0) return -1/zero; /* -inf and divide by zero exception.  */
        if(hx<0) return zero/(zero*x);
        if(ix >= 0x40000000) {  /* |x| >= 2.0 */
        /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))

diff --git a/sysdeps/ieee754/flt-32/e_j1f.c b/sysdeps/ieee754/flt-32/e_j1f.c
index f359a3d9ba6bf55074c000936d2190250d52ae76..805a87d85b732b9d30d3b79464e039df2b8e832b 100644 (file)
--- a/sysdeps/ieee754/flt-32/e_j1f.c
+++ b/sysdeps/ieee754/flt-32/e_j1f.c
@@ -112,7 +112,7 @@ __ieee754_y1f(float x)
     /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
        if(__builtin_expect(ix>=0x7f800000, 0)) return  one/(x+x*x);
        if(__builtin_expect(ix==0, 0))
-               return -HUGE_VALF+x;  /* -inf and overflow exception.  */
+               return -1/zero; /* -inf and divide by zero exception.  */
        if(__builtin_expect(hx<0, 0)) return zero/(zero*x);
        if(ix >= 0x40000000) {  /* |x| >= 2.0 */
                SET_RESTORE_ROUNDF (FE_TONEAREST);

diff --git a/sysdeps/ieee754/ldbl-128/e_j0l.c b/sysdeps/ieee754/ldbl-128/e_j0l.c
index 855b5a578ba7b88aa7abd0a8914c591ebaa5de09..fb8d3518cea2617f59d08599ac03e8ed9dcc32b6 100644 (file)
--- a/sysdeps/ieee754/ldbl-128/e_j0l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j0l.c
@@ -834,7 +834,7 @@ _Float128
     {
       if (x < 0)
        return (zero / (zero * x));
-      return -HUGE_VALL + x;
+      return -1 / zero; /* -inf and divide by zero exception.  */
     }
   xx = fabsl (x);
   if (xx <= 0x1p-57)

diff --git a/sysdeps/ieee754/ldbl-128/e_j1l.c b/sysdeps/ieee754/ldbl-128/e_j1l.c
index db8dca0ab13e6747903b3d039e1e926e85fbe71c..6fc69faa3c5fb5d563504817d307981e4ebc757c 100644 (file)
--- a/sysdeps/ieee754/ldbl-128/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j1l.c
@@ -852,7 +852,7 @@ __ieee754_y1l (_Float128 x)
     {
       if (x < 0)
        return (zero / (zero * x));
-      return -HUGE_VALL + x;
+      return -1 / zero; /* -inf and divide by zero exception.  */
     }
   xx = fabsl (x);
   if (xx <= 0x1p-114)