Mercurial > hg > CbC > CbC_gcc
comparison libquadmath/math/casinhq_kernel.c @ 145:1830386684a0
gcc-9.2.0
author | anatofuz |
---|---|
date | Thu, 13 Feb 2020 11:34:05 +0900 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
131:84e7813d76e9 | 145:1830386684a0 |
---|---|
1 /* Return arc hyperbolic sine for a complex float type, with the | |
2 imaginary part of the result possibly adjusted for use in | |
3 computing other functions. | |
4 Copyright (C) 1997-2018 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 "quadmath-imp.h" | |
22 | |
23 /* Return the complex inverse hyperbolic sine of finite nonzero Z, | |
24 with the imaginary part of the result subtracted from pi/2 if ADJ | |
25 is nonzero. */ | |
26 | |
27 __complex128 | |
28 __quadmath_kernel_casinhq (__complex128 x, int adj) | |
29 { | |
30 __complex128 res; | |
31 __float128 rx, ix; | |
32 __complex128 y; | |
33 | |
34 /* Avoid cancellation by reducing to the first quadrant. */ | |
35 rx = fabsq (__real__ x); | |
36 ix = fabsq (__imag__ x); | |
37 | |
38 if (rx >= 1 / FLT128_EPSILON || ix >= 1 / FLT128_EPSILON) | |
39 { | |
40 /* For large x in the first quadrant, x + csqrt (1 + x * x) | |
41 is sufficiently close to 2 * x to make no significant | |
42 difference to the result; avoid possible overflow from | |
43 the squaring and addition. */ | |
44 __real__ y = rx; | |
45 __imag__ y = ix; | |
46 | |
47 if (adj) | |
48 { | |
49 __float128 t = __real__ y; | |
50 __real__ y = copysignq (__imag__ y, __imag__ x); | |
51 __imag__ y = t; | |
52 } | |
53 | |
54 res = clogq (y); | |
55 __real__ res += (__float128) M_LN2q; | |
56 } | |
57 else if (rx >= 0.5Q && ix < FLT128_EPSILON / 8) | |
58 { | |
59 __float128 s = hypotq (1, rx); | |
60 | |
61 __real__ res = logq (rx + s); | |
62 if (adj) | |
63 __imag__ res = atan2q (s, __imag__ x); | |
64 else | |
65 __imag__ res = atan2q (ix, s); | |
66 } | |
67 else if (rx < FLT128_EPSILON / 8 && ix >= 1.5Q) | |
68 { | |
69 __float128 s = sqrtq ((ix + 1) * (ix - 1)); | |
70 | |
71 __real__ res = logq (ix + s); | |
72 if (adj) | |
73 __imag__ res = atan2q (rx, copysignq (s, __imag__ x)); | |
74 else | |
75 __imag__ res = atan2q (s, rx); | |
76 } | |
77 else if (ix > 1 && ix < 1.5Q && rx < 0.5Q) | |
78 { | |
79 if (rx < FLT128_EPSILON * FLT128_EPSILON) | |
80 { | |
81 __float128 ix2m1 = (ix + 1) * (ix - 1); | |
82 __float128 s = sqrtq (ix2m1); | |
83 | |
84 __real__ res = log1pq (2 * (ix2m1 + ix * s)) / 2; | |
85 if (adj) | |
86 __imag__ res = atan2q (rx, copysignq (s, __imag__ x)); | |
87 else | |
88 __imag__ res = atan2q (s, rx); | |
89 } | |
90 else | |
91 { | |
92 __float128 ix2m1 = (ix + 1) * (ix - 1); | |
93 __float128 rx2 = rx * rx; | |
94 __float128 f = rx2 * (2 + rx2 + 2 * ix * ix); | |
95 __float128 d = sqrtq (ix2m1 * ix2m1 + f); | |
96 __float128 dp = d + ix2m1; | |
97 __float128 dm = f / dp; | |
98 __float128 r1 = sqrtq ((dm + rx2) / 2); | |
99 __float128 r2 = rx * ix / r1; | |
100 | |
101 __real__ res = log1pq (rx2 + dp + 2 * (rx * r1 + ix * r2)) / 2; | |
102 if (adj) | |
103 __imag__ res = atan2q (rx + r1, copysignq (ix + r2, __imag__ x)); | |
104 else | |
105 __imag__ res = atan2q (ix + r2, rx + r1); | |
106 } | |
107 } | |
108 else if (ix == 1 && rx < 0.5Q) | |
109 { | |
110 if (rx < FLT128_EPSILON / 8) | |
111 { | |
112 __real__ res = log1pq (2 * (rx + sqrtq (rx))) / 2; | |
113 if (adj) | |
114 __imag__ res = atan2q (sqrtq (rx), copysignq (1, __imag__ x)); | |
115 else | |
116 __imag__ res = atan2q (1, sqrtq (rx)); | |
117 } | |
118 else | |
119 { | |
120 __float128 d = rx * sqrtq (4 + rx * rx); | |
121 __float128 s1 = sqrtq ((d + rx * rx) / 2); | |
122 __float128 s2 = sqrtq ((d - rx * rx) / 2); | |
123 | |
124 __real__ res = log1pq (rx * rx + d + 2 * (rx * s1 + s2)) / 2; | |
125 if (adj) | |
126 __imag__ res = atan2q (rx + s1, copysignq (1 + s2, __imag__ x)); | |
127 else | |
128 __imag__ res = atan2q (1 + s2, rx + s1); | |
129 } | |
130 } | |
131 else if (ix < 1 && rx < 0.5Q) | |
132 { | |
133 if (ix >= FLT128_EPSILON) | |
134 { | |
135 if (rx < FLT128_EPSILON * FLT128_EPSILON) | |
136 { | |
137 __float128 onemix2 = (1 + ix) * (1 - ix); | |
138 __float128 s = sqrtq (onemix2); | |
139 | |
140 __real__ res = log1pq (2 * rx / s) / 2; | |
141 if (adj) | |
142 __imag__ res = atan2q (s, __imag__ x); | |
143 else | |
144 __imag__ res = atan2q (ix, s); | |
145 } | |
146 else | |
147 { | |
148 __float128 onemix2 = (1 + ix) * (1 - ix); | |
149 __float128 rx2 = rx * rx; | |
150 __float128 f = rx2 * (2 + rx2 + 2 * ix * ix); | |
151 __float128 d = sqrtq (onemix2 * onemix2 + f); | |
152 __float128 dp = d + onemix2; | |
153 __float128 dm = f / dp; | |
154 __float128 r1 = sqrtq ((dp + rx2) / 2); | |
155 __float128 r2 = rx * ix / r1; | |
156 | |
157 __real__ res = log1pq (rx2 + dm + 2 * (rx * r1 + ix * r2)) / 2; | |
158 if (adj) | |
159 __imag__ res = atan2q (rx + r1, copysignq (ix + r2, | |
160 __imag__ x)); | |
161 else | |
162 __imag__ res = atan2q (ix + r2, rx + r1); | |
163 } | |
164 } | |
165 else | |
166 { | |
167 __float128 s = hypotq (1, rx); | |
168 | |
169 __real__ res = log1pq (2 * rx * (rx + s)) / 2; | |
170 if (adj) | |
171 __imag__ res = atan2q (s, __imag__ x); | |
172 else | |
173 __imag__ res = atan2q (ix, s); | |
174 } | |
175 math_check_force_underflow_nonneg (__real__ res); | |
176 } | |
177 else | |
178 { | |
179 __real__ y = (rx - ix) * (rx + ix) + 1; | |
180 __imag__ y = 2 * rx * ix; | |
181 | |
182 y = csqrtq (y); | |
183 | |
184 __real__ y += rx; | |
185 __imag__ y += ix; | |
186 | |
187 if (adj) | |
188 { | |
189 __float128 t = __real__ y; | |
190 __real__ y = copysignq (__imag__ y, __imag__ x); | |
191 __imag__ y = t; | |
192 } | |
193 | |
194 res = clogq (y); | |
195 } | |
196 | |
197 /* Give results the correct sign for the original argument. */ | |
198 __real__ res = copysignq (__real__ res, __real__ x); | |
199 __imag__ res = copysignq (__imag__ res, (adj ? 1 : __imag__ x)); | |
200 | |
201 return res; | |
202 } |