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

This commit is contained in:
Joseph Myers 2013-03-21 13:57:21 +00:00
parent 3775a8bc2d
commit 98c48fe5cc
6 changed files with 62 additions and 21 deletions

View File

@ -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

2
NEWS
View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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);

View File

@ -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);