Use round-to-nearest internally in jn, test with ALL_RM_TEST (bug 18602).

Some existing jn tests, if run in non-default rounding modes, produce
errors above those accepted in glibc, which causes problems for moving
tests of jn to use ALL_RM_TEST.  This patch makes jn set rounding
to-nearest internally, as was done for yn some time ago, then computes
the appropriate underflowing value for results that underflowed to
zero in to-nearest, and moves the tests to ALL_RM_TEST.  It does
nothing about the general inaccuracy of Bessel function
implementations in glibc, though it should make jn more accurate on
average in non-default rounding modes through reduced error
accumulation.  The recomputation of results that underflowed to zero
should as a side-effect fix some cases of bug 16559, where jn just
used an exact zero, but that is *not* the goal of this patch and other
cases of that bug remain unfixed.

(Most of the changes in the patch are reindentation to add new scopes
for SET_RESTORE_ROUND*.)

Tested for x86_64, x86, powerpc and mips64.

	[BZ #16559]
	[BZ #18602]
	* sysdeps/ieee754/dbl-64/e_jn.c (__ieee754_jn): Set
	round-to-nearest internally then recompute results that
	underflowed to zero in the original rounding mode.
	* sysdeps/ieee754/flt-32/e_jnf.c (__ieee754_jnf): Likewise.
	* sysdeps/ieee754/ldbl-128/e_jnl.c (__ieee754_jnl): Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_jnl.c (__ieee754_jnl): Likewise.
	* sysdeps/ieee754/ldbl-96/e_jnl.c (__ieee754_jnl): Likewise
	* math/libm-test.inc (jn_test): Use ALL_RM_TEST.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.

10 files changed, 806 insertions, 714 deletions

diff --git a/ChangeLog b/ChangeLog
index be0be0fa6c..b61ea3cec9 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,18 @@
+2015-06-25  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #16559]
+	[BZ #18602]
+	* sysdeps/ieee754/dbl-64/e_jn.c (__ieee754_jn): Set
+	round-to-nearest internally then recompute results that
+	underflowed to zero in the original rounding mode.
+	* sysdeps/ieee754/flt-32/e_jnf.c (__ieee754_jnf): Likewise.
+	* sysdeps/ieee754/ldbl-128/e_jnl.c (__ieee754_jnl): Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_jnl.c (__ieee754_jnl): Likewise.
+	* sysdeps/ieee754/ldbl-96/e_jnl.c (__ieee754_jnl): Likewise
+	* math/libm-test.inc (jn_test): Use ALL_RM_TEST.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+
 2015-06-25  Andrew Senkevich  <andrew.senkevich@intel.com>
 
 	* NEWS: Fixed description of link with vector math library.
diff --git a/NEWS b/NEWS
index c9be0e44f1..24f8c27138 100644
--- a/NEWS
+++ b/NEWS
@@ -25,7 +25,7 @@ Version 2.22
   18498, 18507, 18512, 18513, 18519, 18520, 18522, 18527, 18528, 18529,
   18530, 18532, 18533, 18534, 18536, 18539, 18540, 18542, 18544, 18545,
   18546, 18547, 18549, 18553, 18558, 18569, 18583, 18585, 18586, 18593,
-  18594.
+  18594, 18602.
 
 * Cache information can be queried via sysconf() function on s390 e.g. with
   _SC_LEVEL1_ICACHE_SIZE as argument.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index da8f8caec4..9e402ab634 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -7486,9 +7486,7 @@ static const struct test_if_f_data jn_test_data[] =
 static void
 jn_test (void)
 {
-  START (jn,, 0);
-  RUN_TEST_LOOP_if_f (jn, jn_test_data, );
-  END;
+  ALL_RM_TEST (jn, 0, jn_test_data, RUN_TEST_LOOP_if_f, END);
 }
 
 
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index c9b565f4fc..5a2af000fc 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -1613,6 +1613,30 @@ ifloat: 3
 ildouble: 4
 ldouble: 4
 
+Function: "jn_downward":
+double: 2
+float: 3
+idouble: 2
+ifloat: 3
+ildouble: 4
+ldouble: 4
+
+Function: "jn_towardzero":
+double: 2
+float: 3
+idouble: 2
+ifloat: 3
+ildouble: 5
+ldouble: 5
+
+Function: "jn_upward":
+double: 2
+float: 3
+idouble: 2
+ifloat: 3
+ildouble: 5
+ldouble: 5
+
 Function: "lgamma":
 double: 1
 float: 1
diff --git a/sysdeps/ieee754/dbl-64/e_jn.c b/sysdeps/ieee754/dbl-64/e_jn.c
index 900737c401..b0ddd5e841 100644
--- a/sysdeps/ieee754/dbl-64/e_jn.c
+++ b/sysdeps/ieee754/dbl-64/e_jn.c
@@ -52,7 +52,7 @@ double
 __ieee754_jn (int n, double x)
 {
   int32_t i, hx, ix, lx, sgn;
-  double a, b, temp, di;
+  double a, b, temp, di, ret;
   double z, w;
 
   /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
@@ -75,14 +75,16 @@ __ieee754_jn (int n, double x)
     return (__ieee754_j1 (x));
   sgn = (n & 1) & (hx >> 31);   /* even n -- 0, odd n -- sign(x) */
   x = fabs (x);
-  if (__glibc_unlikely ((ix | lx) == 0 || ix >= 0x7ff00000))
-    /* if x is 0 or inf */
-    b = zero;
-  else if ((double) n <= x)
-    {
-      /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
-      if (ix >= 0x52D00000)      /* x > 2**302 */
-	{ /* (x >> n**2)
+  {
+    SET_RESTORE_ROUND (FE_TONEAREST);
+    if (__glibc_unlikely ((ix | lx) == 0 || ix >= 0x7ff00000))
+      /* if x is 0 or inf */
+      return sgn == 1 ? -zero : zero;
+    else if ((double) n <= x)
+      {
+	/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
+	if (ix >= 0x52D00000)      /* x > 2**302 */
+	  { /* (x >> n**2)
 			 *	    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 			 *	    Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
 			 *	    Let s=sin(x), c=cos(x),
@@ -95,152 +97,156 @@ __ieee754_jn (int n, double x)
 			 *		   2	-s+c		-c-s
 			 *		   3	 s+c		 c-s
 			 */
-	  double s;
-	  double c;
-	  __sincos (x, &s, &c);
-	  switch (n & 3)
-	    {
-	    case 0: temp = c + s; break;
-	    case 1: temp = -c + s; break;
-	    case 2: temp = -c - s; break;
-	    case 3: temp = c - s; break;
-	    }
-	  b = invsqrtpi * temp / __ieee754_sqrt (x);
-	}
-      else
-	{
-	  a = __ieee754_j0 (x);
-	  b = __ieee754_j1 (x);
-	  for (i = 1; i < n; i++)
-	    {
-	      temp = b;
-	      b = b * ((double) (i + i) / x) - a; /* avoid underflow */
-	      a = temp;
-	    }
-	}
-    }
-  else
-    {
-      if (ix < 0x3e100000)      /* x < 2**-29 */
-	{ /* x is tiny, return the first Taylor expansion of J(n,x)
+	    double s;
+	    double c;
+	    __sincos (x, &s, &c);
+	    switch (n & 3)
+	      {
+	      case 0: temp = c + s; break;
+	      case 1: temp = -c + s; break;
+	      case 2: temp = -c - s; break;
+	      case 3: temp = c - s; break;
+	      }
+	    b = invsqrtpi * temp / __ieee754_sqrt (x);
+	  }
+	else
+	  {
+	    a = __ieee754_j0 (x);
+	    b = __ieee754_j1 (x);
+	    for (i = 1; i < n; i++)
+	      {
+		temp = b;
+		b = b * ((double) (i + i) / x) - a; /* avoid underflow */
+		a = temp;
+	      }
+	  }
+      }
+    else
+      {
+	if (ix < 0x3e100000)      /* x < 2**-29 */
+	  { /* x is tiny, return the first Taylor expansion of J(n,x)
 			 * J(n,x) = 1/n!*(x/2)^n  - ...
 			 */
-	  if (n > 33)           /* underflow */
-	    b = zero;
-	  else
-	    {
-	      temp = x * 0.5; b = temp;
-	      for (a = one, i = 2; i <= n; i++)
-		{
-		  a *= (double) i;              /* a = n! */
-		  b *= temp;                    /* b = (x/2)^n */
-		}
-	      b = b / a;
-	    }
-	}
-      else
-	{
-	  /* use backward recurrence */
-	  /*			x      x^2      x^2
-	   *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
-	   *			2n  - 2(n+1) - 2(n+2)
-	   *
-	   *			1      1        1
-	   *  (for large x)   =  ----  ------   ------   .....
-	   *			2n   2(n+1)   2(n+2)
-	   *			-- - ------ - ------ -
-	   *			 x     x         x
-	   *
-	   * Let w = 2n/x and h=2/x, then the above quotient
-	   * is equal to the continued fraction:
-	   *		    1
-	   *	= -----------------------
-	   *		       1
-	   *	   w - -----------------
-	   *			  1
-	   *		w+h - ---------
-	   *		       w+2h - ...
-	   *
-	   * To determine how many terms needed, let
-	   * Q(0) = w, Q(1) = w(w+h) - 1,
-	   * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
-	   * When Q(k) > 1e4	good for single
-	   * When Q(k) > 1e9	good for double
-	   * When Q(k) > 1e17	good for quadruple
-	   */
-	  /* determine k */
-	  double t, v;
-	  double q0, q1, h, tmp; int32_t k, m;
-	  w = (n + n) / (double) x; h = 2.0 / (double) x;
-	  q0 = w;  z = w + h; q1 = w * z - 1.0; k = 1;
-	  while (q1 < 1.0e9)
-	    {
-	      k += 1; z += h;
-	      tmp = z * q1 - q0;
-	      q0 = q1;
-	      q1 = tmp;
-	    }
-	  m = n + n;
-	  for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
-	    t = one / (i / x - t);
-	  a = t;
-	  b = one;
-	  /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
-	   *  Hence, if n*(log(2n/x)) > ...
-	   *  single 8.8722839355e+01
-	   *  double 7.09782712893383973096e+02
-	   *  long double 1.1356523406294143949491931077970765006170e+04
-	   *  then recurrent value may overflow and the result is
-	   *  likely underflow to zero
-	   */
-	  tmp = n;
-	  v = two / x;
-	  tmp = tmp * __ieee754_log (fabs (v * tmp));
-	  if (tmp < 7.09782712893383973096e+02)
-	    {
-	      for (i = n - 1, di = (double) (i + i); i > 0; i--)
-		{
-		  temp = b;
-		  b *= di;
-		  b = b / x - a;
-		  a = temp;
-		  di -= two;
-		}
-	    }
-	  else
-	    {
-	      for (i = n - 1, di = (double) (i + i); i > 0; i--)
-		{
-		  temp = b;
-		  b *= di;
-		  b = b / x - a;
-		  a = temp;
-		  di -= two;
-		  /* scale b to avoid spurious overflow */
-		  if (b > 1e100)
-		    {
-		      a /= b;
-		      t /= b;
-		      b = one;
-		    }
-		}
-	    }
-	  /* j0() and j1() suffer enormous loss of precision at and
-	   * near zero; however, we know that their zero points never
-	   * coincide, so just choose the one further away from zero.
-	   */
-	  z = __ieee754_j0 (x);
-	  w = __ieee754_j1 (x);
-	  if (fabs (z) >= fabs (w))
-	    b = (t * z / b);
-	  else
-	    b = (t * w / a);
-	}
-    }
-  if (sgn == 1)
-    return -b;
-  else
-    return b;
+	    if (n > 33)           /* underflow */
+	      b = zero;
+	    else
+	      {
+		temp = x * 0.5; b = temp;
+		for (a = one, i = 2; i <= n; i++)
+		  {
+		    a *= (double) i;              /* a = n! */
+		    b *= temp;                    /* b = (x/2)^n */
+		  }
+		b = b / a;
+	      }
+	  }
+	else
+	  {
+	    /* use backward recurrence */
+	    /*			x      x^2      x^2
+	     *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
+	     *			2n  - 2(n+1) - 2(n+2)
+	     *
+	     *			1      1        1
+	     *  (for large x)   =  ----  ------   ------   .....
+	     *			2n   2(n+1)   2(n+2)
+	     *			-- - ------ - ------ -
+	     *			 x     x         x
+	     *
+	     * Let w = 2n/x and h=2/x, then the above quotient
+	     * is equal to the continued fraction:
+	     *		    1
+	     *	= -----------------------
+	     *		       1
+	     *	   w - -----------------
+	     *			  1
+	     *		w+h - ---------
+	     *		       w+2h - ...
+	     *
+	     * To determine how many terms needed, let
+	     * Q(0) = w, Q(1) = w(w+h) - 1,
+	     * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
+	     * When Q(k) > 1e4	good for single
+	     * When Q(k) > 1e9	good for double
+	     * When Q(k) > 1e17	good for quadruple
+	     */
+	    /* determine k */
+	    double t, v;
+	    double q0, q1, h, tmp; int32_t k, m;
+	    w = (n + n) / (double) x; h = 2.0 / (double) x;
+	    q0 = w;  z = w + h; q1 = w * z - 1.0; k = 1;
+	    while (q1 < 1.0e9)
+	      {
+		k += 1; z += h;
+		tmp = z * q1 - q0;
+		q0 = q1;
+		q1 = tmp;
+	      }
+	    m = n + n;
+	    for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
+	      t = one / (i / x - t);
+	    a = t;
+	    b = one;
+	    /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
+	     *  Hence, if n*(log(2n/x)) > ...
+	     *  single 8.8722839355e+01
+	     *  double 7.09782712893383973096e+02
+	     *  long double 1.1356523406294143949491931077970765006170e+04
+	     *  then recurrent value may overflow and the result is
+	     *  likely underflow to zero
+	     */
+	    tmp = n;
+	    v = two / x;
+	    tmp = tmp * __ieee754_log (fabs (v * tmp));
+	    if (tmp < 7.09782712893383973096e+02)
+	      {
+		for (i = n - 1, di = (double) (i + i); i > 0; i--)
+		  {
+		    temp = b;
+		    b *= di;
+		    b = b / x - a;
+		    a = temp;
+		    di -= two;
+		  }
+	      }
+	    else
+	      {
+		for (i = n - 1, di = (double) (i + i); i > 0; i--)
+		  {
+		    temp = b;
+		    b *= di;
+		    b = b / x - a;
+		    a = temp;
+		    di -= two;
+		    /* scale b to avoid spurious overflow */
+		    if (b > 1e100)
+		      {
+			a /= b;
+			t /= b;
+			b = one;
+		      }
+		  }
+	      }
+	    /* j0() and j1() suffer enormous loss of precision at and
+	     * near zero; however, we know that their zero points never
+	     * coincide, so just choose the one further away from zero.
+	     */
+	    z = __ieee754_j0 (x);
+	    w = __ieee754_j1 (x);
+	    if (fabs (z) >= fabs (w))
+	      b = (t * z / b);
+	    else
+	      b = (t * w / a);
+	  }
+      }
+    if (sgn == 1)
+      ret = -b;
+    else
+      ret = b;
+  }
+  if (ret == 0)
+    ret = __copysign (DBL_MIN, ret) * DBL_MIN;
+  return ret;
 }
 strong_alias (__ieee754_jn, __jn_finite)
 
diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c
index dc4b371bc1..ec5a81b653 100644
--- a/sysdeps/ieee754/flt-32/e_jnf.c
+++ b/sysdeps/ieee754/flt-32/e_jnf.c
@@ -27,6 +27,8 @@ static const float zero  =  0.0000000000e+00;
 float
 __ieee754_jnf(int n, float x)
 {
+    float ret;
+    {
 	int32_t i,hx,ix, sgn;
 	float a, b, temp, di;
 	float z, w;
@@ -47,8 +49,9 @@ __ieee754_jnf(int n, float x)
 	if(n==1) return(__ieee754_j1f(x));
 	sgn = (n&1)&(hx>>31);	/* even n -- 0, odd n -- sign(x) */
 	x = fabsf(x);
+	SET_RESTORE_ROUNDF (FE_TONEAREST);
 	if(__builtin_expect(ix==0||ix>=0x7f800000, 0))	/* if x is 0 or inf */
-	    b = zero;
+	    return sgn == 1 ? -zero : zero;
 	else if((float)n<=x) {
 		/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
 	    a = __ieee754_j0f(x);
@@ -163,7 +166,11 @@ __ieee754_jnf(int n, float x)
 		  b = (t * w / a);
 	    }
 	}
-	if(sgn==1) return -b; else return b;
+	if(sgn==1) ret = -b; else ret = b;
+    }
+    if (ret == 0)
+	ret = __copysignf (FLT_MIN, ret) * FLT_MIN;
+    return ret;
 }
 strong_alias (__ieee754_jnf, __jnf_finite)
 
diff --git a/sysdeps/ieee754/ldbl-128/e_jnl.c b/sysdeps/ieee754/ldbl-128/e_jnl.c
index 422623f0dc..14d65ff081 100644
--- a/sysdeps/ieee754/ldbl-128/e_jnl.c
+++ b/sysdeps/ieee754/ldbl-128/e_jnl.c
@@ -73,7 +73,7 @@ __ieee754_jnl (int n, long double x)
 {
   u_int32_t se;
   int32_t i, ix, sgn;
-  long double a, b, temp, di;
+  long double a, b, temp, di, ret;
   long double z, w;
   ieee854_long_double_shape_type u;
 
@@ -106,192 +106,198 @@ __ieee754_jnl (int n, long double x)
   sgn = (n & 1) & (se >> 31);	/* even n -- 0, odd n -- sign(x) */
   x = fabsl (x);
 
-  if (x == 0.0L || ix >= 0x7fff0000)	/* if x is 0 or inf */
-    b = zero;
-  else if ((long double) n <= x)
-    {
-      /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
-      if (ix >= 0x412D0000)
-	{			/* x > 2**302 */
+  {
+    SET_RESTORE_ROUNDL (FE_TONEAREST);
+    if (x == 0.0L || ix >= 0x7fff0000)	/* if x is 0 or inf */
+      return sgn == 1 ? -zero : zero;
+    else if ((long double) n <= x)
+      {
+	/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
+	if (ix >= 0x412D0000)
+	  {			/* x > 2**302 */
 
-	  /* ??? Could use an expansion for large x here.  */
+	    /* ??? Could use an expansion for large x here.  */
 
-	  /* (x >> n**2)
-	   *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-	   *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-	   *      Let s=sin(x), c=cos(x),
-	   *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
-	   *
-	   *             n    sin(xn)*sqt2    cos(xn)*sqt2
-	   *          ----------------------------------
-	   *             0     s-c             c+s
-	   *             1    -s-c            -c+s
-	   *             2    -s+c            -c-s
-	   *             3     s+c             c-s
-	   */
-	  long double s;
-	  long double c;
-	  __sincosl (x, &s, &c);
-	  switch (n & 3)
-	    {
-	    case 0:
-	      temp = c + s;
-	      break;
-	    case 1:
-	      temp = -c + s;
-	      break;
-	    case 2:
-	      temp = -c - s;
-	      break;
-	    case 3:
-	      temp = c - s;
-	      break;
-	    }
-	  b = invsqrtpi * temp / __ieee754_sqrtl (x);
-	}
-      else
-	{
-	  a = __ieee754_j0l (x);
-	  b = __ieee754_j1l (x);
-	  for (i = 1; i < n; i++)
-	    {
-	      temp = b;
-	      b = b * ((long double) (i + i) / x) - a;	/* avoid underflow */
-	      a = temp;
-	    }
-	}
-    }
-  else
-    {
-      if (ix < 0x3fc60000)
-	{			/* x < 2**-57 */
-	  /* x is tiny, return the first Taylor expansion of J(n,x)
-	   * J(n,x) = 1/n!*(x/2)^n  - ...
-	   */
-	  if (n >= 400)		/* underflow, result < 10^-4952 */
-	    b = zero;
-	  else
-	    {
-	      temp = x * 0.5;
-	      b = temp;
-	      for (a = one, i = 2; i <= n; i++)
-		{
-		  a *= (long double) i;	/* a = n! */
-		  b *= temp;	/* b = (x/2)^n */
-		}
-	      b = b / a;
-	    }
-	}
-      else
-	{
-	  /* use backward recurrence */
-	  /*                      x      x^2      x^2
-	   *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
-	   *                      2n  - 2(n+1) - 2(n+2)
-	   *
-	   *                      1      1        1
-	   *  (for large x)   =  ----  ------   ------   .....
-	   *                      2n   2(n+1)   2(n+2)
-	   *                      -- - ------ - ------ -
-	   *                       x     x         x
-	   *
-	   * Let w = 2n/x and h=2/x, then the above quotient
-	   * is equal to the continued fraction:
-	   *                  1
-	   *      = -----------------------
-	   *                     1
-	   *         w - -----------------
-	   *                        1
-	   *              w+h - ---------
-	   *                     w+2h - ...
-	   *
-	   * To determine how many terms needed, let
-	   * Q(0) = w, Q(1) = w(w+h) - 1,
-	   * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
-	   * When Q(k) > 1e4      good for single
-	   * When Q(k) > 1e9      good for double
-	   * When Q(k) > 1e17     good for quadruple
-	   */
-	  /* determine k */
-	  long double t, v;
-	  long double q0, q1, h, tmp;
-	  int32_t k, m;
-	  w = (n + n) / (long double) x;
-	  h = 2.0L / (long double) x;
-	  q0 = w;
-	  z = w + h;
-	  q1 = w * z - 1.0L;
-	  k = 1;
-	  while (q1 < 1.0e17L)
-	    {
-	      k += 1;
-	      z += h;
-	      tmp = z * q1 - q0;
-	      q0 = q1;
-	      q1 = tmp;
-	    }
-	  m = n + n;
-	  for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
-	    t = one / (i / x - t);
-	  a = t;
-	  b = one;
-	  /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
-	   *  Hence, if n*(log(2n/x)) > ...
-	   *  single 8.8722839355e+01
-	   *  double 7.09782712893383973096e+02
-	   *  long double 1.1356523406294143949491931077970765006170e+04
-	   *  then recurrent value may overflow and the result is
-	   *  likely underflow to zero
-	   */
-	  tmp = n;
-	  v = two / x;
-	  tmp = tmp * __ieee754_logl (fabsl (v * tmp));
+	    /* (x >> n**2)
+	     *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+	     *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+	     *      Let s=sin(x), c=cos(x),
+	     *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
+	     *
+	     *             n    sin(xn)*sqt2    cos(xn)*sqt2
+	     *          ----------------------------------
+	     *             0     s-c             c+s
+	     *             1    -s-c            -c+s
+	     *             2    -s+c            -c-s
+	     *             3     s+c             c-s
+	     */
+	    long double s;
+	    long double c;
+	    __sincosl (x, &s, &c);
+	    switch (n & 3)
+	      {
+	      case 0:
+		temp = c + s;
+		break;
+	      case 1:
+		temp = -c + s;
+		break;
+	      case 2:
+		temp = -c - s;
+		break;
+	      case 3:
+		temp = c - s;
+		break;
+	      }
+	    b = invsqrtpi * temp / __ieee754_sqrtl (x);
+	  }
+	else
+	  {
+	    a = __ieee754_j0l (x);
+	    b = __ieee754_j1l (x);
+	    for (i = 1; i < n; i++)
+	      {
+		temp = b;
+		b = b * ((long double) (i + i) / x) - a;	/* avoid underflow */
+		a = temp;
+	      }
+	  }
+      }
+    else
+      {
+	if (ix < 0x3fc60000)
+	  {			/* x < 2**-57 */
+	    /* x is tiny, return the first Taylor expansion of J(n,x)
+	     * J(n,x) = 1/n!*(x/2)^n  - ...
+	     */
+	    if (n >= 400)		/* underflow, result < 10^-4952 */
+	      b = zero;
+	    else
+	      {
+		temp = x * 0.5;
+		b = temp;
+		for (a = one, i = 2; i <= n; i++)
+		  {
+		    a *= (long double) i;	/* a = n! */
+		    b *= temp;	/* b = (x/2)^n */
+		  }
+		b = b / a;
+	      }
+	  }
+	else
+	  {
+	    /* use backward recurrence */
+	    /*                      x      x^2      x^2
+	     *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
+	     *                      2n  - 2(n+1) - 2(n+2)
+	     *
+	     *                      1      1        1
+	     *  (for large x)   =  ----  ------   ------   .....
+	     *                      2n   2(n+1)   2(n+2)
+	     *                      -- - ------ - ------ -
+	     *                       x     x         x
+	     *
+	     * Let w = 2n/x and h=2/x, then the above quotient
+	     * is equal to the continued fraction:
+	     *                  1
+	     *      = -----------------------
+	     *                     1
+	     *         w - -----------------
+	     *                        1
+	     *              w+h - ---------
+	     *                     w+2h - ...
+	     *
+	     * To determine how many terms needed, let
+	     * Q(0) = w, Q(1) = w(w+h) - 1,
+	     * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
+	     * When Q(k) > 1e4      good for single
+	     * When Q(k) > 1e9      good for double
+	     * When Q(k) > 1e17     good for quadruple
+	     */
+	    /* determine k */
+	    long double t, v;
+	    long double q0, q1, h, tmp;
+	    int32_t k, m;
+	    w = (n + n) / (long double) x;
+	    h = 2.0L / (long double) x;
+	    q0 = w;
+	    z = w + h;
+	    q1 = w * z - 1.0L;
+	    k = 1;
+	    while (q1 < 1.0e17L)
+	      {
+		k += 1;
+		z += h;
+		tmp = z * q1 - q0;
+		q0 = q1;
+		q1 = tmp;
+	      }
+	    m = n + n;
+	    for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
+	      t = one / (i / x - t);
+	    a = t;
+	    b = one;
+	    /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
+	     *  Hence, if n*(log(2n/x)) > ...
+	     *  single 8.8722839355e+01
+	     *  double 7.09782712893383973096e+02
+	     *  long double 1.1356523406294143949491931077970765006170e+04
+	     *  then recurrent value may overflow and the result is
+	     *  likely underflow to zero
+	     */
+	    tmp = n;
+	    v = two / x;
+	    tmp = tmp * __ieee754_logl (fabsl (v * tmp));
 
-	  if (tmp < 1.1356523406294143949491931077970765006170e+04L)
-	    {
-	      for (i = n - 1, di = (long double) (i + i); i > 0; i--)
-		{
-		  temp = b;
-		  b *= di;
-		  b = b / x - a;
-		  a = temp;
-		  di -= two;
-		}
-	    }
-	  else
-	    {
-	      for (i = n - 1, di = (long double) (i + i); i > 0; i--)
-		{
-		  temp = b;
-		  b *= di;
-		  b = b / x - a;
-		  a = temp;
-		  di -= two;
-		  /* scale b to avoid spurious overflow */
-		  if (b > 1e100L)
-		    {
-		      a /= b;
-		      t /= b;
-		      b = one;
-		    }
-		}
-	    }
-	  /* j0() and j1() suffer enormous loss of precision at and
-	   * near zero; however, we know that their zero points never
-	   * coincide, so just choose the one further away from zero.
-	   */
-	  z = __ieee754_j0l (x);
-	  w = __ieee754_j1l (x);
-	  if (fabsl (z) >= fabsl (w))
-	    b = (t * z / b);
-	  else
-	    b = (t * w / a);
-	}
-    }
-  if (sgn == 1)
-    return -b;
-  else
-    return b;
+	    if (tmp < 1.1356523406294143949491931077970765006170e+04L)
+	      {
+		for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+		  {
+		    temp = b;
+		    b *= di;
+		    b = b / x - a;
+		    a = temp;
+		    di -= two;
+		  }
+	      }
+	    else
+	      {
+		for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+		  {
+		    temp = b;
+		    b *= di;
+		    b = b / x - a;
+		    a = temp;
+		    di -= two;
+		    /* scale b to avoid spurious overflow */
+		    if (b > 1e100L)
+		      {
+			a /= b;
+			t /= b;
+			b = one;
+		      }
+		  }
+	      }
+	    /* j0() and j1() suffer enormous loss of precision at and
+	     * near zero; however, we know that their zero points never
+	     * coincide, so just choose the one further away from zero.
+	     */
+	    z = __ieee754_j0l (x);
+	    w = __ieee754_j1l (x);
+	    if (fabsl (z) >= fabsl (w))
+	      b = (t * z / b);
+	    else
+	      b = (t * w / a);
+	  }
+      }
+    if (sgn == 1)
+      ret = -b;
+    else
+      ret = b;
+  }
+  if (ret == 0)
+    ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN;
+  return ret;
 }
 strong_alias (__ieee754_jnl, __jnl_finite)
 
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
index d2b9318327..5d0a2b5b6a 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
@@ -73,7 +73,7 @@ __ieee754_jnl (int n, long double x)
 {
   uint32_t se, lx;
   int32_t i, ix, sgn;
-  long double a, b, temp, di;
+  long double a, b, temp, di, ret;
   long double z, w;
   double xhi;
 
@@ -106,192 +106,198 @@ __ieee754_jnl (int n, long double x)
   sgn = (n & 1) & (se >> 31);	/* even n -- 0, odd n -- sign(x) */
   x = fabsl (x);
 
-  if (x == 0.0L || ix >= 0x7ff00000)	/* if x is 0 or inf */
-    b = zero;
-  else if ((long double) n <= x)
-    {
-      /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
-      if (ix >= 0x52d00000)
-	{			/* x > 2**302 */
+  {
+    SET_RESTORE_ROUNDL (FE_TONEAREST);
+    if (x == 0.0L || ix >= 0x7ff00000)	/* if x is 0 or inf */
+      return sgn == 1 ? -zero : zero;
+    else if ((long double) n <= x)
+      {
+	/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
+	if (ix >= 0x52d00000)
+	  {			/* x > 2**302 */
 
-	  /* ??? Could use an expansion for large x here.  */
+	    /* ??? Could use an expansion for large x here.  */
 
-	  /* (x >> n**2)
-	   *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-	   *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-	   *      Let s=sin(x), c=cos(x),
-	   *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
-	   *
-	   *             n    sin(xn)*sqt2    cos(xn)*sqt2
-	   *          ----------------------------------
-	   *             0     s-c             c+s
-	   *             1    -s-c            -c+s
-	   *             2    -s+c            -c-s
-	   *             3     s+c             c-s
-	   */
-	  long double s;
-	  long double c;
-	  __sincosl (x, &s, &c);
-	  switch (n & 3)
-	    {
-	    case 0:
-	      temp = c + s;
-	      break;
-	    case 1:
-	      temp = -c + s;
-	      break;
-	    case 2:
-	      temp = -c - s;
-	      break;
-	    case 3:
-	      temp = c - s;
-	      break;
-	    }
-	  b = invsqrtpi * temp / __ieee754_sqrtl (x);
-	}
-      else
-	{
-	  a = __ieee754_j0l (x);
-	  b = __ieee754_j1l (x);
-	  for (i = 1; i < n; i++)
-	    {
-	      temp = b;
-	      b = b * ((long double) (i + i) / x) - a;	/* avoid underflow */
-	      a = temp;
-	    }
-	}
-    }
-  else
-    {
-      if (ix < 0x3e100000)
-	{			/* x < 2**-29 */
-	  /* x is tiny, return the first Taylor expansion of J(n,x)
-	   * J(n,x) = 1/n!*(x/2)^n  - ...
-	   */
-	  if (n >= 33)		/* underflow, result < 10^-300 */
-	    b = zero;
-	  else
-	    {
-	      temp = x * 0.5;
-	      b = temp;
-	      for (a = one, i = 2; i <= n; i++)
-		{
-		  a *= (long double) i;	/* a = n! */
-		  b *= temp;	/* b = (x/2)^n */
-		}
-	      b = b / a;
-	    }
-	}
-      else
-	{
-	  /* use backward recurrence */
-	  /*                      x      x^2      x^2
-	   *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
-	   *                      2n  - 2(n+1) - 2(n+2)
-	   *
-	   *                      1      1        1
-	   *  (for large x)   =  ----  ------   ------   .....
-	   *                      2n   2(n+1)   2(n+2)
-	   *                      -- - ------ - ------ -
-	   *                       x     x         x
-	   *
-	   * Let w = 2n/x and h=2/x, then the above quotient
-	   * is equal to the continued fraction:
-	   *                  1
-	   *      = -----------------------
-	   *                     1
-	   *         w - -----------------
-	   *                        1
-	   *              w+h - ---------
-	   *                     w+2h - ...
-	   *
-	   * To determine how many terms needed, let
-	   * Q(0) = w, Q(1) = w(w+h) - 1,
-	   * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
-	   * When Q(k) > 1e4      good for single
-	   * When Q(k) > 1e9      good for double
-	   * When Q(k) > 1e17     good for quadruple
-	   */
-	  /* determine k */
-	  long double t, v;
-	  long double q0, q1, h, tmp;
-	  int32_t k, m;
-	  w = (n + n) / (long double) x;
-	  h = 2.0L / (long double) x;
-	  q0 = w;
-	  z = w + h;
-	  q1 = w * z - 1.0L;
-	  k = 1;
-	  while (q1 < 1.0e17L)
-	    {
-	      k += 1;
-	      z += h;
-	      tmp = z * q1 - q0;
-	      q0 = q1;
-	      q1 = tmp;
-	    }
-	  m = n + n;
-	  for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
-	    t = one / (i / x - t);
-	  a = t;
-	  b = one;
-	  /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
-	   *  Hence, if n*(log(2n/x)) > ...
-	   *  single 8.8722839355e+01
-	   *  double 7.09782712893383973096e+02
-	   *  long double 1.1356523406294143949491931077970765006170e+04
-	   *  then recurrent value may overflow and the result is
-	   *  likely underflow to zero
-	   */
-	  tmp = n;
-	  v = two / x;
-	  tmp = tmp * __ieee754_logl (fabsl (v * tmp));
+	    /* (x >> n**2)
+	     *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+	     *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+	     *      Let s=sin(x), c=cos(x),
+	     *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
+	     *
+	     *             n    sin(xn)*sqt2    cos(xn)*sqt2
+	     *          ----------------------------------
+	     *             0     s-c             c+s
+	     *             1    -s-c            -c+s
+	     *             2    -s+c            -c-s
+	     *             3     s+c             c-s
+	     */
+	    long double s;
+	    long double c;
+	    __sincosl (x, &s, &c);
+	    switch (n & 3)
+	      {
+	      case 0:
+		temp = c + s;
+		break;
+	      case 1:
+		temp = -c + s;
+		break;
+	      case 2:
+		temp = -c - s;
+		break;
+	      case 3:
+		temp = c - s;
+		break;
+	      }
+	    b = invsqrtpi * temp / __ieee754_sqrtl (x);
+	  }
+	else
+	  {
+	    a = __ieee754_j0l (x);
+	    b = __ieee754_j1l (x);
+	    for (i = 1; i < n; i++)
+	      {
+		temp = b;
+		b = b * ((long double) (i + i) / x) - a;	/* avoid underflow */
+		a = temp;
+	      }
+	  }
+      }
+    else
+      {
+	if (ix < 0x3e100000)
+	  {			/* x < 2**-29 */
+	    /* x is tiny, return the first Taylor expansion of J(n,x)
+	     * J(n,x) = 1/n!*(x/2)^n  - ...
+	     */
+	    if (n >= 33)		/* underflow, result < 10^-300 */
+	      b = zero;
+	    else
+	      {
+		temp = x * 0.5;
+		b = temp;
+		for (a = one, i = 2; i <= n; i++)
+		  {
+		    a *= (long double) i;	/* a = n! */
+		    b *= temp;	/* b = (x/2)^n */
+		  }
+		b = b / a;
+	      }
+	  }
+	else
+	  {
+	    /* use backward recurrence */
+	    /*                      x      x^2      x^2
+	     *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
+	     *                      2n  - 2(n+1) - 2(n+2)
+	     *
+	     *                      1      1        1
+	     *  (for large x)   =  ----  ------   ------   .....
+	     *                      2n   2(n+1)   2(n+2)
+	     *                      -- - ------ - ------ -
+	     *                       x     x         x
+	     *
+	     * Let w = 2n/x and h=2/x, then the above quotient
+	     * is equal to the continued fraction:
+	     *                  1
+	     *      = -----------------------
+	     *                     1
+	     *         w - -----------------
+	     *                        1
+	     *              w+h - ---------
+	     *                     w+2h - ...
+	     *
+	     * To determine how many terms needed, let
+	     * Q(0) = w, Q(1) = w(w+h) - 1,
+	     * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
+	     * When Q(k) > 1e4      good for single
+	     * When Q(k) > 1e9      good for double
+	     * When Q(k) > 1e17     good for quadruple
+	     */
+	    /* determine k */
+	    long double t, v;
+	    long double q0, q1, h, tmp;
+	    int32_t k, m;
+	    w = (n + n) / (long double) x;
+	    h = 2.0L / (long double) x;
+	    q0 = w;
+	    z = w + h;
+	    q1 = w * z - 1.0L;
+	    k = 1;
+	    while (q1 < 1.0e17L)
+	      {
+		k += 1;
+		z += h;
+		tmp = z * q1 - q0;
+		q0 = q1;
+		q1 = tmp;
+	      }
+	    m = n + n;
+	    for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
+	      t = one / (i / x - t);
+	    a = t;
+	    b = one;
+	    /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
+	     *  Hence, if n*(log(2n/x)) > ...
+	     *  single 8.8722839355e+01
+	     *  double 7.09782712893383973096e+02
+	     *  long double 1.1356523406294143949491931077970765006170e+04
+	     *  then recurrent value may overflow and the result is
+	     *  likely underflow to zero
+	     */
+	    tmp = n;
+	    v = two / x;
+	    tmp = tmp * __ieee754_logl (fabsl (v * tmp));
 
-	  if (tmp < 1.1356523406294143949491931077970765006170e+04L)
-	    {
-	      for (i = n - 1, di = (long double) (i + i); i > 0; i--)
-		{
-		  temp = b;
-		  b *= di;
-		  b = b / x - a;
-		  a = temp;
-		  di -= two;
-		}
-	    }
-	  else
-	    {
-	      for (i = n - 1, di = (long double) (i + i); i > 0; i--)
-		{
-		  temp = b;
-		  b *= di;
-		  b = b / x - a;
-		  a = temp;
-		  di -= two;
-		  /* scale b to avoid spurious overflow */
-		  if (b > 1e100L)
-		    {
-		      a /= b;
-		      t /= b;
-		      b = one;
-		    }
-		}
-	    }
-	  /* j0() and j1() suffer enormous loss of precision at and
-	   * near zero; however, we know that their zero points never
-	   * coincide, so just choose the one further away from zero.
-	   */
-	  z = __ieee754_j0l (x);
-	  w = __ieee754_j1l (x);
-	  if (fabsl (z) >= fabsl (w))
-	    b = (t * z / b);
-	  else
-	    b = (t * w / a);
-	}
-    }
-  if (sgn == 1)
-    return -b;
-  else
-    return b;
+	    if (tmp < 1.1356523406294143949491931077970765006170e+04L)
+	      {
+		for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+		  {
+		    temp = b;
+		    b *= di;
+		    b = b / x - a;
+		    a = temp;
+		    di -= two;
+		  }
+	      }
+	    else
+	      {
+		for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+		  {
+		    temp = b;
+		    b *= di;
+		    b = b / x - a;
+		    a = temp;
+		    di -= two;
+		    /* scale b to avoid spurious overflow */
+		    if (b > 1e100L)
+		      {
+			a /= b;
+			t /= b;
+			b = one;
+		      }
+		  }
+	      }
+	    /* j0() and j1() suffer enormous loss of precision at and
+	     * near zero; however, we know that their zero points never
+	     * coincide, so just choose the one further away from zero.
+	     */
+	    z = __ieee754_j0l (x);
+	    w = __ieee754_j1l (x);
+	    if (fabsl (z) >= fabsl (w))
+	      b = (t * z / b);
+	    else
+	      b = (t * w / a);
+	  }
+      }
+    if (sgn == 1)
+      ret = -b;
+    else
+      ret = b;
+  }
+  if (ret == 0)
+    ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN;
+  return ret;
 }
 strong_alias (__ieee754_jnl, __jnl_finite)
 
diff --git a/sysdeps/ieee754/ldbl-96/e_jnl.c b/sysdeps/ieee754/ldbl-96/e_jnl.c
index a6668089dd..49c9c421b0 100644
--- a/sysdeps/ieee754/ldbl-96/e_jnl.c
+++ b/sysdeps/ieee754/ldbl-96/e_jnl.c
@@ -71,7 +71,7 @@ __ieee754_jnl (int n, long double x)
 {
   u_int32_t se, i0, i1;
   int32_t i, ix, sgn;
-  long double a, b, temp, di;
+  long double a, b, temp, di, ret;
   long double z, w;
 
   /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
@@ -96,195 +96,201 @@ __ieee754_jnl (int n, long double x)
     return (__ieee754_j1l (x));
   sgn = (n & 1) & (se >> 15);	/* even n -- 0, odd n -- sign(x) */
   x = fabsl (x);
-  if (__glibc_unlikely ((ix | i0 | i1) == 0 || ix >= 0x7fff))
-    /* if x is 0 or inf */
-    b = zero;
-  else if ((long double) n <= x)
-    {
-      /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
-      if (ix >= 0x412D)
-	{			/* x > 2**302 */
+  {
+    SET_RESTORE_ROUNDL (FE_TONEAREST);
+    if (__glibc_unlikely ((ix | i0 | i1) == 0 || ix >= 0x7fff))
+      /* if x is 0 or inf */
+      return sgn == 1 ? -zero : zero;
+    else if ((long double) n <= x)
+      {
+	/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
+	if (ix >= 0x412D)
+	  {			/* x > 2**302 */
 
-	  /* ??? This might be a futile gesture.
-	     If x exceeds X_TLOSS anyway, the wrapper function
-	     will set the result to zero. */
+	    /* ??? This might be a futile gesture.
+	       If x exceeds X_TLOSS anyway, the wrapper function
+	       will set the result to zero. */
 
-	  /* (x >> n**2)
-	   *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-	   *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
-	   *      Let s=sin(x), c=cos(x),
-	   *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
-	   *
-	   *             n    sin(xn)*sqt2    cos(xn)*sqt2
-	   *          ----------------------------------
-	   *             0     s-c             c+s
-	   *             1    -s-c            -c+s
-	   *             2    -s+c            -c-s
-	   *             3     s+c             c-s
-	   */
-	  long double s;
-	  long double c;
-	  __sincosl (x, &s, &c);
-	  switch (n & 3)
-	    {
-	    case 0:
-	      temp = c + s;
-	      break;
-	    case 1:
-	      temp = -c + s;
-	      break;
-	    case 2:
-	      temp = -c - s;
-	      break;
-	    case 3:
-	      temp = c - s;
-	      break;
-	    }
-	  b = invsqrtpi * temp / __ieee754_sqrtl (x);
-	}
-      else
-	{
-	  a = __ieee754_j0l (x);
-	  b = __ieee754_j1l (x);
-	  for (i = 1; i < n; i++)
-	    {
-	      temp = b;
-	      b = b * ((long double) (i + i) / x) - a;	/* avoid underflow */
-	      a = temp;
-	    }
-	}
-    }
-  else
-    {
-      if (ix < 0x3fde)
-	{			/* x < 2**-33 */
-	  /* x is tiny, return the first Taylor expansion of J(n,x)
-	   * J(n,x) = 1/n!*(x/2)^n  - ...
-	   */
-	  if (n >= 400)		/* underflow, result < 10^-4952 */
-	    b = zero;
-	  else
-	    {
-	      temp = x * 0.5;
-	      b = temp;
-	      for (a = one, i = 2; i <= n; i++)
-		{
-		  a *= (long double) i;	/* a = n! */
-		  b *= temp;	/* b = (x/2)^n */
-		}
-	      b = b / a;
-	    }
-	}
-      else
-	{
-	  /* use backward recurrence */
-	  /*                      x      x^2      x^2
-	   *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
-	   *                      2n  - 2(n+1) - 2(n+2)
-	   *
-	   *                      1      1        1
-	   *  (for large x)   =  ----  ------   ------   .....
-	   *                      2n   2(n+1)   2(n+2)
-	   *                      -- - ------ - ------ -
-	   *                       x     x         x
-	   *
-	   * Let w = 2n/x and h=2/x, then the above quotient
-	   * is equal to the continued fraction:
-	   *                  1
-	   *      = -----------------------
-	   *                     1
-	   *         w - -----------------
-	   *                        1
-	   *              w+h - ---------
-	   *                     w+2h - ...
-	   *
-	   * To determine how many terms needed, let
-	   * Q(0) = w, Q(1) = w(w+h) - 1,
-	   * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
-	   * When Q(k) > 1e4      good for single
-	   * When Q(k) > 1e9      good for double
-	   * When Q(k) > 1e17     good for quadruple
-	   */
-	  /* determine k */
-	  long double t, v;
-	  long double q0, q1, h, tmp;
-	  int32_t k, m;
-	  w = (n + n) / (long double) x;
-	  h = 2.0L / (long double) x;
-	  q0 = w;
-	  z = w + h;
-	  q1 = w * z - 1.0L;
-	  k = 1;
-	  while (q1 < 1.0e11L)
-	    {
-	      k += 1;
-	      z += h;
-	      tmp = z * q1 - q0;
-	      q0 = q1;
-	      q1 = tmp;
-	    }
-	  m = n + n;
-	  for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
-	    t = one / (i / x - t);
-	  a = t;
-	  b = one;
-	  /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
-	   *  Hence, if n*(log(2n/x)) > ...
-	   *  single 8.8722839355e+01
-	   *  double 7.09782712893383973096e+02
-	   *  long double 1.1356523406294143949491931077970765006170e+04
-	   *  then recurrent value may overflow and the result is
-	   *  likely underflow to zero
-	   */
-	  tmp = n;
-	  v = two / x;
-	  tmp = tmp * __ieee754_logl (fabsl (v * tmp));
+	    /* (x >> n**2)
+	     *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+	     *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+	     *      Let s=sin(x), c=cos(x),
+	     *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
+	     *
+	     *             n    sin(xn)*sqt2    cos(xn)*sqt2
+	     *          ----------------------------------
+	     *             0     s-c             c+s
+	     *             1    -s-c            -c+s
+	     *             2    -s+c            -c-s
+	     *             3     s+c             c-s
+	     */
+	    long double s;
+	    long double c;
+	    __sincosl (x, &s, &c);
+	    switch (n & 3)
+	      {
+	      case 0:
+		temp = c + s;
+		break;
+	      case 1:
+		temp = -c + s;
+		break;
+	      case 2:
+		temp = -c - s;
+		break;
+	      case 3:
+		temp = c - s;
+		break;
+	      }
+	    b = invsqrtpi * temp / __ieee754_sqrtl (x);
+	  }
+	else
+	  {
+	    a = __ieee754_j0l (x);
+	    b = __ieee754_j1l (x);
+	    for (i = 1; i < n; i++)
+	      {
+		temp = b;
+		b = b * ((long double) (i + i) / x) - a;	/* avoid underflow */
+		a = temp;
+	      }
+	  }
+      }
+    else
+      {
+	if (ix < 0x3fde)
+	  {			/* x < 2**-33 */
+	    /* x is tiny, return the first Taylor expansion of J(n,x)
+	     * J(n,x) = 1/n!*(x/2)^n  - ...
+	     */
+	    if (n >= 400)		/* underflow, result < 10^-4952 */
+	      b = zero;
+	    else
+	      {
+		temp = x * 0.5;
+		b = temp;
+		for (a = one, i = 2; i <= n; i++)
+		  {
+		    a *= (long double) i;	/* a = n! */
+		    b *= temp;	/* b = (x/2)^n */
+		  }
+		b = b / a;
+	      }
+	  }
+	else
+	  {
+	    /* use backward recurrence */
+	    /*                      x      x^2      x^2
+	     *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
+	     *                      2n  - 2(n+1) - 2(n+2)
+	     *
+	     *                      1      1        1
+	     *  (for large x)   =  ----  ------   ------   .....
+	     *                      2n   2(n+1)   2(n+2)
+	     *                      -- - ------ - ------ -
+	     *                       x     x         x
+	     *
+	     * Let w = 2n/x and h=2/x, then the above quotient
+	     * is equal to the continued fraction:
+	     *                  1
+	     *      = -----------------------
+	     *                     1
+	     *         w - -----------------
+	     *                        1
+	     *              w+h - ---------
+	     *                     w+2h - ...
+	     *
+	     * To determine how many terms needed, let
+	     * Q(0) = w, Q(1) = w(w+h) - 1,
+	     * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
+	     * When Q(k) > 1e4      good for single
+	     * When Q(k) > 1e9      good for double
+	     * When Q(k) > 1e17     good for quadruple
+	     */
+	    /* determine k */
+	    long double t, v;
+	    long double q0, q1, h, tmp;
+	    int32_t k, m;
+	    w = (n + n) / (long double) x;
+	    h = 2.0L / (long double) x;
+	    q0 = w;
+	    z = w + h;
+	    q1 = w * z - 1.0L;
+	    k = 1;
+	    while (q1 < 1.0e11L)
+	      {
+		k += 1;
+		z += h;
+		tmp = z * q1 - q0;
+		q0 = q1;
+		q1 = tmp;
+	      }
+	    m = n + n;
+	    for (t = zero, i = 2 * (n + k); i >= m; i -= 2)
+	      t = one / (i / x - t);
+	    a = t;
+	    b = one;
+	    /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
+	     *  Hence, if n*(log(2n/x)) > ...
+	     *  single 8.8722839355e+01
+	     *  double 7.09782712893383973096e+02
+	     *  long double 1.1356523406294143949491931077970765006170e+04
+	     *  then recurrent value may overflow and the result is
+	     *  likely underflow to zero
+	     */
+	    tmp = n;
+	    v = two / x;
+	    tmp = tmp * __ieee754_logl (fabsl (v * tmp));
 
-	  if (tmp < 1.1356523406294143949491931077970765006170e+04L)
-	    {
-	      for (i = n - 1, di = (long double) (i + i); i > 0; i--)
-		{
-		  temp = b;
-		  b *= di;
-		  b = b / x - a;
-		  a = temp;
-		  di -= two;
-		}
-	    }
-	  else
-	    {
-	      for (i = n - 1, di = (long double) (i + i); i > 0; i--)
-		{
-		  temp = b;
-		  b *= di;
-		  b = b / x - a;
-		  a = temp;
-		  di -= two;
-		  /* scale b to avoid spurious overflow */
-		  if (b > 1e100L)
-		    {
-		      a /= b;
-		      t /= b;
-		      b = one;
-		    }
-		}
-	    }
-	  /* j0() and j1() suffer enormous loss of precision at and
-	   * near zero; however, we know that their zero points never
-	   * coincide, so just choose the one further away from zero.
-	   */
-	  z = __ieee754_j0l (x);
-	  w = __ieee754_j1l (x);
-	  if (fabsl (z) >= fabsl (w))
-	    b = (t * z / b);
-	  else
-	    b = (t * w / a);
-	}
-    }
-  if (sgn == 1)
-    return -b;
-  else
-    return b;
+	    if (tmp < 1.1356523406294143949491931077970765006170e+04L)
+	      {
+		for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+		  {
+		    temp = b;
+		    b *= di;
+		    b = b / x - a;
+		    a = temp;
+		    di -= two;
+		  }
+	      }
+	    else
+	      {
+		for (i = n - 1, di = (long double) (i + i); i > 0; i--)
+		  {
+		    temp = b;
+		    b *= di;
+		    b = b / x - a;
+		    a = temp;
+		    di -= two;
+		    /* scale b to avoid spurious overflow */
+		    if (b > 1e100L)
+		      {
+			a /= b;
+			t /= b;
+			b = one;
+		      }
+		  }
+	      }
+	    /* j0() and j1() suffer enormous loss of precision at and
+	     * near zero; however, we know that their zero points never
+	     * coincide, so just choose the one further away from zero.
+	     */
+	    z = __ieee754_j0l (x);
+	    w = __ieee754_j1l (x);
+	    if (fabsl (z) >= fabsl (w))
+	      b = (t * z / b);
+	    else
+	      b = (t * w / a);
+	  }
+      }
+    if (sgn == 1)
+      ret = -b;
+    else
+      ret = b;
+  }
+  if (ret == 0)
+    ret = __copysignl (LDBL_MIN, ret) * LDBL_MIN;
+  return ret;
 }
 strong_alias (__ieee754_jnl, __jnl_finite)
 
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 48d11a6f4d..12d0c5a7db 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1767,6 +1767,30 @@ ifloat: 4
 ildouble: 4
 ldouble: 4
 
+Function: "jn_downward":
+double: 5
+float: 5
+idouble: 5
+ifloat: 5
+ildouble: 4
+ldouble: 4
+
+Function: "jn_towardzero":
+double: 5
+float: 5
+idouble: 5
+ifloat: 5
+ildouble: 5
+ldouble: 5
+
+Function: "jn_upward":
+double: 5
+float: 5
+idouble: 5
+ifloat: 5
+ildouble: 5
+ldouble: 5
+
 Function: "lgamma":
 double: 2
 float: 2