[PATCH] ldbl-128: Fix y0 and y1 for -Inf input

[Gabriel F. T. Gomes]

The code that this patch changes is not executed, because the wrappers
(in math/w_j0_compat.c and math/w_j1_compat.c) call __kernel_standard
and return.

However, in the float128 branch (not on master), a new wrapper that
does not use __kernel_standard is in use, so this code path gets
executed.  This patch is tested in that branch.

---8<---
The Bessel functions of the second type (Yn) are not defined for
negative x and should return NAN in these cases.  However, current
code checks for infinity return zero, regardless of the sign.  This
error is not exposed for float, double, and long double, because the
wrappers for these functions, which use __kernel_standard
functionality, return the correct value.  For float128 wrappers, which
must not rely on __kernel_standard functionality, this error is
exposed.  This patch fixes it.

2017-01-06  Gabriel F. T. Gomes  <gftg@linux.vnet.ibm.com>

	* sysdeps/ieee754/ldbl-128/e_j0l.c (__ieee754_y0l): Return NAN
	when x is -Inf.
	* sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_y1l): Likewise.
---
 sysdeps/ieee754/ldbl-128/e_j0l.c | 10 +++++++++-
 sysdeps/ieee754/ldbl-128/e_j1l.c | 10 +++++++++-
 2 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/sysdeps/ieee754/ldbl-128/e_j0l.c b/sysdeps/ieee754/ldbl-128/e_j0l.c
index d711007..3bf3569 100644
--- a/sysdeps/ieee754/ldbl-128/e_j0l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j0l.c
@@ -833,7 +833,15 @@ _Float128
       if (x != x)
 	return x + x;
       else
-	return 0;
+	{
+	  /* x = +/-Inf */
+	  if (!signbit (x))
+	    /* y0 (Inf) = 0.  */
+	    return 0;
+	  else
+	    /* y0 (-Inf) = NaN.  */
+	    return NAN;
+	}
     }
   if (x <= 0)
     {
diff --git a/sysdeps/ieee754/ldbl-128/e_j1l.c b/sysdeps/ieee754/ldbl-128/e_j1l.c
index 9e78230..99b6534 100644
--- a/sysdeps/ieee754/ldbl-128/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j1l.c
@@ -851,7 +851,15 @@ __ieee754_y1l (_Float128 x)
       if (x != x)
 	return x + x;
       else
-	return 0;
+	{
+	  /* x = +/-Inf */
+	  if (!signbit (x))
+	    /* y0 (Inf) = 0.  */
+	    return 0;
+	  else
+	    /* y0 (-Inf) = NaN.  */
+	    return NAN;
+	}
     }
   if (x <= 0)
     {
-- 
2.4.11