This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[PATCH] Fix y0 and y1 exception handling for zero input [BZ #21134]
- From: "Gabriel F. T. Gomes" <gftg at linux dot vnet dot ibm dot com>
- To: libc-alpha at sourceware dot org
- Date: Sun, 12 Feb 2017 23:00:52 -0200
- Subject: [PATCH] Fix y0 and y1 exception handling for zero input [BZ #21134]
- Authentication-results: sourceware.org; auth=none
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.
2017-02-12 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.
---
sysdeps/ieee754/dbl-64/e_j0.c | 5 ++++-
sysdeps/ieee754/dbl-64/e_j1.c | 5 ++++-
sysdeps/ieee754/flt-32/e_j0f.c | 5 ++++-
sysdeps/ieee754/flt-32/e_j1f.c | 4 +++-
sysdeps/ieee754/ldbl-128/e_j0l.c | 1 +
sysdeps/ieee754/ldbl-128/e_j1l.c | 1 +
6 files changed, 17 insertions(+), 4 deletions(-)
diff --git a/sysdeps/ieee754/dbl-64/e_j0.c b/sysdeps/ieee754/dbl-64/e_j0.c
index 9f25aa8..06c7555 100644
--- a/sysdeps/ieee754/dbl-64/e_j0.c
+++ b/sysdeps/ieee754/dbl-64/e_j0.c
@@ -169,7 +169,10 @@ __ieee754_y0 (double x)
if (ix >= 0x7ff00000)
return one / (x + x * x);
if ((ix | lx) == 0)
- return -HUGE_VAL + x; /* -inf and overflow exception. */
+ {
+ __feraiseexcept (FE_DIVBYZERO);
+ return -HUGE_VAL + x; /* -inf and overflow 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 4827fbf..1ff6e86 100644
--- a/sysdeps/ieee754/dbl-64/e_j1.c
+++ b/sysdeps/ieee754/dbl-64/e_j1.c
@@ -174,7 +174,10 @@ __ieee754_y1 (double x)
if (__glibc_unlikely (ix >= 0x7ff00000))
return one / (x + x * x);
if (__glibc_unlikely ((ix | lx) == 0))
- return -HUGE_VAL + x;
+ {
+ __feraiseexcept (FE_DIVBYZERO);
+ return -HUGE_VAL + x;
+ }
/* -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 bd0b80f..1a137fc 100644
--- a/sysdeps/ieee754/flt-32/e_j0f.c
+++ b/sysdeps/ieee754/flt-32/e_j0f.c
@@ -105,7 +105,10 @@ __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) {
+ __feraiseexcept (FE_DIVBYZERO);
+ return -HUGE_VALF+x; /* -inf and overflow 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 f359a3d..10d7edf 100644
--- a/sysdeps/ieee754/flt-32/e_j1f.c
+++ b/sysdeps/ieee754/flt-32/e_j1f.c
@@ -111,8 +111,10 @@ __ieee754_y1f(float x)
ix = 0x7fffffff&hx;
/* 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))
+ if(__builtin_expect(ix==0, 0)) {
+ __feraiseexcept (FE_DIVBYZERO);
return -HUGE_VALF+x; /* -inf and overflow 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 855b5a5..16082c6 100644
--- a/sysdeps/ieee754/ldbl-128/e_j0l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j0l.c
@@ -834,6 +834,7 @@ _Float128
{
if (x < 0)
return (zero / (zero * x));
+ __feraiseexcept (FE_DIVBYZERO);
return -HUGE_VALL + x;
}
xx = fabsl (x);
diff --git a/sysdeps/ieee754/ldbl-128/e_j1l.c b/sysdeps/ieee754/ldbl-128/e_j1l.c
index db8dca0..f2b1036 100644
--- a/sysdeps/ieee754/ldbl-128/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j1l.c
@@ -852,6 +852,7 @@ __ieee754_y1l (_Float128 x)
{
if (x < 0)
return (zero / (zero * x));
+ __feraiseexcept (FE_DIVBYZERO);
return -HUGE_VALL + x;
}
xx = fabsl (x);
--
2.4.11