lh | 9ed821d | 2023-04-07 01:36:19 -0700 | [diff] [blame] | 1 | /* Complex cosine hyperbole function for double. |
| 2 | Copyright (C) 1997-2015 Free Software Foundation, Inc. |
| 3 | This file is part of the GNU C Library. |
| 4 | Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. |
| 5 | |
| 6 | The GNU C Library is free software; you can redistribute it and/or |
| 7 | modify it under the terms of the GNU Lesser General Public |
| 8 | License as published by the Free Software Foundation; either |
| 9 | version 2.1 of the License, or (at your option) any later version. |
| 10 | |
| 11 | The GNU C Library is distributed in the hope that it will be useful, |
| 12 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| 14 | Lesser General Public License for more details. |
| 15 | |
| 16 | You should have received a copy of the GNU Lesser General Public |
| 17 | License along with the GNU C Library; if not, see |
| 18 | <http://www.gnu.org/licenses/>. */ |
| 19 | |
| 20 | #include <complex.h> |
| 21 | #include <fenv.h> |
| 22 | #include <math.h> |
| 23 | #include <math_private.h> |
| 24 | #include <float.h> |
| 25 | |
| 26 | __complex__ double |
| 27 | __ccosh (__complex__ double x) |
| 28 | { |
| 29 | __complex__ double retval; |
| 30 | int rcls = fpclassify (__real__ x); |
| 31 | int icls = fpclassify (__imag__ x); |
| 32 | |
| 33 | if (__glibc_likely (rcls >= FP_ZERO)) |
| 34 | { |
| 35 | /* Real part is finite. */ |
| 36 | if (__glibc_likely (icls >= FP_ZERO)) |
| 37 | { |
| 38 | /* Imaginary part is finite. */ |
| 39 | const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2); |
| 40 | double sinix, cosix; |
| 41 | |
| 42 | if (__glibc_likely (fabs (__imag__ x) > DBL_MIN)) |
| 43 | { |
| 44 | __sincos (__imag__ x, &sinix, &cosix); |
| 45 | } |
| 46 | else |
| 47 | { |
| 48 | sinix = __imag__ x; |
| 49 | cosix = 1.0; |
| 50 | } |
| 51 | |
| 52 | if (fabs (__real__ x) > t) |
| 53 | { |
| 54 | double exp_t = __ieee754_exp (t); |
| 55 | double rx = fabs (__real__ x); |
| 56 | if (signbit (__real__ x)) |
| 57 | sinix = -sinix; |
| 58 | rx -= t; |
| 59 | sinix *= exp_t / 2.0; |
| 60 | cosix *= exp_t / 2.0; |
| 61 | if (rx > t) |
| 62 | { |
| 63 | rx -= t; |
| 64 | sinix *= exp_t; |
| 65 | cosix *= exp_t; |
| 66 | } |
| 67 | if (rx > t) |
| 68 | { |
| 69 | /* Overflow (original real part of x > 3t). */ |
| 70 | __real__ retval = DBL_MAX * cosix; |
| 71 | __imag__ retval = DBL_MAX * sinix; |
| 72 | } |
| 73 | else |
| 74 | { |
| 75 | double exp_val = __ieee754_exp (rx); |
| 76 | __real__ retval = exp_val * cosix; |
| 77 | __imag__ retval = exp_val * sinix; |
| 78 | } |
| 79 | } |
| 80 | else |
| 81 | { |
| 82 | __real__ retval = __ieee754_cosh (__real__ x) * cosix; |
| 83 | __imag__ retval = __ieee754_sinh (__real__ x) * sinix; |
| 84 | } |
| 85 | |
| 86 | if (fabs (__real__ retval) < DBL_MIN) |
| 87 | { |
| 88 | volatile double force_underflow |
| 89 | = __real__ retval * __real__ retval; |
| 90 | (void) force_underflow; |
| 91 | } |
| 92 | if (fabs (__imag__ retval) < DBL_MIN) |
| 93 | { |
| 94 | volatile double force_underflow |
| 95 | = __imag__ retval * __imag__ retval; |
| 96 | (void) force_underflow; |
| 97 | } |
| 98 | } |
| 99 | else |
| 100 | { |
| 101 | __imag__ retval = __real__ x == 0.0 ? 0.0 : __nan (""); |
| 102 | __real__ retval = __nan ("") + __nan (""); |
| 103 | |
| 104 | if (icls == FP_INFINITE) |
| 105 | feraiseexcept (FE_INVALID); |
| 106 | } |
| 107 | } |
| 108 | else if (rcls == FP_INFINITE) |
| 109 | { |
| 110 | /* Real part is infinite. */ |
| 111 | if (__glibc_likely (icls > FP_ZERO)) |
| 112 | { |
| 113 | /* Imaginary part is finite. */ |
| 114 | double sinix, cosix; |
| 115 | |
| 116 | if (__glibc_likely (fabs (__imag__ x) > DBL_MIN)) |
| 117 | { |
| 118 | __sincos (__imag__ x, &sinix, &cosix); |
| 119 | } |
| 120 | else |
| 121 | { |
| 122 | sinix = __imag__ x; |
| 123 | cosix = 1.0; |
| 124 | } |
| 125 | |
| 126 | __real__ retval = __copysign (HUGE_VAL, cosix); |
| 127 | __imag__ retval = (__copysign (HUGE_VAL, sinix) |
| 128 | * __copysign (1.0, __real__ x)); |
| 129 | } |
| 130 | else if (icls == FP_ZERO) |
| 131 | { |
| 132 | /* Imaginary part is 0.0. */ |
| 133 | __real__ retval = HUGE_VAL; |
| 134 | __imag__ retval = __imag__ x * __copysign (1.0, __real__ x); |
| 135 | } |
| 136 | else |
| 137 | { |
| 138 | /* The addition raises the invalid exception. */ |
| 139 | __real__ retval = HUGE_VAL; |
| 140 | __imag__ retval = __nan ("") + __nan (""); |
| 141 | |
| 142 | if (icls == FP_INFINITE) |
| 143 | feraiseexcept (FE_INVALID); |
| 144 | } |
| 145 | } |
| 146 | else |
| 147 | { |
| 148 | __real__ retval = __nan (""); |
| 149 | __imag__ retval = __imag__ x == 0.0 ? __imag__ x : __nan (""); |
| 150 | } |
| 151 | |
| 152 | return retval; |
| 153 | } |
| 154 | weak_alias (__ccosh, ccosh) |
| 155 | #ifdef NO_LONG_DOUBLE |
| 156 | strong_alias (__ccosh, __ccoshl) |
| 157 | weak_alias (__ccosh, ccoshl) |
| 158 | #endif |