Fix jn precision

7 files changed, 198 insertions, 53 deletions

diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 4b1a9e7..ebd46b0 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/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
diff --git a/sysdeps/ieee754/dbl-64/e_jn.c b/sysdeps/ieee754/dbl-64/e_jn.c
index bf4a13d..d9d6f91 100644
--- a/sysdeps/ieee754/dbl-64/e_jn.c
+++ b/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;
diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c
index de2e53d..dd3d551 100644
--- a/sysdeps/ieee754/flt-32/e_jnf.c
+++ 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;
diff --git a/sysdeps/ieee754/ldbl-128/e_jnl.c b/sysdeps/ieee754/ldbl-128/e_jnl.c
index a4a4e24..a7f6772 100644
--- a/sysdeps/ieee754/ldbl-128/e_jnl.c
+++ 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)
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
index 0eea745..372f942 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c
+++ b/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)
diff --git a/sysdeps/ieee754/ldbl-96/e_jnl.c b/sysdeps/ieee754/ldbl-96/e_jnl.c
index 3d715d3..bedff7d 100644
--- a/sysdeps/ieee754/ldbl-96/e_jnl.c
+++ b/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)
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 0ced4be..25135e3 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -103,8 +103,8 @@ double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ldouble: 1
 ildouble: 1
+ldouble: 1
 
 # catan
 Test "Real part of: catan (-2 - 3 i) == -1.4099210495965755225306193844604208 - 0.22907268296853876629588180294200276 i":
@@ -168,10 +168,10 @@ ifloat: 1
 ildouble: 1
 ldouble: 1
 Test "Imaginary part of: ccos (0.75 + 1.25 i) == 1.38173873063425888530729933139078645 - 1.09193013555397466170919531722024128 i":
-ildouble: 1
-ldouble: 1
 float: 1
 ifloat: 1
+ildouble: 1
+ldouble: 1
 
 # ccosh
 Test "Real part of: ccosh (-2 - 3 i) == -3.72454550491532256547397070325597253 + 0.511822569987384608834463849801875634 i":
@@ -304,6 +304,9 @@ idouble: 1
 ifloat: 1
 
 # cos
+Test "cos (0.80190127184058835) == 0.69534156199418473":
+double: 1
+idouble: 1
 Test "cos (M_PI_6l * 2.0) == 0.5":
 double: 1
 float: 1
@@ -323,16 +326,13 @@ idouble: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1
-Test "cos (0.80190127184058835) == 0.69534156199418473":
-double: 1
-idouble: 1
 
 # cpow
 Test "Real part of: cpow (0.75 + 1.25 i, 0.0 + 1.0 i) == 0.331825439177608832276067945276730566 + 0.131338600281188544930936345230903032 i":
 float: 1
 ifloat: 1
-ldouble: 1
 ildouble: 1
+ldouble: 1
 Test "Imaginary part of: cpow (0.75 + 1.25 i, 0.0 + 1.0 i) == 0.331825439177608832276067945276730566 + 0.131338600281188544930936345230903032 i":
 float: 1
 ifloat: 1
@@ -380,15 +380,15 @@ ildouble: 1
 ldouble: 1
 
 # csin
+Test "Imaginary part of: csin (-2 - 3 i) == -9.15449914691142957346729954460983256 + 4.16890695996656435075481305885375484 i":
+double: 1
+idouble: 1
 Test "Real part of: csin (0.75 + 1.25 i) == 1.28722291002649188575873510790565441 + 1.17210635989270256101081285116138863 i":
 ildouble: 1
 ldouble: 1
 Test "Imaginary part of: csin (0.75 + 1.25 i) == 1.28722291002649188575873510790565441 + 1.17210635989270256101081285116138863 i":
 float: 1
 ifloat: 1
-Test "Imaginary part of: csin (-2 - 3 i) == -9.15449914691142957346729954460983256 + 4.16890695996656435075481305885375484 i":
-double: 1
-idouble: 1
 
 # csinh
 Test "Real part of: csinh (-2 - 3 i) == 3.59056458998577995201256544779481679 - 0.530921086248519805267040090660676560 i":
@@ -440,12 +440,12 @@ ldouble: 3
 
 # ctanh
 Test "Real part of: ctanh (-2 - 3 i) == -0.965385879022133124278480269394560686 + 0.988437503832249372031403430350121098e-2 i":
+double: 1
 float: 2
+idouble: 1
 ifloat: 2
 ildouble: 5
 ldouble: 5
-double: 1
-idouble: 1
 Test "Imaginary part of: ctanh (-2 - 3 i) == -0.965385879022133124278480269394560686 + 0.988437503832249372031403430350121098e-2 i":
 ildouble: 25
 ldouble: 25
@@ -456,10 +456,10 @@ Test "Real part of: ctanh (0.75 + 1.25 i) == 1.372607570533783202580486065712268
 double: 1
 idouble: 1
 Test "Imaginary part of: ctanh (0.75 + 1.25 i) == 1.37260757053378320258048606571226857 + 0.385795952609750664177596760720790220 i":
-ildouble: 1
-ldouble: 1
 double: 1
 idouble: 1
+ildouble: 1
+ldouble: 1
 
 # erf
 Test "erf (1.25) == 0.922900128256458230136523481197281140":
@@ -481,26 +481,26 @@ ldouble: 1
 
 # exp10
 Test "exp10 (-1) == 0.1":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
 ildouble: 1
 ldouble: 1
+Test "exp10 (0.75) == 5.62341325190349080394951039776481231":
+double: 1
 float: 1
+idouble: 1
 ifloat: 1
-double: 2
-idouble: 2
-Test "exp10 (0.75) == 5.62341325190349080394951039776481231":
 ildouble: 2
 ldouble: 2
-float: 1
-ifloat: 1
-double: 1
-idouble: 1
 Test "exp10 (3) == 1000":
-ildouble: 8
-ldouble: 8
-float: 2
-ifloat: 2
 double: 6
+float: 2
 idouble: 6
+ifloat: 2
+ildouble: 8
+ldouble: 8
 
 # expm1
 Test "expm1 (0.75) == 1.11700001661267466854536981983709561":
@@ -660,10 +660,19 @@ ifloat: 3
 ildouble: 2
 ldouble: 2
 Test "jn (10, 2.0) == 0.251538628271673670963516093751820639e-6":
+double: 1
 float: 4
+idouble: 1
 ifloat: 4
 ildouble: 1
 ldouble: 1
+Test "jn (2, 2.4048255576957729) == 0.43175480701968038399746111312430703":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+ildouble: 82
+ldouble: 82
 Test "jn (3, -1.0) == -0.0195633539826684059189053216217515083":
 ildouble: 1
 ldouble: 1
@@ -694,6 +703,51 @@ idouble: 1
 ifloat: 2
 ildouble: 1
 ldouble: 1
+Test "jn (3, 2.4048255576957729) == 0.19899990535769083404042146764530813":
+double: 3
+idouble: 3
+ildouble: 186
+ldouble: 186
+Test "jn (4, 2.4048255576957729) == 0.647466661641779720084932282551219891E-1":
+double: 1
+idouble: 1
+ildouble: 185
+ldouble: 185
+Test "jn (5, 2.4048255576957729) == 0.163892432048058525099230549946147698E-1":
+double: 3
+float: 1
+idouble: 3
+ifloat: 1
+ildouble: 249
+ldouble: 249
+Test "jn (6, 2.4048255576957729) == 0.34048184720278336646673682895929161E-2":
+double: 4
+float: 3
+idouble: 4
+ifloat: 3
+ildouble: 511
+ldouble: 511
+Test "jn (7, 2.4048255576957729) == 0.60068836573295394221291569249883076E-3":
+double: 3
+float: 5
+idouble: 3
+ifloat: 5
+ildouble: 428
+ldouble: 428
+Test "jn (8, 2.4048255576957729) == 0.92165786705344923232879022467054148E-4":
+double: 3
+float: 2
+idouble: 3
+ifloat: 2
+ildouble: 609
+ldouble: 609
+Test "jn (9, 2.4048255576957729) == 0.12517270977961513005428966643852564E-4":
+double: 1
+float: 2
+idouble: 1
+ifloat: 2
+ildouble: 750
+ldouble: 750
 
 # lgamma
 Test "lgamma (-0.5) == log(2*sqrt(pi))":
@@ -714,12 +768,12 @@ ldouble: 1
 
 # log10
 Test "log10 (0.75) == -0.124938736608299953132449886193870744":
-ildouble: 1
-ldouble: 1
-float: 2
-ifloat: 2
 double: 1
+float: 2
 idouble: 1
+ifloat: 2
+ildouble: 1
+ldouble: 1
 Test "log10 (e) == log10(e)":
 float: 1
 ifloat: 1
@@ -732,6 +786,9 @@ float: 1
 ifloat: 1
 
 # sincos
+Test "sincos (0.80190127184058835, &sin_res, &cos_res) puts 0.69534156199418473 in cos_res":
+double: 1
+idouble: 1
 Test "sincos (M_PI_6l*2.0, &sin_res, &cos_res) puts 0.5 in cos_res":
 double: 1
 float: 1
@@ -754,9 +811,6 @@ ldouble: 1
 Test "sincos (pi/6, &sin_res, &cos_res) puts 0.86602540378443864676372317075293616 in cos_res":
 float: 1
 ifloat: 1
-Test "sincos (0.80190127184058835, &sin_res, &cos_res) puts 0.69534156199418473 in cos_res":
-double: 1
-idouble: 1
 
 # tan
 Test "tan (pi/4) == 1":
@@ -1178,12 +1232,12 @@ ildouble: 5
 ldouble: 5
 
 Function: Imaginary part of "ctanh":
+double: 1
 float: 1
+idouble: 1
 ifloat: 1
 ildouble: 25
 ldouble: 25
-double: 1
-idouble: 1
 
 Function: "erf":
 double: 1
@@ -1196,12 +1250,12 @@ ildouble: 1
 ldouble: 1
 
 Function: "exp10":
-ildouble: 8
-ldouble: 8
-float: 2
-ifloat: 2
 double: 6
+float: 2
 idouble: 6
+ifloat: 2
+ildouble: 8
+ldouble: 8
 
 Function: "expm1":
 double: 1
@@ -1235,11 +1289,11 @@ ldouble: 1
 
 Function: "jn":
 double: 4
-float: 4
+float: 5
 idouble: 4
-ifloat: 4
-ildouble: 2
-ldouble: 2
+ifloat: 5
+ildouble: 750
+ldouble: 750
 
 Function: "lgamma":
 double: 1
@@ -1250,12 +1304,12 @@ ildouble: 1
 ldouble: 1
 
 Function: "log10":
+double: 1
 float: 2
+idouble: 1
 ifloat: 2
 ildouble: 1
 ldouble: 1
-double: 1
-idouble: 1
 
 Function: "log1p":
 float: 1