AuroraMiddleware/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2025-01-18 06:30:05 +00:00
Code Issues Packages Projects Releases Wiki Activity

Fix yn overflow handling in non-default rounding modes (bug 16561, bug 16562).

Browse Source 
This patch fixes bugs 16561 and 16562, bad results of yn in overflow
cases in non-default rounding modes, both because an intermediate
overflow in the recurrence does not get detected if the result is not
an infinity and because an overflowing result may occur in the wrong
sign.  The fix is to set FE_TONEAREST mode internally for the parts of
the function where such overflows can occur (which includes the call
to y1 - where yn is used to compute a Bessel function of order -1,
negating the result of y1 isn't correct for overflowing results in
directed rounding modes) and then compute an overflowing value in the
original rounding mode if the to-nearest result was an infinity.

Tested x86_64 and x86 and ulps updated accordingly.  Also tested for
mips64 and powerpc32 to test the ldbl-128 and ldbl-128ibm changes.

(The tests for these bugs were added in my previous y1 patch, so the
only thing this patch has to do with the testsuite is enable yn
testing in all rounding modes.)

	[BZ #16561]
	[BZ #16562]
	* sysdeps/ieee754/dbl-64/e_jn.c: Include <float.h>.
	(__ieee754_yn): Set FE_TONEAREST mode internally and then
	recompute overflowing results in original rounding mode.
	* sysdeps/ieee754/flt-32/e_jnf.c: Include <float.h>.
	(__ieee754_ynf): Set FE_TONEAREST mode internally and then
	recompute overflowing results in original rounding mode.
	* sysdeps/ieee754/ldbl-128/e_jnl.c: Include <float.h>.
	(__ieee754_ynl): Set FE_TONEAREST mode internally and then
	recompute overflowing results in original rounding mode.
	* sysdeps/ieee754/ldbl-128ibm/e_jnl.c: Include <float.h>.
	(__ieee754_ynl): Set FE_TONEAREST mode internally and then
	recompute overflowing results in original rounding mode.
	* sysdeps/ieee754/ldbl-96/e_jnl.c: Include <float.h>.
	(__ieee754_ynl): Set FE_TONEAREST mode internally and then
	recompute overflowing results in original rounding mode.
	* sysdeps/i386/fpu/fenv_private.h [!__SSE2_MATH__]
	(libc_feholdsetround_ctx): New macro.
	* math/libm-test.inc (yn_test): Use ALL_RM_TEST.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps : Likewise.
This commit is contained in:
Joseph Myers 2014-06-27 14:52:13 +00:00
parent b351d85aa2
commit be25493251
11 changed files with 389 additions and 261 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
ChangeLog
View File

@ -1,3 +1,28 @@
2014-06-27 Joseph Myers <joseph@codesourcery.com>
[BZ #16561]
[BZ #16562]
* sysdeps/ieee754/dbl-64/e_jn.c: Include <float.h>.
(__ieee754_yn): Set FE_TONEAREST mode internally and then
recompute overflowing results in original rounding mode.
* sysdeps/ieee754/flt-32/e_jnf.c: Include <float.h>.
(__ieee754_ynf): Set FE_TONEAREST mode internally and then
recompute overflowing results in original rounding mode.
* sysdeps/ieee754/ldbl-128/e_jnl.c: Include <float.h>.
(__ieee754_ynl): Set FE_TONEAREST mode internally and then
recompute overflowing results in original rounding mode.
* sysdeps/ieee754/ldbl-128ibm/e_jnl.c: Include <float.h>.
(__ieee754_ynl): Set FE_TONEAREST mode internally and then
recompute overflowing results in original rounding mode.
* sysdeps/ieee754/ldbl-96/e_jnl.c: Include <float.h>.
(__ieee754_ynl): Set FE_TONEAREST mode internally and then
recompute overflowing results in original rounding mode.
* sysdeps/i386/fpu/fenv_private.h [!__SSE2_MATH__]
(libc_feholdsetround_ctx): New macro.
* math/libm-test.inc (yn_test): Use ALL_RM_TEST.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps : Likewise.
2014-06-26 Joseph Myers <joseph@codesourcery.com> 2014-06-26 Joseph Myers <joseph@codesourcery.com>
* sysdeps/unix/sysv/linux/mips/mips64/nptl/sysdep-cancel.h * sysdeps/unix/sysv/linux/mips/mips64/nptl/sysdep-cancel.h

20
NEWS
View File

@ -12,16 +12,16 @@ Version 2.20
6804, 9894, 12994, 13347, 13651, 14308, 14770, 15119, 15132, 15347, 15514, 6804, 9894, 12994, 13347, 13651, 14308, 14770, 15119, 15132, 15347, 15514,
15698, 15804, 15894, 15946, 16002, 16064, 16095, 16198, 16284, 16287, 15698, 15804, 15894, 15946, 16002, 16064, 16095, 16198, 16284, 16287,
16315, 16348, 16349, 16354, 16357, 16362, 16447, 16516, 16532, 16539, 16315, 16348, 16349, 16354, 16357, 16362, 16447, 16516, 16532, 16539,
16545, 16564, 16574, 16599, 16600, 16609, 16610, 16611, 16613, 16619, 16545, 16561, 16562, 16564, 16574, 16599, 16600, 16609, 16610, 16611,
16623, 16629, 16632, 16634, 16639, 16642, 16648, 16649, 16670, 16674, 16613, 16619, 16623, 16629, 16632, 16634, 16639, 16642, 16648, 16649,
16677, 16680, 16681, 16683, 16689, 16695, 16701, 16706, 16707, 16712, 16670, 16674, 16677, 16680, 16681, 16683, 16689, 16695, 16701, 16706,
16713, 16714, 16724, 16731, 16739, 16740, 16743, 16754, 16758, 16759, 16707, 16712, 16713, 16714, 16724, 16731, 16739, 16740, 16743, 16754,
16760, 16770, 16786, 16789, 16791, 16796, 16799, 16800, 16815, 16823, 16758, 16759, 16760, 16770, 16786, 16789, 16791, 16796, 16799, 16800,
16824, 16831, 16838, 16849, 16854, 16876, 16877, 16878, 16882, 16885, 16815, 16823, 16824, 16831, 16838, 16849, 16854, 16876, 16877, 16878,
16888, 16890, 16912, 16915, 16916, 16917, 16918, 16922, 16927, 16928, 16882, 16885, 16888, 16890, 16912, 16915, 16916, 16917, 16918, 16922,
16932, 16943, 16958, 16965, 16966, 16967, 16977, 16978, 16984, 16990, 16927, 16928, 16932, 16943, 16958, 16965, 16966, 16967, 16977, 16978,
16996, 17009, 17022, 17031, 17042, 17048, 17050, 17058, 17061, 17062, 16984, 16990, 16996, 17009, 17022, 17031, 17042, 17048, 17050, 17058,
17069, 17075, 17079, 17084, 17086. 17061, 17062, 17069, 17075, 17079, 17084, 17086.
* Optimized strchr implementation for AArch64. Contributed by ARM Ltd. * Optimized strchr implementation for AArch64. Contributed by ARM Ltd.

4
math/libm-test.inc
View File

@ -9609,9 +9609,7 @@ static const struct test_if_f_data yn_test_data[] =
static void static void
yn_test (void) yn_test (void)
{ {
START (yn, 0); ALL_RM_TEST (yn, 0, yn_test_data, RUN_TEST_LOOP_if_f, END);
RUN_TEST_LOOP_if_f (yn, yn_test_data, );
END;
} }

1
sysdeps/i386/fpu/fenv_private.h
View File

@ -482,6 +482,7 @@ libc_feupdateenv_387_ctx (struct rm_ctx *ctx)
#else #else
# define libc_feholdexcept_setround_ctx libc_feholdexcept_setround_387_ctx # define libc_feholdexcept_setround_ctx libc_feholdexcept_setround_387_ctx
# define libc_feupdateenv_ctx libc_feupdateenv_387_ctx # define libc_feupdateenv_ctx libc_feupdateenv_387_ctx
# define libc_feholdsetround_ctx libc_feholdsetround_387_ctx
# define libc_feresetround_ctx libc_feresetround_387_ctx # define libc_feresetround_ctx libc_feresetround_387_ctx
#endif /* __SSE2_MATH__ */ #endif /* __SSE2_MATH__ */

24
sysdeps/i386/fpu/libm-test-ulps
View File

@ -1931,4 +1931,28 @@ ifloat: 3
ildouble: 4 ildouble: 4
ldouble: 4 ldouble: 4
Function: "yn_downward":
double: 2
float: 2
idouble: 2
ifloat: 2
ildouble: 5
ldouble: 5
Function: "yn_towardzero":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 5
ldouble: 5
Function: "yn_upward":
double: 2
float: 3
idouble: 2
ifloat: 3
ildouble: 4
ldouble: 4
# end of automatic generation # end of automatic generation

117
sysdeps/ieee754/dbl-64/e_jn.c
View File

@ -37,6 +37,7 @@
*/ */
#include <errno.h> #include <errno.h>
#include <float.h>
#include <math.h> #include <math.h>
#include <math_private.h> #include <math_private.h>
@ -248,7 +249,7 @@ __ieee754_yn (int n, double x)
{ {
int32_t i, hx, ix, lx; int32_t i, hx, ix, lx;
int32_t sign; int32_t sign;
double a, b, temp; double a, b, temp, ret;
EXTRACT_WORDS (hx, lx, x); EXTRACT_WORDS (hx, lx, x);
ix = 0x7fffffff & hx; ix = 0x7fffffff & hx;
@ -268,57 +269,67 @@ __ieee754_yn (int n, double x)
} }
if (n == 0) if (n == 0)
return (__ieee754_y0 (x)); return (__ieee754_y0 (x));
if (n == 1) {
return (sign * __ieee754_y1 (x)); SET_RESTORE_ROUND (FE_TONEAREST);
if (__glibc_unlikely (ix == 0x7ff00000)) if (n == 1)
return zero; {
if (ix >= 0x52D00000) /* x > 2**302 */ ret = sign * __ieee754_y1 (x);
{ /* (x >> n**2) goto out;
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) }
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) if (__glibc_unlikely (ix == 0x7ff00000))
* Let s=sin(x), c=cos(x), return zero;
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then if (ix >= 0x52D00000) /* x > 2**302 */
* { /* (x >> n**2)
* n sin(xn)*sqt2 cos(xn)*sqt2 * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* ---------------------------------- * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* 0 s-c c+s * Let s=sin(x), c=cos(x),
* 1 -s-c -c+s * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
* 2 -s+c -c-s *
* 3 s+c c-s * n sin(xn)*sqt2 cos(xn)*sqt2
*/ * ----------------------------------
double c; * 0 s-c c+s
double s; * 1 -s-c -c+s
__sincos (x, &s, &c); * 2 -s+c -c-s
switch (n & 3) * 3 s+c c-s
{ */
case 0: temp = s - c; break; double c;
case 1: temp = -s - c; break; double s;
case 2: temp = -s + c; break; __sincos (x, &s, &c);
case 3: temp = s + c; break; switch (n & 3)
} {
b = invsqrtpi * temp / __ieee754_sqrt (x); case 0: temp = s - c; break;
} case 1: temp = -s - c; break;
else case 2: temp = -s + c; break;
{ case 3: temp = s + c; break;
u_int32_t high; }
a = __ieee754_y0 (x); b = invsqrtpi * temp / __ieee754_sqrt (x);
b = __ieee754_y1 (x); }
/* quit if b is -inf */ else
GET_HIGH_WORD (high, b); {
for (i = 1; i < n && high != 0xfff00000; i++) u_int32_t high;
{ a = __ieee754_y0 (x);
temp = b; b = __ieee754_y1 (x);
b = ((double) (i + i) / x) * b - a; /* quit if b is -inf */
GET_HIGH_WORD (high, b); GET_HIGH_WORD (high, b);
a = temp; for (i = 1; i < n && high != 0xfff00000; i++)
} {
/* If B is +-Inf, set up errno accordingly. */ temp = b;
if (!__finite (b)) b = ((double) (i + i) / x) * b - a;
__set_errno (ERANGE); GET_HIGH_WORD (high, b);
} a = temp;
if (sign > 0) }
return b; /* If B is +-Inf, set up errno accordingly. */
else if (!__finite (b))
return -b; __set_errno (ERANGE);
}
if (sign > 0)
ret = b;
else
ret = -b;
}
out:
if (__isinf (ret))
ret = __copysign (DBL_MAX, ret) * DBL_MAX;
return ret;
} }
strong_alias (__ieee754_yn, __yn_finite) strong_alias (__ieee754_yn, __yn_finite)

16
sysdeps/ieee754/flt-32/e_jnf.c
View File

@ -14,6 +14,7 @@
*/ */
#include <errno.h> #include <errno.h>
#include <float.h>
#include <math.h> #include <math.h>
#include <math_private.h> #include <math_private.h>
@ -169,6 +170,8 @@ strong_alias (__ieee754_jnf, __jnf_finite)
float float
__ieee754_ynf(int n, float x) __ieee754_ynf(int n, float x)
{ {
float ret;
{
int32_t i,hx,ix; int32_t i,hx,ix;
u_int32_t ib; u_int32_t ib;
int32_t sign; int32_t sign;
@ -187,7 +190,11 @@ __ieee754_ynf(int n, float x)
sign = 1 - ((n&1)<<1); sign = 1 - ((n&1)<<1);
} }
if(n==0) return(__ieee754_y0f(x)); if(n==0) return(__ieee754_y0f(x));
if(n==1) return(sign*__ieee754_y1f(x)); SET_RESTORE_ROUNDF (FE_TONEAREST);
if(n==1) {
ret = sign*__ieee754_y1f(x);
goto out;
}
if(__builtin_expect(ix==0x7f800000, 0)) return zero; if(__builtin_expect(ix==0x7f800000, 0)) return zero;
a = __ieee754_y0f(x); a = __ieee754_y0f(x);
@ -203,6 +210,11 @@ __ieee754_ynf(int n, float x)
/* If B is +-Inf, set up errno accordingly. */ /* If B is +-Inf, set up errno accordingly. */
if (! __finitef (b)) if (! __finitef (b))
__set_errno (ERANGE); __set_errno (ERANGE);
if(sign>0) return b; else return -b; if(sign>0) ret = b; else ret = -b;
}
out:
if (__isinff (ret))
ret = __copysignf (FLT_MAX, ret) * FLT_MAX;
return ret;
} }
strong_alias (__ieee754_ynf, __ynf_finite) strong_alias (__ieee754_ynf, __ynf_finite)

139
sysdeps/ieee754/ldbl-128/e_jnl.c
View File

@ -57,6 +57,7 @@
*/ */
#include <errno.h> #include <errno.h>
#include <float.h>
#include <math.h> #include <math.h>
#include <math_private.h> #include <math_private.h>
@ -300,7 +301,7 @@ __ieee754_ynl (int n, long double x)
u_int32_t se; u_int32_t se;
int32_t i, ix; int32_t i, ix;
int32_t sign; int32_t sign;
long double a, b, temp; long double a, b, temp, ret;
ieee854_long_double_shape_type u; ieee854_long_double_shape_type u;
u.value = x; u.value = x;
@ -328,70 +329,80 @@ __ieee754_ynl (int n, long double x)
} }
if (n == 0) if (n == 0)
return (__ieee754_y0l (x)); return (__ieee754_y0l (x));
if (n == 1) {
return (sign * __ieee754_y1l (x)); SET_RESTORE_ROUNDL (FE_TONEAREST);
if (ix >= 0x7fff0000) if (n == 1)
return zero; {
if (ix >= 0x412D0000) ret = sign * __ieee754_y1l (x);
{ /* x > 2**302 */ goto out;
}
if (ix >= 0x7fff0000)
return zero;
if (ix >= 0x412D0000)
{ /* x > 2**302 */
/* ??? See comment above on the possible futility of this. */ /* ??? See comment above on the possible futility of this. */
/* (x >> n**2) /* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x), * Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
* *
* n sin(xn)*sqt2 cos(xn)*sqt2 * n sin(xn)*sqt2 cos(xn)*sqt2
* ---------------------------------- * ----------------------------------
* 0 s-c c+s * 0 s-c c+s
* 1 -s-c -c+s * 1 -s-c -c+s
* 2 -s+c -c-s * 2 -s+c -c-s
* 3 s+c c-s * 3 s+c c-s
*/ */
long double s; long double s;
long double c; long double c;
__sincosl (x, &s, &c); __sincosl (x, &s, &c);
switch (n & 3) switch (n & 3)
{ {
case 0: case 0:
temp = s - c; temp = s - c;
break; break;
case 1: case 1:
temp = -s - c; temp = -s - c;
break; break;
case 2: case 2:
temp = -s + c; temp = -s + c;
break; break;
case 3: case 3:
temp = s + c; temp = s + c;
break; break;
} }
b = invsqrtpi * temp / __ieee754_sqrtl (x); b = invsqrtpi * temp / __ieee754_sqrtl (x);
} }
else else
{ {
a = __ieee754_y0l (x); a = __ieee754_y0l (x);
b = __ieee754_y1l (x); b = __ieee754_y1l (x);
/* quit if b is -inf */ /* quit if b is -inf */
u.value = b; u.value = b;
se = u.parts32.w0 & 0xffff0000; se = u.parts32.w0 & 0xffff0000;
for (i = 1; i < n && se != 0xffff0000; i++) for (i = 1; i < n && se != 0xffff0000; i++)
{ {
temp = b; temp = b;
b = ((long double) (i + i) / x) * b - a; b = ((long double) (i + i) / x) * b - a;
u.value = b; u.value = b;
se = u.parts32.w0 & 0xffff0000; se = u.parts32.w0 & 0xffff0000;
a = temp; a = temp;
} }
} }
/* If B is +-Inf, set up errno accordingly. */ /* If B is +-Inf, set up errno accordingly. */
if (! __finitel (b)) if (! __finitel (b))
__set_errno (ERANGE); __set_errno (ERANGE);
if (sign > 0) if (sign > 0)
return b; ret = b;
else else
return -b; ret = -b;
}
out:
if (__isinfl (ret))
ret = __copysignl (LDBL_MAX, ret) * LDBL_MAX;
return ret;
} }
strong_alias (__ieee754_ynl, __ynl_finite) strong_alias (__ieee754_ynl, __ynl_finite)

143
sysdeps/ieee754/ldbl-128ibm/e_jnl.c
View File

@ -57,6 +57,7 @@
*/ */
#include <errno.h> #include <errno.h>
#include <float.h>
#include <math.h> #include <math.h>
#include <math_private.h> #include <math_private.h>
@ -300,7 +301,7 @@ __ieee754_ynl (int n, long double x)
uint32_t se, lx; uint32_t se, lx;
int32_t i, ix; int32_t i, ix;
int32_t sign; int32_t sign;
long double a, b, temp; long double a, b, temp, ret;
double xhi; double xhi;
xhi = ldbl_high (x); xhi = ldbl_high (x);
@ -328,72 +329,82 @@ __ieee754_ynl (int n, long double x)
} }
if (n == 0) if (n == 0)
return (__ieee754_y0l (x)); return (__ieee754_y0l (x));
if (n == 1) {
return (sign * __ieee754_y1l (x)); SET_RESTORE_ROUNDL (FE_TONEAREST);
if (ix >= 0x7ff00000) if (n == 1)
return zero; {
if (ix >= 0x52D00000) ret = sign * __ieee754_y1l (x);
{ /* x > 2**302 */ goto out;
}
if (ix >= 0x7ff00000)
return zero;
if (ix >= 0x52D00000)
{ /* x > 2**302 */
/* ??? See comment above on the possible futility of this. */ /* ??? See comment above on the possible futility of this. */
/* (x >> n**2) /* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x), * Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
* *
* n sin(xn)*sqt2 cos(xn)*sqt2 * n sin(xn)*sqt2 cos(xn)*sqt2
* ---------------------------------- * ----------------------------------
* 0 s-c c+s * 0 s-c c+s
* 1 -s-c -c+s * 1 -s-c -c+s
* 2 -s+c -c-s * 2 -s+c -c-s
* 3 s+c c-s * 3 s+c c-s
*/ */
long double s; long double s;
long double c; long double c;
__sincosl (x, &s, &c); __sincosl (x, &s, &c);
switch (n & 3) switch (n & 3)
{ {
case 0: case 0:
temp = s - c; temp = s - c;
break; break;
case 1: case 1:
temp = -s - c; temp = -s - c;
break; break;
case 2: case 2:
temp = -s + c; temp = -s + c;
break; break;
case 3: case 3:
temp = s + c; temp = s + c;
break; break;
} }
b = invsqrtpi * temp / __ieee754_sqrtl (x); b = invsqrtpi * temp / __ieee754_sqrtl (x);
} }
else else
{ {
a = __ieee754_y0l (x); a = __ieee754_y0l (x);
b = __ieee754_y1l (x); b = __ieee754_y1l (x);
/* quit if b is -inf */ /* quit if b is -inf */
xhi = ldbl_high (b); xhi = ldbl_high (b);
GET_HIGH_WORD (se, xhi); GET_HIGH_WORD (se, xhi);
se &= 0xfff00000; se &= 0xfff00000;
for (i = 1; i < n && se != 0xfff00000; i++) for (i = 1; i < n && se != 0xfff00000; i++)
{ {
temp = b; temp = b;
b = ((long double) (i + i) / x) * b - a; b = ((long double) (i + i) / x) * b - a;
xhi = ldbl_high (b); xhi = ldbl_high (b);
GET_HIGH_WORD (se, xhi); GET_HIGH_WORD (se, xhi);
se &= 0xfff00000; se &= 0xfff00000;
a = temp; a = temp;
} }
} }
/* If B is +-Inf, set up errno accordingly. */ /* If B is +-Inf, set up errno accordingly. */
if (! __finitel (b)) if (! __finitel (b))
__set_errno (ERANGE); __set_errno (ERANGE);
if (sign > 0) if (sign > 0)
return b; ret = b;
else else
return -b; ret = -b;
}
out:
if (__isinfl (ret))
ret = __copysignl (LDBL_MAX, ret) * LDBL_MAX;
return ret;
} }
strong_alias (__ieee754_ynl, __ynl_finite) strong_alias (__ieee754_ynl, __ynl_finite)

137
sysdeps/ieee754/ldbl-96/e_jnl.c
View File

@ -57,6 +57,7 @@
*/ */
#include <errno.h> #include <errno.h>
#include <float.h>
#include <math.h> #include <math.h>
#include <math_private.h> #include <math_private.h>
@ -293,7 +294,7 @@ __ieee754_ynl (int n, long double x)
u_int32_t se, i0, i1; u_int32_t se, i0, i1;
int32_t i, ix; int32_t i, ix;
int32_t sign; int32_t sign;
long double a, b, temp; long double a, b, temp, ret;
GET_LDOUBLE_WORDS (se, i0, i1, x); GET_LDOUBLE_WORDS (se, i0, i1, x);
@ -314,69 +315,79 @@ __ieee754_ynl (int n, long double x)
} }
if (n == 0) if (n == 0)
return (__ieee754_y0l (x)); return (__ieee754_y0l (x));
if (n == 1) {
return (sign * __ieee754_y1l (x)); SET_RESTORE_ROUNDL (FE_TONEAREST);
if (__glibc_unlikely (ix == 0x7fff)) if (n == 1)
return zero; {
if (ix >= 0x412D) ret = sign * __ieee754_y1l (x);
{ /* x > 2**302 */ goto out;
}
if (__glibc_unlikely (ix == 0x7fff))
return zero;
if (ix >= 0x412D)
{ /* x > 2**302 */
/* ??? See comment above on the possible futility of this. */ /* ??? See comment above on the possible futility of this. */
/* (x >> n**2) /* (x >> n**2)
* Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
* Let s=sin(x), c=cos(x), * Let s=sin(x), c=cos(x),
* xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
* *
* n sin(xn)*sqt2 cos(xn)*sqt2 * n sin(xn)*sqt2 cos(xn)*sqt2
* ---------------------------------- * ----------------------------------
* 0 s-c c+s * 0 s-c c+s
* 1 -s-c -c+s * 1 -s-c -c+s
* 2 -s+c -c-s * 2 -s+c -c-s
* 3 s+c c-s * 3 s+c c-s
*/ */
long double s; long double s;
long double c; long double c;
__sincosl (x, &s, &c); __sincosl (x, &s, &c);
switch (n & 3) switch (n & 3)
{ {
case 0: case 0:
temp = s - c; temp = s - c;
break; break;
case 1: case 1:
temp = -s - c; temp = -s - c;
break; break;
case 2: case 2:
temp = -s + c; temp = -s + c;
break; break;
case 3: case 3:
temp = s + c; temp = s + c;
break; break;
} }
b = invsqrtpi * temp / __ieee754_sqrtl (x); b = invsqrtpi * temp / __ieee754_sqrtl (x);
} }
else else
{ {
a = __ieee754_y0l (x); a = __ieee754_y0l (x);
b = __ieee754_y1l (x); b = __ieee754_y1l (x);
/* quit if b is -inf */ /* quit if b is -inf */
GET_LDOUBLE_WORDS (se, i0, i1, b); GET_LDOUBLE_WORDS (se, i0, i1, b);
/* Use 0xffffffff since GET_LDOUBLE_WORDS sign-extends SE. */ /* Use 0xffffffff since GET_LDOUBLE_WORDS sign-extends SE. */
for (i = 1; i < n && se != 0xffffffff; i++) for (i = 1; i < n && se != 0xffffffff; i++)
{ {
temp = b; temp = b;
b = ((long double) (i + i) / x) * b - a; b = ((long double) (i + i) / x) * b - a;
GET_LDOUBLE_WORDS (se, i0, i1, b); GET_LDOUBLE_WORDS (se, i0, i1, b);
a = temp; a = temp;
} }
} }
/* If B is +-Inf, set up errno accordingly. */ /* If B is +-Inf, set up errno accordingly. */
if (! __finitel (b)) if (! __finitel (b))
__set_errno (ERANGE); __set_errno (ERANGE);
if (sign > 0) if (sign > 0)
return b; ret = b;
else else
return -b; ret = -b;
}
out:
if (__isinfl (ret))
ret = __copysignl (LDBL_MAX, ret) * LDBL_MAX;
return ret;
} }
strong_alias (__ieee754_ynl, __ynl_finite) strong_alias (__ieee754_ynl, __ynl_finite)

24
sysdeps/x86_64/fpu/libm-test-ulps
View File

@ -2033,4 +2033,28 @@ ifloat: 3
ildouble: 4 ildouble: 4
ldouble: 4 ldouble: 4
Function: "yn_downward":
double: 3
float: 4
idouble: 3
ifloat: 4
ildouble: 5
ldouble: 5
Function: "yn_towardzero":
double: 3
float: 3
idouble: 3
ifloat: 3
ildouble: 5
ldouble: 5
Function: "yn_upward":
double: 4
float: 5
idouble: 4
ifloat: 5
ildouble: 4
ldouble: 4
# end of automatic generation # end of automatic generation