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