This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[PATCH v2] [BZ 21130] ldbl-128: Fix y0 and y1 for -Inf input
- From: "Gabriel F. T. Gomes" <gftg at linux dot vnet dot ibm dot com>
- To: libc-alpha at sourceware dot org
- Date: Fri, 10 Feb 2017 17:28:15 -0200
- Subject: [PATCH v2] [BZ 21130] ldbl-128: Fix y0 and y1 for -Inf input
- Authentication-results: sourceware.org; auth=none
The Bessel functions of the second type (Yn) are not defined for
negative input and should return NAN with the "invalid" exception
raised, in these cases. However, current code checks for infinity and
return zero, regardless of the sign. This error is exposed for 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, return the correct value.
Tested for powerpc64le.
2017-01-10 Gabriel F. T. Gomes <gftg@linux.vnet.ibm.com>
[BZ 21130]
* sysdeps/ieee754/ldbl-128/e_j0l.c (__ieee754_y0l): Return NAN
with the "invalid" exception raised when x is -Inf.
* sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_y1l): Likewise.
---
sysdeps/ieee754/ldbl-128/e_j0l.c | 7 +------
sysdeps/ieee754/ldbl-128/e_j1l.c | 7 +------
2 files changed, 2 insertions(+), 12 deletions(-)
diff --git a/sysdeps/ieee754/ldbl-128/e_j0l.c b/sysdeps/ieee754/ldbl-128/e_j0l.c
index d711007..855b5a5 100644
--- a/sysdeps/ieee754/ldbl-128/e_j0l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j0l.c
@@ -829,12 +829,7 @@ _Float128
_Float128 xx, xinv, z, p, q, c, s, cc, ss;
if (! isfinite (x))
- {
- if (x != x)
- return x + x;
- else
- return 0;
- }
+ return 1 / (x + x * x);
if (x <= 0)
{
if (x < 0)
diff --git a/sysdeps/ieee754/ldbl-128/e_j1l.c b/sysdeps/ieee754/ldbl-128/e_j1l.c
index 9e78230..db8dca0 100644
--- a/sysdeps/ieee754/ldbl-128/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j1l.c
@@ -847,12 +847,7 @@ __ieee754_y1l (_Float128 x)
_Float128 xx, xinv, z, p, q, c, s, cc, ss;
if (! isfinite (x))
- {
- if (x != x)
- return x + x;
- else
- return 0;
- }
+ return 1 / (x + x * x);
if (x <= 0)
{
if (x < 0)
--
2.4.11