diff options
Diffstat (limited to 'math')
| -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 |
11 files changed, 366 insertions, 128 deletions
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; |