Fix our <complex.h> support.

We build libm with -fvisibility=hidden, so we weren't exporting any
of the <complex.h> functions.

We also weren't building many of the functions anyway.

We were also missing the complex inverse trigonometric functions.

And because we didn't even have perfunctory "call each function once"
tests, we didn't notice that we weren't exporting any symbols, so this
patch adds at least that level of testing.

Change-Id: Ibcf2843f507126c51d134cc5fc8d67747e033a0d

2 files changed, 1032 insertions, 0 deletions

diff --git a/libm/upstream-freebsd/lib/msun/src/catrig.c b/libm/upstream-freebsd/lib/msun/src/catrig.c
new file mode 100644
index 000000000..200977c0a
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/catrig.c
@@ -0,0 +1,639 @@
+/*-
+ * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <complex.h>
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+#undef isinf
+#define isinf(x)	(fabs(x) == INFINITY)
+#undef isnan
+#define isnan(x)	((x) != (x))
+#define	raise_inexact()	do { volatile float junk = 1 + tiny; } while(0)
+#undef signbit
+#define signbit(x)	(__builtin_signbit(x))
+
+/* We need that DBL_EPSILON^2/128 is larger than FOUR_SQRT_MIN. */
+static const double
+A_crossover =		10, /* Hull et al suggest 1.5, but 10 works better */
+B_crossover =		0.6417,			/* suggested by Hull et al */
+FOUR_SQRT_MIN =		0x1p-509,		/* >= 4 * sqrt(DBL_MIN) */
+QUARTER_SQRT_MAX =	0x1p509,		/* <= sqrt(DBL_MAX) / 4 */
+m_e =			2.7182818284590452e0,	/*  0x15bf0a8b145769.0p-51 */
+m_ln2 =			6.9314718055994531e-1,	/*  0x162e42fefa39ef.0p-53 */
+pio2_hi =		1.5707963267948966e0,	/*  0x1921fb54442d18.0p-52 */
+RECIP_EPSILON =		1 / DBL_EPSILON,
+SQRT_3_EPSILON =	2.5809568279517849e-8,	/*  0x1bb67ae8584caa.0p-78 */
+SQRT_6_EPSILON =	3.6500241499888571e-8,	/*  0x13988e1409212e.0p-77 */
+SQRT_MIN =		0x1p-511;		/* >= sqrt(DBL_MIN) */
+
+static const volatile double
+pio2_lo =		6.1232339957367659e-17;	/*  0x11a62633145c07.0p-106 */
+static const volatile float
+tiny =			0x1p-100; 
+
+static double complex clog_for_large_values(double complex z);
+
+/*
+ * Testing indicates that all these functions are accurate up to 4 ULP.
+ * The functions casin(h) and cacos(h) are about 2.5 times slower than asinh.
+ * The functions catan(h) are a little under 2 times slower than atanh.
+ *
+ * The code for casinh, casin, cacos, and cacosh comes first.  The code is
+ * rather complicated, and the four functions are highly interdependent.
+ *
+ * The code for catanh and catan comes at the end.  It is much simpler than
+ * the other functions, and the code for these can be disconnected from the
+ * rest of the code.
+ */
+
+/*
+ *			================================
+ *			| casinh, casin, cacos, cacosh |
+ *			================================
+ */
+
+/*
+ * The algorithm is very close to that in "Implementing the complex arcsine
+ * and arccosine functions using exception handling" by T. E. Hull, Thomas F.
+ * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on
+ * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
+ * http://dl.acm.org/citation.cfm?id=275324.
+ *
+ * Throughout we use the convention z = x + I*y.
+ *
+ * casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B)
+ * where
+ * A = (|z+I| + |z-I|) / 2
+ * B = (|z+I| - |z-I|) / 2 = y/A
+ *
+ * These formulas become numerically unstable:
+ *   (a) for Re(casinh(z)) when z is close to the line segment [-I, I] (that
+ *       is, Re(casinh(z)) is close to 0);
+ *   (b) for Im(casinh(z)) when z is close to either of the intervals
+ *       [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is
+ *       close to PI/2).
+ *
+ * These numerical problems are overcome by defining
+ * f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2
+ * Then if A < A_crossover, we use
+ *   log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1)))
+ *   A-1 = f(x, 1+y) + f(x, 1-y)
+ * and if B > B_crossover, we use
+ *   asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y)))
+ *   A-y = f(x, y+1) + f(x, y-1)
+ * where without loss of generality we have assumed that x and y are
+ * non-negative.
+ *
+ * Much of the difficulty comes because the intermediate computations may
+ * produce overflows or underflows.  This is dealt with in the paper by Hull
+ * et al by using exception handling.  We do this by detecting when
+ * computations risk underflow or overflow.  The hardest part is handling the
+ * underflows when computing f(a, b).
+ *
+ * Note that the function f(a, b) does not appear explicitly in the paper by
+ * Hull et al, but the idea may be found on pages 308 and 309.  Introducing the
+ * function f(a, b) allows us to concentrate many of the clever tricks in this
+ * paper into one function.
+ */
+
+/*
+ * Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2.
+ * Pass hypot(a, b) as the third argument.
+ */
+static inline double
+f(double a, double b, double hypot_a_b)
+{
+	if (b < 0)
+		return ((hypot_a_b - b) / 2);
+	if (b == 0)
+		return (a / 2);
+	return (a * a / (hypot_a_b + b) / 2);
+}
+
+/*
+ * All the hard work is contained in this function.
+ * x and y are assumed positive or zero, and less than RECIP_EPSILON.
+ * Upon return:
+ * rx = Re(casinh(z)) = -Im(cacos(y + I*x)).
+ * B_is_usable is set to 1 if the value of B is usable.
+ * If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y.
+ * If returning sqrt_A2my2 has potential to result in an underflow, it is
+ * rescaled, and new_y is similarly rescaled.
+ */
+static inline void
+do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
+    double *sqrt_A2my2, double *new_y)
+{
+	double R, S, A; /* A, B, R, and S are as in Hull et al. */
+	double Am1, Amy; /* A-1, A-y. */
+
+	R = hypot(x, y + 1);		/* |z+I| */
+	S = hypot(x, y - 1);		/* |z-I| */
+
+	/* A = (|z+I| + |z-I|) / 2 */
+	A = (R + S) / 2;
+	/*
+	 * Mathematically A >= 1.  There is a small chance that this will not
+	 * be so because of rounding errors.  So we will make certain it is
+	 * so.
+	 */
+	if (A < 1)
+		A = 1;
+
+	if (A < A_crossover) {
+		/*
+		 * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y).
+		 * rx = log1p(Am1 + sqrt(Am1*(A+1)))
+		 */
+		if (y == 1 && x < DBL_EPSILON * DBL_EPSILON / 128) {
+			/*
+			 * fp is of order x^2, and fm = x/2.
+			 * A = 1 (inexactly).
+			 */
+			*rx = sqrt(x);
+		} else if (x >= DBL_EPSILON * fabs(y - 1)) {
+			/*
+			 * Underflow will not occur because
+			 * x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN
+			 */
+			Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
+			*rx = log1p(Am1 + sqrt(Am1 * (A + 1)));
+		} else if (y < 1) {
+			/*
+			 * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and
+			 * A = 1 (inexactly).
+			 */
+			*rx = x / sqrt((1 - y) * (1 + y));
+		} else {		/* if (y > 1) */
+			/*
+			 * A-1 = y-1 (inexactly).
+			 */
+			*rx = log1p((y - 1) + sqrt((y - 1) * (y + 1)));
+		}
+	} else {
+		*rx = log(A + sqrt(A * A - 1));
+	}
+
+	*new_y = y;
+
+	if (y < FOUR_SQRT_MIN) {
+		/*
+		 * Avoid a possible underflow caused by y/A.  For casinh this
+		 * would be legitimate, but will be picked up by invoking atan2
+		 * later on.  For cacos this would not be legitimate.
+		 */
+		*B_is_usable = 0;
+		*sqrt_A2my2 = A * (2 / DBL_EPSILON);
+		*new_y = y * (2 / DBL_EPSILON);
+		return;
+	}
+
+	/* B = (|z+I| - |z-I|) / 2 = y/A */
+	*B = y / A;
+	*B_is_usable = 1;
+
+	if (*B > B_crossover) {
+		*B_is_usable = 0;
+		/*
+		 * Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1).
+		 * sqrt_A2my2 = sqrt(Amy*(A+y))
+		 */
+		if (y == 1 && x < DBL_EPSILON / 128) {
+			/*
+			 * fp is of order x^2, and fm = x/2.
+			 * A = 1 (inexactly).
+			 */
+			*sqrt_A2my2 = sqrt(x) * sqrt((A + y) / 2);
+		} else if (x >= DBL_EPSILON * fabs(y - 1)) {
+			/*
+			 * Underflow will not occur because
+			 * x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN
+			 * and
+			 * x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN
+			 */
+			Amy = f(x, y + 1, R) + f(x, y - 1, S);
+			*sqrt_A2my2 = sqrt(Amy * (A + y));
+		} else if (y > 1) {
+			/*
+			 * fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and
+			 * A = y (inexactly).
+			 *
+			 * y < RECIP_EPSILON.  So the following
+			 * scaling should avoid any underflow problems.
+			 */
+			*sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y /
+			    sqrt((y + 1) * (y - 1));
+			*new_y = y * (4 / DBL_EPSILON / DBL_EPSILON);
+		} else {		/* if (y < 1) */
+			/*
+			 * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and
+			 * A = 1 (inexactly).
+			 */
+			*sqrt_A2my2 = sqrt((1 - y) * (1 + y));
+		}
+	}
+}
+
+/*
+ * casinh(z) = z + O(z^3)   as z -> 0
+ *
+ * casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2)   as z -> infinity
+ * The above formula works for the imaginary part as well, because
+ * Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3)
+ *    as z -> infinity, uniformly in y
+ */
+double complex
+casinh(double complex z)
+{
+	double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;
+	int B_is_usable;
+	double complex w;
+
+	x = creal(z);
+	y = cimag(z);
+	ax = fabs(x);
+	ay = fabs(y);
+
+	if (isnan(x) || isnan(y)) {
+		/* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */
+		if (isinf(x))
+			return (cpack(x, y + y));
+		/* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
+		if (isinf(y))
+			return (cpack(y, x + x));
+		/* casinh(NaN + I*0) = NaN + I*0 */
+		if (y == 0)
+			return (cpack(x + x, y));
+		/*
+		 * All other cases involving NaN return NaN + I*NaN.
+		 * C99 leaves it optional whether to raise invalid if one of
+		 * the arguments is not NaN, so we opt not to raise it.
+		 */
+		return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+		/* clog...() will raise inexact unless x or y is infinite. */
+		if (signbit(x) == 0)
+			w = clog_for_large_values(z) + m_ln2;
+		else
+			w = clog_for_large_values(-z) + m_ln2;
+		return (cpack(copysign(creal(w), x), copysign(cimag(w), y)));
+	}
+
+	/* Avoid spuriously raising inexact for z = 0. */
+	if (x == 0 && y == 0)
+		return (z);
+
+	/* All remaining cases are inexact. */
+	raise_inexact();
+
+	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+		return (z);
+
+	do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
+	if (B_is_usable)
+		ry = asin(B);
+	else
+		ry = atan2(new_y, sqrt_A2my2);
+	return (cpack(copysign(rx, x), copysign(ry, y)));
+}
+
+/*
+ * casin(z) = reverse(casinh(reverse(z)))
+ * where reverse(x + I*y) = y + I*x = I*conj(z).
+ */
+double complex
+casin(double complex z)
+{
+	double complex w = casinh(cpack(cimag(z), creal(z)));
+
+	return (cpack(cimag(w), creal(w)));
+}
+
+/*
+ * cacos(z) = PI/2 - casin(z)
+ * but do the computation carefully so cacos(z) is accurate when z is
+ * close to 1.
+ *
+ * cacos(z) = PI/2 - z + O(z^3)   as z -> 0
+ *
+ * cacos(z) = -sign(y)*I*clog(z) + O(1/z^2)   as z -> infinity
+ * The above formula works for the real part as well, because
+ * Re(cacos(z)) = atan2(fabs(y), x) + O(y/z^3)
+ *    as z -> infinity, uniformly in y
+ */
+double complex
+cacos(double complex z)
+{
+	double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
+	int sx, sy;
+	int B_is_usable;
+	double complex w;
+
+	x = creal(z);
+	y = cimag(z);
+	sx = signbit(x);
+	sy = signbit(y);
+	ax = fabs(x);
+	ay = fabs(y);
+
+	if (isnan(x) || isnan(y)) {
+		/* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
+		if (isinf(x))
+			return (cpack(y + y, -INFINITY));
+		/* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */
+		if (isinf(y))
+			return (cpack(x + x, -y));
+		/* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */
+		if (x == 0)
+			return (cpack(pio2_hi + pio2_lo, y + y));
+		/*
+		 * All other cases involving NaN return NaN + I*NaN.
+		 * C99 leaves it optional whether to raise invalid if one of
+		 * the arguments is not NaN, so we opt not to raise it.
+		 */
+		return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+		/* clog...() will raise inexact unless x or y is infinite. */
+		w = clog_for_large_values(z);
+		rx = fabs(cimag(w));
+		ry = creal(w) + m_ln2;
+		if (sy == 0)
+			ry = -ry;
+		return (cpack(rx, ry));
+	}
+
+	/* Avoid spuriously raising inexact for z = 1. */
+	if (x == 1 && y == 0)
+		return (cpack(0, -y));
+
+	/* All remaining cases are inexact. */
+	raise_inexact();
+
+	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+		return (cpack(pio2_hi - (x - pio2_lo), -y));
+
+	do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
+	if (B_is_usable) {
+		if (sx == 0)
+			rx = acos(B);
+		else
+			rx = acos(-B);
+	} else {
+		if (sx == 0)
+			rx = atan2(sqrt_A2mx2, new_x);
+		else
+			rx = atan2(sqrt_A2mx2, -new_x);
+	}
+	if (sy == 0)
+		ry = -ry;
+	return (cpack(rx, ry));
+}
+
+/*
+ * cacosh(z) = I*cacos(z) or -I*cacos(z)
+ * where the sign is chosen so Re(cacosh(z)) >= 0.
+ */
+double complex
+cacosh(double complex z)
+{
+	double complex w;
+	double rx, ry;
+
+	w = cacos(z);
+	rx = creal(w);
+	ry = cimag(w);
+	/* cacosh(NaN + I*NaN) = NaN + I*NaN */
+	if (isnan(rx) && isnan(ry))
+		return (cpack(ry, rx));
+	/* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */
+	/* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */
+	if (isnan(rx))
+		return (cpack(fabs(ry), rx));
+	/* cacosh(0 + I*NaN) = NaN + I*NaN */
+	if (isnan(ry))
+		return (cpack(ry, ry));
+	return (cpack(fabs(ry), copysign(rx, cimag(z))));
+}
+
+/*
+ * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON.
+ */
+static double complex
+clog_for_large_values(double complex z)
+{
+	double x, y;
+	double ax, ay, t;
+
+	x = creal(z);
+	y = cimag(z);
+	ax = fabs(x);
+	ay = fabs(y);
+	if (ax < ay) {
+		t = ax;
+		ax = ay;
+		ay = t;
+	}
+
+	/*
+	 * Avoid overflow in hypot() when x and y are both very large.
+	 * Divide x and y by E, and then add 1 to the logarithm.  This depends
+	 * on E being larger than sqrt(2).
+	 * Dividing by E causes an insignificant loss of accuracy; however
+	 * this method is still poor since it is uneccessarily slow.
+	 */
+	if (ax > DBL_MAX / 2)
+		return (cpack(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x)));
+
+	/*
+	 * Avoid overflow when x or y is large.  Avoid underflow when x or
+	 * y is small.
+	 */
+	if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
+		return (cpack(log(hypot(x, y)), atan2(y, x)));
+
+	return (cpack(log(ax * ax + ay * ay) / 2, atan2(y, x)));
+}
+
+/*
+ *				=================
+ *				| catanh, catan |
+ *				=================
+ */
+
+/*
+ * sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow).
+ * Assumes x*x and y*y will not overflow.
+ * Assumes x and y are finite.
+ * Assumes y is non-negative.
+ * Assumes fabs(x) >= DBL_EPSILON.
+ */
+static inline double
+sum_squares(double x, double y)
+{
+
+	/* Avoid underflow when y is small. */
+	if (y < SQRT_MIN)
+		return (x * x);
+
+	return (x * x + y * y);
+}
+
+/*
+ * real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y).
+ * Assumes x and y are not NaN, and one of x and y is larger than
+ * RECIP_EPSILON.  We avoid unwarranted underflow.  It is important to not use
+ * the code creal(1/z), because the imaginary part may produce an unwanted
+ * underflow.
+ * This is only called in a context where inexact is always raised before
+ * the call, so no effort is made to avoid or force inexact.
+ */
+static inline double
+real_part_reciprocal(double x, double y)
+{
+	double scale;
+	uint32_t hx, hy;
+	int32_t ix, iy;
+
+	/*
+	 * This code is inspired by the C99 document n1124.pdf, Section G.5.1,
+	 * example 2.
+	 */
+	GET_HIGH_WORD(hx, x);
+	ix = hx & 0x7ff00000;
+	GET_HIGH_WORD(hy, y);
+	iy = hy & 0x7ff00000;
+#define	BIAS	(DBL_MAX_EXP - 1)
+/* XXX more guard digits are useful iff there is extra precision. */
+#define	CUTOFF	(DBL_MANT_DIG / 2 + 1)	/* just half or 1 guard digit */
+	if (ix - iy >= CUTOFF << 20 || isinf(x))
+		return (1 / x);		/* +-Inf -> +-0 is special */
+	if (iy - ix >= CUTOFF << 20)
+		return (x / y / y);	/* should avoid double div, but hard */
+	if (ix <= (BIAS + DBL_MAX_EXP / 2 - CUTOFF) << 20)
+		return (x / (x * x + y * y));
+	scale = 1;
+	SET_HIGH_WORD(scale, 0x7ff00000 - ix);	/* 2**(1-ilogb(x)) */
+	x *= scale;
+	y *= scale;
+	return (x / (x * x + y * y) * scale);
+}
+
+/*
+ * catanh(z) = log((1+z)/(1-z)) / 2
+ *           = log1p(4*x / |z-1|^2) / 4
+ *             + I * atan2(2*y, (1-x)*(1+x)-y*y) / 2
+ *
+ * catanh(z) = z + O(z^3)   as z -> 0
+ *
+ * catanh(z) = 1/z + sign(y)*I*PI/2 + O(1/z^3)   as z -> infinity
+ * The above formula works for the real part as well, because
+ * Re(catanh(z)) = x/|z|^2 + O(x/z^4)
+ *    as z -> infinity, uniformly in x
+ */
+double complex
+catanh(double complex z)
+{
+	double x, y, ax, ay, rx, ry;
+
+	x = creal(z);
+	y = cimag(z);
+	ax = fabs(x);
+	ay = fabs(y);
+
+	/* This helps handle many cases. */
+	if (y == 0 && ax <= 1)
+		return (cpack(atanh(x), y));
+
+	/* To ensure the same accuracy as atan(), and to filter out z = 0. */
+	if (x == 0)
+		return (cpack(x, atan(y)));
+
+	if (isnan(x) || isnan(y)) {
+		/* catanh(+-Inf + I*NaN) = +-0 + I*NaN */
+		if (isinf(x))
+			return (cpack(copysign(0, x), y + y));
+		/* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
+		if (isinf(y))
+			return (cpack(copysign(0, x),
+			    copysign(pio2_hi + pio2_lo, y)));
+		/*
+		 * All other cases involving NaN return NaN + I*NaN.
+		 * C99 leaves it optional whether to raise invalid if one of
+		 * the arguments is not NaN, so we opt not to raise it.
+		 */
+		return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
+		return (cpack(real_part_reciprocal(x, y),
+		    copysign(pio2_hi + pio2_lo, y)));
+
+	if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
+		/*
+		 * z = 0 was filtered out above.  All other cases must raise
+		 * inexact, but this is the only only that needs to do it
+		 * explicitly.
+		 */
+		raise_inexact();
+		return (z);
+	}
+
+	if (ax == 1 && ay < DBL_EPSILON)
+		rx = (m_ln2 - log(ay)) / 2;
+	else
+		rx = log1p(4 * ax / sum_squares(ax - 1, ay)) / 4;
+
+	if (ax == 1)
+		ry = atan2(2, -ay) / 2;
+	else if (ay < DBL_EPSILON)
+		ry = atan2(2 * ay, (1 - ax) * (1 + ax)) / 2;
+	else
+		ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
+
+	return (cpack(copysign(rx, x), copysign(ry, y)));
+}
+
+/*
+ * catan(z) = reverse(catanh(reverse(z)))
+ * where reverse(x + I*y) = y + I*x = I*conj(z).
+ */
+double complex
+catan(double complex z)
+{
+	double complex w = catanh(cpack(cimag(z), creal(z)));
+
+	return (cpack(cimag(w), creal(w)));
+}
diff --git a/libm/upstream-freebsd/lib/msun/src/catrigf.c b/libm/upstream-freebsd/lib/msun/src/catrigf.c
new file mode 100644
index 000000000..08ebef783
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/catrigf.c
@@ -0,0 +1,393 @@
+/*-
+ * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+/*
+ * The algorithm is very close to that in "Implementing the complex arcsine
+ * and arccosine functions using exception handling" by T. E. Hull, Thomas F.
+ * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on
+ * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
+ * http://dl.acm.org/citation.cfm?id=275324.
+ *
+ * See catrig.c for complete comments.
+ *
+ * XXX comments were removed automatically, and even short ones on the right
+ * of statements were removed (all of them), contrary to normal style.  Only
+ * a few comments on the right of declarations remain.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <complex.h>
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+#undef isinf
+#define isinf(x)	(fabsf(x) == INFINITY)
+#undef isnan
+#define isnan(x)	((x) != (x))
+#define	raise_inexact()	do { volatile float junk = 1 + tiny; } while(0)
+#undef signbit
+#define signbit(x)	(__builtin_signbitf(x))
+
+static const float
+A_crossover =		10,
+B_crossover =		0.6417,
+FOUR_SQRT_MIN =		0x1p-61,
+QUARTER_SQRT_MAX =	0x1p61,
+m_e =			2.7182818285e0,		/*  0xadf854.0p-22 */
+m_ln2 =			6.9314718056e-1,	/*  0xb17218.0p-24 */
+pio2_hi =		1.5707962513e0,		/*  0xc90fda.0p-23 */
+RECIP_EPSILON =		1 / FLT_EPSILON,
+SQRT_3_EPSILON =	5.9801995673e-4,	/*  0x9cc471.0p-34 */
+SQRT_6_EPSILON =	8.4572793338e-4,	/*  0xddb3d7.0p-34 */
+SQRT_MIN =		0x1p-63;
+
+static const volatile float
+pio2_lo =		7.5497899549e-8,	/*  0xa22169.0p-47 */
+tiny =			0x1p-100;
+
+static float complex clog_for_large_values(float complex z);
+
+static inline float
+f(float a, float b, float hypot_a_b)
+{
+	if (b < 0)
+		return ((hypot_a_b - b) / 2);
+	if (b == 0)
+		return (a / 2);
+	return (a * a / (hypot_a_b + b) / 2);
+}
+
+static inline void
+do_hard_work(float x, float y, float *rx, int *B_is_usable, float *B,
+    float *sqrt_A2my2, float *new_y)
+{
+	float R, S, A;
+	float Am1, Amy;
+
+	R = hypotf(x, y + 1);
+	S = hypotf(x, y - 1);
+
+	A = (R + S) / 2;
+	if (A < 1)
+		A = 1;
+
+	if (A < A_crossover) {
+		if (y == 1 && x < FLT_EPSILON * FLT_EPSILON / 128) {
+			*rx = sqrtf(x);
+		} else if (x >= FLT_EPSILON * fabsf(y - 1)) {
+			Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
+			*rx = log1pf(Am1 + sqrtf(Am1 * (A + 1)));
+		} else if (y < 1) {
+			*rx = x / sqrtf((1 - y) * (1 + y));
+		} else {
+			*rx = log1pf((y - 1) + sqrtf((y - 1) * (y + 1)));
+		}
+	} else {
+		*rx = logf(A + sqrtf(A * A - 1));
+	}
+
+	*new_y = y;
+
+	if (y < FOUR_SQRT_MIN) {
+		*B_is_usable = 0;
+		*sqrt_A2my2 = A * (2 / FLT_EPSILON);
+		*new_y = y * (2 / FLT_EPSILON);
+		return;
+	}
+
+	*B = y / A;
+	*B_is_usable = 1;
+
+	if (*B > B_crossover) {
+		*B_is_usable = 0;
+		if (y == 1 && x < FLT_EPSILON / 128) {
+			*sqrt_A2my2 = sqrtf(x) * sqrtf((A + y) / 2);
+		} else if (x >= FLT_EPSILON * fabsf(y - 1)) {
+			Amy = f(x, y + 1, R) + f(x, y - 1, S);
+			*sqrt_A2my2 = sqrtf(Amy * (A + y));
+		} else if (y > 1) {
+			*sqrt_A2my2 = x * (4 / FLT_EPSILON / FLT_EPSILON) * y /
+			    sqrtf((y + 1) * (y - 1));
+			*new_y = y * (4 / FLT_EPSILON / FLT_EPSILON);
+		} else {
+			*sqrt_A2my2 = sqrtf((1 - y) * (1 + y));
+		}
+	}
+}
+
+float complex
+casinhf(float complex z)
+{
+	float x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;
+	int B_is_usable;
+	float complex w;
+
+	x = crealf(z);
+	y = cimagf(z);
+	ax = fabsf(x);
+	ay = fabsf(y);
+
+	if (isnan(x) || isnan(y)) {
+		if (isinf(x))
+			return (cpackf(x, y + y));
+		if (isinf(y))
+			return (cpackf(y, x + x));
+		if (y == 0)
+			return (cpackf(x + x, y));
+		return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+		if (signbit(x) == 0)
+			w = clog_for_large_values(z) + m_ln2;
+		else
+			w = clog_for_large_values(-z) + m_ln2;
+		return (cpackf(copysignf(crealf(w), x),
+		    copysignf(cimagf(w), y)));
+	}
+
+	if (x == 0 && y == 0)
+		return (z);
+
+	raise_inexact();
+
+	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+		return (z);
+
+	do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
+	if (B_is_usable)
+		ry = asinf(B);
+	else
+		ry = atan2f(new_y, sqrt_A2my2);
+	return (cpackf(copysignf(rx, x), copysignf(ry, y)));
+}
+
+float complex
+casinf(float complex z)
+{
+	float complex w = casinhf(cpackf(cimagf(z), crealf(z)));
+
+	return (cpackf(cimagf(w), crealf(w)));
+}
+
+float complex
+cacosf(float complex z)
+{
+	float x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
+	int sx, sy;
+	int B_is_usable;
+	float complex w;
+
+	x = crealf(z);
+	y = cimagf(z);
+	sx = signbit(x);
+	sy = signbit(y);
+	ax = fabsf(x);
+	ay = fabsf(y);
+
+	if (isnan(x) || isnan(y)) {
+		if (isinf(x))
+			return (cpackf(y + y, -INFINITY));
+		if (isinf(y))
+			return (cpackf(x + x, -y));
+		if (x == 0)
+			return (cpackf(pio2_hi + pio2_lo, y + y));
+		return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+		w = clog_for_large_values(z);
+		rx = fabsf(cimagf(w));
+		ry = crealf(w) + m_ln2;
+		if (sy == 0)
+			ry = -ry;
+		return (cpackf(rx, ry));
+	}
+
+	if (x == 1 && y == 0)
+		return (cpackf(0, -y));
+
+	raise_inexact();
+
+	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+		return (cpackf(pio2_hi - (x - pio2_lo), -y));
+
+	do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
+	if (B_is_usable) {
+		if (sx == 0)
+			rx = acosf(B);
+		else
+			rx = acosf(-B);
+	} else {
+		if (sx == 0)
+			rx = atan2f(sqrt_A2mx2, new_x);
+		else
+			rx = atan2f(sqrt_A2mx2, -new_x);
+	}
+	if (sy == 0)
+		ry = -ry;
+	return (cpackf(rx, ry));
+}
+
+float complex
+cacoshf(float complex z)
+{
+	float complex w;
+	float rx, ry;
+
+	w = cacosf(z);
+	rx = crealf(w);
+	ry = cimagf(w);
+	if (isnan(rx) && isnan(ry))
+		return (cpackf(ry, rx));
+	if (isnan(rx))
+		return (cpackf(fabsf(ry), rx));
+	if (isnan(ry))
+		return (cpackf(ry, ry));
+	return (cpackf(fabsf(ry), copysignf(rx, cimagf(z))));
+}
+
+static float complex
+clog_for_large_values(float complex z)
+{
+	float x, y;
+	float ax, ay, t;
+
+	x = crealf(z);
+	y = cimagf(z);
+	ax = fabsf(x);
+	ay = fabsf(y);
+	if (ax < ay) {
+		t = ax;
+		ax = ay;
+		ay = t;
+	}
+
+	if (ax > FLT_MAX / 2)
+		return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1,
+		    atan2f(y, x)));
+
+	if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
+		return (cpackf(logf(hypotf(x, y)), atan2f(y, x)));
+
+	return (cpackf(logf(ax * ax + ay * ay) / 2, atan2f(y, x)));
+}
+
+static inline float
+sum_squares(float x, float y)
+{
+
+	if (y < SQRT_MIN)
+		return (x * x);
+
+	return (x * x + y * y);
+}
+
+static inline float
+real_part_reciprocal(float x, float y)
+{
+	float scale;
+	uint32_t hx, hy;
+	int32_t ix, iy;
+
+	GET_FLOAT_WORD(hx, x);
+	ix = hx & 0x7f800000;
+	GET_FLOAT_WORD(hy, y);
+	iy = hy & 0x7f800000;
+#define	BIAS	(FLT_MAX_EXP - 1)
+#define	CUTOFF	(FLT_MANT_DIG / 2 + 1)
+	if (ix - iy >= CUTOFF << 23 || isinf(x))
+		return (1 / x);
+	if (iy - ix >= CUTOFF << 23)
+		return (x / y / y);
+	if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23)
+		return (x / (x * x + y * y));
+	SET_FLOAT_WORD(scale, 0x7f800000 - ix);
+	x *= scale;
+	y *= scale;
+	return (x / (x * x + y * y) * scale);
+}
+
+float complex
+catanhf(float complex z)
+{
+	float x, y, ax, ay, rx, ry;
+
+	x = crealf(z);
+	y = cimagf(z);
+	ax = fabsf(x);
+	ay = fabsf(y);
+
+	if (y == 0 && ax <= 1)
+		return (cpackf(atanhf(x), y));
+
+	if (x == 0)
+		return (cpackf(x, atanf(y)));
+
+	if (isnan(x) || isnan(y)) {
+		if (isinf(x))
+			return (cpackf(copysignf(0, x), y + y));
+		if (isinf(y))
+			return (cpackf(copysignf(0, x),
+			    copysignf(pio2_hi + pio2_lo, y)));
+		return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
+		return (cpackf(real_part_reciprocal(x, y),
+		    copysignf(pio2_hi + pio2_lo, y)));
+
+	if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
+		raise_inexact();
+		return (z);
+	}
+
+	if (ax == 1 && ay < FLT_EPSILON)
+		rx = (m_ln2 - logf(ay)) / 2;
+	else
+		rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4;
+
+	if (ax == 1)
+		ry = atan2f(2, -ay) / 2;
+	else if (ay < FLT_EPSILON)
+		ry = atan2f(2 * ay, (1 - ax) * (1 + ax)) / 2;
+	else
+		ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
+
+	return (cpackf(copysignf(rx, x), copysignf(ry, y)));
+}
+
+float complex
+catanf(float complex z)
+{
+	float complex w = catanhf(cpackf(cimagf(z), crealf(z)));
+
+	return (cpackf(cimagf(w), crealf(w)));
+}