]>
Commit | Line | Data |
---|---|---|
728d7b43 JM |
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. | |
f7a9f785 | 4 | Copyright (C) 1997-2016 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 | /* 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 | } | |
8cf28c5e JM |
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 | } | |
3a7182a1 JM |
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 | } | |
0a1b2ae6 JM |
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 | } | |
ccc8cadf JM |
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 | } | |
d96164c3 | 192 | math_check_force_underflow_nonneg (__real__ res); |
ccc8cadf | 193 | } |
728d7b43 JM |
194 | else |
195 | { | |
6b18bea6 JM |
196 | __real__ y = (rx - ix) * (rx + ix) + 1.0L; |
197 | __imag__ y = 2.0L * rx * ix; | |
728d7b43 JM |
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 | } |