about summary refs log tree commit diff
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2013-03-21 13:57:21 +0000
committerJoseph Myers <joseph@codesourcery.com>2013-03-21 13:57:21 +0000
commit98c48fe5cc4317123a168490a8fb37540e23f104 (patch)
treebef07c691890f1a2f0a339df2722c64a1b7c08a1
parent3775a8bc2d2e0c29c8a7e673f5f42537ced2b3c7 (diff)
downloadglibc-98c48fe5cc4317123a168490a8fb37540e23f104.tar.gz
glibc-98c48fe5cc4317123a168490a8fb37540e23f104.tar.xz
glibc-98c48fe5cc4317123a168490a8fb37540e23f104.zip
Fix Bessel function spurious overflows for ldbl-128 / ldbl-128ibm (bug 15285).
-rw-r--r--ChangeLog17
-rw-r--r--NEWS2
-rw-r--r--math/libm-test.inc4
-rw-r--r--sysdeps/i386/fpu/libm-test-ulps6
-rw-r--r--sysdeps/ieee754/ldbl-128/e_j0l.c27
-rw-r--r--sysdeps/ieee754/ldbl-128/e_j1l.c27
6 files changed, 62 insertions, 21 deletions
diff --git a/ChangeLog b/ChangeLog
index 3ffaa38b5b..dffd409d5b 100644
--- 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 55f491aeed..fe3f0b7911 100644
--- 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 aafab40aac..1b70c35f39 100644
--- 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 560e27c936..27b96ce5f1 100644
--- 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 9e7880c49d..108eff4435 100644
--- 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 95e01a39cc..70a1c86fd2 100644
--- 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);