Fix Bessel function spurious overflows for ldbl-128 / ldbl-128ibm (bug 15285).

diff --git a/ChangeLog b/ChangeLog
index 3ffaa38b5bac27f0e7d30106c0e992a60af5e90e..dffd409d5b5fd36d138aeb8ff4e0b2660c9b1198 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,20 @@
+2013-03-21  Joseph Myers  <joseph@codesourcery.com>
+
+       [BZ #15285]
+       * sysdeps/ieee754/ldbl-128/e_j0l.c: Include <float.h>.
+       (__ieee754_j0l): Do not improve calculations using cos of twice
+       input for inputs above LDBL_MAX / 2.0L.
+       (__ieee754_y0l): Likewise.
+       * sysdeps/ieee754/ldbl-128/e_j1l.c: Include <float.h>.
+       (__ieee754_j1l): Do not improve calculations using cos of twice
+       input for inputs above LDBL_MAX / 2.0L.
+       (__ieee754_y1l): Likewise.
+       * math/libm-test.inc (j0_test): Add another test.
+       (j1_test): Likewise.
+       (y0_test): Likewise.
+       (y1_test): Likewise.
+       * sysdeps/i386/fpu/libm-test-ulps: Update.
+
 2013-03-21  Siddhesh Poyarekar  <siddhesh@redhat.com>
 
        * Rules ($(objpfx)bench-%.c): Include code from a C source

diff --git a/NEWS b/NEWS
index 55f491aeed0f23d883022f5c2f0c0c7c61d61b6c..fe3f0b79115878ec47aed2d1b0775b624591409b 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -12,7 +12,7 @@ Version 2.18
   11561, 12723, 13550, 13951, 14142, 14176, 14200, 14317, 14327, 14496,
   14812, 14920, 14964, 14981, 14982, 14985, 14994, 14996, 15003, 15006,
   15020, 15023, 15036, 15054, 15055, 15062, 15078, 15160, 15232, 15234,
-  15283, 15287.
+  15283, 15285, 15287.
 
 * Add support for calling C++11 thread_local object destructors on thread
   and program exit.  This needs compiler support for offloading C++11

diff --git a/math/libm-test.inc b/math/libm-test.inc
index aafab40aac0d06be9b732cee0052f4998ce544f8..1b70c35f3923927e33aa99d2f38209c030287a1f 100644 (file)
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -6499,6 +6499,7 @@ j0_test (void)
 
 #ifndef TEST_FLOAT
   TEST_f_f (j0, -0x1.001000001p+593L, -3.927269966354206207832593635798954916263e-90L);
+  TEST_f_f (j0, 0x1p1023L, -1.5665258060609012834424478437196679802783e-155L);
 #endif
 
 #if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
@@ -6545,6 +6546,7 @@ j1_test (void)
 
 #ifndef TEST_FLOAT
   TEST_f_f (j1, 0x1.ff00000000002p+840L, 1.846591691699331493194965158699937660696e-127L);
+  TEST_f_f (j1, 0x1p1023L, 8.2687542933709649327986678723012001545638e-155L);
 #endif
 
 #if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
@@ -10721,6 +10723,7 @@ y0_test (void)
 
 #ifndef TEST_FLOAT
   TEST_f_f (y0, 0x1.ff00000000002p+840L, 1.846591691699331493194965158699937660696e-127L);
+  TEST_f_f (y0, 0x1p1023L, 8.2687542933709649327986678723012001545638e-155L);
 #endif
 
 #if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
@@ -10779,6 +10782,7 @@ y1_test (void)
 
 #ifndef TEST_FLOAT
   TEST_f_f (y1, 0x1.001000001p+593L, 3.927269966354206207832593635798954916263e-90L);
+  TEST_f_f (y1, 0x1p1023L, 1.5665258060609012834424478437196679802783e-155L);
 #endif
 
 #if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384

diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 560e27c936ba94863d656e4876a631cc54f742f4..27b96ce5f17d82caf853d76f9e0cece11044f29a 100644 (file)
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -3153,6 +3153,9 @@ ldouble: 2
 Test "j0 (0x1.d7ce3ap+107) == 2.775523647291230802651040996274861694514e-17":
 float: 1
 ifloat: 1
+Test "j0 (0x1p1023) == -1.5665258060609012834424478437196679802783e-155":
+double: 1
+idouble: 1
 Test "j0 (0x1p16382) == -1.2193782500509000574176799046642541129387e-2466":
 ildouble: 1
 ldouble: 1
@@ -4016,6 +4019,9 @@ ldouble: 1
 Test "y1 (0x1p-10) == -6.5190099301063115047395187618929589514382e+02":
 float: 1
 ifloat: 1
+Test "y1 (0x1p1023) == 1.5665258060609012834424478437196679802783e-155":
+double: 1
+idouble: 1
 Test "y1 (0x1p16382) == 1.2193782500509000574176799046642541129387e-2466":
 ildouble: 1
 ldouble: 1

diff --git a/sysdeps/ieee754/ldbl-128/e_j0l.c b/sysdeps/ieee754/ldbl-128/e_j0l.c
index 9e7880c49d1de96c6aabc035bf25d812aacaee09..108eff443500c817c29a60074cb0cba9066a5b30 100644 (file)
--- a/sysdeps/ieee754/ldbl-128/e_j0l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j0l.c
@@ -93,6 +93,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <float.h>
 
 /* 1 / sqrt(pi) */
 static const long double ONEOSQPI = 5.6418958354775628694807945156077258584405E-1L;
@@ -710,11 +711,14 @@ __ieee754_j0l (long double x)
   __sincosl (xx, &s, &c);
   ss = s - c;
   cc = s + c;
-  z = -__cosl (xx + xx);
-  if ((s * c) < 0)
-    cc = z / ss;
-  else
-    ss = z / cc;
+  if (xx <= LDBL_MAX / 2.0L)
+    {
+      z = -__cosl (xx + xx);
+      if ((s * c) < 0)
+       cc = z / ss;
+      else
+       ss = z / cc;
+    }
 
   if (xx > 0x1p256L)
     return ONEOSQPI * cc / __ieee754_sqrtl (xx);
@@ -857,11 +861,14 @@ long double
   __sincosl (x, &s, &c);
   ss = s - c;
   cc = s + c;
-  z = -__cosl (x + x);
-  if ((s * c) < 0)
-    cc = z / ss;
-  else
-    ss = z / cc;
+  if (xx <= LDBL_MAX / 2.0L)
+    {
+      z = -__cosl (x + x);
+      if ((s * c) < 0)
+       cc = z / ss;
+      else
+       ss = z / cc;
+    }
 
   if (xx > 0x1p256L)
     return ONEOSQPI * ss / __ieee754_sqrtl (x);

diff --git a/sysdeps/ieee754/ldbl-128/e_j1l.c b/sysdeps/ieee754/ldbl-128/e_j1l.c
index 95e01a39cc9e2db5652f66a29fdf67932edda487..70a1c86fd2d55f2fa89447bf6a6eb027452ed455 100644 (file)
--- a/sysdeps/ieee754/ldbl-128/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j1l.c
@@ -97,6 +97,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <float.h>
 
 /* 1 / sqrt(pi) */
 static const long double ONEOSQPI = 5.6418958354775628694807945156077258584405E-1L;
@@ -715,11 +716,14 @@ __ieee754_j1l (long double x)
   __sincosl (xx, &s, &c);
   ss = -s - c;
   cc = s - c;
-  z = __cosl (xx + xx);
-  if ((s * c) > 0)
-    cc = z / ss;
-  else
-    ss = z / cc;
+  if (xx <= LDBL_MAX / 2.0L)
+    {
+      z = __cosl (xx + xx);
+      if ((s * c) > 0)
+       cc = z / ss;
+      else
+       ss = z / cc;
+    }
 
   if (xx > 0x1p256L)
     {
@@ -868,11 +872,14 @@ __ieee754_y1l (long double x)
   __sincosl (xx, &s, &c);
   ss = -s - c;
   cc = s - c;
-  z = __cosl (xx + xx);
-  if ((s * c) > 0)
-    cc = z / ss;
-  else
-    ss = z / cc;
+  if (xx <= LDBL_MAX / 2.0L)
+    {
+      z = __cosl (xx + xx);
+      if ((s * c) > 0)
+       cc = z / ss;
+      else
+       ss = z / cc;
+    }
 
   if (xx > 0x1p256L)
     return ONEOSQPI * ss / __ieee754_sqrtl (xx);