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-2016 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 = __ieee754_hypotf (1.0f, rx); |
63 | |
64 | __real__ res = __ieee754_logf (rx + s); |
65 | if (adj) |
66 | __imag__ res = __ieee754_atan2f (s, __imag__ x); |
67 | else |
68 | __imag__ res = __ieee754_atan2f (ix, s); |
69 | } |
70 | else if (rx < FLT_EPSILON / 8.0f && ix >= 1.5f) |
71 | { |
72 | float s = __ieee754_sqrtf ((ix + 1.0f) * (ix - 1.0f)); |
73 | |
74 | __real__ res = __ieee754_logf (ix + s); |
75 | if (adj) |
76 | __imag__ res = __ieee754_atan2f (rx, __copysignf (s, __imag__ x)); |
77 | else |
78 | __imag__ res = __ieee754_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 = __ieee754_sqrtf (ix2m1); |
86 | |
87 | __real__ res = __log1pf (2.0f * (ix2m1 + ix * s)) / 2.0f; |
88 | if (adj) |
89 | __imag__ res = __ieee754_atan2f (rx, __copysignf (s, __imag__ x)); |
90 | else |
91 | __imag__ res = __ieee754_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 = __ieee754_sqrtf (ix2m1 * ix2m1 + f); |
99 | float dp = d + ix2m1; |
100 | float dm = f / dp; |
101 | float r1 = __ieee754_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 = __ieee754_atan2f (rx + r1, __copysignf (ix + r2, |
108 | __imag__ x)); |
109 | else |
110 | __imag__ res = __ieee754_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 + __ieee754_sqrtf (rx))) / 2.0f; |
118 | if (adj) |
119 | __imag__ res = __ieee754_atan2f (__ieee754_sqrtf (rx), |
120 | __copysignf (1.0f, __imag__ x)); |
121 | else |
122 | __imag__ res = __ieee754_atan2f (1.0f, __ieee754_sqrtf (rx)); |
123 | } |
124 | else |
125 | { |
126 | float d = rx * __ieee754_sqrtf (4.0f + rx * rx); |
127 | float s1 = __ieee754_sqrtf ((d + rx * rx) / 2.0f); |
128 | float s2 = __ieee754_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 = __ieee754_atan2f (rx + s1, |
133 | __copysignf (1.0f + s2, |
134 | __imag__ x)); |
135 | else |
136 | __imag__ res = __ieee754_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 = __ieee754_sqrtf (onemix2); |
147 | |
148 | __real__ res = __log1pf (2.0f * rx / s) / 2.0f; |
149 | if (adj) |
150 | __imag__ res = __ieee754_atan2f (s, __imag__ x); |
151 | else |
152 | __imag__ res = __ieee754_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 = __ieee754_sqrtf (onemix2 * onemix2 + f); |
160 | float dp = d + onemix2; |
161 | float dm = f / dp; |
162 | float r1 = __ieee754_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 = __ieee754_atan2f (rx + r1, |
169 | __copysignf (ix + r2, |
170 | __imag__ x)); |
171 | else |
172 | __imag__ res = __ieee754_atan2f (ix + r2, rx + r1); |
173 | } |
174 | } |
175 | else |
176 | { |
177 | float s = __ieee754_hypotf (1.0f, rx); |
178 | |
179 | __real__ res = __log1pf (2.0f * rx * (rx + s)) / 2.0f; |
180 | if (adj) |
181 | __imag__ res = __ieee754_atan2f (s, __imag__ x); |
182 | else |
183 | __imag__ res = __ieee754_atan2f (ix, s); |
184 | } |
185 | math_check_force_underflow_nonneg (__real__ res); |
186 | } |
187 | else |
188 | { |
189 | __real__ y = (rx - ix) * (rx + ix) + 1.0f; |
190 | __imag__ y = 2.0f * rx * ix; |
191 | |
192 | y = __csqrtf (y); |
193 | |
194 | __real__ y += rx; |
195 | __imag__ y += ix; |
196 | |
197 | if (adj) |
198 | { |
199 | float t = __real__ y; |
200 | __real__ y = __copysignf (__imag__ y, __imag__ x); |
201 | __imag__ y = t; |
202 | } |
203 | |
204 | res = __clogf (y); |
205 | } |
206 | |
207 | /* Give results the correct sign for the original argument. */ |
208 | __real__ res = __copysignf (__real__ res, __real__ x); |
209 | __imag__ res = __copysignf (__imag__ res, (adj ? 1.0f : __imag__ x)); |
210 | |
211 | return res; |
212 | } |
213 | |