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-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__ 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 = __ieee754_hypot (1.0, rx); |
63 | |
64 | __real__ res = __ieee754_log (rx + s); |
65 | if (adj) |
66 | __imag__ res = __ieee754_atan2 (s, __imag__ x); |
67 | else |
68 | __imag__ res = __ieee754_atan2 (ix, s); |
69 | } |
70 | else if (rx < DBL_EPSILON / 8.0 && ix >= 1.5) |
71 | { |
72 | double s = __ieee754_sqrt ((ix + 1.0) * (ix - 1.0)); |
73 | |
74 | __real__ res = __ieee754_log (ix + s); |
75 | if (adj) |
76 | __imag__ res = __ieee754_atan2 (rx, __copysign (s, __imag__ x)); |
77 | else |
78 | __imag__ res = __ieee754_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 = __ieee754_sqrt (ix2m1); |
86 | |
87 | __real__ res = __log1p (2.0 * (ix2m1 + ix * s)) / 2.0; |
88 | if (adj) |
89 | __imag__ res = __ieee754_atan2 (rx, __copysign (s, __imag__ x)); |
90 | else |
91 | __imag__ res = __ieee754_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 = __ieee754_sqrt (ix2m1 * ix2m1 + f); |
99 | double dp = d + ix2m1; |
100 | double dm = f / dp; |
101 | double r1 = __ieee754_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 = __ieee754_atan2 (rx + r1, __copysign (ix + r2, |
107 | __imag__ x)); |
108 | else |
109 | __imag__ res = __ieee754_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 + __ieee754_sqrt (rx))) / 2.0; |
117 | if (adj) |
118 | __imag__ res = __ieee754_atan2 (__ieee754_sqrt (rx), |
119 | __copysign (1.0, __imag__ x)); |
120 | else |
121 | __imag__ res = __ieee754_atan2 (1.0, __ieee754_sqrt (rx)); |
122 | } |
123 | else |
124 | { |
125 | double d = rx * __ieee754_sqrt (4.0 + rx * rx); |
126 | double s1 = __ieee754_sqrt ((d + rx * rx) / 2.0); |
127 | double s2 = __ieee754_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 = __ieee754_atan2 (rx + s1, __copysign (1.0 + s2, |
132 | __imag__ x)); |
133 | else |
134 | __imag__ res = __ieee754_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 = __ieee754_sqrt (onemix2); |
145 | |
146 | __real__ res = __log1p (2.0 * rx / s) / 2.0; |
147 | if (adj) |
148 | __imag__ res = __ieee754_atan2 (s, __imag__ x); |
149 | else |
150 | __imag__ res = __ieee754_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 = __ieee754_sqrt (onemix2 * onemix2 + f); |
158 | double dp = d + onemix2; |
159 | double dm = f / dp; |
160 | double r1 = __ieee754_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 = __ieee754_atan2 (rx + r1, |
167 | __copysign (ix + r2, |
168 | __imag__ x)); |
169 | else |
170 | __imag__ res = __ieee754_atan2 (ix + r2, rx + r1); |
171 | } |
172 | } |
173 | else |
174 | { |
175 | double s = __ieee754_hypot (1.0, rx); |
176 | |
177 | __real__ res = __log1p (2.0 * rx * (rx + s)) / 2.0; |
178 | if (adj) |
179 | __imag__ res = __ieee754_atan2 (s, __imag__ x); |
180 | else |
181 | __imag__ res = __ieee754_atan2 (ix, s); |
182 | } |
183 | math_check_force_underflow_nonneg (__real__ res); |
184 | } |
185 | else |
186 | { |
187 | __real__ y = (rx - ix) * (rx + ix) + 1.0; |
188 | __imag__ y = 2.0 * rx * ix; |
189 | |
190 | y = __csqrt (y); |
191 | |
192 | __real__ y += rx; |
193 | __imag__ y += ix; |
194 | |
195 | if (adj) |
196 | { |
197 | double t = __real__ y; |
198 | __real__ y = __copysign (__imag__ y, __imag__ x); |
199 | __imag__ y = t; |
200 | } |
201 | |
202 | res = __clog (y); |
203 | } |
204 | |
205 | /* Give results the correct sign for the original argument. */ |
206 | __real__ res = __copysign (__real__ res, __real__ x); |
207 | __imag__ res = __copysign (__imag__ res, (adj ? 1.0 : __imag__ x)); |
208 | |
209 | return res; |
210 | } |
211 | |