diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index 69011345a0d6..b2502742ac83 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -1,227 +1,227 @@
 #  @(#)Makefile 5.1beta 93/09/24
 # $FreeBSD$
 #
 #  ====================================================
 #  Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 #
 #  Developed at SunPro, a Sun Microsystems, Inc. business.
 #  Permission to use, copy, modify, and distribute this
 #  software is freely granted, provided that this notice
 #  is preserved.
 #  ====================================================
 #
 #
 
 PACKAGE=	clibs
 ARCH_SUBDIR=	${MACHINE_CPUARCH:S/i386/i387/}
 
 .include "${ARCH_SUBDIR}/Makefile.inc"
 
 .PATH:	${.CURDIR}/${ARCH_SUBDIR}
 .if ${MACHINE_CPUARCH} == "i386" || ${MACHINE_CPUARCH} == "amd64"
 .PATH:	${.CURDIR}/x86
 CFLAGS+=	-I${.CURDIR}/x86
 .endif
 
 # long double format
 .if ${LDBL_PREC} == 64
 .PATH:  ${.CURDIR}/ld80
 CFLAGS+=	-I${.CURDIR}/ld80
 .elif ${LDBL_PREC} == 113
 .PATH:  ${.CURDIR}/ld128
 CFLAGS+=	-I${.CURDIR}/ld128
 .endif
 
 CFLAGS+=	-I${.CURDIR}/${ARCH_SUBDIR}
 
 .PATH:	${.CURDIR}/bsdsrc
 .PATH:	${.CURDIR}/src
 .PATH:	${.CURDIR}/man
 
 LIB=	m
 SHLIBDIR?= /lib
 SHLIB_MAJOR= 5
 WARNS?=	1
 IGNORE_PRAGMA=
 COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
 	e_acos.c e_acosf.c e_acosh.c e_acoshf.c e_asin.c e_asinf.c \
 	e_atan2.c e_atan2f.c e_atanh.c e_atanhf.c e_cosh.c e_coshf.c e_exp.c \
 	e_expf.c e_fmod.c e_fmodf.c e_gamma.c e_gamma_r.c e_gammaf.c \
 	e_gammaf_r.c e_hypot.c e_hypotf.c e_j0.c e_j0f.c e_j1.c e_j1f.c \
 	e_jn.c e_jnf.c e_lgamma.c e_lgamma_r.c e_lgammaf.c e_lgammaf_r.c \
 	e_log.c e_log10.c e_log10f.c e_log2.c e_log2f.c e_logf.c \
 	e_pow.c e_powf.c e_rem_pio2.c \
 	e_rem_pio2f.c e_remainder.c e_remainderf.c e_scalb.c e_scalbf.c \
 	e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \
 	imprecise.c \
 	k_cos.c k_cosf.c k_exp.c k_expf.c k_rem_pio2.c k_sin.c k_sinf.c \
 	k_tan.c k_tanf.c \
 	s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c s_cargl.c \
 	s_cbrt.c s_cbrtf.c s_ceil.c s_ceilf.c \
 	s_copysign.c s_copysignf.c s_cos.c s_cosf.c \
 	s_csqrt.c s_csqrtf.c s_erf.c s_erff.c \
 	s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c \
 	s_finite.c s_finitef.c \
 	s_floor.c s_floorf.c s_fma.c s_fmaf.c \
 	s_fmax.c s_fmaxf.c s_fmin.c \
 	s_fminf.c s_frexp.c s_frexpf.c s_ilogb.c s_ilogbf.c \
 	s_ilogbl.c s_isfinite.c s_isnan.c s_isnormal.c \
 	s_llrint.c s_llrintf.c s_llround.c s_llroundf.c s_llroundl.c \
 	s_log1p.c s_log1pf.c s_logb.c s_logbf.c s_lrint.c s_lrintf.c \
 	s_lround.c s_lroundf.c s_lroundl.c s_modff.c \
 	s_nan.c s_nearbyint.c s_nextafter.c s_nextafterf.c \
 	s_nexttowardf.c s_remquo.c s_remquof.c \
 	s_rint.c s_rintf.c s_round.c s_roundf.c \
 	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
 	s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \
 	s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \
 	w_cabs.c w_cabsf.c w_drem.c w_dremf.c
 
 # Location of fpmath.h and _fpmath.h
 .if exists(${LIBCSRCDIR}/${MACHINE_ARCH})
 LIBC_ARCH=${MACHINE_ARCH}
 .else
 LIBC_ARCH=${MACHINE_CPUARCH}
 .endif
 CFLAGS+=	-I${.CURDIR}/src -I${LIBCSRCDIR}/include \
 	-I${LIBCSRCDIR}/${LIBC_ARCH}
 SYM_MAPS+=	${.CURDIR}/Symbol.map
 
 VERSION_DEF=	${LIBCSRCDIR}/Versions.def
 SYMBOL_MAPS=	${SYM_MAPS}
 
 # C99 long double functions
 COMMON_SRCS+=	s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c
 .if ${LDBL_PREC} != 53
 # If long double != double use these; otherwise, we alias the double versions.
 COMMON_SRCS+=	e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \
 	e_coshl.c e_fmodl.c e_hypotl.c \
 	e_lgammal.c e_lgammal_r.c \
 	e_remainderl.c e_sinhl.c e_sqrtl.c \
 	invtrig.c k_cosl.c k_sinl.c k_tanl.c \
 	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \
 	s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
 	s_fmaxl.c s_fminl.c s_frexpl.c s_logbl.c s_logl.c s_nanl.c \
 	s_nextafterl.c s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c \
 	s_scalbnl.c s_sinl.c s_tanhl.c s_tanl.c s_truncl.c w_cabsl.c
 .endif
 
 # C99 complex functions
-COMMON_SRCS+=	catrig.c catrigf.c \
+COMMON_SRCS+=	catrig.c catrigf.c catrigl.c \
 	s_ccosh.c s_ccoshf.c s_cexp.c s_cexpf.c \
 	s_cimag.c s_cimagf.c s_cimagl.c \
 	s_conj.c s_conjf.c s_conjl.c \
 	s_cproj.c s_cprojf.c s_creal.c s_crealf.c s_creall.c \
 	s_csinh.c s_csinhf.c s_ctanh.c s_ctanhf.c
 
 # FreeBSD's C library supplies these functions:
 #COMMON_SRCS+=	s_fabs.c s_frexp.c s_isnan.c s_ldexp.c s_modf.c
 
 # Exclude the generic versions of what we provide in the MD area.
 .if defined(ARCH_SRCS)
 .for i in ${ARCH_SRCS}
 COMMON_SRCS:=  ${COMMON_SRCS:N${i:R}.c}
 .endfor
 .endif
 
 SRCS=	${COMMON_SRCS} ${ARCH_SRCS}
 
 INCS+=	fenv.h math.h
 
 MAN=	acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
 	ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \
 	cimag.3 copysign.3 cos.3 cosh.3 csqrt.3 erf.3 exp.3 fabs.3 fdim.3 \
 	feclearexcept.3 feenableexcept.3 fegetenv.3 \
 	fegetround.3 fenv.3 floor.3 \
 	fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \
 	lgamma.3 log.3 lrint.3 lround.3 math.3 nan.3 \
 	nextafter.3 remainder.3 rint.3 \
 	round.3 scalbn.3 signbit.3 sin.3 sinh.3 sqrt.3 tan.3 tanh.3 trunc.3 \
 	complex.3
 
 MLINKS+=acos.3 acosf.3 acos.3 acosl.3
 MLINKS+=acosh.3 acoshf.3 acosh.3 acoshl.3
 MLINKS+=asin.3 asinf.3 asin.3 asinl.3
 MLINKS+=asinh.3 asinhf.3 asinh.3 asinhl.3
 MLINKS+=atan.3 atanf.3 atan.3 atanl.3
 MLINKS+=atanh.3 atanhf.3 atanh.3 atanhl.3
 MLINKS+=atan2.3 atan2f.3 atan2.3 atan2l.3 \
 	atan2.3 carg.3 atan2.3 cargf.3 atan2.3 cargl.3
 MLINKS+=cacos.3 cacosf.3 cacos.3 cacosh.3 cacos.3 cacoshf.3 \
 	cacos.3 casin.3 cacos.3 casinf.3 cacos.3 casinh.3 cacos.3 casinhf.3 \
 	cacos.3 catan.3 cacos.3 catanf.3 cacos.3 catanh.3 cacos.3 catanhf.3
 MLINKS+=ccos.3 ccosf.3 ccos.3 csin.3 ccos.3 csinf.3 ccos.3 ctan.3 ccos.3 ctanf.3
 MLINKS+=ccosh.3 ccoshf.3 ccosh.3 csinh.3 ccosh.3 csinhf.3 \
 	ccosh.3 ctanh.3 ccosh.3 ctanhf.3
 MLINKS+=ceil.3 ceilf.3 ceil.3 ceill.3
 MLINKS+=cexp.3 cexpf.3
 MLINKS+=cimag.3 cimagf.3 cimag.3 cimagl.3 \
 	cimag.3 conj.3 cimag.3 conjf.3 cimag.3 conjl.3 \
 	cimag.3 cproj.3 cimag.3 cprojf.3 cimag.3 cprojl.3 \
 	cimag.3 creal.3 cimag.3 crealf.3 cimag.3 creall.3
 MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3
 MLINKS+=cos.3 cosf.3 cos.3 cosl.3
 MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3
 MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3
 MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 erf.3 erfl.3 erf.3 erfcl.3
 MLINKS+=exp.3 expm1.3 exp.3 expm1f.3 exp.3 expm1l.3 exp.3 pow.3 exp.3 powf.3 \
 	exp.3 powl.3 exp.3 exp2.3 exp.3 exp2f.3 exp.3 exp2l.3 exp.3 expf.3 \
 	exp.3 expl.3
 MLINKS+=fabs.3 fabsf.3 fabs.3 fabsl.3
 MLINKS+=fdim.3 fdimf.3 fdim.3 fdiml.3
 MLINKS+=feclearexcept.3 fegetexceptflag.3 feclearexcept.3 feraiseexcept.3 \
 	feclearexcept.3 fesetexceptflag.3 feclearexcept.3 fetestexcept.3
 MLINKS+=feenableexcept.3 fedisableexcept.3 feenableexcept.3 fegetexcept.3
 MLINKS+=fegetenv.3 feholdexcept.3 fegetenv.3 fesetenv.3 \
 	fegetenv.3 feupdateenv.3
 MLINKS+=fegetround.3 fesetround.3
 MLINKS+=floor.3 floorf.3 floor.3 floorl.3
 MLINKS+=fma.3 fmaf.3 fma.3 fmal.3
 MLINKS+=fmax.3 fmaxf.3 fmax.3 fmaxl.3 \
 	fmax.3 fmin.3 fmax.3 fminf.3 fmax.3 fminl.3
 MLINKS+=fmod.3 fmodf.3 fmod.3 fmodl.3
 MLINKS+=hypot.3 cabs.3 hypot.3 cabsf.3 hypot.3 cabsl.3 \
 	hypot.3 hypotf.3 hypot.3 hypotl.3
 MLINKS+=ieee_test.3 scalb.3 ieee_test.3 scalbf.3
 MLINKS+=ieee_test.3 significand.3 ieee_test.3 significandf.3
 MLINKS+=ilogb.3 ilogbf.3 ilogb.3 ilogbl.3 \
 	ilogb.3 logb.3 ilogb.3 logbf.3 ilogb.3 logbl.3
 MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3
 MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3
 MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 \
 	lgamma.3 lgammaf.3 lgamma.3 lgammal.3 \
 	lgamma.3 tgamma.3 lgamma.3 tgammaf.3
 MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log10l.3 \
 	log.3 log1p.3 log.3 log1pf.3 log.3 log1pl.3 \
 	log.3 logf.3 log.3 logl.3 \
 	log.3 log2.3 log.3 log2f.3 log.3 log2l.3
 MLINKS+=lrint.3 llrint.3 lrint.3 llrintf.3 lrint.3 llrintl.3 \
 	lrint.3 lrintf.3 lrint.3 lrintl.3
 MLINKS+=lround.3 llround.3 lround.3 llroundf.3 lround.3 llroundl.3 \
 	lround.3 lroundf.3 lround.3 lroundl.3
 MLINKS+=nan.3 nanf.3 nan.3 nanl.3
 MLINKS+=nextafter.3 nextafterf.3 nextafter.3 nextafterl.3
 MLINKS+=nextafter.3 nexttoward.3 nextafter.3 nexttowardf.3
 MLINKS+=nextafter.3 nexttowardl.3
 MLINKS+=remainder.3 remainderf.3 remainder.3 remainderl.3 \
 	remainder.3 remquo.3 remainder.3 remquof.3 remainder.3 remquol.3
 MLINKS+=rint.3 rintf.3 rint.3 rintl.3 \
 	rint.3 nearbyint.3 rint.3 nearbyintf.3 rint.3 nearbyintl.3
 MLINKS+=round.3 roundf.3 round.3 roundl.3
 MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3
 MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
 MLINKS+=sin.3 sinf.3 sin.3 sinl.3
 MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3
 MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
 	sqrt.3 sqrtl.3
 MLINKS+=tan.3 tanf.3 tan.3 tanl.3
 MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3
 MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3
 
 .include <src.opts.mk>
 
 .if ${MK_TESTS} != "no"
 SUBDIR+=	tests
 .endif
 
 .include <bsd.lib.mk>
diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map
index d5b7f4699f10..5ba7c49c1413 100644
--- a/lib/msun/Symbol.map
+++ b/lib/msun/Symbol.map
@@ -1,287 +1,297 @@
 /*
  * $FreeBSD$
  */
 
 /* 7.0-CURRENT */
 FBSD_1.0 {
 	__fe_dfl_env;
 	tgamma;
 	acos;
 	acosf;
 	acosh;
 	acoshf;
 	asin;
 	asinf;
 	atan2;
 	atan2f;
 	atanh;
 	atanhf;
 	cosh;
 	coshf;
 	exp;
 	expf;
 	fmod;
 	fmodf;
 	gamma;
 	gamma_r;
 	gammaf;
 	gammaf_r;
 	hypot;
 	hypotf;
 	j0;
 	y0;
 	j0f;
 	y0f;
 	j1;
 	y1;
 	j1f;
 	y1f;
 	jn;
 	yn;
 	jnf;
 	ynf;
 	lgamma;
 	lgamma_r;
 	lgammaf;
 	lgammaf_r;
 	log;
 	log10;
 	log10f;
 	logf;
 	pow;
 	powf;
 	remainder;
 	remainderf;
 	scalb;
 	scalbf;
 	sinh;
 	sinhf;
 	sqrt;
 	sqrtf;
 	asinh;
 	asinhf;
 	atan;
 	atanf;
 	cbrt;
 	cbrtf;
 	ceil;
 	ceilf;
 	ceill;
 	cimag;
 	cimagf;
 	cimagl;
 	conj;
 	conjf;
 	conjl;
 	copysign;
 	copysignf;
 	copysignl;
 	cos;
 	cosf;
 	creal;
 	crealf;
 	creall;
 	erf;
 	erfc;
 	erff;
 	erfcf;
 	exp2;
 	exp2f;
 	expm1;
 	expm1f;
 	fabs;
 	fabsf;
 	fabsl;
 	fdim;
 	fdimf;
 	fdiml;
 	finite;
 	finitef;
 	floor;
 	floorf;
 	floorl;
 	fma;
 	fmaf;
 	fmal;
 	fmax;
 	fmaxf;
 	fmaxl;
 	fmin;
 	fminf;
 	fminl;
 	frexp;
 	frexpf;
 	frexpl;
 	ilogb;
 	ilogbf;
 	ilogbl;
 	__isfinite;
 	__isfinitef;
 	__isfinitel;
 	isnanf;
 	__isnanl;
 	__isnormal;
 	__isnormalf;
 	__isnormall;
 	llrint;
 	llrintf;
 	llround;
 	llroundf;
 	llroundl;
 	log1p;
 	log1pf;
 	logb;
 	logbf;
 	lrint;
 	lrintf;
 	lround;
 	lroundf;
 	lroundl;
 	modff;
 	modfl;
 	nearbyint;
 	nearbyintf;
 	nextafter;
 	nexttoward;
 	nexttowardl;
 	nextafterl;
 	nextafterf;
 	nexttowardf;
 	remquo;
 	remquof;
 	rint;
 	rintf;
 	round;
 	roundf;
 	roundl;
 	scalbln;
 	scalblnf;
 	scalblnl;
 	scalbn;
 	scalbnl;
 	scalbnf;
 	ldexpf;
 	ldexpl;
 	__signbit;
 	__signbitf;
 	__signbitl;
 	signgam;
 	significand;
 	significandf;
 	sin;
 	sinf;
 	tan;
 	tanf;
 	tanh;
 	tanhf;
 	trunc;
 	truncf;
 	truncl;
 	cabs;
 	cabsf;
 	drem;
 	dremf;
 };
 
 /* First added in 8.0-CURRENT */
 FBSD_1.1 {
 	carg;
 	cargf;
 	csqrt;
 	csqrtf;
 	logbl;
 	nan;
 	nanf;
 	nanl;
 	llrintl;
 	lrintl;
 	nearbyintl;
 	rintl;
 	exp2l;
 	sinl;
 	cosl;
 	tanl;
 	tgammaf;
 	sqrtl;
 	hypotl;
 	cabsl;
 	csqrtl;
 	remquol;
 	remainderl;
 	fmodl;
 	acosl;
 	asinl;
 	atan2l;
 	atanl;
 	cargl;
 	cproj;
 	cprojf;
 	cprojl;
 };
 
 /* First added in 9.0-CURRENT */
 FBSD_1.2 {
 	__isnanf;
 	cbrtl;
 	cexp;
 	cexpf;
 	log2;
 	log2f;
 };
 
 /* First added in 10.0-CURRENT */
 FBSD_1.3 {
 	feclearexcept;
 	fegetexceptflag;
 	fetestexcept;
 	fegetround;
 	fesetround;
 	fesetenv;
 	acoshl;
 	asinhl;
 	atanhl;
 	cacos;
 	cacosf;
 	cacosh;
 	cacoshf;
 	casin;
 	casinf;
 	casinh;
 	casinhf;
 	catan;
 	catanf;
 	catanh;
 	catanhf;
 	csin;
 	csinf;
 	csinh;
 	csinhf;
 	ccos;
 	ccosf;
 	ccosh;
 	ccoshf;
 	coshl;
 	ctan;
 	ctanf;
 	ctanh;
 	ctanhf;
 	erfcl;
 	erfl;
 	expl;
 	expm1l;
 	lgammal;
 	log10l;
 	log1pl;
 	log2l;
 	logl;
 	sinhl;
 	tanhl;
 	/* Implemented as weak aliases for imprecise versions */
 	powl;
 	tgammal;
 };
 
 /* First added in 11.0-CURRENT */
 FBSD_1.4 {
 	lgammal_r;
 };
+
+/* First added in 12.0-CURRENT */
+FBSD_1.5 {
+	cacoshl;
+	cacosl;
+	casinhl;
+	casinl;
+	catanl;
+	catanhl;
+};
diff --git a/lib/msun/src/catrigl.c b/lib/msun/src/catrigl.c
new file mode 100644
index 000000000000..f71a3118de7e
--- /dev/null
+++ b/lib/msun/src/catrigl.c
@@ -0,0 +1,412 @@
+/*-
+ * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
+ * Copyright (c) 2017 Mahdi Mokhtari <mmokhi@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 "invtrig.h"
+#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 = 1 + tiny; } while(0)
+#undef signbit
+#define signbit(x)	(__builtin_signbitl(x))
+
+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;
+
+#if LDBL_MANT_DIG == 64
+static const union IEEEl2bits
+um_e =		LD80C(0xadf85458a2bb4a9b,  1, 2.71828182845904523536e+0L),
+um_ln2 =	LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
+#define		m_e	um_e.e
+#define		m_ln2	um_ln2.e
+static const long double
+/* The next 2 literals for non-i386.  Misrounding them on i386 is harmless. */
+SQRT_3_EPSILON = 5.70316273435758915310e-10,	/*  0x9cc470a0490973e8.0p-94 */
+SQRT_6_EPSILON = 8.06549008734932771664e-10;	/*  0xddb3d742c265539e.0p-94 */
+#elif LDBL_MANT_DIG == 113
+static const long double
+m_e =		2.71828182845904523536028747135266250e0L,	/* 0x15bf0a8b1457695355fb8ac404e7a.0p-111 */
+m_ln2 =		6.93147180559945309417232121458176568e-1L,	/* 0x162e42fefa39ef35793c7673007e6.0p-113 */
+SQRT_3_EPSILON = 2.40370335797945490975336727199878124e-17,	/*  0x1bb67ae8584caa73b25742d7078b8.0p-168 */
+SQRT_6_EPSILON = 3.39934988877629587239082586223300391e-17;	/*  0x13988e1409212e7d0321914321a55.0p-167 */
+#else
+#error "Unsupported long double format"
+#endif
+
+static const volatile float
+tiny =			0x1p-100;
+
+static long double complex clog_for_large_values(long double complex z);
+
+static inline 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);
+}
+
+static inline 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;
+
+	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)));
+}
+
+static inline long double
+sum_squares(long double x, long double y)
+{
+
+	if (y < SQRT_MIN)
+		return (x * x);
+
+	return (x * x + y * y);
+}
+
+static inline 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));
+
+	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)
+		rx = (m_ln2 - logl(ay)) / 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;
+
+	w = catanhl(CMPLXL(cimagl(z), creall(z)));
+	return (CMPLXL(cimagl(w), creall(w)));
+}