]>
Commit | Line | Data |
---|---|---|
728d7b43 JM |
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. | |
d4697bc9 | 4 | Copyright (C) 1997-2014 Free Software Foundation, Inc. |
728d7b43 JM |
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 | } | |
8cf28c5e JM |
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 | } | |
3a7182a1 JM |
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 | } | |
0a1b2ae6 JM |
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 | } | |
ccc8cadf JM |
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 | if (__real__ res < FLT_MIN) | |
186 | { | |
187 | volatile float force_underflow = __real__ res * __real__ res; | |
188 | (void) force_underflow; | |
189 | } | |
190 | } | |
728d7b43 JM |
191 | else |
192 | { | |
6b18bea6 JM |
193 | __real__ y = (rx - ix) * (rx + ix) + 1.0f; |
194 | __imag__ y = 2.0f * rx * ix; | |
728d7b43 JM |
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 | } |