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.

diff --git a/ChangeLog b/ChangeLog
index cd03015ad804cb0c02b9873fc505f14cf0c121d1..9f91f183b58712fc680938f22ea2758268ef4ff3 100644 (file)
--- 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 c454e912bebe54b59baebe3b337a1b05ab962c97..08899a8162be90c8ad5d1b5a2f1f6e0f2aafde94 100644 (file)
--- 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 34b02c91f7ff1d69eb0001ddef7ed2b8ccde9ff1..f4313a2b99b72ebf3f49bc16aeef427fee06d3ce 100644 (file)
--- 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 8ae005d081560c7a8f669929a5a80c55727d89c2..b71429b4ac2fe279467ca5f7262b9f8502d5c0b4 100644 (file)
--- 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 958077d92773c3debd1bd979fc4c0545e74b5ef9..591c38efd0f0e38d45d3a833ebdf010cd1bb65fb 100644 (file)
--- 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 */