Fix jn precision

ChangeLog

@@ -1,3 +1,23 @@
+2011-09-09  Ulrich Drepper  <[email protected]>
+
+	* sysdeps/i386/fpu/libm-test-ulps: Adjust ULPs for jn tests.
+
+2011-07-03  Andreas Jaeger  <[email protected]>
+
+	* math/libm-test.inc (jn_test): Add tests for BZ#11589.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Add new ULPs for jn_test,
+	regenerate with gen-libm-tests.pl.
+
+2010-05-12  Petr Baudis  <[email protected]>
+
+	[BZ #11589]
+	* sysdeps/ieee754/dbl-64/e_jn.c: Compensate major precision loss
+	around j0() zero points by switching to j1().
+	* sysdeps/ieee754/flt-32/e_jnf.c: Likewise.
+	* sysdeps/ieee754/ldbl-128/e_jnl.c: Likewise.
+	* sysdeps/ieee754/ldbl-128ibm/e_jnl.c: Likewise.
+	* sysdeps/ieee754/ldbl-96/e_jnl.c: Likewise.
+
 2011-09-09  Ulrich Drepper  <[email protected]>
 
 	* sysdeps/unix/bsd/bsd4.4/bits/socket.h (__cmsg_nxthdr): Use NULL

NEWS

@@ -9,8 +9,9 @@ Version 2.15
 
 * The following bugs are resolved with this release:
 
-  9696, 12403, 12847, 12868, 12852, 12874, 12885, 12907, 12922, 12935,
-  13007, 13021, 13068, 13092, 13114, 13118, 13123, 13134, 13138, 13150
+  9696, 11589, 12403, 12847, 12868, 12852, 12874, 12885, 12907, 12922,
+  12935, 13007, 13021, 13068, 13092, 13114, 13118, 13123, 13134, 13138,
+  13150
 
 * New program pldd to list loaded object of a process
   Implemented by Ulrich Drepper.

math/libm-test.inc

@@ -1,4 +1,4 @@
-/* Copyright (C) 1997-2006, 2007, 2009, 2010 Free Software Foundation, Inc.
+/* Copyright (C) 1997-2006, 2007, 2009, 2010, 2011 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Andreas Jaeger <[email protected]>, 1997.
 
@@ -3285,6 +3285,16 @@ jn_test (void)
   TEST_ff_f (jn, 10, 2.0, 0.251538628271673670963516093751820639e-6L);
   TEST_ff_f (jn, 10, 10.0, 0.207486106633358857697278723518753428L);
 
+  /* BZ #11589 .*/
+  TEST_ff_f (jn, 2, 2.4048255576957729, 0.43175480701968038399746111312430703L);
+  TEST_ff_f (jn, 3, 2.4048255576957729, 0.19899990535769083404042146764530813L);
+  TEST_ff_f (jn, 4, 2.4048255576957729, 0.647466661641779720084932282551219891E-1L);
+  TEST_ff_f (jn, 5, 2.4048255576957729, 0.163892432048058525099230549946147698E-1L);
+  TEST_ff_f (jn, 6, 2.4048255576957729, 0.34048184720278336646673682895929161E-2L);
+  TEST_ff_f (jn, 7, 2.4048255576957729, 0.60068836573295394221291569249883076E-3L);
+  TEST_ff_f (jn, 8, 2.4048255576957729, 0.92165786705344923232879022467054148E-4L);
+  TEST_ff_f (jn, 9, 2.4048255576957729, 0.12517270977961513005428966643852564E-4L)
+
   END (jn);
 }
 

sysdeps/i386/fpu/libm-test-ulps

@@ -640,6 +640,52 @@ double: 1
 idouble: 1
 ildouble: 1
 ldouble: 1
+Test "jn (2, 2.4048255576957729) == 0.43175480701968038399746111312430703":
+float: 1
+ifloat: 1
+double: 1
+idouble: 1
+ldouble: 82
+ildouble: 82
+Test "jn (3, 2.4048255576957729) == 0.19899990535769083404042146764530813":
+ldouble: 186
+ildouble: 186
+Test "jn (4, 2.4048255576957729) == 0.647466661641779720084932282551219891E-1":
+ldouble: 185
+ildouble: 185
+Test: "jn (5, 2.4048255576957729) == 0.163892432048058525099230549946147698E-1":
+float: 1
+ifloat: 1
+double: 1
+idouble: 1
+ldouble: 249
+ildouble: 249
+Test "jn (6, 2.4048255576957729) == 0.34048184720278336646673682895929161E-2":
+float: 2
+ifloat: 2
+double: 1
+idouble: 1
+ldouble: 511
+ildouble: 511
+Test "jn (7, 2.4048255576957729) == 0.60068836573295394221291569249883076E-3":
+float: 2
+ifloat: 2
+double: 1
+idouble: 1
+ldouble: 428
+ildouble: 428
+Test "jn (8, 2.4048255576957729) == 0.92165786705344923232879022467054148E-4":
+float: 3
+ifloat: 3
+double: 1
+idouble: 1
+ldouble: 609
+ildouble: 609
+Test "jn (9, 2.4048255576957729) == 0.12517270977961513005428966643852564E-4":
+float: 4
+ifloat: 4
+ldouble: 750
+ildouble: 750
 
 # lgamma
 Test "lgamma (-0.5) == log(2*sqrt(pi))":
@@ -1168,11 +1214,11 @@ ldouble: 1
 
 Function: "jn":
 double: 5
-float: 2
+float: 4
 idouble: 5
-ifloat: 2
-ildouble: 2
-ldouble: 2
+ifloat: 4
+ildouble: 750
+ldouble: 750
 
 Function: "lgamma":
 double: 1

sysdeps/ieee754/dbl-64/e_jn.c

@@ -215,7 +215,16 @@ static double zero  =  0.00000000000000000000e+00;
 			}
 	     	    }
 		}
-	    	b = (t*__ieee754_j0(x)/b);
+		/* 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;

sysdeps/ieee754/flt-32/e_jnf.c

@@ -165,7 +165,16 @@ static float zero  =  0.0000000000e+00;
 			}
 		    }
 		}
-		b = (t*__ieee754_j0f(x)/b);
+		/* 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_j0f (x);
+		w = __ieee754_j1f (x);
+		if (fabsf (z) >= fabsf (w))
+		  b = (t * z / b);
+		else
+		  b = (t * w / a);
 	    }
 	}
 	if(sgn==1) return -b; else return b;

sysdeps/ieee754/ldbl-128/e_jnl.c

@@ -285,7 +285,16 @@ __ieee754_jnl (n, x)
 		    }
 		}
 	    }
-	  b = (t * __ieee754_j0l (x) / b);
+	  /* 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)

sysdeps/ieee754/ldbl-128ibm/e_jnl.c

@@ -286,7 +286,16 @@ __ieee754_jnl (n, x)
 		    }
 		}
 	    }
-	  b = (t * __ieee754_j0l (x) / b);
+	  /* 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)

sysdeps/ieee754/ldbl-96/e_jnl.c

@@ -281,7 +281,16 @@ __ieee754_jnl (n, x)
 		    }
 		}
 	    }
-	  b = (t * __ieee754_j0l (x) / b);
+	  /* 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)