From 4918e5f4cd63290bb0b2c614f52092ca6a779126 Mon Sep 17 00:00:00 2001 From: "Gabriel F. T. Gomes" Date: Sun, 12 Feb 2017 22:36:27 -0200 Subject: 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. --- sysdeps/ieee754/dbl-64/e_j0.c | 2 +- sysdeps/ieee754/dbl-64/e_j1.c | 2 +- sysdeps/ieee754/flt-32/e_j0f.c | 2 +- sysdeps/ieee754/flt-32/e_j1f.c | 2 +- sysdeps/ieee754/ldbl-128/e_j0l.c | 2 +- sysdeps/ieee754/ldbl-128/e_j1l.c | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) (limited to 'sysdeps/ieee754') diff --git a/sysdeps/ieee754/dbl-64/e_j0.c b/sysdeps/ieee754/dbl-64/e_j0.c index 9f25aa8..4b440cf 100644 --- 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 4827fbf..eb446fd1 100644 --- 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 bd0b80f..b783dd0 100644 --- 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 f359a3d..805a87d 100644 --- 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 855b5a5..fb8d351 100644 --- 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 db8dca0..6fc69fa 100644 --- 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) -- cgit v1.1