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.
This commit is contained in:
Gabriel F. T. Gomes 2017-02-12 22:36:27 -02:00
parent 10303eb74b
commit 4918e5f4cd
7 changed files with 17 additions and 6 deletions

View File

@ -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> 2017-02-15 Joseph Myers <joseph@codesourcery.com>
* sysdeps/x86_64/fpu/test-double-vlen2.c: Move most contents to, * sysdeps/x86_64/fpu/test-double-vlen2.c: Move most contents to,

View File

@ -169,7 +169,7 @@ __ieee754_y0 (double x)
if (ix >= 0x7ff00000) if (ix >= 0x7ff00000)
return one / (x + x * x); return one / (x + x * x);
if ((ix | lx) == 0) if ((ix | lx) == 0)
return -HUGE_VAL + x; /* -inf and overflow exception. */ return -1 / zero; /* -inf and divide by zero exception. */
if (hx < 0) if (hx < 0)
return zero / (zero * x); return zero / (zero * x);
if (ix >= 0x40000000) /* |x| >= 2.0 */ if (ix >= 0x40000000) /* |x| >= 2.0 */

View File

@ -174,7 +174,7 @@ __ieee754_y1 (double x)
if (__glibc_unlikely (ix >= 0x7ff00000)) if (__glibc_unlikely (ix >= 0x7ff00000))
return one / (x + x * x); return one / (x + x * x);
if (__glibc_unlikely ((ix | lx) == 0)) if (__glibc_unlikely ((ix | lx) == 0))
return -HUGE_VAL + x; return -1 / zero; /* -inf and divide by zero exception. */
/* -inf and overflow exception. */; /* -inf and overflow exception. */;
if (__glibc_unlikely (hx < 0)) if (__glibc_unlikely (hx < 0))
return zero / (zero * x); return zero / (zero * x);

View File

@ -105,7 +105,7 @@ __ieee754_y0f(float x)
ix = 0x7fffffff&hx; ix = 0x7fffffff&hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */ /* 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>=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(hx<0) return zero/(zero*x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */ if(ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))

View File

@ -112,7 +112,7 @@ __ieee754_y1f(float x)
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ /* 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>=0x7f800000, 0)) return one/(x+x*x);
if(__builtin_expect(ix==0, 0)) 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(__builtin_expect(hx<0, 0)) return zero/(zero*x);
if(ix >= 0x40000000) { /* |x| >= 2.0 */ if(ix >= 0x40000000) { /* |x| >= 2.0 */
SET_RESTORE_ROUNDF (FE_TONEAREST); SET_RESTORE_ROUNDF (FE_TONEAREST);

View File

@ -834,7 +834,7 @@ _Float128
{ {
if (x < 0) if (x < 0)
return (zero / (zero * x)); return (zero / (zero * x));
return -HUGE_VALL + x; return -1 / zero; /* -inf and divide by zero exception. */
} }
xx = fabsl (x); xx = fabsl (x);
if (xx <= 0x1p-57) if (xx <= 0x1p-57)

View File

@ -852,7 +852,7 @@ __ieee754_y1l (_Float128 x)
{ {
if (x < 0) if (x < 0)
return (zero / (zero * x)); return (zero / (zero * x));
return -HUGE_VALL + x; return -1 / zero; /* -inf and divide by zero exception. */
} }
xx = fabsl (x); xx = fabsl (x);
if (xx <= 0x1p-114) if (xx <= 0x1p-114)