1 /*-
2 * SPDX-License-Identifier: BSD-2-Clause
3 *
4 * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
9 * are met:
10 * 1. Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 *
16 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
17 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
20 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
21 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
22 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
23 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
24 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
25 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
26 * SUCH DAMAGE.
27 */
28
29 #include <sys/cdefs.h>
30 __FBSDID("$FreeBSD$");
31
32 #include <complex.h>
33 #include <float.h>
34
35 #include "math.h"
36 #include "math_private.h"
37
38 #undef isinf
39 #define isinf(x) (fabs(x) == INFINITY)
40 #undef isnan
41 #define isnan(x) ((x) != (x))
42 #define raise_inexact() do { volatile float junk __unused = 1 + tiny; } while(0)
43 #undef signbit
44 #define signbit(x) (__builtin_signbit(x))
45
46 /* We need that DBL_EPSILON^2/128 is larger than FOUR_SQRT_MIN. */
47 static const double
48 A_crossover = 10, /* Hull et al suggest 1.5, but 10 works better */
49 B_crossover = 0.6417, /* suggested by Hull et al */
50 FOUR_SQRT_MIN = 0x1p-509, /* >= 4 * sqrt(DBL_MIN) */
51 QUARTER_SQRT_MAX = 0x1p509, /* <= sqrt(DBL_MAX) / 4 */
52 m_e = 2.7182818284590452e0, /* 0x15bf0a8b145769.0p-51 */
53 m_ln2 = 6.9314718055994531e-1, /* 0x162e42fefa39ef.0p-53 */
54 pio2_hi = 1.5707963267948966e0, /* 0x1921fb54442d18.0p-52 */
55 RECIP_EPSILON = 1 / DBL_EPSILON,
56 SQRT_3_EPSILON = 2.5809568279517849e-8, /* 0x1bb67ae8584caa.0p-78 */
57 SQRT_6_EPSILON = 3.6500241499888571e-8, /* 0x13988e1409212e.0p-77 */
58 SQRT_MIN = 0x1p-511; /* >= sqrt(DBL_MIN) */
59
60 static const volatile double
61 pio2_lo = 6.1232339957367659e-17; /* 0x11a62633145c07.0p-106 */
62 static const volatile float
63 tiny = 0x1p-100;
64
65 static double complex clog_for_large_values(double complex z);
66
67 /*
68 * Testing indicates that all these functions are accurate up to 4 ULP.
69 * The functions casin(h) and cacos(h) are about 2.5 times slower than asinh.
70 * The functions catan(h) are a little under 2 times slower than atanh.
71 *
72 * The code for casinh, casin, cacos, and cacosh comes first. The code is
73 * rather complicated, and the four functions are highly interdependent.
74 *
75 * The code for catanh and catan comes at the end. It is much simpler than
76 * the other functions, and the code for these can be disconnected from the
77 * rest of the code.
78 */
79
80 /*
81 * ================================
82 * | casinh, casin, cacos, cacosh |
83 * ================================
84 */
85
86 /*
87 * The algorithm is very close to that in "Implementing the complex arcsine
88 * and arccosine functions using exception handling" by T. E. Hull, Thomas F.
89 * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on
90 * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
91 * http://dl.acm.org/citation.cfm?id=275324.
92 *
93 * Throughout we use the convention z = x + I*y.
94 *
95 * casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B)
96 * where
97 * A = (|z+I| + |z-I|) / 2
98 * B = (|z+I| - |z-I|) / 2 = y/A
99 *
100 * These formulas become numerically unstable:
101 * (a) for Re(casinh(z)) when z is close to the line segment [-I, I] (that
102 * is, Re(casinh(z)) is close to 0);
103 * (b) for Im(casinh(z)) when z is close to either of the intervals
104 * [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is
105 * close to PI/2).
106 *
107 * These numerical problems are overcome by defining
108 * f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2
109 * Then if A < A_crossover, we use
110 * log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1)))
111 * A-1 = f(x, 1+y) + f(x, 1-y)
112 * and if B > B_crossover, we use
113 * asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y)))
114 * A-y = f(x, y+1) + f(x, y-1)
115 * where without loss of generality we have assumed that x and y are
116 * non-negative.
117 *
118 * Much of the difficulty comes because the intermediate computations may
119 * produce overflows or underflows. This is dealt with in the paper by Hull
120 * et al by using exception handling. We do this by detecting when
121 * computations risk underflow or overflow. The hardest part is handling the
122 * underflows when computing f(a, b).
123 *
124 * Note that the function f(a, b) does not appear explicitly in the paper by
125 * Hull et al, but the idea may be found on pages 308 and 309. Introducing the
126 * function f(a, b) allows us to concentrate many of the clever tricks in this
127 * paper into one function.
128 */
129
130 /*
131 * Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2.
132 * Pass hypot(a, b) as the third argument.
133 */
134 static inline double
f(double a,double b,double hypot_a_b)135 f(double a, double b, double hypot_a_b)
136 {
137 if (b < 0)
138 return ((hypot_a_b - b) / 2);
139 if (b == 0)
140 return (a / 2);
141 return (a * a / (hypot_a_b + b) / 2);
142 }
143
144 /*
145 * All the hard work is contained in this function.
146 * x and y are assumed positive or zero, and less than RECIP_EPSILON.
147 * Upon return:
148 * rx = Re(casinh(z)) = -Im(cacos(y + I*x)).
149 * B_is_usable is set to 1 if the value of B is usable.
150 * If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y.
151 * If returning sqrt_A2my2 has potential to result in an underflow, it is
152 * rescaled, and new_y is similarly rescaled.
153 */
154 static inline void
do_hard_work(double x,double y,double * rx,int * B_is_usable,double * B,double * sqrt_A2my2,double * new_y)155 do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
156 double *sqrt_A2my2, double *new_y)
157 {
158 double R, S, A; /* A, B, R, and S are as in Hull et al. */
159 double Am1, Amy; /* A-1, A-y. */
160
161 R = hypot(x, y + 1); /* |z+I| */
162 S = hypot(x, y - 1); /* |z-I| */
163
164 /* A = (|z+I| + |z-I|) / 2 */
165 A = (R + S) / 2;
166 /*
167 * Mathematically A >= 1. There is a small chance that this will not
168 * be so because of rounding errors. So we will make certain it is
169 * so.
170 */
171 if (A < 1)
172 A = 1;
173
174 if (A < A_crossover) {
175 /*
176 * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y).
177 * rx = log1p(Am1 + sqrt(Am1*(A+1)))
178 */
179 if (y == 1 && x < DBL_EPSILON * DBL_EPSILON / 128) {
180 /*
181 * fp is of order x^2, and fm = x/2.
182 * A = 1 (inexactly).
183 */
184 *rx = sqrt(x);
185 } else if (x >= DBL_EPSILON * fabs(y - 1)) {
186 /*
187 * Underflow will not occur because
188 * x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN
189 */
190 Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
191 *rx = log1p(Am1 + sqrt(Am1 * (A + 1)));
192 } else if (y < 1) {
193 /*
194 * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and
195 * A = 1 (inexactly).
196 */
197 *rx = x / sqrt((1 - y) * (1 + y));
198 } else { /* if (y > 1) */
199 /*
200 * A-1 = y-1 (inexactly).
201 */
202 *rx = log1p((y - 1) + sqrt((y - 1) * (y + 1)));
203 }
204 } else {
205 *rx = log(A + sqrt(A * A - 1));
206 }
207
208 *new_y = y;
209
210 if (y < FOUR_SQRT_MIN) {
211 /*
212 * Avoid a possible underflow caused by y/A. For casinh this
213 * would be legitimate, but will be picked up by invoking atan2
214 * later on. For cacos this would not be legitimate.
215 */
216 *B_is_usable = 0;
217 *sqrt_A2my2 = A * (2 / DBL_EPSILON);
218 *new_y = y * (2 / DBL_EPSILON);
219 return;
220 }
221
222 /* B = (|z+I| - |z-I|) / 2 = y/A */
223 *B = y / A;
224 *B_is_usable = 1;
225
226 if (*B > B_crossover) {
227 *B_is_usable = 0;
228 /*
229 * Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1).
230 * sqrt_A2my2 = sqrt(Amy*(A+y))
231 */
232 if (y == 1 && x < DBL_EPSILON / 128) {
233 /*
234 * fp is of order x^2, and fm = x/2.
235 * A = 1 (inexactly).
236 */
237 *sqrt_A2my2 = sqrt(x) * sqrt((A + y) / 2);
238 } else if (x >= DBL_EPSILON * fabs(y - 1)) {
239 /*
240 * Underflow will not occur because
241 * x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN
242 * and
243 * x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN
244 */
245 Amy = f(x, y + 1, R) + f(x, y - 1, S);
246 *sqrt_A2my2 = sqrt(Amy * (A + y));
247 } else if (y > 1) {
248 /*
249 * fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and
250 * A = y (inexactly).
251 *
252 * y < RECIP_EPSILON. So the following
253 * scaling should avoid any underflow problems.
254 */
255 *sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y /
256 sqrt((y + 1) * (y - 1));
257 *new_y = y * (4 / DBL_EPSILON / DBL_EPSILON);
258 } else { /* if (y < 1) */
259 /*
260 * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and
261 * A = 1 (inexactly).
262 */
263 *sqrt_A2my2 = sqrt((1 - y) * (1 + y));
264 }
265 }
266 }
267
268 /*
269 * casinh(z) = z + O(z^3) as z -> 0
270 *
271 * casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2) as z -> infinity
272 * The above formula works for the imaginary part as well, because
273 * Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3)
274 * as z -> infinity, uniformly in y
275 */
276 double complex
casinh(double complex z)277 casinh(double complex z)
278 {
279 double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;
280 int B_is_usable;
281 double complex w;
282
283 x = creal(z);
284 y = cimag(z);
285 ax = fabs(x);
286 ay = fabs(y);
287
288 if (isnan(x) || isnan(y)) {
289 /* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */
290 if (isinf(x))
291 return (CMPLX(x, y + y));
292 /* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
293 if (isinf(y))
294 return (CMPLX(y, x + x));
295 /* casinh(NaN + I*0) = NaN + I*0 */
296 if (y == 0)
297 return (CMPLX(x + x, y));
298 /*
299 * All other cases involving NaN return NaN + I*NaN.
300 * C99 leaves it optional whether to raise invalid if one of
301 * the arguments is not NaN, so we opt not to raise it.
302 */
303 return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
304 }
305
306 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
307 /* clog...() will raise inexact unless x or y is infinite. */
308 if (signbit(x) == 0)
309 w = clog_for_large_values(z) + m_ln2;
310 else
311 w = clog_for_large_values(-z) + m_ln2;
312 return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y)));
313 }
314
315 /* Avoid spuriously raising inexact for z = 0. */
316 if (x == 0 && y == 0)
317 return (z);
318
319 /* All remaining cases are inexact. */
320 raise_inexact();
321
322 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
323 return (z);
324
325 do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
326 if (B_is_usable)
327 ry = asin(B);
328 else
329 ry = atan2(new_y, sqrt_A2my2);
330 return (CMPLX(copysign(rx, x), copysign(ry, y)));
331 }
332
333 /*
334 * casin(z) = reverse(casinh(reverse(z)))
335 * where reverse(x + I*y) = y + I*x = I*conj(z).
336 */
337 double complex
casin(double complex z)338 casin(double complex z)
339 {
340 double complex w = casinh(CMPLX(cimag(z), creal(z)));
341
342 return (CMPLX(cimag(w), creal(w)));
343 }
344
345 /*
346 * cacos(z) = PI/2 - casin(z)
347 * but do the computation carefully so cacos(z) is accurate when z is
348 * close to 1.
349 *
350 * cacos(z) = PI/2 - z + O(z^3) as z -> 0
351 *
352 * cacos(z) = -sign(y)*I*clog(z) + O(1/z^2) as z -> infinity
353 * The above formula works for the real part as well, because
354 * Re(cacos(z)) = atan2(fabs(y), x) + O(y/z^3)
355 * as z -> infinity, uniformly in y
356 */
357 double complex
cacos(double complex z)358 cacos(double complex z)
359 {
360 double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
361 int sx, sy;
362 int B_is_usable;
363 double complex w;
364
365 x = creal(z);
366 y = cimag(z);
367 sx = signbit(x);
368 sy = signbit(y);
369 ax = fabs(x);
370 ay = fabs(y);
371
372 if (isnan(x) || isnan(y)) {
373 /* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
374 if (isinf(x))
375 return (CMPLX(y + y, -INFINITY));
376 /* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */
377 if (isinf(y))
378 return (CMPLX(x + x, -y));
379 /* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */
380 if (x == 0)
381 return (CMPLX(pio2_hi + pio2_lo, y + y));
382 /*
383 * All other cases involving NaN return NaN + I*NaN.
384 * C99 leaves it optional whether to raise invalid if one of
385 * the arguments is not NaN, so we opt not to raise it.
386 */
387 return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
388 }
389
390 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
391 /* clog...() will raise inexact unless x or y is infinite. */
392 w = clog_for_large_values(z);
393 rx = fabs(cimag(w));
394 ry = creal(w) + m_ln2;
395 if (sy == 0)
396 ry = -ry;
397 return (CMPLX(rx, ry));
398 }
399
400 /* Avoid spuriously raising inexact for z = 1. */
401 if (x == 1 && y == 0)
402 return (CMPLX(0, -y));
403
404 /* All remaining cases are inexact. */
405 raise_inexact();
406
407 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
408 return (CMPLX(pio2_hi - (x - pio2_lo), -y));
409
410 do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
411 if (B_is_usable) {
412 if (sx == 0)
413 rx = acos(B);
414 else
415 rx = acos(-B);
416 } else {
417 if (sx == 0)
418 rx = atan2(sqrt_A2mx2, new_x);
419 else
420 rx = atan2(sqrt_A2mx2, -new_x);
421 }
422 if (sy == 0)
423 ry = -ry;
424 return (CMPLX(rx, ry));
425 }
426
427 /*
428 * cacosh(z) = I*cacos(z) or -I*cacos(z)
429 * where the sign is chosen so Re(cacosh(z)) >= 0.
430 */
431 double complex
cacosh(double complex z)432 cacosh(double complex z)
433 {
434 double complex w;
435 double rx, ry;
436
437 w = cacos(z);
438 rx = creal(w);
439 ry = cimag(w);
440 /* cacosh(NaN + I*NaN) = NaN + I*NaN */
441 if (isnan(rx) && isnan(ry))
442 return (CMPLX(ry, rx));
443 /* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */
444 /* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */
445 if (isnan(rx))
446 return (CMPLX(fabs(ry), rx));
447 /* cacosh(0 + I*NaN) = NaN + I*NaN */
448 if (isnan(ry))
449 return (CMPLX(ry, ry));
450 return (CMPLX(fabs(ry), copysign(rx, cimag(z))));
451 }
452
453 /*
454 * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON.
455 */
456 static double complex
clog_for_large_values(double complex z)457 clog_for_large_values(double complex z)
458 {
459 double x, y;
460 double ax, ay, t;
461
462 x = creal(z);
463 y = cimag(z);
464 ax = fabs(x);
465 ay = fabs(y);
466 if (ax < ay) {
467 t = ax;
468 ax = ay;
469 ay = t;
470 }
471
472 /*
473 * Avoid overflow in hypot() when x and y are both very large.
474 * Divide x and y by E, and then add 1 to the logarithm. This
475 * depends on E being larger than sqrt(2), since the return value of
476 * hypot cannot overflow if neither argument is greater in magnitude
477 * than 1/sqrt(2) of the maximum value of the return type. Likewise
478 * this determines the necessary threshold for using this method
479 * (however, actually use 1/2 instead as it is simpler).
480 *
481 * Dividing by E causes an insignificant loss of accuracy; however
482 * this method is still poor since it is uneccessarily slow.
483 */
484 if (ax > DBL_MAX / 2)
485 return (CMPLX(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x)));
486
487 /*
488 * Avoid overflow when x or y is large. Avoid underflow when x or
489 * y is small.
490 */
491 if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
492 return (CMPLX(log(hypot(x, y)), atan2(y, x)));
493
494 return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x)));
495 }
496
497 /*
498 * =================
499 * | catanh, catan |
500 * =================
501 */
502
503 /*
504 * sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow).
505 * Assumes x*x and y*y will not overflow.
506 * Assumes x and y are finite.
507 * Assumes y is non-negative.
508 * Assumes fabs(x) >= DBL_EPSILON.
509 */
510 static inline double
sum_squares(double x,double y)511 sum_squares(double x, double y)
512 {
513
514 /* Avoid underflow when y is small. */
515 if (y < SQRT_MIN)
516 return (x * x);
517
518 return (x * x + y * y);
519 }
520
521 /*
522 * real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y).
523 * Assumes x and y are not NaN, and one of x and y is larger than
524 * RECIP_EPSILON. We avoid unwarranted underflow. It is important to not use
525 * the code creal(1/z), because the imaginary part may produce an unwanted
526 * underflow.
527 * This is only called in a context where inexact is always raised before
528 * the call, so no effort is made to avoid or force inexact.
529 */
530 static inline double
real_part_reciprocal(double x,double y)531 real_part_reciprocal(double x, double y)
532 {
533 double scale;
534 uint32_t hx, hy;
535 int32_t ix, iy;
536
537 /*
538 * This code is inspired by the C99 document n1124.pdf, Section G.5.1,
539 * example 2.
540 */
541 GET_HIGH_WORD(hx, x);
542 ix = hx & 0x7ff00000;
543 GET_HIGH_WORD(hy, y);
544 iy = hy & 0x7ff00000;
545 #define BIAS (DBL_MAX_EXP - 1)
546 /* XXX more guard digits are useful iff there is extra precision. */
547 #define CUTOFF (DBL_MANT_DIG / 2 + 1) /* just half or 1 guard digit */
548 if (ix - iy >= CUTOFF << 20 || isinf(x))
549 return (1 / x); /* +-Inf -> +-0 is special */
550 if (iy - ix >= CUTOFF << 20)
551 return (x / y / y); /* should avoid double div, but hard */
552 if (ix <= (BIAS + DBL_MAX_EXP / 2 - CUTOFF) << 20)
553 return (x / (x * x + y * y));
554 scale = 1;
555 SET_HIGH_WORD(scale, 0x7ff00000 - ix); /* 2**(1-ilogb(x)) */
556 x *= scale;
557 y *= scale;
558 return (x / (x * x + y * y) * scale);
559 }
560
561 /*
562 * catanh(z) = log((1+z)/(1-z)) / 2
563 * = log1p(4*x / |z-1|^2) / 4
564 * + I * atan2(2*y, (1-x)*(1+x)-y*y) / 2
565 *
566 * catanh(z) = z + O(z^3) as z -> 0
567 *
568 * catanh(z) = 1/z + sign(y)*I*PI/2 + O(1/z^3) as z -> infinity
569 * The above formula works for the real part as well, because
570 * Re(catanh(z)) = x/|z|^2 + O(x/z^4)
571 * as z -> infinity, uniformly in x
572 */
573 double complex
catanh(double complex z)574 catanh(double complex z)
575 {
576 double x, y, ax, ay, rx, ry;
577
578 x = creal(z);
579 y = cimag(z);
580 ax = fabs(x);
581 ay = fabs(y);
582
583 /* This helps handle many cases. */
584 if (y == 0 && ax <= 1)
585 return (CMPLX(atanh(x), y));
586
587 /* To ensure the same accuracy as atan(), and to filter out z = 0. */
588 if (x == 0)
589 return (CMPLX(x, atan(y)));
590
591 if (isnan(x) || isnan(y)) {
592 /* catanh(+-Inf + I*NaN) = +-0 + I*NaN */
593 if (isinf(x))
594 return (CMPLX(copysign(0, x), y + y));
595 /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
596 if (isinf(y))
597 return (CMPLX(copysign(0, x),
598 copysign(pio2_hi + pio2_lo, y)));
599 /*
600 * All other cases involving NaN return NaN + I*NaN.
601 * C99 leaves it optional whether to raise invalid if one of
602 * the arguments is not NaN, so we opt not to raise it.
603 */
604 return (CMPLX(nan_mix(x, y), nan_mix(x, y)));
605 }
606
607 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
608 return (CMPLX(real_part_reciprocal(x, y),
609 copysign(pio2_hi + pio2_lo, y)));
610
611 if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
612 /*
613 * z = 0 was filtered out above. All other cases must raise
614 * inexact, but this is the only case that needs to do it
615 * explicitly.
616 */
617 raise_inexact();
618 return (z);
619 }
620
621 if (ax == 1 && ay < DBL_EPSILON)
622 rx = (m_ln2 - log(ay)) / 2;
623 else
624 rx = log1p(4 * ax / sum_squares(ax - 1, ay)) / 4;
625
626 if (ax == 1)
627 ry = atan2(2, -ay) / 2;
628 else if (ay < DBL_EPSILON)
629 ry = atan2(2 * ay, (1 - ax) * (1 + ax)) / 2;
630 else
631 ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
632
633 return (CMPLX(copysign(rx, x), copysign(ry, y)));
634 }
635
636 /*
637 * catan(z) = reverse(catanh(reverse(z)))
638 * where reverse(x + I*y) = y + I*x = I*conj(z).
639 */
640 double complex
catan(double complex z)641 catan(double complex z)
642 {
643 double complex w = catanh(CMPLX(cimag(z), creal(z)));
644
645 return (CMPLX(cimag(w), creal(w)));
646 }
647
648 #if LDBL_MANT_DIG == 53
649 __weak_reference(cacosh, cacoshl);
650 __weak_reference(cacos, cacosl);
651 __weak_reference(casinh, casinhl);
652 __weak_reference(casin, casinl);
653 __weak_reference(catanh, catanhl);
654 __weak_reference(catan, catanl);
655 #endif
656