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