Module Name:    src
Committed By:   christos
Date:           Mon Sep 19 22:05:05 UTC 2016

Modified Files:
        src/lib/libm/complex: Makefile.inc
        src/lib/libm/src: math_private.h
Added Files:
        src/lib/libm/complex: catrig.c catrigf.c catrigl.c

Log Message:
Add the complex trig functions from FreeBSD


To generate a diff of this commit:
cvs rdiff -u -r1.8 -r1.9 src/lib/libm/complex/Makefile.inc
cvs rdiff -u -r0 -r1.1 src/lib/libm/complex/catrig.c \
    src/lib/libm/complex/catrigf.c src/lib/libm/complex/catrigl.c
cvs rdiff -u -r1.22 -r1.23 src/lib/libm/src/math_private.h

Please note that diffs are not public domain; they are subject to the
copyright notices on the relevant files.

Modified files:

Index: src/lib/libm/complex/Makefile.inc
diff -u src/lib/libm/complex/Makefile.inc:1.8 src/lib/libm/complex/Makefile.inc:1.9
--- src/lib/libm/complex/Makefile.inc:1.8	Fri Oct 10 08:43:07 2014
+++ src/lib/libm/complex/Makefile.inc	Mon Sep 19 18:05:05 2016
@@ -1,14 +1,25 @@
-# $NetBSD: Makefile.inc,v 1.8 2014/10/10 12:43:07 christos Exp $
+# $NetBSD: Makefile.inc,v 1.9 2016/09/19 22:05:05 christos Exp $
 
 .PATH: ${.CURDIR}/complex
 
-COMPLEX_SRCS = cabs.c cacos.c cacosh.c carg.c casin.c casinh.c catan.c \
+COMPLEX_SRCS = cabs.c carg.c \
 	ccos.c ccosh.c cephes_subr.c cexp.c clog.c conj.c cpow.c cproj.c \
-	cimag.c creal.c csin.c csinh.c csqrt.c ctan.c ctanh.c catanh.c
+	cimag.c creal.c csin.c csinh.c csqrt.c ctan.c ctanh.c \
+	catrig.c
+CATRIG_SRCS = cacos.c cacosh.c casin.c casinh.c catan.c catanh.c
 
+CPPFLAGS+=-I${.CURDIR}/src
 .for i in ${COMPLEX_SRCS}
 SRCS+=	$i ${i:S/.c/f.c/} ${i:S/.c/l.c/}
-MAN+= ${i:Ncephes_*:S/.c/.3/}
-MLINKS+= ${i:Ncephes_*:S/.c/.3/} ${i:Ncephes_*:S/.c/f.3/}
-MLINKS+= ${i:Ncephes_*:S/.c/.3/} ${i:Ncephes_*:S/.c/l.3/}
+MAN+= ${i:Ncatrig*:Ncephes_*:S/.c/.3/}
+MLINKS+= ${i:Ncatrig*:Ncephes_*:S/.c/.3/} ${i:Ncatrig*:Ncephes_*:S/.c/f.3/}
+MLINKS+= ${i:Ncatrig*:Ncephes_*:S/.c/.3/} ${i:Ncatrig*:Ncephes_*:S/.c/l.3/}
 .endfor
+
+.for i in ${CATRIG_SRCS}
+MAN+= ${i:S/.c/.3/}
+MLINKS+= ${i:S/.c/.3/} ${i:S/.c/f.3/}
+MLINKS+= ${i:S/.c/.3/} ${i:S/.c/l.3/}
+.endfor
+
+

Index: src/lib/libm/src/math_private.h
diff -u src/lib/libm/src/math_private.h:1.22 src/lib/libm/src/math_private.h:1.23
--- src/lib/libm/src/math_private.h:1.22	Thu Mar 26 07:59:38 2015
+++ src/lib/libm/src/math_private.h	Mon Sep 19 18:05:05 2016
@@ -11,7 +11,7 @@
 
 /*
  * from: @(#)fdlibm.h 5.1 93/09/24
- * $NetBSD: math_private.h,v 1.22 2015/03/26 11:59:38 justin Exp $
+ * $NetBSD: math_private.h,v 1.23 2016/09/19 22:05:05 christos Exp $
  */
 
 #ifndef _MATH_PRIVATE_H_
@@ -223,6 +223,57 @@ typedef union {
 #define	REAL_PART(z)	((z).parts[0])
 #define	IMAG_PART(z)	((z).parts[1])
 
+/*
+ * Inline functions that can be used to construct complex values.
+ *
+ * The C99 standard intends x+I*y to be used for this, but x+I*y is
+ * currently unusable in general since gcc introduces many overflow,
+ * underflow, sign and efficiency bugs by rewriting I*y as
+ * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
+ * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
+ * to -0.0+I*0.0.
+ *
+ * The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL()
+ * to construct complex values.  Compilers that conform to the C99
+ * standard require the following functions to avoid the above issues.
+ */
+
+#ifndef CMPLXF
+static __inline float complex
+CMPLXF(float x, float y)
+{
+	float_complex z;
+
+	REAL_PART(z) = x;
+	IMAG_PART(z) = y;
+	return (z.z);
+}
+#endif
+
+#ifndef CMPLX
+static __inline double complex
+CMPLX(double x, double y)
+{
+	double_complex z;
+
+	REAL_PART(z) = x;
+	IMAG_PART(z) = y;
+	return (z.z);
+}
+#endif
+
+#ifndef CMPLXL
+static __inline long double complex
+CMPLXL(long double x, long double y)
+{
+	long_double_complex z;
+
+	REAL_PART(z) = x;
+	IMAG_PART(z) = y;
+	return (z.z);
+}
+#endif
+
 #endif	/* _COMPLEX_H */
 
 /* ieee style elementary functions */

Added files:

Index: src/lib/libm/complex/catrig.c
diff -u /dev/null src/lib/libm/complex/catrig.c:1.1
--- /dev/null	Mon Sep 19 18:05:05 2016
+++ src/lib/libm/complex/catrig.c	Mon Sep 19 18:05:05 2016
@@ -0,0 +1,653 @@
+/*	$NetBSD: catrig.c,v 1.1 2016/09/19 22:05:05 christos Exp $	*/
+/*-
+ * Copyright (c) 2012 Stephen Montgomery-Smith <step...@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>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/catrig.c 275819 2014-12-16 09:21:56Z ed $");
+#endif
+__RCSID("$NetBSD: catrig.c,v 1.1 2016/09/19 22:05:05 christos Exp $");
+
+#include "namespace.h"
+#ifdef __weak_alias
+__weak_alias(casin, _casin)
+#endif
+#ifdef __weak_alias
+__weak_alias(catan, _catan)
+#endif
+
+#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 __unused = /*LINTED*/1 + tiny; } while(/*CONSTCOND*/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 (CMPLX(x, y + y));
+		/* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
+		if (isinf(y))
+			return (CMPLX(y, x + x));
+		/* casinh(NaN + I*0) = NaN + I*0 */
+		if (y == 0)
+			return (CMPLX(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 (CMPLX(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 (CMPLX(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 (CMPLX(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(CMPLX(cimag(z), creal(z)));
+
+	return (CMPLX(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 (CMPLX(y + y, -INFINITY));
+		/* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */
+		if (isinf(y))
+			return (CMPLX(x + x, -y));
+		/* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */
+		if (x == 0)
+			return (CMPLX(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 (CMPLX(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 (CMPLX(rx, ry));
+	}
+
+	/* Avoid spuriously raising inexact for z = 1. */
+	if (x == 1 && y == 0)
+		return (CMPLX(0, -y));
+
+	/* All remaining cases are inexact. */
+	raise_inexact();
+
+	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+		return (CMPLX(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 (CMPLX(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 (CMPLX(ry, rx));
+	/* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */
+	/* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */
+	if (isnan(rx))
+		return (CMPLX(fabs(ry), rx));
+	/* cacosh(0 + I*NaN) = NaN + I*NaN */
+	if (isnan(ry))
+		return (CMPLX(ry, ry));
+	return (CMPLX(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 (CMPLX(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 (CMPLX(log(hypot(x, y)), atan2(y, x)));
+
+	return (CMPLX(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 (CMPLX(atanh(x), y));
+
+	/* To ensure the same accuracy as atan(), and to filter out z = 0. */
+	if (x == 0)
+		return (CMPLX(x, atan(y)));
+
+	if (isnan(x) || isnan(y)) {
+		/* catanh(+-Inf + I*NaN) = +-0 + I*NaN */
+		if (isinf(x))
+			return (CMPLX(copysign(0, x), y + y));
+		/* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
+		if (isinf(y))
+			return (CMPLX(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 (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
+		return (CMPLX(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 (CMPLX(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(CMPLX(cimag(z), creal(z)));
+
+	return (CMPLX(cimag(w), creal(w)));
+}
Index: src/lib/libm/complex/catrigf.c
diff -u /dev/null src/lib/libm/complex/catrigf.c:1.1
--- /dev/null	Mon Sep 19 18:05:05 2016
+++ src/lib/libm/complex/catrigf.c	Mon Sep 19 18:05:05 2016
@@ -0,0 +1,406 @@
+/*	$NetBSD: catrigf.c,v 1.1 2016/09/19 22:05:05 christos Exp $	*/
+/*-
+ * Copyright (c) 2012 Stephen Montgomery-Smith <step...@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>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/catrigf.c 275819 2014-12-16 09:21:56Z ed $");
+#endif
+__RCSID("$NetBSD: catrigf.c,v 1.1 2016/09/19 22:05:05 christos Exp $");
+
+#include "namespace.h"
+#ifdef __weak_alias
+__weak_alias(casinf, _casinf)
+#endif
+#ifdef __weak_alias
+__weak_alias(catanf, _catanf)
+#endif
+
+
+#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 __unused = /*LINTED*/1 + tiny; } while(/*CONSTCOND*/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 (CMPLXF(x, y + y));
+		if (isinf(y))
+			return (CMPLXF(y, x + x));
+		if (y == 0)
+			return (CMPLXF(x + x, y));
+		return (CMPLXF(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 (CMPLXF(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 (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
+}
+
+float complex
+casinf(float complex z)
+{
+	float complex w = casinhf(CMPLXF(cimagf(z), crealf(z)));
+
+	return (CMPLXF(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 (CMPLXF(y + y, -INFINITY));
+		if (isinf(y))
+			return (CMPLXF(x + x, -y));
+		if (x == 0)
+			return (CMPLXF(pio2_hi + pio2_lo, y + y));
+		return (CMPLXF(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 (CMPLXF(rx, ry));
+	}
+
+	if (x == 1 && y == 0)
+		return (CMPLXF(0, -y));
+
+	raise_inexact();
+
+	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
+		return (CMPLXF(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 (CMPLXF(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 (CMPLXF(ry, rx));
+	if (isnan(rx))
+		return (CMPLXF(fabsf(ry), rx));
+	if (isnan(ry))
+		return (CMPLXF(ry, ry));
+	return (CMPLXF(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 (CMPLXF(logf(hypotf(x / m_e, y / m_e)) + 1,
+		    atan2f(y, x)));
+
+	if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
+		return (CMPLXF(logf(hypotf(x, y)), atan2f(y, x)));
+
+	return (CMPLXF(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 (CMPLXF(atanhf(x), y));
+
+	if (x == 0)
+		return (CMPLXF(x, atanf(y)));
+
+	if (isnan(x) || isnan(y)) {
+		if (isinf(x))
+			return (CMPLXF(copysignf(0, x), y + y));
+		if (isinf(y))
+			return (CMPLXF(copysignf(0, x),
+			    copysignf(pio2_hi + pio2_lo, y)));
+		return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
+		return (CMPLXF(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 (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
+}
+
+float complex
+catanf(float complex z)
+{
+	float complex w = catanhf(CMPLXF(cimagf(z), crealf(z)));
+
+	return (CMPLXF(cimagf(w), crealf(w)));
+}
Index: src/lib/libm/complex/catrigl.c
diff -u /dev/null src/lib/libm/complex/catrigl.c:1.1
--- /dev/null	Mon Sep 19 18:05:05 2016
+++ src/lib/libm/complex/catrigl.c	Mon Sep 19 18:05:05 2016
@@ -0,0 +1,443 @@
+/*	$NetBSD: catrigl.c,v 1.1 2016/09/19 22:05:05 christos Exp $	*/
+/*-
+ * Copyright (c) 2012 Stephen Montgomery-Smith <step...@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.
+ *
+ * The code for catrig.c contains complete comments.
+ */
+#include <sys/cdefs.h>
+__RCSID("$NetBSD: catrigl.c,v 1.1 2016/09/19 22:05:05 christos Exp $");
+
+#include "namespace.h"
+#ifdef __weak_alias
+__weak_alias(casinl, _casinl)
+#endif
+#ifdef __weak_alias
+__weak_alias(catanl, _catanl)
+#endif
+
+
+#include <complex.h>
+#include <float.h>
+#ifdef __HAVE_LONG_DOUBLE
+
+#include "math.h"
+#include "math_private.h"
+
+#undef isinf
+#define isinf(x)	(fabsl(x) == INFINITY)
+#undef isnan
+#define isnan(x)	((x) != (x))
+#define	raise_inexact()	do { volatile float junk __unused = /*LINTED*/1 + tiny; } while(/*CONSTCOND*/0)
+#undef signbit
+#define signbit(x)	(__builtin_signbitl(x)) 
+
+#if __HAVE_LONG_DOUBLE + 0 == 128
+// Ok
+#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
+// XXX: Byte order
+struct ieee_ext {
+	uint64_t ext_frac;
+	uint16_t ext_exp:15;
+	uint16_t ext_sign:1;
+	uint16_t ext_pad;
+};
+#define extu_exp	extu_ext.ext_exp
+#define extu_sign	extu_ext.ext_sign
+#define extu_frac	extu_ext.ext_frac
+union ieee_ext_u {
+	long double extu_ld;
+	struct ieee_ext extu_ext;
+};
+#else
+	#error "unsupported long double format"
+#endif
+
+#define GET_LDBL_EXPSIGN(r, s) \
+    do { \
+	    union ieee_ext_u u; \
+	    u.extu_ld = s; \
+	    r = u.extu_sign; \
+	    r >>= EXT_EXPBITS - 1;
+    } while (/*CONSTCOND*/0)
+#define SET_LDBL_EXPSIGN(r, s) \
+    do { \
+	    union ieee_ext_u u; \
+	    u.extu_ld = s; \
+	    u.extu_exp &= __BITS(0, EXT_EXPBITS - 1); \
+	    u.extu_exp |= r << (EXT_EXPBITS - 1); \
+	    s = u.extu_ld; \
+    } while (/*CONSTCOND*/0)
+
+static const long double
+A_crossover =		10,
+B_crossover =		0.6417,
+FOUR_SQRT_MIN =		0x1p-8189L,
+QUARTER_SQRT_MAX =	0x1p8189L,
+RECIP_EPSILON =		1/LDBL_EPSILON,
+SQRT_MIN =		0x1p-8191L;
+
+static const long double
+m_e =		2.71828182845904523536028747135266250e0L,	/* 0x15bf0a8b1457695355fb8ac404e7a.0p-111 */
+m_ln2 =		6.93147180559945309417232121458176568e-1L,	/* 0x162e42fefa39ef35793c7673007e6.0p-113 */
+pio2_hi =      1.5707963267948966192313216916397514L, /* pi/2 */
+SQRT_3_EPSILON = 2.40370335797945490975336727199878124e-17L,	/*  0x1bb67ae8584caa73b25742d7078b8.0p-168 */
+SQRT_6_EPSILON = 3.39934988877629587239082586223300391e-17L;	/*  0x13988e1409212e7d0321914321a55.0p-167 */
+
+static const volatile double
+pio2_lo =               6.1232339957367659e-17; /*  0x11a62633145c07.0p-106 */
+static const volatile float
+tiny =			0x1p-100;
+
+static long double complex clog_for_large_values(long double complex z);
+
+inline static long double
+f(long double a, long double b, long 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);
+}
+
+inline static void
+do_hard_work(long double x, long double y, long double *rx, int *B_is_usable, long double *B, long double *sqrt_A2my2, long double *new_y)
+{
+	long double R, S, A;
+	long double Am1, Amy;
+
+	R = hypotl(x, y+1);
+	S = hypotl(x, y-1);
+
+	A = (R + S) / 2;
+	if (A < 1)
+		A = 1;
+
+	if (A < A_crossover) {
+		if (y == 1 && x < LDBL_EPSILON*LDBL_EPSILON/128) {
+			*rx = sqrtl(x);
+		} else if (x >= LDBL_EPSILON * fabsl(y-1)) {
+			Am1 = f(x, 1+y, R) + f(x, 1-y, S);
+			*rx = log1pl(Am1 + sqrtl(Am1*(A+1)));
+		} else if (y < 1) {
+			*rx = x/sqrtl((1-y)*(1+y));
+		} else {
+			*rx = log1pl((y-1) + sqrtl((y-1)*(y+1)));
+		}
+	} else
+		*rx = logl(A + sqrtl(A*A-1));
+
+	*new_y = y;
+
+	if (y < FOUR_SQRT_MIN) {
+		*B_is_usable = 0;
+		*sqrt_A2my2 = A * (2 / LDBL_EPSILON);
+		*new_y= y * (2 / LDBL_EPSILON);
+		return;
+	}
+
+	*B = y/A;
+	*B_is_usable = 1;
+
+	if (*B > B_crossover) {
+		*B_is_usable = 0;
+		if (y == 1 && x < LDBL_EPSILON/128) {
+			*sqrt_A2my2 = sqrtl(x)*sqrtl((A+y)/2);
+		} else if (x >= LDBL_EPSILON * fabsl(y-1)) {
+			Amy = f(x, y+1, R) + f(x, y-1, S);
+			*sqrt_A2my2 = sqrtl(Amy*(A+y));
+		} else if (y > 1) {
+			*sqrt_A2my2 = x * (4/LDBL_EPSILON/LDBL_EPSILON) * y /
+			    sqrtl((y+1)*(y-1));
+			*new_y = y * (4/LDBL_EPSILON/LDBL_EPSILON);
+		} else {
+			*sqrt_A2my2 = sqrtl((1-y)*(1+y));
+		}
+	}
+}
+
+long double complex
+casinhl(long double complex z)
+{
+	long double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;
+	int B_is_usable;
+	long double complex w;
+
+	x = creall(z);
+	y = cimagl(z);
+	ax = fabsl(x);
+	ay = fabsl(y);
+
+	if (isnan(x) || isnan(y)) {
+		if (isinf(x))
+			return (CMPLXL(x, y+y));
+		if (isinf(y))
+			return (CMPLXL(y, x+x));
+		if (y == 0) return (CMPLXL(x+x, y));
+		return (CMPLXL(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 (CMPLXL(copysignl(creall(w), x), copysignl(cimagl(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 = asinl(B);
+	else
+		ry = atan2l(new_y, sqrt_A2my2);
+	return (CMPLXL(copysignl(rx, x), copysignl(ry, y)));
+}
+
+long double complex
+casinl(long double complex z)
+{
+	long double complex w = casinhl(CMPLXL(cimagl(z), creall(z)));
+	return (CMPLXL(cimagl(w), creall(w)));
+}
+
+long double complex
+cacosl(long double complex z)
+{
+	long double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
+	int sx, sy;
+	int B_is_usable;
+	long double complex w;
+
+	x = creall(z);
+	y = cimagl(z);
+	sx = signbit(x);
+	sy = signbit(y);
+	ax = fabsl(x);
+	ay = fabsl(y);
+
+	if (isnan(x) || isnan(y)) {
+		if (isinf(x))
+			return (CMPLXL(y+y, -INFINITY));
+		if (isinf(y))
+			return (CMPLXL(x+x, -y));
+		if (x == 0) return (CMPLXL(pio2_hi + pio2_lo, y+y));
+		return (CMPLXL(x+0.0L+(y+0), x+0.0L+(y+0)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+		w = clog_for_large_values(z);
+		rx = fabsl(cimagl(w));
+		ry = creall(w) + m_ln2;
+		if (sy == 0)
+			ry = -ry;
+		return (CMPLXL(rx, ry));
+	}
+
+	if (x == 1 && y == 0)
+		return (CMPLXL(0, -y));
+
+	raise_inexact();
+
+	if (ax < SQRT_6_EPSILON/4 && ay < SQRT_6_EPSILON/4)
+		return (CMPLXL(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 = acosl(B);
+		else
+			rx = acosl(-B);
+	} else {
+		if (sx==0)
+			rx = atan2l(sqrt_A2mx2, new_x);
+		else
+			rx = atan2l(sqrt_A2mx2, -new_x);
+	}
+	if (sy==0)
+		ry = -ry;
+	return (CMPLXL(rx, ry));
+}
+
+long double complex
+cacoshl(long double complex z)
+{
+	long double complex w;
+	long double rx, ry;
+
+	w = cacosl(z);
+	rx = creall(w);
+	ry = cimagl(w);
+	if (isnan(rx) && isnan(ry))
+		return (CMPLXL(ry, rx));
+	if (isnan(rx))
+		return (CMPLXL(fabsl(ry), rx));
+	if (isnan(ry))
+		return (CMPLXL(ry, ry));
+	return (CMPLXL(fabsl(ry), copysignl(rx, cimagl(z))));
+}
+
+static long double complex
+clog_for_large_values(long double complex z)
+{
+	long double x, y;
+	long double ax, ay, t;
+
+	x = creall(z);
+	y = cimagl(z);
+	ax = fabsl(x);
+	ay = fabsl(y);
+	if (ax < ay) {
+		t = ax;
+		ax = ay;
+		ay = t;
+	}
+
+	if (ax > LDBL_MAX / 2)
+		return (CMPLXL(logl(hypotl(x / m_e, y / m_e)) + 1, atan2l(y, x)));
+
+	if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
+		return (CMPLXL(logl(hypotl(x, y)), atan2l(y, x)));
+
+	return (CMPLXL(logl(ax*ax + ay*ay) / 2, atan2l(y, x)));
+}
+
+inline static long double
+sum_squares(long double x, long double y)
+{
+	if (y < SQRT_MIN)
+		return (x*x);
+
+	return (x*x + y*y);
+}
+
+inline static long double
+real_part_reciprocal(long double x, long double y)
+{
+	long double scale;
+	uint16_t hx, hy;
+	int16_t ix, iy;
+
+	GET_LDBL_EXPSIGN(hx, x);
+	ix = hx & 0x7fff;
+	GET_LDBL_EXPSIGN(hy, y);
+	iy = hy & 0x7fff;
+#define	BIAS	(LDBL_MAX_EXP - 1)
+#define	CUTOFF	(LDBL_MANT_DIG / 2 + 1)
+	if (ix - iy >= CUTOFF || isinf(x))
+		return (1/x);
+	if (iy - ix >= CUTOFF)
+		return (x/y/y);
+	if (ix <= BIAS + LDBL_MAX_EXP / 2 - CUTOFF)
+		return (x/(x*x + y*y));
+	scale = 1;
+	SET_LDBL_EXPSIGN(scale, 0x7fff - ix);
+	x *= scale;
+	y *= scale;
+	return (x/(x*x + y*y) * scale);
+}
+
+long double complex
+catanhl(long double complex z)
+{
+	long double x, y, ax, ay, rx, ry;
+
+	x = creall(z);
+	y = cimagl(z);
+	ax = fabsl(x);
+	ay = fabsl(y);
+
+	if (y == 0 && ax <= 1)
+		return (CMPLXL(atanhl(x), y)); 	/* XXX need atanhl() */
+
+	if (x == 0)
+		return (CMPLXL(x, atanl(y)));
+
+	if (isnan(x) || isnan(y)) {
+		if (isinf(x))
+			return (CMPLXL(copysignl(0, x), y+y));
+		if (isinf(y))
+			return (CMPLXL(copysignl(0, x), copysignl(pio2_hi + pio2_lo, y)));
+		return (CMPLXL(x+0.0L+(y+0), x+0.0L+(y+0)));
+	}
+
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
+		return (CMPLXL(real_part_reciprocal(x, y), copysignl(pio2_hi + pio2_lo, y)));
+
+	if (ax < SQRT_3_EPSILON/2 && ay < SQRT_3_EPSILON/2) {
+		raise_inexact();
+		return (z);
+	}
+
+	if (ax == 1 && ay < LDBL_EPSILON) {
+#if 0
+		if (ay > 2*LDBL_MIN)
+			rx = - logl(ay/2) / 2;
+		else
+#endif
+			rx = - (logl(ay) - m_ln2) / 2;
+	} else
+		rx = log1pl(4*ax / sum_squares(ax-1, ay)) / 4;
+
+	if (ax == 1)
+		ry = atan2l(2, -ay) / 2;
+	else if (ay < LDBL_EPSILON)
+		ry = atan2l(2*ay, (1-ax)*(1+ax)) / 2;
+	else
+		ry = atan2l(2*ay, (1-ax)*(1+ax) - ay*ay) / 2;
+
+	return (CMPLXL(copysignl(rx, x), copysignl(ry, y)));
+}
+
+long double complex
+catanl(long double complex z)
+{
+	long double complex w = catanhl(CMPLXL(cimagl(z), creall(z)));
+	return (CMPLXL(cimagl(w), creall(w)));
+}
+
+#else
+__strong_alias(_casinl, casin)
+__strong_alias(_catanl, catan)
+__strong_alias(cacoshl, cacosh)
+__strong_alias(cacosl, cacos)
+__strong_alias(casinhl, casinh)
+__strong_alias(catanhl, catanh)
+#endif

Reply via email to