From 98c48fe5cc4317123a168490a8fb37540e23f104 Mon Sep 17 00:00:00 2001
From: Joseph Myers <joseph@codesourcery.com>
Date: Thu, 21 Mar 2013 13:57:21 +0000
Subject: Fix Bessel function spurious overflows for ldbl-128 / ldbl-128ibm
 (bug 15285).

---
 sysdeps/i386/fpu/libm-test-ulps  |  6 ++++++
 sysdeps/ieee754/ldbl-128/e_j0l.c | 27 +++++++++++++++++----------
 sysdeps/ieee754/ldbl-128/e_j1l.c | 27 +++++++++++++++++----------
 3 files changed, 40 insertions(+), 20 deletions(-)

(limited to 'sysdeps')

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);
-- 
cgit 1.4.1