1 /* Return arc hyperbole sine for float 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__ float 31 __kernel_casinhf (__complex__ float x, int adj) 32 { 33 __complex__ float res; 34 float rx, ix; 35 __complex__ float y; 36 37 /* Avoid cancellation by reducing to the first quadrant. */ 38 rx = fabsf (__real__ x); 39 ix = fabsf (__imag__ x); 40 41 if (rx >= 1.0f / FLT_EPSILON || ix >= 1.0f / FLT_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 float t = __real__ y; 53 __real__ y = copysignf (__imag__ y, __imag__ x); 54 __imag__ y = t; 55 } 56 57 res = clogf (y); 58 __real__ res += (float) M_LN2; 59 } 60 else if (rx >= 0.5f && ix < FLT_EPSILON / 8.0f) 61 { 62 float s = hypotf (1.0f, rx); 63 64 __real__ res = logf (rx + s); 65 if (adj) 66 __imag__ res = atan2f (s, __imag__ x); 67 else 68 __imag__ res = atan2f (ix, s); 69 } 70 else if (rx < FLT_EPSILON / 8.0f && ix >= 1.5f) 71 { 72 float s = sqrtf ((ix + 1.0f) * (ix - 1.0f)); 73 74 __real__ res = logf (ix + s); 75 if (adj) 76 __imag__ res = atan2f (rx, copysignf (s, __imag__ x)); 77 else 78 __imag__ res = atan2f (s, rx); 79 } 80 else if (ix > 1.0f && ix < 1.5f && rx < 0.5f) 81 { 82 if (rx < FLT_EPSILON * FLT_EPSILON) 83 { 84 float ix2m1 = (ix + 1.0f) * (ix - 1.0f); 85 float s = sqrtf (ix2m1); 86 87 __real__ res = log1pf (2.0f * (ix2m1 + ix * s)) / 2.0f; 88 if (adj) 89 __imag__ res = atan2f (rx, copysignf (s, __imag__ x)); 90 else 91 __imag__ res = atan2f (s, rx); 92 } 93 else 94 { 95 float ix2m1 = (ix + 1.0f) * (ix - 1.0f); 96 float rx2 = rx * rx; 97 float f = rx2 * (2.0f + rx2 + 2.0f * ix * ix); 98 float d = sqrtf (ix2m1 * ix2m1 + f); 99 float dp = d + ix2m1; 100 float dm = f / dp; 101 float r1 = sqrtf ((dm + rx2) / 2.0f); 102 float r2 = rx * ix / r1; 103 104 __real__ res 105 = log1pf (rx2 + dp + 2.0f * (rx * r1 + ix * r2)) / 2.0f; 106 if (adj) 107 __imag__ res = atan2f (rx + r1, copysignf (ix + r2, 108 __imag__ x)); 109 else 110 __imag__ res = atan2f (ix + r2, rx + r1); 111 } 112 } 113 else if (ix == 1.0f && rx < 0.5f) 114 { 115 if (rx < FLT_EPSILON / 8.0f) 116 { 117 __real__ res = log1pf (2.0f * (rx + sqrtf (rx))) / 2.0f; 118 if (adj) 119 __imag__ res = atan2f (sqrtf (rx), 120 copysignf (1.0f, __imag__ x)); 121 else 122 __imag__ res = atan2f (1.0f, sqrtf (rx)); 123 } 124 else 125 { 126 float d = rx * sqrtf (4.0f + rx * rx); 127 float s1 = sqrtf ((d + rx * rx) / 2.0f); 128 float s2 = sqrtf ((d - rx * rx) / 2.0f); 129 130 __real__ res = log1pf (rx * rx + d + 2.0f * (rx * s1 + s2)) / 2.0f; 131 if (adj) 132 __imag__ res = atan2f (rx + s1, 133 copysignf (1.0f + s2, 134 __imag__ x)); 135 else 136 __imag__ res = atan2f (1.0f + s2, rx + s1); 137 } 138 } 139 else if (ix < 1.0f && rx < 0.5f) 140 { 141 if (ix >= FLT_EPSILON) 142 { 143 if (rx < FLT_EPSILON * FLT_EPSILON) 144 { 145 float onemix2 = (1.0f + ix) * (1.0f - ix); 146 float s = sqrtf (onemix2); 147 148 __real__ res = log1pf (2.0f * rx / s) / 2.0f; 149 if (adj) 150 __imag__ res = atan2f (s, __imag__ x); 151 else 152 __imag__ res = atan2f (ix, s); 153 } 154 else 155 { 156 float onemix2 = (1.0f + ix) * (1.0f - ix); 157 float rx2 = rx * rx; 158 float f = rx2 * (2.0f + rx2 + 2.0f * ix * ix); 159 float d = sqrtf (onemix2 * onemix2 + f); 160 float dp = d + onemix2; 161 float dm = f / dp; 162 float r1 = sqrtf ((dp + rx2) / 2.0f); 163 float r2 = rx * ix / r1; 164 165 __real__ res 166 = log1pf (rx2 + dm + 2.0f * (rx * r1 + ix * r2)) / 2.0f; 167 if (adj) 168 __imag__ res = atan2f (rx + r1, 169 copysignf (ix + r2, 170 __imag__ x)); 171 else 172 __imag__ res = atan2f (ix + r2, rx + r1); 173 } 174 } 175 else 176 { 177 float s = hypotf (1.0f, rx); 178 179 __real__ res = log1pf (2.0f * rx * (rx + s)) / 2.0f; 180 if (adj) 181 __imag__ res = atan2f (s, __imag__ x); 182 else 183 __imag__ res = atan2f (ix, s); 184 } 185 if (__real__ res < FLT_MIN) 186 { 187 volatile float force_underflow = __real__ res * __real__ res; 188 (void) force_underflow; 189 } 190 } 191 else 192 { 193 __real__ y = (rx - ix) * (rx + ix) + 1.0f; 194 __imag__ y = 2.0f * rx * ix; 195 196 y = csqrtf (y); 197 198 __real__ y += rx; 199 __imag__ y += ix; 200 201 if (adj) 202 { 203 float t = __real__ y; 204 __real__ y = copysignf (__imag__ y, __imag__ x); 205 __imag__ y = t; 206 } 207 208 res = clogf (y); 209 } 210 211 /* Give results the correct sign for the original argument. */ 212 __real__ res = copysignf (__real__ res, __real__ x); 213 __imag__ res = copysignf (__imag__ res, (adj ? 1.0f : __imag__ x)); 214 215 return res; 216 } 217