90  * Throughout we use the convention z = x + I*y.
95  * B = (|z+I| - |z-I|) / 2 = y/A
108  *   A-1 = f(x, 1+y) + f(x, 1-y)
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
143  * x and y are assumed positive or zero, and less than RECIP_EPSILON.
145  * rx = Re(casinh(z)) = -Im(cacos(y + I*x)).
147  * If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y.
152 do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,  in do_hard_work()  argument
156 	double Am1, Amy; /* A-1, A-y. */  in do_hard_work()
158 	R = hypot(x, y + 1);		/* |z+I| */  in do_hard_work()
159 	S = hypot(x, y - 1);		/* |z-I| */  in do_hard_work()
173 		 * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y).  in do_hard_work()
176 		if (y == 1 && x < DBL_EPSILON * DBL_EPSILON / 128) {  in do_hard_work()
182 		} else if (x >= DBL_EPSILON * fabs(y - 1)) {  in do_hard_work()
187 			Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);  in do_hard_work()
189 		} else if (y < 1) {  in do_hard_work()
191 			 * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and  in do_hard_work()
194 			*rx = x / sqrt((1 - y) * (1 + y));  in do_hard_work()
195 		} else {		/* if (y > 1) */  in do_hard_work()
197 			 * A-1 = y-1 (inexactly).  in do_hard_work()
199 			*rx = log1p((y - 1) + sqrt((y - 1) * (y + 1)));  in do_hard_work()
205 	*new_y = y;  in do_hard_work()
207 	if (y < FOUR_SQRT_MIN) {  in do_hard_work()
209 		 * Avoid a possible underflow caused by y/A.  For casinh this  in do_hard_work()
215 		*new_y = y * (2 / DBL_EPSILON);  in do_hard_work()
219 	/* B = (|z+I| - |z-I|) / 2 = y/A */  in do_hard_work()
220 	*B = y / A;  in do_hard_work()
226 		 * Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1).  in do_hard_work()
227 		 * sqrt_A2my2 = sqrt(Amy*(A+y))  in do_hard_work()
229 		if (y == 1 && x < DBL_EPSILON / 128) {  in do_hard_work()
234 			*sqrt_A2my2 = sqrt(x) * sqrt((A + y) / 2);  in do_hard_work()
235 		} else if (x >= DBL_EPSILON * fabs(y - 1)) {  in do_hard_work()
242 			Amy = f(x, y + 1, R) + f(x, y - 1, S);  in do_hard_work()
243 			*sqrt_A2my2 = sqrt(Amy * (A + y));  in do_hard_work()
244 		} else if (y > 1) {  in do_hard_work()
246 			 * fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and  in do_hard_work()
247 			 * A = y (inexactly).  in do_hard_work()
249 			 * y < RECIP_EPSILON.  So the following  in do_hard_work()
252 			*sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y /  in do_hard_work()
253 			    sqrt((y + 1) * (y - 1));  in do_hard_work()
254 			*new_y = y * (4 / DBL_EPSILON / DBL_EPSILON);  in do_hard_work()
255 		} else {		/* if (y < 1) */  in do_hard_work()
257 			 * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and  in do_hard_work()
260 			*sqrt_A2my2 = sqrt((1 - y) * (1 + y));  in do_hard_work()
270  * Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3)
271  *    as z -> infinity, uniformly in y
276 	double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;  in casinh()  local
281 	y = cimag(z);  in casinh()
283 	ay = fabs(y);  in casinh()
285 	if (isnan(x) || isnan(y)) {  in casinh()
288 			return (CMPLX(x, y + y));  in casinh()
290 		if (isinf(y))  in casinh()
291 			return (CMPLX(y, x + x));  in casinh()
293 		if (y == 0)  in casinh()
294 			return (CMPLX(x + x, y));  in casinh()
300 		return (CMPLX(nan_mix(x, y), nan_mix(x, y)));  in casinh()
304 		/* clog...() will raise inexact unless x or y is infinite. */  in casinh()
309 		return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y)));  in casinh()
313 	if (x == 0 && y == 0)  in casinh()
327 	return (CMPLX(copysign(rx, x), copysign(ry, y)));  in casinh()
332  * where reverse(x + I*y) = y + I*x = I*conj(z).
349  * cacos(z) = -sign(y)*I*clog(z) + O(1/z^2)   as z -> infinity
351  * Re(cacos(z)) = atan2(fabs(y), x) + O(y/z^3)
352  *    as z -> infinity, uniformly in y
357 	double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;  in cacos()  local
363 	y = cimag(z);  in cacos()
365 	sy = signbit(y);  in cacos()
367 	ay = fabs(y);  in cacos()
369 	if (isnan(x) || isnan(y)) {  in cacos()
372 			return (CMPLX(y + y, -INFINITY));  in cacos()
374 		if (isinf(y))  in cacos()
375 			return (CMPLX(x + x, -y));  in cacos()
378 			return (CMPLX(pio2_hi + pio2_lo, y + y));  in cacos()
384 		return (CMPLX(nan_mix(x, y), nan_mix(x, y)));  in cacos()
388 		/* clog...() will raise inexact unless x or y is infinite. */  in cacos()
398 	if (x == 1 && y == 0)  in cacos()
399 		return (CMPLX(0, -y));  in cacos()
405 		return (CMPLX(pio2_hi - (x - pio2_lo), -y));  in cacos()
456 	double x, y;  in clog_for_large_values()  local
460 	y = cimag(z);  in clog_for_large_values()
462 	ay = fabs(y);  in clog_for_large_values()
470 	 * Avoid overflow in hypot() when x and y are both very large.  in clog_for_large_values()
471 	 * Divide x and y by E, and then add 1 to the logarithm.  This  in clog_for_large_values()
482 		return (CMPLX(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x)));  in clog_for_large_values()
485 	 * Avoid overflow when x or y is large.  Avoid underflow when x or  in clog_for_large_values()
486 	 * y is small.  in clog_for_large_values()
489 		return (CMPLX(log(hypot(x, y)), atan2(y, x)));  in clog_for_large_values()
491 	return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x)));  in clog_for_large_values()
501  * sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow).
502  * Assumes x*x and y*y will not overflow.
503  * Assumes x and y are finite.
504  * Assumes y is non-negative.
508 sum_squares(double x, double y)  in sum_squares()  argument
511 	/* Avoid underflow when y is small. */  in sum_squares()
512 	if (y < SQRT_MIN)  in sum_squares()
515 	return (x * x + y * y);  in sum_squares()
519  * real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y).
520  * Assumes x and y are not NaN, and one of x and y is larger than
528 real_part_reciprocal(double x, double y)  in real_part_reciprocal()  argument
540 	GET_HIGH_WORD(hy, y);  in real_part_reciprocal()
548 		return (x / y / y);	/* should avoid double div, but hard */  in real_part_reciprocal()
550 		return (x / (x * x + y * y));  in real_part_reciprocal()
554 	y *= scale;  in real_part_reciprocal()
555 	return (x / (x * x + y * y) * scale);  in real_part_reciprocal()
561  *             + I * atan2(2*y, (1-x)*(1+x)-y*y) / 2
565  * catanh(z) = 1/z + sign(y)*I*PI/2 + O(1/z^3)   as z -> infinity
573 	double x, y, ax, ay, rx, ry;  in catanh()  local
576 	y = cimag(z);  in catanh()
578 	ay = fabs(y);  in catanh()
581 	if (y == 0 && ax <= 1)  in catanh()
582 		return (CMPLX(atanh(x), y));  in catanh()
586 		return (CMPLX(x, atan(y)));  in catanh()
588 	if (isnan(x) || isnan(y)) {  in catanh()
591 			return (CMPLX(copysign(0, x), y + y));  in catanh()
593 		if (isinf(y))  in catanh()
595 			    copysign(pio2_hi + pio2_lo, y)));  in catanh()
601 		return (CMPLX(nan_mix(x, y), nan_mix(x, y)));  in catanh()
605 		return (CMPLX(real_part_reciprocal(x, y),  in catanh()
606 		    copysign(pio2_hi + pio2_lo, y)));  in catanh()
630 	return (CMPLX(copysign(rx, x), copysign(ry, y)));  in catanh()
635  * where reverse(x + I*y) = y + I*x = I*conj(z).