This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix Bessel function spurious overflows for ldbl-128 / ldbl-128ibm(bug 15285)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Sun, 17 Mar 2013 20:49:28 +0000
- Subject: Fix Bessel function spurious overflows for ldbl-128 / ldbl-128ibm(bug 15285)
Where the Bessel function implementations improve the accuracy of
their calculations by using cos of twice the input, this requires a
cut-off to avoid internal overflow when twice the input would
overflow. Bug 15285 is the absence of that cut-off in the ldbl-128 /
ldbl-128ibm implementation. This patch adds the missing cut-off
checks and adds testcases that confirm such large inputs are handled
OK for double and IBM long double.
Tested x86_64, x86 and powerpc64.
2013-03-17 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.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 914aab3..a732ea6 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -6244,6 +6244,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
@@ -6290,6 +6291,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
@@ -10466,6 +10468,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
@@ -10524,6 +10527,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 6186c99..7e77567 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -2475,6 +2475,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
@@ -3338,6 +3341,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 9e7880c..108eff4 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 95e01a3..70a1c86 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);
--
Joseph S. Myers
joseph@codesourcery.com