1/* Return arc hyperbolic sine for a complex float type, with the 2 imaginary part of the result possibly adjusted for use in 3 computing other functions. 4 Copyright (C) 1997-2018 Free Software Foundation, Inc. 5 This file is part of the GNU C Library. 6 7 The GNU C Library is free software; you can redistribute it and/or 8 modify it under the terms of the GNU Lesser General Public 9 License as published by the Free Software Foundation; either 10 version 2.1 of the License, or (at your option) any later version. 11 12 The GNU C Library is distributed in the hope that it will be useful, 13 but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15 Lesser General Public License for more details. 16 17 You should have received a copy of the GNU Lesser General Public 18 License along with the GNU C Library; if not, see 19 <http://www.gnu.org/licenses/>. */ 20 21#include "quadmath-imp.h" 22 23/* Return the complex inverse hyperbolic sine of finite nonzero Z, 24 with the imaginary part of the result subtracted from pi/2 if ADJ 25 is nonzero. */ 26 27__complex128 28__quadmath_kernel_casinhq (__complex128 x, int adj) 29{ 30 __complex128 res; 31 __float128 rx, ix; 32 __complex128 y; 33 34 /* Avoid cancellation by reducing to the first quadrant. */ 35 rx = fabsq (__real__ x); 36 ix = fabsq (__imag__ x); 37 38 if (rx >= 1 / FLT128_EPSILON || ix >= 1 / FLT128_EPSILON) 39 { 40 /* For large x in the first quadrant, x + csqrt (1 + x * x) 41 is sufficiently close to 2 * x to make no significant 42 difference to the result; avoid possible overflow from 43 the squaring and addition. */ 44 __real__ y = rx; 45 __imag__ y = ix; 46 47 if (adj) 48 { 49 __float128 t = __real__ y; 50 __real__ y = copysignq (__imag__ y, __imag__ x); 51 __imag__ y = t; 52 } 53 54 res = clogq (y); 55 __real__ res += (__float128) M_LN2q; 56 } 57 else if (rx >= 0.5Q && ix < FLT128_EPSILON / 8) 58 { 59 __float128 s = hypotq (1, rx); 60 61 __real__ res = logq (rx + s); 62 if (adj) 63 __imag__ res = atan2q (s, __imag__ x); 64 else 65 __imag__ res = atan2q (ix, s); 66 } 67 else if (rx < FLT128_EPSILON / 8 && ix >= 1.5Q) 68 { 69 __float128 s = sqrtq ((ix + 1) * (ix - 1)); 70 71 __real__ res = logq (ix + s); 72 if (adj) 73 __imag__ res = atan2q (rx, copysignq (s, __imag__ x)); 74 else 75 __imag__ res = atan2q (s, rx); 76 } 77 else if (ix > 1 && ix < 1.5Q && rx < 0.5Q) 78 { 79 if (rx < FLT128_EPSILON * FLT128_EPSILON) 80 { 81 __float128 ix2m1 = (ix + 1) * (ix - 1); 82 __float128 s = sqrtq (ix2m1); 83 84 __real__ res = log1pq (2 * (ix2m1 + ix * s)) / 2; 85 if (adj) 86 __imag__ res = atan2q (rx, copysignq (s, __imag__ x)); 87 else 88 __imag__ res = atan2q (s, rx); 89 } 90 else 91 { 92 __float128 ix2m1 = (ix + 1) * (ix - 1); 93 __float128 rx2 = rx * rx; 94 __float128 f = rx2 * (2 + rx2 + 2 * ix * ix); 95 __float128 d = sqrtq (ix2m1 * ix2m1 + f); 96 __float128 dp = d + ix2m1; 97 __float128 dm = f / dp; 98 __float128 r1 = sqrtq ((dm + rx2) / 2); 99 __float128 r2 = rx * ix / r1; 100 101 __real__ res = log1pq (rx2 + dp + 2 * (rx * r1 + ix * r2)) / 2; 102 if (adj) 103 __imag__ res = atan2q (rx + r1, copysignq (ix + r2, __imag__ x)); 104 else 105 __imag__ res = atan2q (ix + r2, rx + r1); 106 } 107 } 108 else if (ix == 1 && rx < 0.5Q) 109 { 110 if (rx < FLT128_EPSILON / 8) 111 { 112 __real__ res = log1pq (2 * (rx + sqrtq (rx))) / 2; 113 if (adj) 114 __imag__ res = atan2q (sqrtq (rx), copysignq (1, __imag__ x)); 115 else 116 __imag__ res = atan2q (1, sqrtq (rx)); 117 } 118 else 119 { 120 __float128 d = rx * sqrtq (4 + rx * rx); 121 __float128 s1 = sqrtq ((d + rx * rx) / 2); 122 __float128 s2 = sqrtq ((d - rx * rx) / 2); 123 124 __real__ res = log1pq (rx * rx + d + 2 * (rx * s1 + s2)) / 2; 125 if (adj) 126 __imag__ res = atan2q (rx + s1, copysignq (1 + s2, __imag__ x)); 127 else 128 __imag__ res = atan2q (1 + s2, rx + s1); 129 } 130 } 131 else if (ix < 1 && rx < 0.5Q) 132 { 133 if (ix >= FLT128_EPSILON) 134 { 135 if (rx < FLT128_EPSILON * FLT128_EPSILON) 136 { 137 __float128 onemix2 = (1 + ix) * (1 - ix); 138 __float128 s = sqrtq (onemix2); 139 140 __real__ res = log1pq (2 * rx / s) / 2; 141 if (adj) 142 __imag__ res = atan2q (s, __imag__ x); 143 else 144 __imag__ res = atan2q (ix, s); 145 } 146 else 147 { 148 __float128 onemix2 = (1 + ix) * (1 - ix); 149 __float128 rx2 = rx * rx; 150 __float128 f = rx2 * (2 + rx2 + 2 * ix * ix); 151 __float128 d = sqrtq (onemix2 * onemix2 + f); 152 __float128 dp = d + onemix2; 153 __float128 dm = f / dp; 154 __float128 r1 = sqrtq ((dp + rx2) / 2); 155 __float128 r2 = rx * ix / r1; 156 157 __real__ res = log1pq (rx2 + dm + 2 * (rx * r1 + ix * r2)) / 2; 158 if (adj) 159 __imag__ res = atan2q (rx + r1, copysignq (ix + r2, 160 __imag__ x)); 161 else 162 __imag__ res = atan2q (ix + r2, rx + r1); 163 } 164 } 165 else 166 { 167 __float128 s = hypotq (1, rx); 168 169 __real__ res = log1pq (2 * rx * (rx + s)) / 2; 170 if (adj) 171 __imag__ res = atan2q (s, __imag__ x); 172 else 173 __imag__ res = atan2q (ix, s); 174 } 175 math_check_force_underflow_nonneg (__real__ res); 176 } 177 else 178 { 179 __real__ y = (rx - ix) * (rx + ix) + 1; 180 __imag__ y = 2 * rx * ix; 181 182 y = csqrtq (y); 183 184 __real__ y += rx; 185 __imag__ y += ix; 186 187 if (adj) 188 { 189 __float128 t = __real__ y; 190 __real__ y = copysignq (__imag__ y, __imag__ x); 191 __imag__ y = t; 192 } 193 194 res = clogq (y); 195 } 196 197 /* Give results the correct sign for the original argument. */ 198 __real__ res = copysignq (__real__ res, __real__ x); 199 __imag__ res = copysignq (__imag__ res, (adj ? 1 : __imag__ x)); 200 201 return res; 202} 203