1 /* Return arc hyperbole sine for long double value, with the imaginary
2 part of the result possibly adjusted for use in computing other
4 Copyright (C) 1997-2013 Free Software Foundation, Inc.
5 This file is part of the GNU C Library.
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.
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.
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/>. */
23 #include <math_private.h>
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
30 # define LDBL_EPSILON 0x1p-106L
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
37 __complex__
long double
38 __kernel_casinhl (__complex__
long double x
, int adj
)
40 __complex__
long double res
;
42 __complex__
long double y
;
44 /* Avoid cancellation by reducing to the first quadrant. */
45 rx
= fabsl (__real__ x
);
46 ix
= fabsl (__imag__ x
);
48 if (rx
>= 1.0L / LDBL_EPSILON
|| ix
>= 1.0L / LDBL_EPSILON
)
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. */
59 long double t
= __real__ y
;
60 __real__ y
= __copysignl (__imag__ y
, __imag__ x
);
65 __real__ res
+= M_LN2l
;
67 else if (rx
>= 0.5L && ix
< LDBL_EPSILON
/ 8.0L)
69 long double s
= __ieee754_hypotl (1.0L, rx
);
71 __real__ res
= __ieee754_logl (rx
+ s
);
73 __imag__ res
= __ieee754_atan2l (s
, __imag__ x
);
75 __imag__ res
= __ieee754_atan2l (ix
, s
);
77 else if (rx
< LDBL_EPSILON
/ 8.0L && ix
>= 1.5L)
79 long double s
= __ieee754_sqrtl ((ix
+ 1.0L) * (ix
- 1.0L));
81 __real__ res
= __ieee754_logl (ix
+ s
);
83 __imag__ res
= __ieee754_atan2l (rx
, __copysignl (s
, __imag__ x
));
85 __imag__ res
= __ieee754_atan2l (s
, rx
);
87 else if (ix
> 1.0L && ix
< 1.5L && rx
< 0.5L)
89 if (rx
< LDBL_EPSILON
* LDBL_EPSILON
)
91 long double ix2m1
= (ix
+ 1.0L) * (ix
- 1.0L);
92 long double s
= __ieee754_sqrtl (ix2m1
);
94 __real__ res
= __log1pl (2.0L * (ix2m1
+ ix
* s
)) / 2.0L;
96 __imag__ res
= __ieee754_atan2l (rx
, __copysignl (s
, __imag__ x
));
98 __imag__ res
= __ieee754_atan2l (s
, rx
);
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
;
112 = __log1pl (rx2
+ dp
+ 2.0L * (rx
* r1
+ ix
* r2
)) / 2.0L;
114 __imag__ res
= __ieee754_atan2l (rx
+ r1
, __copysignl (ix
+ r2
,
117 __imag__ res
= __ieee754_atan2l (ix
+ r2
, rx
+ r1
);
120 else if (ix
== 1.0L && rx
< 0.5L)
122 if (rx
< LDBL_EPSILON
/ 8.0L)
124 __real__ res
= __log1pl (2.0L * (rx
+ __ieee754_sqrtl (rx
))) / 2.0L;
126 __imag__ res
= __ieee754_atan2l (__ieee754_sqrtl (rx
),
127 __copysignl (1.0L, __imag__ x
));
129 __imag__ res
= __ieee754_atan2l (1.0L, __ieee754_sqrtl (rx
));
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);
137 __real__ res
= __log1pl (rx
* rx
+ d
+ 2.0L * (rx
* s1
+ s2
)) / 2.0L;
139 __imag__ res
= __ieee754_atan2l (rx
+ s1
,
140 __copysignl (1.0L + s2
,
143 __imag__ res
= __ieee754_atan2l (1.0L + s2
, rx
+ s1
);
146 else if (ix
< 1.0L && rx
< 0.5L)
148 if (ix
>= LDBL_EPSILON
)
150 if (rx
< LDBL_EPSILON
* LDBL_EPSILON
)
152 long double onemix2
= (1.0L + ix
) * (1.0L - ix
);
153 long double s
= __ieee754_sqrtl (onemix2
);
155 __real__ res
= __log1pl (2.0L * rx
/ s
) / 2.0L;
157 __imag__ res
= __ieee754_atan2l (s
, __imag__ x
);
159 __imag__ res
= __ieee754_atan2l (ix
, s
);
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
;
173 = __log1pl (rx2
+ dm
+ 2.0L * (rx
* r1
+ ix
* r2
)) / 2.0L;
175 __imag__ res
= __ieee754_atan2l (rx
+ r1
,
176 __copysignl (ix
+ r2
,
179 __imag__ res
= __ieee754_atan2l (ix
+ r2
, rx
+ r1
);
184 long double s
= __ieee754_hypotl (1.0L, rx
);
186 __real__ res
= __log1pl (2.0L * rx
* (rx
+ s
)) / 2.0L;
188 __imag__ res
= __ieee754_atan2l (s
, __imag__ x
);
190 __imag__ res
= __ieee754_atan2l (ix
, s
);
192 if (__real__ res
< LDBL_MIN
)
194 volatile long double force_underflow
= __real__ res
* __real__ res
;
195 (void) force_underflow
;
200 __real__ y
= (rx
- ix
) * (rx
+ ix
) + 1.0L;
201 __imag__ y
= 2.0L * rx
* ix
;
210 long double t
= __real__ y
;
211 __real__ y
= __copysignl (__imag__ y
, __imag__ x
);
218 /* Give results the correct sign for the original argument. */
219 __real__ res
= __copysignl (__real__ res
, __real__ x
);
220 __imag__ res
= __copysignl (__imag__ res
, (adj
? 1.0L : __imag__ x
));