Fix ldbl-128 j1l spurious underflows (bug 18612).

The ldbl-128 implementation of j1l produces spurious underflow
exceptions for some small arguments, as a result of squaring the
argument.  This patch fixes it just to use a linear approximation for
sufficiently small arguments, and then to force an underflow exception
only in the cases where it is required.

Tested for mips64.

	[BZ #18612]
	* sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_j1l): For small
	arguments, just return 0.5 times the argument, with underflow
	forced as needed.
	* math/auto-libm-test-in: Add more tests of j1.
	* math/auto-libm-test-out: Regenerated.

5 files changed, 196 insertions, 1 deletions

diff --git a/ChangeLog b/ChangeLog
index cd03015ad8..9f91f183b5 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,12 @@
 2015-06-29  Joseph Myers  <joseph@codesourcery.com>
 
+	[BZ #18612]
+	* sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_j1l): For small
+	arguments, just return 0.5 times the argument, with underflow
+	forced as needed.
+	* math/auto-libm-test-in: Add more tests of j1.
+	* math/auto-libm-test-out: Regenerated.
+
 	[BZ #16559]
 	* sysdeps/ieee754/dbl-64/e_j1.c: Include <float.h>.
 	(__ieee754_j1): Force underflow exception for small results.
diff --git a/NEWS b/NEWS
index c454e912be..08899a8162 100644
--- a/NEWS
+++ b/NEWS
@@ -25,7 +25,7 @@ Version 2.22
   18497, 18498, 18502, 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, 18602.
+  18586, 18593, 18594, 18602, 18612.
 
 * Cache information can be queried via sysconf() function on s390 e.g. with
   _SC_LEVEL1_ICACHE_SIZE as argument.
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 34b02c91f7..f4313a2b99 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -1812,6 +1812,9 @@ j1 0x1.ff00000000002p+840
 j1 0x1p1023
 j1 0x1p16382
 j1 0x1p16383
+j1 0x1p-100
+j1 0x1p-600
+j1 0x1p-10000
 # Bug 18611: errno setting may be missing.
 j1 min missing-errno
 j1 -min missing-errno
diff --git a/math/auto-libm-test-out b/math/auto-libm-test-out
index 8ae005d081..b71429b4ac 100644
--- a/math/auto-libm-test-out
+++ b/math/auto-libm-test-out
@@ -142437,6 +142437,181 @@ j1 0x1p16383
 = j1 tonearest ldbl-128ibm 0xf.ffffffffffffbffffffffffffcp+1020L : -0x5.34b8be56c9cb044a0ef191b026p-516L : inexact-ok
 = j1 towardzero ldbl-128ibm 0xf.ffffffffffffbffffffffffffcp+1020L : -0x5.34b8be56c9cb044a0ef191b024p-516L : inexact-ok
 = j1 upward ldbl-128ibm 0xf.ffffffffffffbffffffffffffcp+1020L : -0x5.34b8be56c9cb044a0ef191b024p-516L : inexact-ok
+j1 0x1p-100
+= j1 downward flt-32 0x1p-100f : 0x7.fffff8p-104f : inexact-ok
+= j1 tonearest flt-32 0x1p-100f : 0x8p-104f : inexact-ok
+= j1 towardzero flt-32 0x1p-100f : 0x7.fffff8p-104f : inexact-ok
+= j1 upward flt-32 0x1p-100f : 0x8p-104f : inexact-ok
+= j1 downward dbl-64 0x1p-100 : 0x7.ffffffffffffcp-104 : inexact-ok
+= j1 tonearest dbl-64 0x1p-100 : 0x8p-104 : inexact-ok
+= j1 towardzero dbl-64 0x1p-100 : 0x7.ffffffffffffcp-104 : inexact-ok
+= j1 upward dbl-64 0x1p-100 : 0x8p-104 : inexact-ok
+= j1 downward ldbl-96-intel 0x1p-100L : 0x7.fffffffffffffff8p-104L : inexact-ok
+= j1 tonearest ldbl-96-intel 0x1p-100L : 0x8p-104L : inexact-ok
+= j1 towardzero ldbl-96-intel 0x1p-100L : 0x7.fffffffffffffff8p-104L : inexact-ok
+= j1 upward ldbl-96-intel 0x1p-100L : 0x8p-104L : inexact-ok
+= j1 downward ldbl-96-m68k 0x1p-100L : 0x7.fffffffffffffff8p-104L : inexact-ok
+= j1 tonearest ldbl-96-m68k 0x1p-100L : 0x8p-104L : inexact-ok
+= j1 towardzero ldbl-96-m68k 0x1p-100L : 0x7.fffffffffffffff8p-104L : inexact-ok
+= j1 upward ldbl-96-m68k 0x1p-100L : 0x8p-104L : inexact-ok
+= j1 downward ldbl-128 0x1p-100L : 0x7.fffffffffffffffffffffffffffcp-104L : inexact-ok
+= j1 tonearest ldbl-128 0x1p-100L : 0x8p-104L : inexact-ok
+= j1 towardzero ldbl-128 0x1p-100L : 0x7.fffffffffffffffffffffffffffcp-104L : inexact-ok
+= j1 upward ldbl-128 0x1p-100L : 0x8p-104L : inexact-ok
+= j1 downward ldbl-128ibm 0x1p-100L : 0x7.fffffffffffffffffffffffffep-104L : inexact-ok
+= j1 tonearest ldbl-128ibm 0x1p-100L : 0x8p-104L : inexact-ok
+= j1 towardzero ldbl-128ibm 0x1p-100L : 0x7.fffffffffffffffffffffffffep-104L : inexact-ok
+= j1 upward ldbl-128ibm 0x1p-100L : 0x8p-104L : inexact-ok
+j1 0x1p-600
+= j1 downward flt-32 0x8p-152f : 0x0p+0f : inexact-ok underflow errno-erange
+= j1 tonearest flt-32 0x8p-152f : 0x0p+0f : inexact-ok underflow errno-erange
+= j1 towardzero flt-32 0x8p-152f : 0x0p+0f : inexact-ok underflow errno-erange
+= j1 upward flt-32 0x8p-152f : 0x8p-152f : inexact-ok underflow errno-erange-ok
+= j1 downward dbl-64 0x8p-152 : 0x3.ffffffffffffep-152 : inexact-ok
+= j1 tonearest dbl-64 0x8p-152 : 0x4p-152 : inexact-ok
+= j1 towardzero dbl-64 0x8p-152 : 0x3.ffffffffffffep-152 : inexact-ok
+= j1 upward dbl-64 0x8p-152 : 0x4p-152 : inexact-ok
+= j1 downward ldbl-96-intel 0x8p-152L : 0x3.fffffffffffffffcp-152L : inexact-ok
+= j1 tonearest ldbl-96-intel 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 towardzero ldbl-96-intel 0x8p-152L : 0x3.fffffffffffffffcp-152L : inexact-ok
+= j1 upward ldbl-96-intel 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 downward ldbl-96-m68k 0x8p-152L : 0x3.fffffffffffffffcp-152L : inexact-ok
+= j1 tonearest ldbl-96-m68k 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 towardzero ldbl-96-m68k 0x8p-152L : 0x3.fffffffffffffffcp-152L : inexact-ok
+= j1 upward ldbl-96-m68k 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 downward ldbl-128 0x8p-152L : 0x3.fffffffffffffffffffffffffffep-152L : inexact-ok
+= j1 tonearest ldbl-128 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 towardzero ldbl-128 0x8p-152L : 0x3.fffffffffffffffffffffffffffep-152L : inexact-ok
+= j1 upward ldbl-128 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 downward ldbl-128ibm 0x8p-152L : 0x3.ffffffffffffffffffffffffffp-152L : inexact-ok
+= j1 tonearest ldbl-128ibm 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 towardzero ldbl-128ibm 0x8p-152L : 0x3.ffffffffffffffffffffffffffp-152L : inexact-ok
+= j1 upward ldbl-128ibm 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 downward flt-32 0x0p+0f : 0x0p+0f : inexact-ok
+= j1 tonearest flt-32 0x0p+0f : 0x0p+0f : inexact-ok
+= j1 towardzero flt-32 0x0p+0f : 0x0p+0f : inexact-ok
+= j1 upward flt-32 0x0p+0f : 0x0p+0f : inexact-ok
+= j1 downward dbl-64 0x0p+0 : 0x0p+0 : inexact-ok
+= j1 tonearest dbl-64 0x0p+0 : 0x0p+0 : inexact-ok
+= j1 towardzero dbl-64 0x0p+0 : 0x0p+0 : inexact-ok
+= j1 upward dbl-64 0x0p+0 : 0x0p+0 : inexact-ok
+= j1 downward ldbl-96-intel 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 tonearest ldbl-96-intel 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 towardzero ldbl-96-intel 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 upward ldbl-96-intel 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 downward ldbl-96-m68k 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 tonearest ldbl-96-m68k 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 towardzero ldbl-96-m68k 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 upward ldbl-96-m68k 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 downward ldbl-128 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 tonearest ldbl-128 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 towardzero ldbl-128 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 upward ldbl-128 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 downward ldbl-128ibm 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 tonearest ldbl-128ibm 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 towardzero ldbl-128ibm 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 upward ldbl-128ibm 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 downward dbl-64 0x1p-600 : 0x7.ffffffffffffcp-604 : inexact-ok
+= j1 tonearest dbl-64 0x1p-600 : 0x8p-604 : inexact-ok
+= j1 towardzero dbl-64 0x1p-600 : 0x7.ffffffffffffcp-604 : inexact-ok
+= j1 upward dbl-64 0x1p-600 : 0x8p-604 : inexact-ok
+= j1 downward ldbl-96-intel 0x1p-600L : 0x7.fffffffffffffff8p-604L : inexact-ok
+= j1 tonearest ldbl-96-intel 0x1p-600L : 0x8p-604L : inexact-ok
+= j1 towardzero ldbl-96-intel 0x1p-600L : 0x7.fffffffffffffff8p-604L : inexact-ok
+= j1 upward ldbl-96-intel 0x1p-600L : 0x8p-604L : inexact-ok
+= j1 downward ldbl-96-m68k 0x1p-600L : 0x7.fffffffffffffff8p-604L : inexact-ok
+= j1 tonearest ldbl-96-m68k 0x1p-600L : 0x8p-604L : inexact-ok
+= j1 towardzero ldbl-96-m68k 0x1p-600L : 0x7.fffffffffffffff8p-604L : inexact-ok
+= j1 upward ldbl-96-m68k 0x1p-600L : 0x8p-604L : inexact-ok
+= j1 downward ldbl-128 0x1p-600L : 0x7.fffffffffffffffffffffffffffcp-604L : inexact-ok
+= j1 tonearest ldbl-128 0x1p-600L : 0x8p-604L : inexact-ok
+= j1 towardzero ldbl-128 0x1p-600L : 0x7.fffffffffffffffffffffffffffcp-604L : inexact-ok
+= j1 upward ldbl-128 0x1p-600L : 0x8p-604L : inexact-ok
+= j1 downward ldbl-128ibm 0x1p-600L : 0x7.fffffffffffffffffffffffffep-604L : inexact-ok
+= j1 tonearest ldbl-128ibm 0x1p-600L : 0x8p-604L : inexact-ok
+= j1 towardzero ldbl-128ibm 0x1p-600L : 0x7.fffffffffffffffffffffffffep-604L : inexact-ok
+= j1 upward ldbl-128ibm 0x1p-600L : 0x8p-604L : inexact-ok
+j1 0x1p-10000
+= j1 downward flt-32 0x8p-152f : 0x0p+0f : inexact-ok underflow errno-erange
+= j1 tonearest flt-32 0x8p-152f : 0x0p+0f : inexact-ok underflow errno-erange
+= j1 towardzero flt-32 0x8p-152f : 0x0p+0f : inexact-ok underflow errno-erange
+= j1 upward flt-32 0x8p-152f : 0x8p-152f : inexact-ok underflow errno-erange-ok
+= j1 downward dbl-64 0x8p-152 : 0x3.ffffffffffffep-152 : inexact-ok
+= j1 tonearest dbl-64 0x8p-152 : 0x4p-152 : inexact-ok
+= j1 towardzero dbl-64 0x8p-152 : 0x3.ffffffffffffep-152 : inexact-ok
+= j1 upward dbl-64 0x8p-152 : 0x4p-152 : inexact-ok
+= j1 downward ldbl-96-intel 0x8p-152L : 0x3.fffffffffffffffcp-152L : inexact-ok
+= j1 tonearest ldbl-96-intel 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 towardzero ldbl-96-intel 0x8p-152L : 0x3.fffffffffffffffcp-152L : inexact-ok
+= j1 upward ldbl-96-intel 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 downward ldbl-96-m68k 0x8p-152L : 0x3.fffffffffffffffcp-152L : inexact-ok
+= j1 tonearest ldbl-96-m68k 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 towardzero ldbl-96-m68k 0x8p-152L : 0x3.fffffffffffffffcp-152L : inexact-ok
+= j1 upward ldbl-96-m68k 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 downward ldbl-128 0x8p-152L : 0x3.fffffffffffffffffffffffffffep-152L : inexact-ok
+= j1 tonearest ldbl-128 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 towardzero ldbl-128 0x8p-152L : 0x3.fffffffffffffffffffffffffffep-152L : inexact-ok
+= j1 upward ldbl-128 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 downward ldbl-128ibm 0x8p-152L : 0x3.ffffffffffffffffffffffffffp-152L : inexact-ok
+= j1 tonearest ldbl-128ibm 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 towardzero ldbl-128ibm 0x8p-152L : 0x3.ffffffffffffffffffffffffffp-152L : inexact-ok
+= j1 upward ldbl-128ibm 0x8p-152L : 0x4p-152L : inexact-ok
+= j1 downward flt-32 0x0p+0f : 0x0p+0f : inexact-ok
+= j1 tonearest flt-32 0x0p+0f : 0x0p+0f : inexact-ok
+= j1 towardzero flt-32 0x0p+0f : 0x0p+0f : inexact-ok
+= j1 upward flt-32 0x0p+0f : 0x0p+0f : inexact-ok
+= j1 downward dbl-64 0x0p+0 : 0x0p+0 : inexact-ok
+= j1 tonearest dbl-64 0x0p+0 : 0x0p+0 : inexact-ok
+= j1 towardzero dbl-64 0x0p+0 : 0x0p+0 : inexact-ok
+= j1 upward dbl-64 0x0p+0 : 0x0p+0 : inexact-ok
+= j1 downward ldbl-96-intel 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 tonearest ldbl-96-intel 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 towardzero ldbl-96-intel 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 upward ldbl-96-intel 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 downward ldbl-96-m68k 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 tonearest ldbl-96-m68k 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 towardzero ldbl-96-m68k 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 upward ldbl-96-m68k 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 downward ldbl-128 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 tonearest ldbl-128 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 towardzero ldbl-128 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 upward ldbl-128 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 downward ldbl-128ibm 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 tonearest ldbl-128ibm 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 towardzero ldbl-128ibm 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 upward ldbl-128ibm 0x0p+0L : 0x0p+0L : inexact-ok
+= j1 downward dbl-64 0x4p-1076 : 0x0p+0 : inexact-ok underflow errno-erange
+= j1 tonearest dbl-64 0x4p-1076 : 0x0p+0 : inexact-ok underflow errno-erange
+= j1 towardzero dbl-64 0x4p-1076 : 0x0p+0 : inexact-ok underflow errno-erange
+= j1 upward dbl-64 0x4p-1076 : 0x4p-1076 : inexact-ok underflow errno-erange-ok
+= j1 downward ldbl-96-intel 0x4p-1076L : 0x1.fffffffffffffffep-1076L : inexact-ok
+= j1 tonearest ldbl-96-intel 0x4p-1076L : 0x2p-1076L : inexact-ok
+= j1 towardzero ldbl-96-intel 0x4p-1076L : 0x1.fffffffffffffffep-1076L : inexact-ok
+= j1 upward ldbl-96-intel 0x4p-1076L : 0x2p-1076L : inexact-ok
+= j1 downward ldbl-96-m68k 0x4p-1076L : 0x1.fffffffffffffffep-1076L : inexact-ok
+= j1 tonearest ldbl-96-m68k 0x4p-1076L : 0x2p-1076L : inexact-ok
+= j1 towardzero ldbl-96-m68k 0x4p-1076L : 0x1.fffffffffffffffep-1076L : inexact-ok
+= j1 upward ldbl-96-m68k 0x4p-1076L : 0x2p-1076L : inexact-ok
+= j1 downward ldbl-128 0x4p-1076L : 0x1.ffffffffffffffffffffffffffffp-1076L : inexact-ok
+= j1 tonearest ldbl-128 0x4p-1076L : 0x2p-1076L : inexact-ok
+= j1 towardzero ldbl-128 0x4p-1076L : 0x1.ffffffffffffffffffffffffffffp-1076L : inexact-ok
+= j1 upward ldbl-128 0x4p-1076L : 0x2p-1076L : inexact-ok
+= j1 downward ldbl-128ibm 0x4p-1076L : 0x0p+0L : inexact-ok underflow errno-erange
+= j1 tonearest ldbl-128ibm 0x4p-1076L : 0x0p+0L : inexact-ok underflow errno-erange
+= j1 towardzero ldbl-128ibm 0x4p-1076L : 0x0p+0L : inexact-ok underflow errno-erange
+= j1 upward ldbl-128ibm 0x4p-1076L : 0x4p-1076L : inexact-ok underflow errno-erange-ok
+= j1 downward ldbl-96-intel 0x1p-10000L : 0x7.fffffffffffffff8p-10004L : inexact-ok
+= j1 tonearest ldbl-96-intel 0x1p-10000L : 0x8p-10004L : inexact-ok
+= j1 towardzero ldbl-96-intel 0x1p-10000L : 0x7.fffffffffffffff8p-10004L : inexact-ok
+= j1 upward ldbl-96-intel 0x1p-10000L : 0x8p-10004L : inexact-ok
+= j1 downward ldbl-96-m68k 0x1p-10000L : 0x7.fffffffffffffff8p-10004L : inexact-ok
+= j1 tonearest ldbl-96-m68k 0x1p-10000L : 0x8p-10004L : inexact-ok
+= j1 towardzero ldbl-96-m68k 0x1p-10000L : 0x7.fffffffffffffff8p-10004L : inexact-ok
+= j1 upward ldbl-96-m68k 0x1p-10000L : 0x8p-10004L : inexact-ok
+= j1 downward ldbl-128 0x1p-10000L : 0x7.fffffffffffffffffffffffffffcp-10004L : inexact-ok
+= j1 tonearest ldbl-128 0x1p-10000L : 0x8p-10004L : inexact-ok
+= j1 towardzero ldbl-128 0x1p-10000L : 0x7.fffffffffffffffffffffffffffcp-10004L : inexact-ok
+= j1 upward ldbl-128 0x1p-10000L : 0x8p-10004L : inexact-ok
 j1 min missing-errno
 = j1 downward flt-32 0x4p-128f : 0x1.fffff8p-128f : inexact-ok underflow errno-erange-ok
 = j1 tonearest flt-32 0x4p-128f : 0x2p-128f : inexact-ok underflow errno-erange-ok
diff --git a/sysdeps/ieee754/ldbl-128/e_j1l.c b/sysdeps/ieee754/ldbl-128/e_j1l.c
index 958077d927..591c38efd0 100644
--- a/sysdeps/ieee754/ldbl-128/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j1l.c
@@ -697,6 +697,16 @@ __ieee754_j1l (long double x)
   if (x == 0.0L)
     return x;
   xx = fabsl (x);
+  if (xx <= 0x1p-58L)
+    {
+      long double ret = x * 0.5L;
+      if (fabsl (ret) < LDBL_MIN)
+	{
+	  long double force_underflow = ret * ret;
+	  math_force_eval (force_underflow);
+	}
+      return ret;
+    }
   if (xx <= 2.0L)
     {
       /* 0 <= x <= 2 */