diff options
| -rw-r--r-- | ChangeLog | 31 | ||||
| -rw-r--r-- | NEWS | 2 | ||||
| -rw-r--r-- | include/complex.h | 12 | ||||
| -rw-r--r-- | math/Makefile | 2 | ||||
| -rw-r--r-- | math/k_casinh.c | 85 | ||||
| -rw-r--r-- | math/k_casinhf.c | 85 | ||||
| -rw-r--r-- | math/k_casinhl.c | 92 | ||||
| -rw-r--r-- | math/libm-test.inc | 37 | ||||
| -rw-r--r-- | math/s_cacos.c | 26 | ||||
| -rw-r--r-- | math/s_cacosf.c | 26 | ||||
| -rw-r--r-- | math/s_cacosl.c | 26 | ||||
| -rw-r--r-- | math/s_casinh.c | 36 | ||||
| -rw-r--r-- | math/s_casinhf.c | 36 | ||||
| -rw-r--r-- | math/s_casinhl.c | 43 | ||||
| -rw-r--r-- | sysdeps/i386/fpu/libm-test-ulps | 6 | ||||
| -rw-r--r-- | sysdeps/x86_64/fpu/libm-test-ulps | 15 |
16 files changed, 430 insertions, 130 deletions
diff --git a/ChangeLog b/ChangeLog index e7cfb64231..890e3d4930 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,34 @@ +2013-01-17 Joseph Myers <joseph@codesourcery.com> + + [BZ #15023] + * include/complex.h: Condition contents on [!_COMPLEX_H]. + (__kernel_casinhf): New prototype. + (__kernel_casinh): Likewise. + (__kernel_casinhl): Likewise. + * math/Makefile (libm_calls): Add k_casinh. + * math/k_casinh.c: New file. + * math/k_casinhf.c: Likewise. + * math/k_casinhl.c: Likewise. + * math/s_cacos.c (__cacos): Implement using __kernel_casinh for + finite nonzero arguments. + * math/s_cacosf.c (__cacosf): Implement using __kernel_casinhf for + finite nonzero arguments. + * math/s_cacosl.c (__cacosl): Implement using __kernel_casinhl for + finite nonzero arguments. + * math/s_casinh.c: Do not include <float.h>. + (__casinh): Move code for finite nonzero arguments to k_casinh.c. + * math/s_casinhf.c: Do not include <float.h>. + (__casinhf): Move code for finite nonzero arguments to + k_casinhf.c. + * math/s_casinhl.c: Do not include <float.h>. + [LDBL_MANT_DIG == 106] (LDBL_EPSILON): Do not undefine and + redefine. + (__casinhl): Move code for finite nonzero arguments to + k_casinhl.c. + * math/libm-test.inc (cacos_test): Add more tests. + * sysdeps/i386/fpu/libm-test-ulps: Update. + * sysdeps/x86_64/fpu/libm-test-ulps: Likewise. + 2013-01-17 Pino Toscano <toscano.pino@tiscali.it> * sysdeps/unix/sysv/linux/malloc-sysdep.h (HAVE_MREMAP): New define. diff --git a/NEWS b/NEWS index 55f65fd686..99fb54226f 100644 --- a/NEWS +++ b/NEWS @@ -10,7 +10,7 @@ Version 2.18 * The following bugs are resolved with this release: 13951, 14200, 14317, 14327, 14964, 14981, 14982, 14985, 14994, 14996, - 15003. + 15003, 15023. Version 2.17 diff --git a/include/complex.h b/include/complex.h index acf8cf14ba..e173f1f6a3 100644 --- a/include/complex.h +++ b/include/complex.h @@ -1 +1,11 @@ -#include <math/complex.h> +#ifndef _COMPLEX_H +# include <math/complex.h> + +/* Return the complex inverse hyperbolic sine of finite nonzero Z, + with the imaginary part of the result subtracted from pi/2 if ADJ + is nonzero. */ +extern complex float __kernel_casinhf (complex float z, int adj); +extern complex double __kernel_casinh (complex double z, int adj); +extern complex long double __kernel_casinhl (complex long double z, int adj); + +#endif diff --git a/math/Makefile b/math/Makefile index b9519cfc24..da18b56d4b 100644 --- a/math/Makefile +++ b/math/Makefile @@ -58,7 +58,7 @@ libm-calls = e_acos e_acosh e_asin e_atan2 e_atanh e_cosh e_exp e_fmod \ s_catan s_casin s_ccos s_csin s_ctan s_ctanh s_cacos \ s_casinh s_cacosh s_catanh s_csqrt s_cpow s_cproj s_clog10 \ s_fma s_lrint s_llrint s_lround s_llround e_exp10 w_log2 \ - s_isinf_ns $(calls:s_%=m_%) x2y2m1 + s_isinf_ns $(calls:s_%=m_%) x2y2m1 k_casinh include ../Makeconfig diff --git a/math/k_casinh.c b/math/k_casinh.c new file mode 100644 index 0000000000..7f98f24a80 --- /dev/null +++ b/math/k_casinh.c @@ -0,0 +1,85 @@ +/* Return arc hyperbole sine for double value, with the imaginary part + of the result possibly adjusted for use in computing other + functions. + Copyright (C) 1997-2013 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <complex.h> +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Return the complex inverse hyperbolic sine of finite nonzero Z, + with the imaginary part of the result subtracted from pi/2 if ADJ + is nonzero. */ + +__complex__ double +__kernel_casinh (__complex__ double x, int adj) +{ + __complex__ double res; + double rx, ix; + __complex__ double y; + + /* Avoid cancellation by reducing to the first quadrant. */ + rx = fabs (__real__ x); + ix = fabs (__imag__ x); + + if (rx >= 1.0 / DBL_EPSILON || ix >= 1.0 / DBL_EPSILON) + { + /* For large x in the first quadrant, x + csqrt (1 + x * x) + is sufficiently close to 2 * x to make no significant + difference to the result; avoid possible overflow from + the squaring and addition. */ + __real__ y = rx; + __imag__ y = ix; + + if (adj) + { + double t = __real__ y; + __real__ y = __copysign (__imag__ y, __imag__ x); + __imag__ y = t; + } + + res = __clog (y); + __real__ res += M_LN2; + } + else + { + __real__ y = (rx - ix) * (rx + ix) + 1.0; + __imag__ y = 2.0 * rx * ix; + + y = __csqrt (y); + + __real__ y += rx; + __imag__ y += ix; + + if (adj) + { + double t = __real__ y; + __real__ y = copysign (__imag__ y, __imag__ x); + __imag__ y = t; + } + + res = __clog (y); + } + + /* Give results the correct sign for the original argument. */ + __real__ res = __copysign (__real__ res, __real__ x); + __imag__ res = __copysign (__imag__ res, (adj ? 1.0 : __imag__ x)); + + return res; +} diff --git a/math/k_casinhf.c b/math/k_casinhf.c new file mode 100644 index 0000000000..9401636348 --- /dev/null +++ b/math/k_casinhf.c @@ -0,0 +1,85 @@ +/* Return arc hyperbole sine for float value, with the imaginary part + of the result possibly adjusted for use in computing other + functions. + Copyright (C) 1997-2013 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <complex.h> +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* Return the complex inverse hyperbolic sine of finite nonzero Z, + with the imaginary part of the result subtracted from pi/2 if ADJ + is nonzero. */ + +__complex__ float +__kernel_casinhf (__complex__ float x, int adj) +{ + __complex__ float res; + float rx, ix; + __complex__ float y; + + /* Avoid cancellation by reducing to the first quadrant. */ + rx = fabsf (__real__ x); + ix = fabsf (__imag__ x); + + if (rx >= 1.0f / FLT_EPSILON || ix >= 1.0f / FLT_EPSILON) + { + /* For large x in the first quadrant, x + csqrt (1 + x * x) + is sufficiently close to 2 * x to make no significant + difference to the result; avoid possible overflow from + the squaring and addition. */ + __real__ y = rx; + __imag__ y = ix; + + if (adj) + { + float t = __real__ y; + __real__ y = __copysignf (__imag__ y, __imag__ x); + __imag__ y = t; + } + + res = __clogf (y); + __real__ res += (float) M_LN2; + } + else + { + __real__ y = (rx - ix) * (rx + ix) + 1.0; + __imag__ y = 2.0 * rx * ix; + + y = __csqrtf (y); + + __real__ y += rx; + __imag__ y += ix; + + if (adj) + { + float t = __real__ y; + __real__ y = __copysignf (__imag__ y, __imag__ x); + __imag__ y = t; + } + + res = __clogf (y); + } + + /* Give results the correct sign for the original argument. */ + __real__ res = __copysignf (__real__ res, __real__ x); + __imag__ res = __copysignf (__imag__ res, (adj ? 1.0f : __imag__ x)); + + return res; +} diff --git a/math/k_casinhl.c b/math/k_casinhl.c new file mode 100644 index 0000000000..6412979755 --- /dev/null +++ b/math/k_casinhl.c @@ -0,0 +1,92 @@ +/* Return arc hyperbole sine for long double value, with the imaginary + part of the result possibly adjusted for use in computing other + functions. + Copyright (C) 1997-2013 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <complex.h> +#include <math.h> +#include <math_private.h> +#include <float.h> + +/* To avoid spurious overflows, use this definition to treat IBM long + double as approximating an IEEE-style format. */ +#if LDBL_MANT_DIG == 106 +# undef LDBL_EPSILON +# define LDBL_EPSILON 0x1p-106L +#endif + +/* Return the complex inverse hyperbolic sine of finite nonzero Z, + with the imaginary part of the result subtracted from pi/2 if ADJ + is nonzero. */ + +__complex__ long double +__kernel_casinhl (__complex__ long double x, int adj) +{ + __complex__ long double res; + long double rx, ix; + __complex__ long double y; + + /* Avoid cancellation by reducing to the first quadrant. */ + rx = fabsl (__real__ x); + ix = fabsl (__imag__ x); + + if (rx >= 1.0L / LDBL_EPSILON || ix >= 1.0L / LDBL_EPSILON) + { + /* For large x in the first quadrant, x + csqrt (1 + x * x) + is sufficiently close to 2 * x to make no significant + difference to the result; avoid possible overflow from + the squaring and addition. */ + __real__ y = rx; + __imag__ y = ix; + + if (adj) + { + long double t = __real__ y; + __real__ y = __copysignl (__imag__ y, __imag__ x); + __imag__ y = t; + } + + res = __clogl (y); + __real__ res += M_LN2l; + } + else + { + __real__ y = (rx - ix) * (rx + ix) + 1.0; + __imag__ y = 2.0 * rx * ix; + + y = __csqrtl (y); + + __real__ y += rx; + __imag__ y += ix; + + if (adj) + { + long double t = __real__ y; + __real__ y = __copysignl (__imag__ y, __imag__ x); + __imag__ y = t; + } + + res = __clogl (y); + } + + /* Give results the correct sign for the original argument. */ + __real__ res = __copysignl (__real__ res, __real__ x); + __imag__ res = __copysignl (__imag__ res, (adj ? 1.0L : __imag__ x)); + + return res; +} diff --git a/math/libm-test.inc b/math/libm-test.inc index 56e321781e..1c2970fc79 100644 --- a/math/libm-test.inc +++ b/math/libm-test.inc @@ -1453,6 +1453,43 @@ cacos_test (void) TEST_c_c (cacos, 1.5L, plus_zero, plus_zero, -0.9624236501192068949955178268487368462704L); TEST_c_c (cacos, 1.5L, minus_zero, plus_zero, 0.9624236501192068949955178268487368462704L); + TEST_c_c (cacos, 0x1p50L, 1.0L, 8.881784197001252323389053344727730248720e-16L, -3.535050620855721078027883819436720218708e1L); + TEST_c_c (cacos, 0x1p50L, -1.0L, 8.881784197001252323389053344727730248720e-16L, 3.535050620855721078027883819436720218708e1L); + TEST_c_c (cacos, -0x1p50L, 1.0L, 3.141592653589792350284223683154270545292L, -3.535050620855721078027883819436720218708e1L); + TEST_c_c (cacos, -0x1p50L, -1.0L, 3.141592653589792350284223683154270545292L, 3.535050620855721078027883819436720218708e1L); + TEST_c_c (cacos, 1.0L, 0x1p50L, 1.570796326794895731052901991514519103193L, -3.535050620855721078027883819436759661753e1L); + TEST_c_c (cacos, -1.0L, 0x1p50L, 1.570796326794897507409741391764983781004L, -3.535050620855721078027883819436759661753e1L); + TEST_c_c (cacos, 1.0L, -0x1p50L, 1.570796326794895731052901991514519103193L, 3.535050620855721078027883819436759661753e1L); + TEST_c_c (cacos, -1.0L, -0x1p50L, 1.570796326794897507409741391764983781004L, 3.535050620855721078027883819436759661753e1L); +#ifndef TEST_FLOAT + TEST_c_c (cacos, 0x1p500L, 1.0L, 3.054936363499604682051979393213617699789e-151L, -3.472667374605326000180332928505464606058e2L); + TEST_c_c (cacos, 0x1p500L, -1.0L, 3.054936363499604682051979393213617699789e-151L, 3.472667374605326000180332928505464606058e2L); + TEST_c_c (cacos, -0x1p500L, 1.0L, 3.141592653589793238462643383279502884197L, -3.472667374605326000180332928505464606058e2L); + TEST_c_c (cacos, -0x1p500L, -1.0L, 3.141592653589793238462643383279502884197L, 3.472667374605326000180332928505464606058e2L); + TEST_c_c (cacos, 1.0L, 0x1p500L, 1.570796326794896619231321691639751442099L, -3.472667374605326000180332928505464606058e2L); + TEST_c_c (cacos, -1.0L, 0x1p500L, 1.570796326794896619231321691639751442099L, -3.472667374605326000180332928505464606058e2L); + TEST_c_c (cacos, 1.0L, -0x1p500L, 1.570796326794896619231321691639751442099L, 3.472667374605326000180332928505464606058e2L); + TEST_c_c (cacos, -1.0L, -0x1p500L, 1.570796326794896619231321691639751442099L, 3.472667374605326000180332928505464606058e2L); +#endif +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 + TEST_c_c (cacos, 0x1p5000L, 1.0L, 7.079811261048172892385615158694057552948e-1506L, -3.466429049980286492395577839412341016946e3L); + TEST_c_c (cacos, 0x1p5000L, -1.0L, 7.079811261048172892385615158694057552948e-1506L, 3.466429049980286492395577839412341016946e3L); + TEST_c_c (cacos, -0x1p5000L, 1.0L, 3.141592653589793238462643383279502884197L, -3.466429049980286492395577839412341016946e3L); + TEST_c_c (cacos, -0x1p5000L, -1.0L, 3.141592653589793238462643383279502884197L, 3.466429049980286492395577839412341016946e3L); + TEST_c_c (cacos, 1.0L, 0x1p5000L, 1.570796326794896619231321691639751442099L, -3.466429049980286492395577839412341016946e3L); + TEST_c_c (cacos, -1.0L, 0x1p5000L, 1.570796326794896619231321691639751442099L, -3.466429049980286492395577839412341016946e3L); + TEST_c_c (cacos, 1.0L, -0x1p5000L, 1.570796326794896619231321691639751442099L, 3.466429049980286492395577839412341016946e3L); + TEST_c_c (cacos, -1.0L, -0x1p5000L, 1.570796326794896619231321691639751442099L, 3.466429049980286492395577839412341016946e3L); +#endif + + TEST_c_c (cacos, 0x1.fp127L, 0x1.fp127L, 7.853981633974483096156608458198757210493e-1L, -8.973081118419833726837456344608533993585e1L); +#ifndef TEST_FLOAT + TEST_c_c (cacos, 0x1.fp1023L, 0x1.fp1023L, 7.853981633974483096156608458198757210493e-1L, -7.107906849659093345062145442726115449315e2L); +#endif +#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384 + TEST_c_c (cacos, 0x1.fp16383L, 0x1.fp16383L, 7.853981633974483096156608458198757210493e-1L, -1.135753137836666928715489992987020363057e4L); +#endif + TEST_c_c (cacos, 0.75L, 1.25L, 1.11752014915610270578240049553777969L, -1.13239363160530819522266333696834467L); TEST_c_c (cacos, -2, -3, 2.1414491111159960199416055713254211L, 1.9833870299165354323470769028940395L); diff --git a/math/s_cacos.c b/math/s_cacos.c index 6604b5aec6..acd9b2462a 100644 --- a/math/s_cacos.c +++ b/math/s_cacos.c @@ -25,11 +25,27 @@ __cacos (__complex__ double x) { __complex__ double y; __complex__ double res; - - y = __casin (x); - - __real__ res = (double) M_PI_2 - __real__ y; - __imag__ res = -__imag__ y; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE + || (rcls == FP_ZERO && icls == FP_ZERO)) + { + y = __casin (x); + + __real__ res = (double) M_PI_2 - __real__ y; + __imag__ res = -__imag__ y; + } + else + { + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + y = __kernel_casinh (y, 1); + + __real__ res = __imag__ y; + __imag__ res = __real__ y; + } return res; } diff --git a/math/s_cacosf.c b/math/s_cacosf.c index 04c13e4fa5..df2bf218a3 100644 --- a/math/s_cacosf.c +++ b/math/s_cacosf.c @@ -25,11 +25,27 @@ __cacosf (__complex__ float x) { __complex__ float y; __complex__ float res; - - y = __casinf (x); - - __real__ res = (float) M_PI_2 - __real__ y; - __imag__ res = -__imag__ y; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE + || (rcls == FP_ZERO && icls == FP_ZERO)) + { + y = __casinf (x); + + __real__ res = (float) M_PI_2 - __real__ y; + __imag__ res = -__imag__ y; + } + else + { + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + y = __kernel_casinhf (y, 1); + + __real__ res = __imag__ y; + __imag__ res = __real__ y; + } return res; } diff --git a/math/s_cacosl.c b/math/s_cacosl.c index 304076ddfe..8eab1f0004 100644 --- a/math/s_cacosl.c +++ b/math/s_cacosl.c @@ -25,11 +25,27 @@ __cacosl (__complex__ long double x) { __complex__ long double y; __complex__ long double res; - - y = __casinl (x); - - __real__ res = M_PI_2l - __real__ y; - __imag__ res = -__imag__ y; + int rcls = fpclassify (__real__ x); + int icls = fpclassify (__imag__ x); + + if (rcls <= FP_INFINITE || icls <= FP_INFINITE + || (rcls == FP_ZERO && icls == FP_ZERO)) + { + y = __casinl (x); + + __real__ res = M_PI_2l - __real__ y; + __imag__ res = -__imag__ y; + } + else + { + __real__ y = -__imag__ x; + __imag__ y = __real__ x; + + y = __kernel_casinhl (y, 1); + + __real__ res = __imag__ y; + __imag__ res = __real__ y; + } return res; } diff --git a/math/s_casinh.c b/math/s_casinh.c index b493982d88..657e269ac1 100644 --- a/math/s_casinh.c +++ b/math/s_casinh.c @@ -20,7 +20,6 @@ #include <complex.h> #include <math.h> #include <math_private.h> -#include <float.h> __complex__ double __casinh (__complex__ double x) @@ -62,40 +61,7 @@ __casinh (__complex__ double x) } else { - double rx, ix; - __complex__ double y; - - /* Avoid cancellation by reducing to the first quadrant. */ - rx = fabs (__real__ x); - ix = fabs (__imag__ x); - - if (rx >= 1.0 / DBL_EPSILON || ix >= 1.0 / DBL_EPSILON) - { - /* For large x in the first quadrant, x + csqrt (1 + x * x) - is sufficiently close to 2 * x to make no significant - difference to the result; avoid possible overflow from - the squaring and addition. */ - __real__ y = rx; - __imag__ y = ix; - res = __clog (y); - __real__ res += M_LN2; - } - else - { - __real__ y = (rx - ix) * (rx + ix) + 1.0; - __imag__ y = 2.0 * rx * ix; - - y = __csqrt (y); - - __real__ y += rx; - __imag__ y += ix; - - res = __clog (y); - } - - /* Give results the correct sign for the original argument. */ - __real__ res = __copysign (__real__ res, __real__ x); - __imag__ res = __copysign (__imag__ res, __imag__ x); + res = __kernel_casinh (x, 0); } return res; diff --git a/math/s_casinhf.c b/math/s_casinhf.c index f865e14490..8663c2e7cc 100644 --- a/math/s_casinhf.c +++ b/math/s_casinhf.c @@ -20,7 +20,6 @@ #include <complex.h> #include <math.h> #include <math_private.h> -#include <float.h> __complex__ float __casinhf (__complex__ float x) @@ -62,40 +61,7 @@ __casinhf (__complex__ float x) } else { - float rx, ix; - __complex__ float y; - - /* Avoid cancellation by reducing to the first quadrant. */ - rx = fabsf (__real__ x); - ix = fabsf (__imag__ x); - - if (rx >= 1.0f / FLT_EPSILON || ix >= 1.0f / FLT_EPSILON) - { - /* For large x in the first quadrant, x + csqrt (1 + x * x) - is sufficiently close to 2 * x to make no significant - difference to the result; avoid possible overflow from - the squaring and addition. */ - __real__ y = rx; - __imag__ y = ix; - res = __clogf (y); - __real__ res += (float) M_LN2; - } - else - { - __real__ y = (rx - ix) * (rx + ix) + 1.0; - __imag__ y = 2.0 * rx * ix; - - y = __csqrtf (y); - - __real__ y += rx; - __imag__ y += ix; - - res = __clogf (y); - } - - /* Give results the correct sign for the original argument. */ - __real__ res = __copysignf (__real__ res, __real__ x); - __imag__ res = __copysignf (__imag__ res, __imag__ x); + res = __kernel_casinhf (x, 0); } return res; diff --git a/math/s_casinhl.c b/math/s_casinhl.c index d7c74593e4..2afc52714e 100644 --- a/math/s_casinhl.c +++ b/math/s_casinhl.c @@ -20,14 +20,6 @@ #include <complex.h> #include <math.h> #include <math_private.h> -#include <float.h> - -/* To avoid spurious overflows, use this definition to treat IBM long - double as approximating an IEEE-style format. */ -#if LDBL_MANT_DIG == 106 -# undef LDBL_EPSILON -# define LDBL_EPSILON 0x1p-106L -#endif __complex__ long double __casinhl (__complex__ long double x) @@ -69,40 +61,7 @@ __casinhl (__complex__ long double x) } else { - long double rx, ix; - __complex__ long double y; - - /* Avoid cancellation by reducing to the first quadrant. */ - rx = fabsl (__real__ x); - ix = fabsl (__imag__ x); - - if (rx >= 1.0L / LDBL_EPSILON || ix >= 1.0L / LDBL_EPSILON) - { - /* For large x in the first quadrant, x + csqrt (1 + x * x) - is sufficiently close to 2 * x to make no significant - difference to the result; avoid possible overflow from - the squaring and addition. */ - __real__ y = rx; - __imag__ y = ix; - res = __clogl (y); - __real__ res += M_LN2l; - } - else - { - __real__ y = (rx - ix) * (rx + ix) + 1.0; - __imag__ y = 2.0 * rx * ix; - - y = __csqrtl (y); - - __real__ y += rx; - __imag__ y += ix; - - res = __clogl (y); - } - - /* Give results the correct sign for the original argument. */ - __real__ res = __copysignl (__real__ res, __real__ x); - __imag__ res = __copysignl (__imag__ res, __imag__ x); + res = __kernel_casinhl (x, 0); } return res; diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps index 3fc30de462..1525b16f37 100644 --- a/sysdeps/i386/fpu/libm-test-ulps +++ b/sysdeps/i386/fpu/libm-test-ulps @@ -303,6 +303,12 @@ float: 1 ifloat: 1 ildouble: 2 ldouble: 2 +Test "Imaginary part of: cacos (0x1.fp1023 + 0x1.fp1023 i) == 7.853981633974483096156608458198757210493e-1 - 7.107906849659093345062145442726115449315e2 i": +double: 1 +idouble: 1 +Test "Imaginary part of: cacos (0x1.fp127 + 0x1.fp127 i) == 7.853981633974483096156608458198757210493e-1 - 8.973081118419833726837456344608533993585e1 i": +double: 1 +idouble: 1 Test "Imaginary part of: cacos (1.5 + +0 i) == +0 - 0.9624236501192068949955178268487368462704 i": double: 1 float: 1 diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps index 95b6aec81e..63c6aed2a5 100644 --- a/sysdeps/x86_64/fpu/libm-test-ulps +++ b/sysdeps/x86_64/fpu/libm-test-ulps @@ -244,6 +244,12 @@ ifloat: 1 Test "Imaginary part of: cacos (-0 - 1.5 i) == pi/2 + 1.194763217287109304111930828519090523536 i": double: 1 idouble: 1 +Test "Real part of: cacos (-1.0 + 0x1p50 i) == 1.570796326794897507409741391764983781004 - 3.535050620855721078027883819436759661753e1 i": +float: 1 +ifloat: 1 +Test "Real part of: cacos (-1.0 - 0x1p50 i) == 1.570796326794897507409741391764983781004 + 3.535050620855721078027883819436759661753e1 i": +float: 1 +ifloat: 1 Test "Imaginary part of: cacos (-1.5 + +0 i) == pi - 0.9624236501192068949955178268487368462704 i": double: 1 float: 1 @@ -254,6 +260,9 @@ ldouble: 1 Test "Imaginary part of: cacos (-1.5 - 0 i) == pi + 0.9624236501192068949955178268487368462704 i": ildouble: 1 ldouble: 1 +Test "Real part of: cacos (-2 - 3 i) == 2.1414491111159960199416055713254211 + 1.9833870299165354323470769028940395 i": +float: 1 +ifloat: 1 Test "Real part of: cacos (0.5 + +0 i) == 1.047197551196597746154214461093167628066 - 0 i": double: 1 idouble: 1 @@ -272,6 +281,12 @@ float: 1 ifloat: 1 ildouble: 2 ldouble: 2 +Test "Imaginary part of: cacos (0x1.fp1023 + 0x1.fp1023 i) == 7.853981633974483096156608458198757210493e-1 - 7.107906849659093345062145442726115449315e2 i": +double: 1 +idouble: 1 +Test "Imaginary part of: cacos (0x1.fp127 + 0x1.fp127 i) == 7.853981633974483096156608458198757210493e-1 - 8.973081118419833726837456344608533993585e1 i": +double: 1 +idouble: 1 Test "Imaginary part of: cacos (1.5 + +0 i) == +0 - 0.9624236501192068949955178268487368462704 i": double: 1 float: 1 |