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