1 /* Return arc hyperbole sine for long double value, with the imaginary 2 part of the result possibly adjusted for use in computing other 3 functions. 4 Copyright (C) 1997-2015 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 <complex.h> 22 #include <math.h> 23 #include <math_private.h> 24 #include <float.h> 25 26 /* To avoid spurious overflows, use this definition to treat IBM long 27 double as approximating an IEEE-style format. */ 28 #if LDBL_MANT_DIG == 106 29 # undef LDBL_EPSILON 30 # define LDBL_EPSILON 0x1p-106L 31 #endif 32 33 /* Return the complex inverse hyperbolic sine of finite nonzero Z, 34 with the imaginary part of the result subtracted from pi/2 if ADJ 35 is nonzero. */ 36 37 __complex__ long double 38 __kernel_casinhl (__complex__ long double x, int adj) 39 { 40 __complex__ long double res; 41 long double rx, ix; 42 __complex__ long double y; 43 44 /* Avoid cancellation by reducing to the first quadrant. */ 45 rx = fabsl (__real__ x); 46 ix = fabsl (__imag__ x); 47 48 if (rx >= 1.0L / LDBL_EPSILON || ix >= 1.0L / LDBL_EPSILON) 49 { 50 /* For large x in the first quadrant, x + csqrt (1 + x * x) 51 is sufficiently close to 2 * x to make no significant 52 difference to the result; avoid possible overflow from 53 the squaring and addition. */ 54 __real__ y = rx; 55 __imag__ y = ix; 56 57 if (adj) 58 { 59 long double t = __real__ y; 60 __real__ y = copysignl (__imag__ y, __imag__ x); 61 __imag__ y = t; 62 } 63 64 res = clogl (y); 65 __real__ res += M_LN2l; 66 } 67 else if (rx >= 0.5L && ix < LDBL_EPSILON / 8.0L) 68 { 69 long double s = hypotl (1.0L, rx); 70 71 __real__ res = logl (rx + s); 72 if (adj) 73 __imag__ res = atan2l (s, __imag__ x); 74 else 75 __imag__ res = atan2l (ix, s); 76 } 77 else if (rx < LDBL_EPSILON / 8.0L && ix >= 1.5L) 78 { 79 long double s = sqrtl ((ix + 1.0L) * (ix - 1.0L)); 80 81 __real__ res = logl (ix + s); 82 if (adj) 83 __imag__ res = atan2l (rx, copysignl (s, __imag__ x)); 84 else 85 __imag__ res = atan2l (s, rx); 86 } 87 else if (ix > 1.0L && ix < 1.5L && rx < 0.5L) 88 { 89 if (rx < LDBL_EPSILON * LDBL_EPSILON) 90 { 91 long double ix2m1 = (ix + 1.0L) * (ix - 1.0L); 92 long double s = sqrtl (ix2m1); 93 94 __real__ res = log1pl (2.0L * (ix2m1 + ix * s)) / 2.0L; 95 if (adj) 96 __imag__ res = atan2l (rx, copysignl (s, __imag__ x)); 97 else 98 __imag__ res = atan2l (s, rx); 99 } 100 else 101 { 102 long double ix2m1 = (ix + 1.0L) * (ix - 1.0L); 103 long double rx2 = rx * rx; 104 long double f = rx2 * (2.0L + rx2 + 2.0L * ix * ix); 105 long double d = sqrtl (ix2m1 * ix2m1 + f); 106 long double dp = d + ix2m1; 107 long double dm = f / dp; 108 long double r1 = sqrtl ((dm + rx2) / 2.0L); 109 long double r2 = rx * ix / r1; 110 111 __real__ res 112 = log1pl (rx2 + dp + 2.0L * (rx * r1 + ix * r2)) / 2.0L; 113 if (adj) 114 __imag__ res = atan2l (rx + r1, copysignl (ix + r2, 115 __imag__ x)); 116 else 117 __imag__ res = atan2l (ix + r2, rx + r1); 118 } 119 } 120 else if (ix == 1.0L && rx < 0.5L) 121 { 122 if (rx < LDBL_EPSILON / 8.0L) 123 { 124 __real__ res = log1pl (2.0L * (rx + sqrtl (rx))) / 2.0L; 125 if (adj) 126 __imag__ res = atan2l (sqrtl (rx), 127 copysignl (1.0L, __imag__ x)); 128 else 129 __imag__ res = atan2l (1.0L, sqrtl (rx)); 130 } 131 else 132 { 133 long double d = rx * sqrtl (4.0L + rx * rx); 134 long double s1 = sqrtl ((d + rx * rx) / 2.0L); 135 long double s2 = sqrtl ((d - rx * rx) / 2.0L); 136 137 __real__ res = log1pl (rx * rx + d + 2.0L * (rx * s1 + s2)) / 2.0L; 138 if (adj) 139 __imag__ res = atan2l (rx + s1, 140 copysignl (1.0L + s2, 141 __imag__ x)); 142 else 143 __imag__ res = atan2l (1.0L + s2, rx + s1); 144 } 145 } 146 else if (ix < 1.0L && rx < 0.5L) 147 { 148 if (ix >= LDBL_EPSILON) 149 { 150 if (rx < LDBL_EPSILON * LDBL_EPSILON) 151 { 152 long double onemix2 = (1.0L + ix) * (1.0L - ix); 153 long double s = sqrtl (onemix2); 154 155 __real__ res = log1pl (2.0L * rx / s) / 2.0L; 156 if (adj) 157 __imag__ res = atan2l (s, __imag__ x); 158 else 159 __imag__ res = atan2l (ix, s); 160 } 161 else 162 { 163 long double onemix2 = (1.0L + ix) * (1.0L - ix); 164 long double rx2 = rx * rx; 165 long double f = rx2 * (2.0L + rx2 + 2.0L * ix * ix); 166 long double d = sqrtl (onemix2 * onemix2 + f); 167 long double dp = d + onemix2; 168 long double dm = f / dp; 169 long double r1 = sqrtl ((dp + rx2) / 2.0L); 170 long double r2 = rx * ix / r1; 171 172 __real__ res 173 = log1pl (rx2 + dm + 2.0L * (rx * r1 + ix * r2)) / 2.0L; 174 if (adj) 175 __imag__ res = atan2l (rx + r1, 176 copysignl (ix + r2, 177 __imag__ x)); 178 else 179 __imag__ res = atan2l (ix + r2, rx + r1); 180 } 181 } 182 else 183 { 184 long double s = hypotl (1.0L, rx); 185 186 __real__ res = log1pl (2.0L * rx * (rx + s)) / 2.0L; 187 if (adj) 188 __imag__ res = atan2l (s, __imag__ x); 189 else 190 __imag__ res = atan2l (ix, s); 191 } 192 if (__real__ res < LDBL_MIN) 193 { 194 volatile long double force_underflow = __real__ res * __real__ res; 195 (void) force_underflow; 196 } 197 } 198 else 199 { 200 __real__ y = (rx - ix) * (rx + ix) + 1.0L; 201 __imag__ y = 2.0L * rx * ix; 202 203 y = csqrtl (y); 204 205 __real__ y += rx; 206 __imag__ y += ix; 207 208 if (adj) 209 { 210 long double t = __real__ y; 211 __real__ y = copysignl (__imag__ y, __imag__ x); 212 __imag__ y = t; 213 } 214 215 res = clogl (y); 216 } 217 218 /* Give results the correct sign for the original argument. */ 219 __real__ res = copysignl (__real__ res, __real__ x); 220 __imag__ res = copysignl (__imag__ res, (adj ? 1.0L : __imag__ x)); 221 222 return res; 223 } 224