1 /*-
2  * SPDX-License-Identifier: BSD-2-Clause
4  * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
18  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
47 FOUR_SQRT_MIN =		0x1p-509,		/* >= 4 * sqrt(DBL_MIN) */
49 m_e =			2.7182818284590452e0,	/*  0x15bf0a8b145769.0p-51 */
50 m_ln2 =			6.9314718055994531e-1,	/*  0x162e42fefa39ef.0p-53 */
51 pio2_hi =		1.5707963267948966e0,	/*  0x1921fb54442d18.0p-52 */
53 SQRT_3_EPSILON =	2.5809568279517849e-8,	/*  0x1bb67ae8584caa.0p-78 */
54 SQRT_6_EPSILON =	3.6500241499888571e-8,	/*  0x13988e1409212e.0p-77 */
55 SQRT_MIN =		0x1p-511;		/* >= sqrt(DBL_MIN) */
58 pio2_lo =		6.1232339957367659e-17;	/*  0x11a62633145c07.0p-106 */
60 tiny =			0x1p-100; 
62 static double complex clog_for_large_values(double complex z);
67  * The functions catan(h) are a little under 2 times slower than atanh.
87  * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
90  * Throughout we use the convention z = x + I*y.
92  * casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B)
94  * A = (|z+I| + |z-I|) / 2
95  * B = (|z+I| - |z-I|) / 2 = y/A
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
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)
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)
113  * non-negative.
119  * underflows when computing f(a, b).
121  * Note that the function f(a, b) does not appear explicitly in the paper by
123  * function f(a, b) allows us to concentrate many of the clever tricks in this
128  * Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2.
129  * Pass hypot(a, b) as the third argument.
132 f(double a, double b, double hypot_a_b)  in f()  argument
135 		return ((hypot_a_b - b) / 2);  in f()
137 		return (a / 2);  in f()
138 	return (a * a / (hypot_a_b + b) / 2);  in f()
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.
155 	double R, S, A; /* A, B, R, and S are as in Hull et al. */  in do_hard_work()  local
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()
161 	/* A = (|z+I| + |z-I|) / 2 */  in do_hard_work()
162 	A = (R + S) / 2;  in do_hard_work()
164 	 * Mathematically A >= 1.  There is a small chance that this will not  in do_hard_work()
168 	if (A < 1)  in do_hard_work()
169 		A = 1;  in do_hard_work()
171 	if (A < A_crossover) {  in do_hard_work()
173 		 * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y).  in do_hard_work()
174 		 * rx = log1p(Am1 + sqrt(Am1*(A+1)))  in do_hard_work()
179 			 * A = 1 (inexactly).  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()
188 			*rx = log1p(Am1 + sqrt(Am1 * (A + 1)));  in do_hard_work()
191 			 * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and  in do_hard_work()
192 			 * A = 1 (inexactly).  in do_hard_work()
194 			*rx = x / sqrt((1 - y) * (1 + y));  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()
202 		*rx = log(A + sqrt(A * A - 1));  in do_hard_work()
209 		 * Avoid a possible underflow caused by y/A.  For casinh this  in do_hard_work()
214 		*sqrt_A2my2 = A * (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()
232 			 * A = 1 (inexactly).  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()
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()
253 			    sqrt((y + 1) * (y - 1));  in do_hard_work()
257 			 * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and  in do_hard_work()
258 			 * A = 1 (inexactly).  in do_hard_work()
260 			*sqrt_A2my2 = sqrt((1 - y) * (1 + y));  in do_hard_work()
266  * casinh(z) = z + O(z^3)   as z -> 0
268  * casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2)   as z -> infinity
270  * Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3)
271  *    as z -> infinity, uniformly in y
274 casinh(double complex z)  in casinh()  argument
280 	x = creal(z);  in casinh()
281 	y = cimag(z);  in casinh()
286 		/* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */  in casinh()
289 		/* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */  in casinh()
306 			w = clog_for_large_values(z) + m_ln2;  in casinh()
308 			w = clog_for_large_values(-z) + m_ln2;  in casinh()
312 	/* Avoid spuriously raising inexact for z = 0. */  in casinh()
314 		return (z);  in casinh()
320 		return (z);  in casinh()
331  * casin(z) = reverse(casinh(reverse(z)))
332  * where reverse(x + I*y) = y + I*x = I*conj(z).
335 casin(double complex z)  in casin()  argument
337 	double complex w = casinh(CMPLX(cimag(z), creal(z)));  in casin()
343  * cacos(z) = PI/2 - casin(z)
344  * but do the computation carefully so cacos(z) is accurate when z is
347  * cacos(z) = PI/2 - z + O(z^3)   as z -> 0
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
355 cacos(double complex z)  in cacos()  argument
362 	x = creal(z);  in cacos()
363 	y = cimag(z);  in cacos()
370 		/* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */  in cacos()
372 			return (CMPLX(y + y, -INFINITY));  in cacos()
373 		/* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */  in cacos()
375 			return (CMPLX(x + x, -y));  in cacos()
389 		w = clog_for_large_values(z);  in cacos()
393 			ry = -ry;  in cacos()
397 	/* Avoid spuriously raising inexact for z = 1. */  in cacos()
399 		return (CMPLX(0, -y));  in cacos()
405 		return (CMPLX(pio2_hi - (x - pio2_lo), -y));  in cacos()
412 			rx = acos(-B);  in cacos()
417 			rx = atan2(sqrt_A2mx2, -new_x);  in cacos()
420 		ry = -ry;  in cacos()
425  * cacosh(z) = I*cacos(z) or -I*cacos(z)
426  * where the sign is chosen so Re(cacosh(z)) >= 0.
429 cacosh(double complex z)  in cacosh()  argument
434 	w = cacos(z);  in cacosh()
440 	/* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */  in cacosh()
441 	/* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */  in cacosh()
447 	return (CMPLX(fabs(ry), copysign(rx, cimag(z))));  in cacosh()
451  * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON.
454 clog_for_large_values(double complex z)  in clog_for_large_values()  argument
459 	x = creal(z);  in clog_for_large_values()
460 	y = cimag(z);  in clog_for_large_values()
504  * Assumes y is non-negative.
522  * the code creal(1/z), because the imaginary part may produce an unwanted
524  * This is only called in a context where inexact is always raised before
542 #define	BIAS	(DBL_MAX_EXP - 1)  in real_part_reciprocal()
545 	if (ix - iy >= CUTOFF << 20 || isinf(x))  in real_part_reciprocal()
546 		return (1 / x);		/* +-Inf -> +-0 is special */  in real_part_reciprocal()
547 	if (iy - ix >= CUTOFF << 20)  in real_part_reciprocal()
549 	if (ix <= (BIAS + DBL_MAX_EXP / 2 - CUTOFF) << 20)  in real_part_reciprocal()
552 	SET_HIGH_WORD(scale, 0x7ff00000 - ix);	/* 2**(1-ilogb(x)) */  in real_part_reciprocal()
559  * catanh(z) = log((1+z)/(1-z)) / 2
560  *           = log1p(4*x / |z-1|^2) / 4
561  *             + I * atan2(2*y, (1-x)*(1+x)-y*y) / 2
563  * catanh(z) = z + O(z^3)   as z -> 0
565  * catanh(z) = 1/z + sign(y)*I*PI/2 + O(1/z^3)   as z -> infinity
567  * Re(catanh(z)) = x/|z|^2 + O(x/z^4)
568  *    as z -> infinity, uniformly in x
571 catanh(double complex z)  in catanh()  argument
575 	x = creal(z);  in catanh()
576 	y = cimag(z);  in catanh()
584 	/* To ensure the same accuracy as atan(), and to filter out z = 0. */  in catanh()
589 		/* catanh(+-Inf + I*NaN) = +-0 + I*NaN */  in catanh()
592 		/* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */  in catanh()
610 		 * z = 0 was filtered out above.  All other cases must raise  in catanh()
615 		return (z);  in catanh()
619 		rx = (m_ln2 - log(ay)) / 2;  in catanh()
621 		rx = log1p(4 * ax / sum_squares(ax - 1, ay)) / 4;  in catanh()
624 		ry = atan2(2, -ay) / 2;  in catanh()
626 		ry = atan2(2 * ay, (1 - ax) * (1 + ax)) / 2;  in catanh()
628 		ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;  in catanh()
634  * catan(z) = reverse(catanh(reverse(z)))
635  * where reverse(x + I*y) = y + I*x = I*conj(z).
638 catan(double complex z)  in catan()  argument
640 	double complex w = catanh(CMPLX(cimag(z), creal(z)));  in catan()