diff options
Diffstat (limited to 'sysdeps')
| -rw-r--r-- | sysdeps/i386/fpu/fenv_private.h | 1 | ||||
| -rw-r--r-- | sysdeps/i386/fpu/libm-test-ulps | 24 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_jn.c | 117 | ||||
| -rw-r--r-- | sysdeps/ieee754/flt-32/e_jnf.c | 16 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128/e_jnl.c | 139 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-128ibm/e_jnl.c | 143 | ||||
| -rw-r--r-- | sysdeps/ieee754/ldbl-96/e_jnl.c | 137 | ||||
| -rw-r--r-- | sysdeps/x86_64/fpu/libm-test-ulps | 24 |
8 files changed, 353 insertions, 248 deletions
diff --git a/sysdeps/i386/fpu/fenv_private.h b/sysdeps/i386/fpu/fenv_private.h index 3998387c57..e20e1f1662 100644 --- a/sysdeps/i386/fpu/fenv_private.h +++ b/sysdeps/i386/fpu/fenv_private.h @@ -482,6 +482,7 @@ libc_feupdateenv_387_ctx (struct rm_ctx *ctx) #else # define libc_feholdexcept_setround_ctx libc_feholdexcept_setround_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 #endif /* __SSE2_MATH__ */ diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps index 7cf958608d..be0e7c3668 100644 --- a/sysdeps/i386/fpu/libm-test-ulps +++ b/sysdeps/i386/fpu/libm-test-ulps @@ -1931,4 +1931,28 @@ ifloat: 3 ildouble: 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 diff --git a/sysdeps/ieee754/dbl-64/e_jn.c b/sysdeps/ieee754/dbl-64/e_jn.c index 236878ba1d..bbe04260b8 100644 --- a/sysdeps/ieee754/dbl-64/e_jn.c +++ b/sysdeps/ieee754/dbl-64/e_jn.c @@ -37,6 +37,7 @@ */ #include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> @@ -248,7 +249,7 @@ __ieee754_yn (int n, double x) { int32_t i, hx, ix, lx; int32_t sign; - double a, b, temp; + double a, b, temp, ret; EXTRACT_WORDS (hx, lx, x); ix = 0x7fffffff & hx; @@ -268,57 +269,67 @@ __ieee754_yn (int n, double x) } if (n == 0) return (__ieee754_y0 (x)); - if (n == 1) - return (sign * __ieee754_y1 (x)); - if (__glibc_unlikely (ix == 0x7ff00000)) - return zero; - if (ix >= 0x52D00000) /* x > 2**302 */ - { /* (x >> n**2) - * Jn(x) = cos(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), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then - * - * n sin(xn)*sqt2 cos(xn)*sqt2 - * ---------------------------------- - * 0 s-c c+s - * 1 -s-c -c+s - * 2 -s+c -c-s - * 3 s+c c-s - */ - double c; - double s; - __sincos (x, &s, &c); - switch (n & 3) - { - case 0: temp = s - c; break; - case 1: temp = -s - c; break; - case 2: temp = -s + c; break; - case 3: temp = s + c; break; - } - b = invsqrtpi * temp / __ieee754_sqrt (x); - } - else - { - u_int32_t high; - a = __ieee754_y0 (x); - b = __ieee754_y1 (x); - /* quit if b is -inf */ - GET_HIGH_WORD (high, b); - for (i = 1; i < n && high != 0xfff00000; i++) - { - temp = b; - b = ((double) (i + i) / x) * b - a; - GET_HIGH_WORD (high, b); - a = temp; - } - /* If B is +-Inf, set up errno accordingly. */ - if (!__finite (b)) - __set_errno (ERANGE); - } - if (sign > 0) - return b; - else - return -b; + { + SET_RESTORE_ROUND (FE_TONEAREST); + if (n == 1) + { + ret = sign * __ieee754_y1 (x); + goto out; + } + if (__glibc_unlikely (ix == 0x7ff00000)) + return zero; + if (ix >= 0x52D00000) /* x > 2**302 */ + { /* (x >> n**2) + * Jn(x) = cos(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), + * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then + * + * n sin(xn)*sqt2 cos(xn)*sqt2 + * ---------------------------------- + * 0 s-c c+s + * 1 -s-c -c+s + * 2 -s+c -c-s + * 3 s+c c-s + */ + double c; + double s; + __sincos (x, &s, &c); + switch (n & 3) + { + case 0: temp = s - c; break; + case 1: temp = -s - c; break; + case 2: temp = -s + c; break; + case 3: temp = s + c; break; + } + b = invsqrtpi * temp / __ieee754_sqrt (x); + } + else + { + u_int32_t high; + a = __ieee754_y0 (x); + b = __ieee754_y1 (x); + /* quit if b is -inf */ + GET_HIGH_WORD (high, b); + for (i = 1; i < n && high != 0xfff00000; i++) + { + temp = b; + b = ((double) (i + i) / x) * b - a; + GET_HIGH_WORD (high, b); + a = temp; + } + /* If B is +-Inf, set up errno accordingly. */ + if (!__finite (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) diff --git a/sysdeps/ieee754/flt-32/e_jnf.c b/sysdeps/ieee754/flt-32/e_jnf.c index 5984d94a3c..86085cc635 100644 --- a/sysdeps/ieee754/flt-32/e_jnf.c +++ b/sysdeps/ieee754/flt-32/e_jnf.c @@ -14,6 +14,7 @@ */ #include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> @@ -169,6 +170,8 @@ strong_alias (__ieee754_jnf, __jnf_finite) float __ieee754_ynf(int n, float x) { + float ret; + { int32_t i,hx,ix; u_int32_t ib; int32_t sign; @@ -187,7 +190,11 @@ __ieee754_ynf(int n, float x) sign = 1 - ((n&1)<<1); } 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; a = __ieee754_y0f(x); @@ -203,6 +210,11 @@ __ieee754_ynf(int n, float x) /* If B is +-Inf, set up errno accordingly. */ if (! __finitef (b)) __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) diff --git a/sysdeps/ieee754/ldbl-128/e_jnl.c b/sysdeps/ieee754/ldbl-128/e_jnl.c index c2a49235c3..4e32d38581 100644 --- a/sysdeps/ieee754/ldbl-128/e_jnl.c +++ b/sysdeps/ieee754/ldbl-128/e_jnl.c @@ -57,6 +57,7 @@ */ #include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> @@ -300,7 +301,7 @@ __ieee754_ynl (int n, long double x) u_int32_t se; int32_t i, ix; int32_t sign; - long double a, b, temp; + long double a, b, temp, ret; ieee854_long_double_shape_type u; u.value = x; @@ -328,70 +329,80 @@ __ieee754_ynl (int n, long double x) } if (n == 0) return (__ieee754_y0l (x)); - if (n == 1) - return (sign * __ieee754_y1l (x)); - if (ix >= 0x7fff0000) - return zero; - if (ix >= 0x412D0000) - { /* x > 2**302 */ + { + SET_RESTORE_ROUNDL (FE_TONEAREST); + if (n == 1) + { + ret = sign * __ieee754_y1l (x); + 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) - * Jn(x) = cos(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), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then - * - * n sin(xn)*sqt2 cos(xn)*sqt2 - * ---------------------------------- - * 0 s-c c+s - * 1 -s-c -c+s - * 2 -s+c -c-s - * 3 s+c c-s - */ - long double s; - long double c; - __sincosl (x, &s, &c); - switch (n & 3) - { - case 0: - temp = s - c; - break; - case 1: - temp = -s - c; - break; - case 2: - temp = -s + c; - break; - case 3: - temp = s + c; - break; - } - b = invsqrtpi * temp / __ieee754_sqrtl (x); - } - else - { - a = __ieee754_y0l (x); - b = __ieee754_y1l (x); - /* quit if b is -inf */ - u.value = b; - se = u.parts32.w0 & 0xffff0000; - for (i = 1; i < n && se != 0xffff0000; i++) - { - temp = b; - b = ((long double) (i + i) / x) * b - a; - u.value = b; - se = u.parts32.w0 & 0xffff0000; - a = temp; - } - } - /* If B is +-Inf, set up errno accordingly. */ - if (! __finitel (b)) - __set_errno (ERANGE); - if (sign > 0) - return b; - else - return -b; + /* (x >> n**2) + * Jn(x) = cos(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), + * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then + * + * n sin(xn)*sqt2 cos(xn)*sqt2 + * ---------------------------------- + * 0 s-c c+s + * 1 -s-c -c+s + * 2 -s+c -c-s + * 3 s+c c-s + */ + long double s; + long double c; + __sincosl (x, &s, &c); + switch (n & 3) + { + case 0: + temp = s - c; + break; + case 1: + temp = -s - c; + break; + case 2: + temp = -s + c; + break; + case 3: + temp = s + c; + break; + } + b = invsqrtpi * temp / __ieee754_sqrtl (x); + } + else + { + a = __ieee754_y0l (x); + b = __ieee754_y1l (x); + /* quit if b is -inf */ + u.value = b; + se = u.parts32.w0 & 0xffff0000; + for (i = 1; i < n && se != 0xffff0000; i++) + { + temp = b; + b = ((long double) (i + i) / x) * b - a; + u.value = b; + se = u.parts32.w0 & 0xffff0000; + a = temp; + } + } + /* If B is +-Inf, set up errno accordingly. */ + if (! __finitel (b)) + __set_errno (ERANGE); + if (sign > 0) + ret = b; + else + ret = -b; + } + out: + if (__isinfl (ret)) + ret = __copysignl (LDBL_MAX, ret) * LDBL_MAX; + return ret; } strong_alias (__ieee754_ynl, __ynl_finite) diff --git a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c index 6761a0d26f..589f1f822a 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_jnl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_jnl.c @@ -57,6 +57,7 @@ */ #include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> @@ -300,7 +301,7 @@ __ieee754_ynl (int n, long double x) uint32_t se, lx; int32_t i, ix; int32_t sign; - long double a, b, temp; + long double a, b, temp, ret; double xhi; xhi = ldbl_high (x); @@ -328,72 +329,82 @@ __ieee754_ynl (int n, long double x) } if (n == 0) return (__ieee754_y0l (x)); - if (n == 1) - return (sign * __ieee754_y1l (x)); - if (ix >= 0x7ff00000) - return zero; - if (ix >= 0x52D00000) - { /* x > 2**302 */ + { + SET_RESTORE_ROUNDL (FE_TONEAREST); + if (n == 1) + { + ret = sign * __ieee754_y1l (x); + 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) - * Jn(x) = cos(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), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then - * - * n sin(xn)*sqt2 cos(xn)*sqt2 - * ---------------------------------- - * 0 s-c c+s - * 1 -s-c -c+s - * 2 -s+c -c-s - * 3 s+c c-s - */ - long double s; - long double c; - __sincosl (x, &s, &c); - switch (n & 3) - { - case 0: - temp = s - c; - break; - case 1: - temp = -s - c; - break; - case 2: - temp = -s + c; - break; - case 3: - temp = s + c; - break; - } - b = invsqrtpi * temp / __ieee754_sqrtl (x); - } - else - { - a = __ieee754_y0l (x); - b = __ieee754_y1l (x); - /* quit if b is -inf */ - xhi = ldbl_high (b); - GET_HIGH_WORD (se, xhi); - se &= 0xfff00000; - for (i = 1; i < n && se != 0xfff00000; i++) - { - temp = b; - b = ((long double) (i + i) / x) * b - a; - xhi = ldbl_high (b); - GET_HIGH_WORD (se, xhi); - se &= 0xfff00000; - a = temp; - } - } - /* If B is +-Inf, set up errno accordingly. */ - if (! __finitel (b)) - __set_errno (ERANGE); - if (sign > 0) - return b; - else - return -b; + /* (x >> n**2) + * Jn(x) = cos(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), + * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then + * + * n sin(xn)*sqt2 cos(xn)*sqt2 + * ---------------------------------- + * 0 s-c c+s + * 1 -s-c -c+s + * 2 -s+c -c-s + * 3 s+c c-s + */ + long double s; + long double c; + __sincosl (x, &s, &c); + switch (n & 3) + { + case 0: + temp = s - c; + break; + case 1: + temp = -s - c; + break; + case 2: + temp = -s + c; + break; + case 3: + temp = s + c; + break; + } + b = invsqrtpi * temp / __ieee754_sqrtl (x); + } + else + { + a = __ieee754_y0l (x); + b = __ieee754_y1l (x); + /* quit if b is -inf */ + xhi = ldbl_high (b); + GET_HIGH_WORD (se, xhi); + se &= 0xfff00000; + for (i = 1; i < n && se != 0xfff00000; i++) + { + temp = b; + b = ((long double) (i + i) / x) * b - a; + xhi = ldbl_high (b); + GET_HIGH_WORD (se, xhi); + se &= 0xfff00000; + a = temp; + } + } + /* If B is +-Inf, set up errno accordingly. */ + if (! __finitel (b)) + __set_errno (ERANGE); + if (sign > 0) + ret = b; + else + ret = -b; + } + out: + if (__isinfl (ret)) + ret = __copysignl (LDBL_MAX, ret) * LDBL_MAX; + return ret; } strong_alias (__ieee754_ynl, __ynl_finite) diff --git a/sysdeps/ieee754/ldbl-96/e_jnl.c b/sysdeps/ieee754/ldbl-96/e_jnl.c index 11d097c271..95ff24201b 100644 --- a/sysdeps/ieee754/ldbl-96/e_jnl.c +++ b/sysdeps/ieee754/ldbl-96/e_jnl.c @@ -57,6 +57,7 @@ */ #include <errno.h> +#include <float.h> #include <math.h> #include <math_private.h> @@ -293,7 +294,7 @@ __ieee754_ynl (int n, long double x) u_int32_t se, i0, i1; int32_t i, ix; int32_t sign; - long double a, b, temp; + long double a, b, temp, ret; GET_LDOUBLE_WORDS (se, i0, i1, x); @@ -314,69 +315,79 @@ __ieee754_ynl (int n, long double x) } if (n == 0) return (__ieee754_y0l (x)); - if (n == 1) - return (sign * __ieee754_y1l (x)); - if (__glibc_unlikely (ix == 0x7fff)) - return zero; - if (ix >= 0x412D) - { /* x > 2**302 */ + { + SET_RESTORE_ROUNDL (FE_TONEAREST); + if (n == 1) + { + ret = sign * __ieee754_y1l (x); + 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) - * Jn(x) = cos(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), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then - * - * n sin(xn)*sqt2 cos(xn)*sqt2 - * ---------------------------------- - * 0 s-c c+s - * 1 -s-c -c+s - * 2 -s+c -c-s - * 3 s+c c-s - */ - long double s; - long double c; - __sincosl (x, &s, &c); - switch (n & 3) - { - case 0: - temp = s - c; - break; - case 1: - temp = -s - c; - break; - case 2: - temp = -s + c; - break; - case 3: - temp = s + c; - break; - } - b = invsqrtpi * temp / __ieee754_sqrtl (x); - } - else - { - a = __ieee754_y0l (x); - b = __ieee754_y1l (x); - /* quit if b is -inf */ - GET_LDOUBLE_WORDS (se, i0, i1, b); - /* Use 0xffffffff since GET_LDOUBLE_WORDS sign-extends SE. */ - for (i = 1; i < n && se != 0xffffffff; i++) - { - temp = b; - b = ((long double) (i + i) / x) * b - a; - GET_LDOUBLE_WORDS (se, i0, i1, b); - a = temp; - } - } - /* If B is +-Inf, set up errno accordingly. */ - if (! __finitel (b)) - __set_errno (ERANGE); - if (sign > 0) - return b; - else - return -b; + /* (x >> n**2) + * Jn(x) = cos(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), + * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then + * + * n sin(xn)*sqt2 cos(xn)*sqt2 + * ---------------------------------- + * 0 s-c c+s + * 1 -s-c -c+s + * 2 -s+c -c-s + * 3 s+c c-s + */ + long double s; + long double c; + __sincosl (x, &s, &c); + switch (n & 3) + { + case 0: + temp = s - c; + break; + case 1: + temp = -s - c; + break; + case 2: + temp = -s + c; + break; + case 3: + temp = s + c; + break; + } + b = invsqrtpi * temp / __ieee754_sqrtl (x); + } + else + { + a = __ieee754_y0l (x); + b = __ieee754_y1l (x); + /* quit if b is -inf */ + GET_LDOUBLE_WORDS (se, i0, i1, b); + /* Use 0xffffffff since GET_LDOUBLE_WORDS sign-extends SE. */ + for (i = 1; i < n && se != 0xffffffff; i++) + { + temp = b; + b = ((long double) (i + i) / x) * b - a; + GET_LDOUBLE_WORDS (se, i0, i1, b); + a = temp; + } + } + /* If B is +-Inf, set up errno accordingly. */ + if (! __finitel (b)) + __set_errno (ERANGE); + if (sign > 0) + ret = b; + else + ret = -b; + } + out: + if (__isinfl (ret)) + ret = __copysignl (LDBL_MAX, ret) * LDBL_MAX; + return ret; } strong_alias (__ieee754_ynl, __ynl_finite) diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps index f68650d36d..36e1b76811 100644 --- a/sysdeps/x86_64/fpu/libm-test-ulps +++ b/sysdeps/x86_64/fpu/libm-test-ulps @@ -2033,4 +2033,28 @@ ifloat: 3 ildouble: 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 |