1/* Return arc hyperbole sine for double value, with the imaginary part 2 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/* Return the complex inverse hyperbolic sine of finite nonzero Z, 27 with the imaginary part of the result subtracted from pi/2 if ADJ 28 is nonzero. */ 29 30__complex__ double 31__kernel_casinh (__complex__ double x, int adj) 32{ 33 __complex__ double res; 34 double rx, ix; 35 __complex__ double y; 36 37 /* Avoid cancellation by reducing to the first quadrant. */ 38 rx = fabs (__real__ x); 39 ix = fabs (__imag__ x); 40 41 if (rx >= 1.0 / DBL_EPSILON || ix >= 1.0 / DBL_EPSILON) 42 { 43 /* For large x in the first quadrant, x + csqrt (1 + x * x) 44 is sufficiently close to 2 * x to make no significant 45 difference to the result; avoid possible overflow from 46 the squaring and addition. */ 47 __real__ y = rx; 48 __imag__ y = ix; 49 50 if (adj) 51 { 52 double t = __real__ y; 53 __real__ y = copysign (__imag__ y, __imag__ x); 54 __imag__ y = t; 55 } 56 57 res = clog (y); 58 __real__ res += M_LN2; 59 } 60 else if (rx >= 0.5 && ix < DBL_EPSILON / 8.0) 61 { 62 double s = hypot (1.0, rx); 63 64 __real__ res = log (rx + s); 65 if (adj) 66 __imag__ res = atan2 (s, __imag__ x); 67 else 68 __imag__ res = atan2 (ix, s); 69 } 70 else if (rx < DBL_EPSILON / 8.0 && ix >= 1.5) 71 { 72 double s = sqrt ((ix + 1.0) * (ix - 1.0)); 73 74 __real__ res = log (ix + s); 75 if (adj) 76 __imag__ res = atan2 (rx, copysign (s, __imag__ x)); 77 else 78 __imag__ res = atan2 (s, rx); 79 } 80 else if (ix > 1.0 && ix < 1.5 && rx < 0.5) 81 { 82 if (rx < DBL_EPSILON * DBL_EPSILON) 83 { 84 double ix2m1 = (ix + 1.0) * (ix - 1.0); 85 double s = sqrt (ix2m1); 86 87 __real__ res = log1p (2.0 * (ix2m1 + ix * s)) / 2.0; 88 if (adj) 89 __imag__ res = atan2 (rx, copysign (s, __imag__ x)); 90 else 91 __imag__ res = atan2 (s, rx); 92 } 93 else 94 { 95 double ix2m1 = (ix + 1.0) * (ix - 1.0); 96 double rx2 = rx * rx; 97 double f = rx2 * (2.0 + rx2 + 2.0 * ix * ix); 98 double d = sqrt (ix2m1 * ix2m1 + f); 99 double dp = d + ix2m1; 100 double dm = f / dp; 101 double r1 = sqrt ((dm + rx2) / 2.0); 102 double r2 = rx * ix / r1; 103 104 __real__ res = log1p (rx2 + dp + 2.0 * (rx * r1 + ix * r2)) / 2.0; 105 if (adj) 106 __imag__ res = atan2 (rx + r1, copysign (ix + r2, 107 __imag__ x)); 108 else 109 __imag__ res = atan2 (ix + r2, rx + r1); 110 } 111 } 112 else if (ix == 1.0 && rx < 0.5) 113 { 114 if (rx < DBL_EPSILON / 8.0) 115 { 116 __real__ res = log1p (2.0 * (rx + sqrt (rx))) / 2.0; 117 if (adj) 118 __imag__ res = atan2 (sqrt (rx), 119 copysign (1.0, __imag__ x)); 120 else 121 __imag__ res = atan2 (1.0, sqrt (rx)); 122 } 123 else 124 { 125 double d = rx * sqrt (4.0 + rx * rx); 126 double s1 = sqrt ((d + rx * rx) / 2.0); 127 double s2 = sqrt ((d - rx * rx) / 2.0); 128 129 __real__ res = log1p (rx * rx + d + 2.0 * (rx * s1 + s2)) / 2.0; 130 if (adj) 131 __imag__ res = atan2 (rx + s1, copysign (1.0 + s2, 132 __imag__ x)); 133 else 134 __imag__ res = atan2 (1.0 + s2, rx + s1); 135 } 136 } 137 else if (ix < 1.0 && rx < 0.5) 138 { 139 if (ix >= DBL_EPSILON) 140 { 141 if (rx < DBL_EPSILON * DBL_EPSILON) 142 { 143 double onemix2 = (1.0 + ix) * (1.0 - ix); 144 double s = sqrt (onemix2); 145 146 __real__ res = log1p (2.0 * rx / s) / 2.0; 147 if (adj) 148 __imag__ res = atan2 (s, __imag__ x); 149 else 150 __imag__ res = atan2 (ix, s); 151 } 152 else 153 { 154 double onemix2 = (1.0 + ix) * (1.0 - ix); 155 double rx2 = rx * rx; 156 double f = rx2 * (2.0 + rx2 + 2.0 * ix * ix); 157 double d = sqrt (onemix2 * onemix2 + f); 158 double dp = d + onemix2; 159 double dm = f / dp; 160 double r1 = sqrt ((dp + rx2) / 2.0); 161 double r2 = rx * ix / r1; 162 163 __real__ res 164 = log1p (rx2 + dm + 2.0 * (rx * r1 + ix * r2)) / 2.0; 165 if (adj) 166 __imag__ res = atan2 (rx + r1, 167 copysign (ix + r2, 168 __imag__ x)); 169 else 170 __imag__ res = atan2 (ix + r2, rx + r1); 171 } 172 } 173 else 174 { 175 double s = hypot (1.0, rx); 176 177 __real__ res = log1p (2.0 * rx * (rx + s)) / 2.0; 178 if (adj) 179 __imag__ res = atan2 (s, __imag__ x); 180 else 181 __imag__ res = atan2 (ix, s); 182 } 183 if (__real__ res < DBL_MIN) 184 { 185 volatile double force_underflow = __real__ res * __real__ res; 186 (void) force_underflow; 187 } 188 } 189 else 190 { 191 __real__ y = (rx - ix) * (rx + ix) + 1.0; 192 __imag__ y = 2.0 * rx * ix; 193 194 y = csqrt (y); 195 196 __real__ y += rx; 197 __imag__ y += ix; 198 199 if (adj) 200 { 201 double t = __real__ y; 202 __real__ y = copysign (__imag__ y, __imag__ x); 203 __imag__ y = t; 204 } 205 206 res = clog (y); 207 } 208 209 /* Give results the correct sign for the original argument. */ 210 __real__ res = copysign (__real__ res, __real__ x); 211 __imag__ res = copysign (__imag__ res, (adj ? 1.0 : __imag__ x)); 212 213 return res; 214} 215