Fix spurious jnf underflows (bug 14155)

["Joseph S. Myers" <joseph at codesourcery dot com>]

This patch fixes the remaining piece of bug 14155, spurious underflow
exceptions from jnf for large arguments.

The exceptions arise from computing Bessel functions of order n using
the recurrence from orders n-1 and n-2, where the computation

                b = b*((float)(i+i)/x) - a; /* avoid underflow */

does, despite the comment, have a spurious underflow roughly when
x^-1.5 underflows (the Bessel functions being on the order of x^-0.5
multiplied by a trig function).  The same issue applies for ynf.  It
does not apply for the functions for other types, because they use a
cut-off of 2^302, above which first-order approximations to magnitude
and phase are used.  Those calculations can suffer from cancellation
when computing sin +/- cos - they don't do the

 *         (To avoid cancellation, use
 *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 *          to compute the worse one.)

used in j0/j1 to avoid that problem - but not spurious underflow.

The 2^302 cut-off appears to come from IEEE binary128.  Suppose a
Bessel function is being computed for a type with M mantissa bits, and
a total of N bits in its representation.  Heuristically (and the
bounds could be made precise for any given floating-point type by
examining continued-fraction approximations to appropriate multiples
of pi), large values representable in that type are going to be no
closer than about 2^-N to a multiple of pi/4 (heuristically treating
the distances from multiples of pi/4 as random).  The first-order
approximation to the phase of a Bessel function of order n is x - (n/2
+ 1/4)pi, so if the error is less than 2^-(N+M) it can be neglected.
Now the first-order error term is (4n^2 - 1)/8x, and n is at most
2^31, so we want x > 2^(61+N+M).  That gives 2^302 in the case where M
= 113 and N = 128.  Unfortunately, for binary32 it gives 2^117, when a
cut-off around 2^84 or below would be needed to avoid underflow even
when the trig functions are approximated by 1 (and a bit smaller when
you allow for the trig functions being as small as 2^-32).

(For the magnitude, you have sqrt((2/(pi*x))(1 + (4n^2 - 1)/8x^2 +
...)), so 2^61/x^2 < 2^-M suffices to neglect the error terms there.)

So the cut-off approach used to avoid underflow for other types cannot
safely be used for float.  Instead, this patch uses the approach of
doing the intermediate calculation for each step of the recurrence in
type double, where the relevant division will not underflow (and
neither will the final result of each recurrence step, which still
gets converted back to float).

Tested x86_64 and x86 and ulps updated accordingly.

2013-08-30  Joseph Myers  <joseph@codesourcery.com>

	[BZ #14155]
	* sysdeps/ieee754/flt-32/e_jnf.c (__ieee754_jnf): Use double for
	intermediate calculations in recurrence.
	(__ieee754_ynf): Likewise.
	* math/libm-test.inc (jn_test_data): Do not allow spurious
	underflow exception.  Add more tests.
	(yn_test_data): Add more tests.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.

diff --git a/math/libm-test.inc b/math/libm-test.inc
index e534fc0..fb5e977 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -9913,8 +9913,14 @@ static const struct test_if_f_data jn_test_data[] =
     TEST_if_f (jn, 8, 2.4048255576957729L, 0.92165786705344923232879022467054148E-4L),
     TEST_if_f (jn, 9, 2.4048255576957729L, 0.12517270977961513005428966643852564E-4L),
 
-    /* Bug 14155: spurious exception may occur.  */
-    TEST_if_f (jn, 2, 0x1.ffff62p+99L, -4.43860668048170034334926693188979974489e-16L, UNDERFLOW_EXCEPTION_OK),
+    TEST_if_f (jn, 2, 0x1.ffff62p+99L, -4.43860668048170034334926693188979974489e-16L),
+    TEST_if_f (jn, 2, 0x1p127L, -6.0784021821505059176832624052765568656702e-20L),
+#ifndef TEST_FLOAT
+    TEST_if_f (jn, 2, 0x1p1023L, 1.5665258060609012834424478437196679802783e-155L),
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+    TEST_if_f (jn, 2, 0x1p16383L, -9.5859502826270374691362975419147645151233e-2467L),
+#endif
   };
 
 static void
@@ -14531,6 +14537,15 @@ static const struct test_if_f_data yn_test_data[] =
     /* Check whether yn returns correct value for LDBL_MIN, DBL_MIN,
        and FLT_MIN.  See Bug 14173.  */
     TEST_if_f (yn, 10, min_value, minus_infty, OVERFLOW_EXCEPTION|ERRNO_ERANGE),
+
+    TEST_if_f (yn, 2, 0x1.ffff62p+99L, -5.5244413477397111790415387179517953221757e-16L),
+    TEST_if_f (yn, 2, 0x1p127L, 6.8569250690166637098111268958532649249771e-21L),
+#ifndef TEST_FLOAT
+    TEST_if_f (yn, 2, 0x1p1023L, -8.2687542933709649327986678723012001545638e-155L),
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+    TEST_if_f (yn, 2, 0x1p16383L, 3.8895531955766020648617743624167352352217e-2467L),
+#endif
   };
 
 static void
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 530dbd7..8244863 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -5635,9 +5635,9 @@ ildouble: 1
 ldouble: 1
 Test "jn (10, 10.0)":
 double: 1
-float: 1
+float: 2
 idouble: 1
-ifloat: 1
+ifloat: 2
 ildouble: 2
 ldouble: 2
 Test "jn (10, 2.0)":
@@ -5648,6 +5648,14 @@ float: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1
+Test "jn (2, 0x1p1023)":
+double: 1
+idouble: 1
+Test "jn (2, 0x1p127)":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
 Test "jn (2, 2.4048255576957729)":
 double: 1
 float: 1
@@ -6844,6 +6852,14 @@ ifloat: 1
 Test "yn (10, 2.0)":
 float: 3
 ifloat: 3
+Test "yn (2, 0x1.ffff62p+99)":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "yn (2, 0x1p127)":
+float: 2
+ifloat: 2
 Test "yn (3, 0.125)":
 ildouble: 1
 ldouble: 1
diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c
index ad26d7e..5984d94 100644
--- a/sysdeps/ieee754/flt-32/e_jnf.c
+++ b/sysdeps/ieee754/flt-32/e_jnf.c
@@ -54,7 +54,7 @@ __ieee754_jnf(int n, float x)
 	    b = __ieee754_j1f(x);
 	    for(i=1;i<n;i++){
 		temp = b;
-		b = b*((float)(i+i)/x) - a; /* avoid underflow */
+		b = b*((double)(i+i)/x) - a; /* avoid underflow */
 		a = temp;
 	    }
 	} else {
@@ -196,7 +196,7 @@ __ieee754_ynf(int n, float x)
 	GET_FLOAT_WORD(ib,b);
 	for(i=1;i<n&&ib!=0xff800000;i++){
 	    temp = b;
-	    b = ((float)(i+i)/x)*b - a;
+	    b = ((double)(i+i)/x)*b - a;
 	    GET_FLOAT_WORD(ib,b);
 	    a = temp;
 	}
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index d02618a..477eedc 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -6403,6 +6403,11 @@ idouble: 2
 ifloat: 2
 ildouble: 1
 ldouble: 1
+Test "jn (2, 0x1p127)":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
 Test "jn (2, 2.4048255576957729)":
 double: 2
 float: 1
@@ -7728,6 +7733,16 @@ double: 3
 float: 1
 idouble: 3
 ifloat: 1
+Test "yn (2, 0x1.ffff62p+99)":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "yn (2, 0x1p127)":
+double: 1
+float: 3
+idouble: 1
+ifloat: 3
 Test "yn (3, 0.125)":
 double: 1
 idouble: 1
@@ -8428,9 +8443,9 @@ ldouble: 2
 
 Function: "yn":
 double: 3
-float: 2
+float: 3
 idouble: 3
-ifloat: 2
+ifloat: 3
 ildouble: 4
 ldouble: 4
 

-- 
Joseph S. Myers
joseph@codesourcery.com