Merge sync of head

45 files changed, 649 insertions, 328 deletions

diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index afb7067..094fcba 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -223,6 +223,8 @@ MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3
 
 .include <src.opts.mk>
 
-.include <bsd.arch.inc.mk>
+.if ${MK_TESTS} != "no"
+SUBDIR+=	tests
+.endif
 
 .include <bsd.lib.mk>
diff --git a/lib/msun/Makefile.amd64 b/lib/msun/Makefile.amd64
deleted file mode 100644
index dd0f5b0..0000000
--- a/lib/msun/Makefile.amd64
+++ /dev/null
@@ -1,6 +0,0 @@
-# $FreeBSD$
-
-.if ${MK_TESTS} != "no"
-SUBDIR+=	tests
-.endif
-
diff --git a/lib/msun/Makefile.i386 b/lib/msun/Makefile.i386
deleted file mode 100644
index dd0f5b0..0000000
--- a/lib/msun/Makefile.i386
+++ /dev/null
@@ -1,6 +0,0 @@
-# $FreeBSD$
-
-.if ${MK_TESTS} != "no"
-SUBDIR+=	tests
-.endif
-
diff --git a/lib/msun/aarch64/Makefile.inc b/lib/msun/aarch64/Makefile.inc
new file mode 100644
index 0000000..286a608
--- /dev/null
+++ b/lib/msun/aarch64/Makefile.inc
@@ -0,0 +1,4 @@
+# $FreeBSD$
+
+LDBL_PREC = 113
+
diff --git a/lib/msun/aarch64/fenv.c b/lib/msun/aarch64/fenv.c
new file mode 100644
index 0000000..4ec362f
--- /dev/null
+++ b/lib/msun/aarch64/fenv.c
@@ -0,0 +1,56 @@
+/*-
+ * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
+ * Copyright (c) 2013 Andrew Turner <andrew@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.
+ *
+ * $FreeBSD$
+ */
+
+#define	__fenv_static
+#include "fenv.h"
+
+/*
+ * Hopefully the system ID byte is immutable, so it's valid to use
+ * this as a default environment.
+ */
+const fenv_t __fe_dfl_env = 0;
+
+#ifdef __GNUC_GNU_INLINE__
+#error "This file must be compiled with C99 'inline' semantics"
+#endif
+
+extern inline int feclearexcept(int __excepts);
+extern inline int fegetexceptflag(fexcept_t *__flagp, int __excepts);
+extern inline int fesetexceptflag(const fexcept_t *__flagp, int __excepts);
+extern inline int feraiseexcept(int __excepts);
+extern inline int fetestexcept(int __excepts);
+extern inline int fegetround(void);
+extern inline int fesetround(int __round);
+extern inline int fegetenv(fenv_t *__envp);
+extern inline int feholdexcept(fenv_t *__envp);
+extern inline int fesetenv(const fenv_t *__envp);
+extern inline int feupdateenv(const fenv_t *__envp);
+extern inline int feenableexcept(int __mask);
+extern inline int fedisableexcept(int __mask);
+extern inline int fegetexcept(void);
diff --git a/lib/msun/aarch64/fenv.h b/lib/msun/aarch64/fenv.h
new file mode 100644
index 0000000..2bd2987
--- /dev/null
+++ b/lib/msun/aarch64/fenv.h
@@ -0,0 +1,246 @@
+/*-
+ * Copyright (c) 2004-2005 David Schultz <das@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.
+ *
+ * $FreeBSD$
+ */
+
+#ifndef	_FENV_H_
+#define	_FENV_H_
+
+#include <sys/_types.h>
+
+#ifndef	__fenv_static
+#define	__fenv_static	static
+#endif
+
+typedef	__uint64_t	fenv_t;
+typedef	__uint64_t	fexcept_t;
+
+/* Exception flags */
+#define	FE_INVALID	0x00000001
+#define	FE_DIVBYZERO	0x00000002
+#define	FE_OVERFLOW	0x00000004
+#define	FE_UNDERFLOW	0x00000008
+#define	FE_INEXACT	0x00000010
+#define	FE_ALL_EXCEPT	(FE_DIVBYZERO | FE_INEXACT | \
+			 FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW)
+
+/*
+ * Rounding modes
+ *
+ * We can't just use the hardware bit values here, because that would
+ * make FE_UPWARD and FE_DOWNWARD negative, which is not allowed.
+ */
+#define	FE_TONEAREST	0x0
+#define	FE_UPWARD	0x1
+#define	FE_DOWNWARD	0x2
+#define	FE_TOWARDZERO	0x3
+#define	_ROUND_MASK	(FE_TONEAREST | FE_DOWNWARD | \
+			 FE_UPWARD | FE_TOWARDZERO)
+#define	_ROUND_SHIFT	22
+
+__BEGIN_DECLS
+
+/* Default floating-point environment */
+extern const fenv_t	__fe_dfl_env;
+#define	FE_DFL_ENV	(&__fe_dfl_env)
+
+/* We need to be able to map status flag positions to mask flag positions */
+#define _FPUSW_SHIFT	8
+#define	_ENABLE_MASK	(FE_ALL_EXCEPT << _FPUSW_SHIFT)
+
+#define	__mrs_fpcr(__r)	__asm __volatile("mrs %0, fpcr" : "=r" (__r))
+#define	__msr_fpcr(__r)	__asm __volatile("msr fpcr, %0" : : "r" (__r))
+
+#define	__mrs_fpsr(__r)	__asm __volatile("mrs %0, fpsr" : "=r" (__r))
+#define	__msr_fpsr(__r)	__asm __volatile("msr fpsr, %0" : : "r" (__r))
+
+__fenv_static __inline int
+feclearexcept(int __excepts)
+{
+	fexcept_t __r;
+
+	__mrs_fpsr(__r);
+	__r &= ~__excepts;
+	__msr_fpsr(__r);
+	return (0);
+}
+
+__fenv_static inline int
+fegetexceptflag(fexcept_t *__flagp, int __excepts)
+{
+	fexcept_t __r;
+
+	__mrs_fpsr(__r);
+	*__flagp = __r & __excepts;
+	return (0);
+}
+
+__fenv_static inline int
+fesetexceptflag(const fexcept_t *__flagp, int __excepts)
+{
+	fexcept_t __r;
+
+	__mrs_fpsr(__r);
+	__r &= ~__excepts;
+	__r |= *__flagp & __excepts;
+	__msr_fpsr(__r);
+	return (0);
+}
+
+__fenv_static inline int
+feraiseexcept(int __excepts)
+{
+	fexcept_t __r;
+
+	__mrs_fpsr(__r);
+	__r |= __excepts;
+	__msr_fpsr(__r);
+	return (0);
+}
+
+__fenv_static inline int
+fetestexcept(int __excepts)
+{
+	fexcept_t __r;
+
+	__mrs_fpsr(__r);
+	return (__r & __excepts);
+}
+
+__fenv_static inline int
+fegetround(void)
+{
+	fenv_t __r;
+
+	__mrs_fpcr(__r);
+	return ((__r >> _ROUND_SHIFT) & _ROUND_MASK);
+}
+
+__fenv_static inline int
+fesetround(int __round)
+{
+	fenv_t __r;
+
+	if (__round & ~_ROUND_MASK)
+		return (-1);
+	__mrs_fpcr(__r);
+	__r &= ~(_ROUND_MASK << _ROUND_SHIFT);
+	__r |= __round << _ROUND_SHIFT;
+	__msr_fpcr(__r);
+	return (0);
+}
+
+__fenv_static inline int
+fegetenv(fenv_t *__envp)
+{
+	fenv_t __r;
+
+	__mrs_fpcr(__r);
+	*__envp = __r & _ENABLE_MASK;
+
+	__mrs_fpsr(__r);
+	*__envp |= __r & (FE_ALL_EXCEPT | (_ROUND_MASK << _ROUND_SHIFT));
+
+	return (0);
+}
+
+__fenv_static inline int
+feholdexcept(fenv_t *__envp)
+{
+	fenv_t __r;
+
+	__mrs_fpcr(__r);
+	*__envp = __r & _ENABLE_MASK;
+	__r &= ~(_ENABLE_MASK);
+	__msr_fpcr(__r);
+
+	__mrs_fpsr(__r);
+	*__envp |= __r & (FE_ALL_EXCEPT | (_ROUND_MASK << _ROUND_SHIFT));
+	__r &= ~(_ENABLE_MASK);
+	__msr_fpsr(__r);
+	return (0);
+}
+
+__fenv_static inline int
+fesetenv(const fenv_t *__envp)
+{
+
+	__msr_fpcr((*__envp) & _ENABLE_MASK);
+	__msr_fpsr((*__envp) & (FE_ALL_EXCEPT | (_ROUND_MASK << _ROUND_SHIFT)));
+	return (0);
+}
+
+__fenv_static inline int
+feupdateenv(const fenv_t *__envp)
+{
+	fexcept_t __r;
+
+	__mrs_fpsr(__r);
+	fesetenv(__envp);
+	feraiseexcept(__r & FE_ALL_EXCEPT);
+	return (0);
+}
+
+#if __BSD_VISIBLE
+
+/* We currently provide no external definitions of the functions below. */
+
+static inline int
+feenableexcept(int __mask)
+{
+	fenv_t __old_r, __new_r;
+
+	__mrs_fpcr(__old_r);
+	__new_r = __old_r | ((__mask & FE_ALL_EXCEPT) << _FPUSW_SHIFT);
+	__msr_fpcr(__new_r);
+	return ((__old_r >> _FPUSW_SHIFT) & FE_ALL_EXCEPT);
+}
+
+static inline int
+fedisableexcept(int __mask)
+{
+	fenv_t __old_r, __new_r;
+
+	__mrs_fpcr(__old_r);
+	__new_r = __old_r & ~((__mask & FE_ALL_EXCEPT) << _FPUSW_SHIFT);
+	__msr_fpcr(__new_r);
+	return ((__old_r >> _FPUSW_SHIFT) & FE_ALL_EXCEPT);
+}
+
+static inline int
+fegetexcept(void)
+{
+	fenv_t __r;
+
+	__mrs_fpcr(__r);
+	return ((__r & _ENABLE_MASK) >> _FPUSW_SHIFT);
+}
+
+#endif /* __BSD_VISIBLE */
+
+__END_DECLS
+
+#endif	/* !_FENV_H_ */
diff --git a/lib/msun/ld128/k_expl.h b/lib/msun/ld128/k_expl.h
index a5668fd..e0a48fc 100644
--- a/lib/msun/ld128/k_expl.h
+++ b/lib/msun/ld128/k_expl.h
@@ -322,7 +322,7 @@ __ldexp_cexpl(long double complex z, int expt)
 	scale2 = 1;
 	SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
 
-	return (cpackl(cos(y) * exp_x * scale1 * scale2,
+	return (CMPLXL(cos(y) * exp_x * scale1 * scale2,
 	    sinl(y) * exp_x * scale1 * scale2));
 }
 #endif /* _COMPLEX_H */
diff --git a/lib/msun/ld80/k_expl.h b/lib/msun/ld80/k_expl.h
index ebfb9a8..9b081fa 100644
--- a/lib/msun/ld80/k_expl.h
+++ b/lib/msun/ld80/k_expl.h
@@ -299,7 +299,7 @@ __ldexp_cexpl(long double complex z, int expt)
 	scale2 = 1;
 	SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
 
-	return (cpackl(cos(y) * exp_x * scale1 * scale2,
+	return (CMPLXL(cos(y) * exp_x * scale1 * scale2,
 	    sinl(y) * exp_x * scale1 * scale2));
 }
 #endif /* _COMPLEX_H */
diff --git a/lib/msun/man/cexp.3 b/lib/msun/man/cexp.3
index 97e36c1..776e6ce 100644
--- a/lib/msun/man/cexp.3
+++ b/lib/msun/man/cexp.3
@@ -103,7 +103,7 @@ is not finite, the sign of the result is indeterminate.
 .Sh SEE ALSO
 .Xr complex 3 ,
 .Xr exp 3 ,
-.Xr math 3 ,
+.Xr math 3
 .Sh STANDARDS
 The
 .Fn cexp
diff --git a/lib/msun/man/complex.3 b/lib/msun/man/complex.3
index 34eb03e..d286494 100644
--- a/lib/msun/man/complex.3
+++ b/lib/msun/man/complex.3
@@ -103,9 +103,9 @@ ctan	tangent
 ctanh	hyperbolic tangent
 .El
 .Sh SEE ALSO
-.Xr math 3 ,
 .Xr fenv 3 ,
 .Xr ieee 3 ,
+.Xr math 3 ,
 .Xr tgmath 3
 .Rs
 .%T "ISO/IEC 9899:TC3"
diff --git a/lib/msun/man/csqrt.3 b/lib/msun/man/csqrt.3
index 1a1dfa0..d8066ae 100644
--- a/lib/msun/man/csqrt.3
+++ b/lib/msun/man/csqrt.3
@@ -85,7 +85,7 @@ an \*(Na is generated, an invalid exception will be thrown.
 .Sh SEE ALSO
 .Xr cabs 3 ,
 .Xr fenv 3 ,
-.Xr math 3 ,
+.Xr math 3
 .Sh STANDARDS
 The
 .Fn csqrt ,
diff --git a/lib/msun/man/j0.3 b/lib/msun/man/j0.3
index 91849da..587b72e 100644
--- a/lib/msun/man/j0.3
+++ b/lib/msun/man/j0.3
@@ -28,7 +28,7 @@
 .\"     from: @(#)j0.3	6.7 (Berkeley) 4/19/91
 .\" $FreeBSD$
 .\"
-.Dd February 18, 2008
+.Dd March 10, 2015
 .Dt J0 3
 .Os
 .Sh NAME
@@ -77,24 +77,17 @@
 The functions
 .Fn j0 ,
 .Fn j0f ,
-.Fn j1
+.Fn j1 ,
 and
 .Fn j1f
-compute the
-.Em Bessel function of the first kind of the order
-0 and the
-.Em order
-1, respectively,
-for the
-real value
+compute the Bessel function of the first kind of orders
+0 and 1 for the real value
 .Fa x ;
 the functions
 .Fn jn
 and
 .Fn jnf
-compute the
-.Em Bessel function of the first kind of the integer
-.Em order
+compute the Bessel function of the first kind of the integer order
 .Fa n
 for the real value
 .Fa x .
@@ -105,13 +98,8 @@ The functions
 .Fn y1 ,
 and
 .Fn y1f
-compute the linearly independent
-.Em Bessel function of the second kind of the order
-0 and the
-.Em order
-1, respectively,
-for the
-positive
+compute the linearly independent Bessel function of the second kind
+of orders 0 and 1 for the positive
 .Em real
 value
 .Fa x ;
@@ -119,9 +107,7 @@ the functions
 .Fn yn
 and
 .Fn ynf
-compute the
-.Em Bessel function of the second kind for the integer
-.Em order
+compute the Bessel function of the second kind for the integer order
 .Fa n
 for the positive
 .Em real
@@ -141,11 +127,20 @@ and
 .Fn ynf .
 If
 .Fa x
-is negative, these routines will generate an invalid exception and
-return \*(Na.
+is negative, including -\*(If, these routines will generate an invalid
+exception and return \*(Na.
+If
+.Fa x
+is \*(Pm0, these routines
+will generate a divide-by-zero exception and return -\*(If.
 If
 .Fa x
-is 0 or a sufficiently small positive number, these routines
+is a sufficiently small positive number, then
+.Fn y1 ,
+.Fn y1f ,
+.Fn yn ,
+and
+.Fn ynf
 will generate an overflow exception and return -\*(If.
 .Sh SEE ALSO
 .Xr math 3
diff --git a/lib/msun/man/lgamma.3 b/lib/msun/man/lgamma.3
index d244722..c8a22a2 100644
--- a/lib/msun/man/lgamma.3
+++ b/lib/msun/man/lgamma.3
@@ -99,7 +99,7 @@ returns the sign of \(*G(x).
 and
 .Fn lgammal_r x signgamp
 provide the same functionality as
-.Fn lgamma x , 
+.Fn lgamma x ,
 .Fn lgammaf x ,
 and
 .Fn lgammal x ,
@@ -124,7 +124,6 @@ are deprecated aliases for
 and
 .Fn lgammaf_r ,
 respectively.
-
 .Sh IDIOSYNCRASIES
 Do not use the expression
 .Dq Li signgam\(**exp(lgamma(x))
diff --git a/lib/msun/man/nextafter.3 b/lib/msun/man/nextafter.3
index c8c4aa3..01c472d 100644
--- a/lib/msun/man/nextafter.3
+++ b/lib/msun/man/nextafter.3
@@ -78,11 +78,11 @@ routines conform to
 They implement the Nextafter function recommended by
 .St -ieee754 ,
 with the extension that
-.Fn nextafter +0.0, -0.0
+.Fn nextafter "+0.0" "-0.0"
 returns
 .Li -0.0 ,
 and
-.Fn nextafter -0.0, +0.0
+.Fn nextafter "-0.0" "+0.0"
 returns
 .Li +0.0 .
 .Sh HISTORY
diff --git a/lib/msun/man/sin.3 b/lib/msun/man/sin.3
index c7daf09..4a5352e 100644
--- a/lib/msun/man/sin.3
+++ b/lib/msun/man/sin.3
@@ -70,9 +70,9 @@ functions return the sine value.
 .Xr asin 3 ,
 .Xr atan 3 ,
 .Xr atan2 3 ,
-.Xr csin 3 ,
 .Xr cos 3 ,
 .Xr cosh 3 ,
+.Xr csin 3 ,
 .Xr math 3 ,
 .Xr sinh 3 ,
 .Xr tan 3 ,
diff --git a/lib/msun/src/catrig.c b/lib/msun/src/catrig.c
index 200977c..c0f5f55 100644
--- a/lib/msun/src/catrig.c
+++ b/lib/msun/src/catrig.c
@@ -286,19 +286,19 @@ casinh(double complex z)
 	if (isnan(x) || isnan(y)) {
 		/* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */
 		if (isinf(x))
-			return (cpack(x, y + y));
+			return (CMPLX(x, y + y));
 		/* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
 		if (isinf(y))
-			return (cpack(y, x + x));
+			return (CMPLX(y, x + x));
 		/* casinh(NaN + I*0) = NaN + I*0 */
 		if (y == 0)
-			return (cpack(x + x, y));
+			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 (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -307,7 +307,7 @@ casinh(double complex z)
 			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)));
+		return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y)));
 	}
 
 	/* Avoid spuriously raising inexact for z = 0. */
@@ -325,7 +325,7 @@ casinh(double complex z)
 		ry = asin(B);
 	else
 		ry = atan2(new_y, sqrt_A2my2);
-	return (cpack(copysign(rx, x), copysign(ry, y)));
+	return (CMPLX(copysign(rx, x), copysign(ry, y)));
 }
 
 /*
@@ -335,9 +335,9 @@ casinh(double complex z)
 double complex
 casin(double complex z)
 {
-	double complex w = casinh(cpack(cimag(z), creal(z)));
+	double complex w = casinh(CMPLX(cimag(z), creal(z)));
 
-	return (cpack(cimag(w), creal(w)));
+	return (CMPLX(cimag(w), creal(w)));
 }
 
 /*
@@ -370,19 +370,19 @@ cacos(double complex z)
 	if (isnan(x) || isnan(y)) {
 		/* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
 		if (isinf(x))
-			return (cpack(y + y, -INFINITY));
+			return (CMPLX(y + y, -INFINITY));
 		/* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */
 		if (isinf(y))
-			return (cpack(x + x, -y));
+			return (CMPLX(x + x, -y));
 		/* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */
 		if (x == 0)
-			return (cpack(pio2_hi + pio2_lo, y + y));
+			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 (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@@ -392,18 +392,18 @@ cacos(double complex z)
 		ry = creal(w) + m_ln2;
 		if (sy == 0)
 			ry = -ry;
-		return (cpack(rx, ry));
+		return (CMPLX(rx, ry));
 	}
 
 	/* Avoid spuriously raising inexact for z = 1. */
 	if (x == 1 && y == 0)
-		return (cpack(0, -y));
+		return (CMPLX(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));
+		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) {
@@ -419,7 +419,7 @@ cacos(double complex z)
 	}
 	if (sy == 0)
 		ry = -ry;
-	return (cpack(rx, ry));
+	return (CMPLX(rx, ry));
 }
 
 /*
@@ -437,15 +437,15 @@ cacosh(double complex z)
 	ry = cimag(w);
 	/* cacosh(NaN + I*NaN) = NaN + I*NaN */
 	if (isnan(rx) && isnan(ry))
-		return (cpack(ry, rx));
+		return (CMPLX(ry, rx));
 	/* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */
 	/* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */
 	if (isnan(rx))
-		return (cpack(fabs(ry), rx));
+		return (CMPLX(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))));
+		return (CMPLX(ry, ry));
+	return (CMPLX(fabs(ry), copysign(rx, cimag(z))));
 }
 
 /*
@@ -475,16 +475,16 @@ clog_for_large_values(double complex z)
 	 * 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)));
+		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 (cpack(log(hypot(x, y)), atan2(y, x)));
+		return (CMPLX(log(hypot(x, y)), atan2(y, x)));
 
-	return (cpack(log(ax * ax + ay * ay) / 2, atan2(y, x)));
+	return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x)));
 }
 
 /*
@@ -575,30 +575,30 @@ catanh(double complex z)
 
 	/* This helps handle many cases. */
 	if (y == 0 && ax <= 1)
-		return (cpack(atanh(x), y));
+		return (CMPLX(atanh(x), y));
 
 	/* To ensure the same accuracy as atan(), and to filter out z = 0. */
 	if (x == 0)
-		return (cpack(x, atan(y)));
+		return (CMPLX(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));
+			return (CMPLX(copysign(0, x), y + y));
 		/* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
 		if (isinf(y))
-			return (cpack(copysign(0, x),
+			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 (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
-		return (cpack(real_part_reciprocal(x, y),
+		return (CMPLX(real_part_reciprocal(x, y),
 		    copysign(pio2_hi + pio2_lo, y)));
 
 	if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
@@ -623,7 +623,7 @@ catanh(double complex z)
 	else
 		ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
 
-	return (cpack(copysign(rx, x), copysign(ry, y)));
+	return (CMPLX(copysign(rx, x), copysign(ry, y)));
 }
 
 /*
@@ -633,7 +633,7 @@ catanh(double complex z)
 double complex
 catan(double complex z)
 {
-	double complex w = catanh(cpack(cimag(z), creal(z)));
+	double complex w = catanh(CMPLX(cimag(z), creal(z)));
 
-	return (cpack(cimag(w), creal(w)));
+	return (CMPLX(cimag(w), creal(w)));
 }
diff --git a/lib/msun/src/catrigf.c b/lib/msun/src/catrigf.c
index 08ebef7..ca84ce2 100644
--- a/lib/msun/src/catrigf.c
+++ b/lib/msun/src/catrigf.c
@@ -156,12 +156,12 @@ casinhf(float complex z)
 
 	if (isnan(x) || isnan(y)) {
 		if (isinf(x))
-			return (cpackf(x, y + y));
+			return (CMPLXF(x, y + y));
 		if (isinf(y))
-			return (cpackf(y, x + x));
+			return (CMPLXF(y, x + x));
 		if (y == 0)
-			return (cpackf(x + x, y));
-		return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (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) {
@@ -169,7 +169,7 @@ casinhf(float complex z)
 			w = clog_for_large_values(z) + m_ln2;
 		else
 			w = clog_for_large_values(-z) + m_ln2;
-		return (cpackf(copysignf(crealf(w), x),
+		return (CMPLXF(copysignf(crealf(w), x),
 		    copysignf(cimagf(w), y)));
 	}
 
@@ -186,15 +186,15 @@ casinhf(float complex z)
 		ry = asinf(B);
 	else
 		ry = atan2f(new_y, sqrt_A2my2);
-	return (cpackf(copysignf(rx, x), copysignf(ry, y)));
+	return (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
 }
 
 float complex
 casinf(float complex z)
 {
-	float complex w = casinhf(cpackf(cimagf(z), crealf(z)));
+	float complex w = casinhf(CMPLXF(cimagf(z), crealf(z)));
 
-	return (cpackf(cimagf(w), crealf(w)));
+	return (CMPLXF(cimagf(w), crealf(w)));
 }
 
 float complex
@@ -214,12 +214,12 @@ cacosf(float complex z)
 
 	if (isnan(x) || isnan(y)) {
 		if (isinf(x))
-			return (cpackf(y + y, -INFINITY));
+			return (CMPLXF(y + y, -INFINITY));
 		if (isinf(y))
-			return (cpackf(x + x, -y));
+			return (CMPLXF(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)));
+			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) {
@@ -228,16 +228,16 @@ cacosf(float complex z)
 		ry = crealf(w) + m_ln2;
 		if (sy == 0)
 			ry = -ry;
-		return (cpackf(rx, ry));
+		return (CMPLXF(rx, ry));
 	}
 
 	if (x == 1 && y == 0)
-		return (cpackf(0, -y));
+		return (CMPLXF(0, -y));
 
 	raise_inexact();
 
 	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
-		return (cpackf(pio2_hi - (x - pio2_lo), -y));
+		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) {
@@ -253,7 +253,7 @@ cacosf(float complex z)
 	}
 	if (sy == 0)
 		ry = -ry;
-	return (cpackf(rx, ry));
+	return (CMPLXF(rx, ry));
 }
 
 float complex
@@ -266,12 +266,12 @@ cacoshf(float complex z)
 	rx = crealf(w);
 	ry = cimagf(w);
 	if (isnan(rx) && isnan(ry))
-		return (cpackf(ry, rx));
+		return (CMPLXF(ry, rx));
 	if (isnan(rx))
-		return (cpackf(fabsf(ry), rx));
+		return (CMPLXF(fabsf(ry), rx));
 	if (isnan(ry))
-		return (cpackf(ry, ry));
-	return (cpackf(fabsf(ry), copysignf(rx, cimagf(z))));
+		return (CMPLXF(ry, ry));
+	return (CMPLXF(fabsf(ry), copysignf(rx, cimagf(z))));
 }
 
 static float complex
@@ -291,13 +291,13 @@ clog_for_large_values(float complex z)
 	}
 
 	if (ax > FLT_MAX / 2)
-		return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1,
+		return (CMPLXF(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 (CMPLXF(logf(hypotf(x, y)), atan2f(y, x)));
 
-	return (cpackf(logf(ax * ax + ay * ay) / 2, atan2f(y, x)));
+	return (CMPLXF(logf(ax * ax + ay * ay) / 2, atan2f(y, x)));
 }
 
 static inline float
@@ -346,22 +346,22 @@ catanhf(float complex z)
 	ay = fabsf(y);
 
 	if (y == 0 && ax <= 1)
-		return (cpackf(atanhf(x), y));
+		return (CMPLXF(atanhf(x), y));
 
 	if (x == 0)
-		return (cpackf(x, atanf(y)));
+		return (CMPLXF(x, atanf(y)));
 
 	if (isnan(x) || isnan(y)) {
 		if (isinf(x))
-			return (cpackf(copysignf(0, x), y + y));
+			return (CMPLXF(copysignf(0, x), y + y));
 		if (isinf(y))
-			return (cpackf(copysignf(0, x),
+			return (CMPLXF(copysignf(0, x),
 			    copysignf(pio2_hi + pio2_lo, y)));
-		return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
+		return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
 	}
 
 	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
-		return (cpackf(real_part_reciprocal(x, y),
+		return (CMPLXF(real_part_reciprocal(x, y),
 		    copysignf(pio2_hi + pio2_lo, y)));
 
 	if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
@@ -381,13 +381,13 @@ catanhf(float complex z)
 	else
 		ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
 
-	return (cpackf(copysignf(rx, x), copysignf(ry, y)));
+	return (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
 }
 
 float complex
 catanf(float complex z)
 {
-	float complex w = catanhf(cpackf(cimagf(z), crealf(z)));
+	float complex w = catanhf(CMPLXF(cimagf(z), crealf(z)));
 
-	return (cpackf(cimagf(w), crealf(w)));
+	return (CMPLXF(cimagf(w), crealf(w)));
 }
diff --git a/lib/msun/src/e_j0.c b/lib/msun/src/e_j0.c
index 8320f25..a253ed1 100644
--- a/lib/msun/src/e_j0.c
+++ b/lib/msun/src/e_j0.c
@@ -62,7 +62,9 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
-static double pzero(double), qzero(double);
+static __inline double pzero(double), qzero(double);
+
+static const volatile double vone = 1, vzero = 0;
 
 static const double
 huge 	= 1e300,
@@ -115,7 +117,7 @@ __ieee754_j0(double x)
 	if(ix<0x3f200000) {	/* |x| < 2**-13 */
 	    if(huge+x>one) {	/* raise inexact if x != 0 */
 	        if(ix<0x3e400000) return one;	/* |x|<2**-27 */
-	        else 	      return one - 0.25*x*x;
+	        else 	      return one - x*x/4;
 	    }
 	}
 	z = x*x;
@@ -150,10 +152,16 @@ __ieee754_y0(double x)
 
 	EXTRACT_WORDS(hx,lx,x);
         ix = 0x7fffffff&hx;
-    /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
-	if(ix>=0x7ff00000) return  one/(x+x*x); 
-        if((ix|lx)==0) return -one/zero;
-        if(hx<0) return zero/zero;
+	/*
+	 * y0(NaN) = NaN.
+	 * y0(Inf) = 0.
+	 * y0(-Inf) = NaN and raise invalid exception.
+	 */
+	if(ix>=0x7ff00000) return vone/(x+x*x); 
+	/* y0(+-0) = -inf and raise divide-by-zero exception. */
+	if((ix|lx)==0) return -one/vzero;
+	/* y0(x<0) = NaN and raise invalid exception. */
+	if(hx<0) return vzero/vzero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
         /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
          * where x0 = x-pi/4
@@ -268,7 +276,8 @@ static const double pS2[5] = {
   1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */
 };
 
-	static double pzero(double x)
+static __inline double
+pzero(double x)
 {
 	const double *p,*q;
 	double z,r,s;
@@ -278,7 +287,7 @@ static const double pS2[5] = {
 	if(ix>=0x40200000)     {p = pR8; q= pS8;}
 	else if(ix>=0x40122E8B){p = pR5; q= pS5;}
 	else if(ix>=0x4006DB6D){p = pR3; q= pS3;}
-	else if(ix>=0x40000000){p = pR2; q= pS2;}
+	else                   {p = pR2; q= pS2;}	/* ix>=0x40000000 */
 	z = one/(x*x);
 	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@@ -363,7 +372,8 @@ static const double qS2[6] = {
  -5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
 };
 
-	static double qzero(double x)
+static __inline double
+qzero(double x)
 {
 	const double *p,*q;
 	double s,r,z;
@@ -373,7 +383,7 @@ static const double qS2[6] = {
 	if(ix>=0x40200000)     {p = qR8; q= qS8;}
 	else if(ix>=0x40122E8B){p = qR5; q= qS5;}
 	else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
-	else if(ix>=0x40000000){p = qR2; q= qS2;}
+	else                   {p = qR2; q= qS2;}	/* ix>=0x40000000 */
 	z = one/(x*x);
 	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
diff --git a/lib/msun/src/e_j0f.c b/lib/msun/src/e_j0f.c
index c45faf3..3e2f7e8 100644
--- a/lib/msun/src/e_j0f.c
+++ b/lib/msun/src/e_j0f.c
@@ -16,10 +16,16 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
+/*
+ * See e_j0.c for complete comments.
+ */
+
 #include "math.h"
 #include "math_private.h"
 
-static float pzerof(float), qzerof(float);
+static __inline float pzerof(float), qzerof(float);
+
+static const volatile float vone = 1,  vzero = 0;
 
 static const float
 huge 	= 1e30,
@@ -62,17 +68,17 @@ __ieee754_j0f(float x)
 	 * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 	 * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 	 */
-		if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);
+		if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(x); /* |x|>2**49 */
 		else {
 		    u = pzerof(x); v = qzerof(x);
 		    z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);
 		}
 		return z;
 	}
-	if(ix<0x39000000) {	/* |x| < 2**-13 */
+	if(ix<0x3b000000) {	/* |x| < 2**-9 */
 	    if(huge+x>one) {	/* raise inexact if x != 0 */
-	        if(ix<0x32000000) return one;	/* |x|<2**-27 */
-	        else 	      return one - (float)0.25*x*x;
+	        if(ix<0x39800000) return one;	/* |x|<2**-12 */
+	        else 	      return one - x*x/4;
 	    }
 	}
 	z = x*x;
@@ -107,10 +113,9 @@ __ieee754_y0f(float x)
 
 	GET_FLOAT_WORD(hx,x);
         ix = 0x7fffffff&hx;
-    /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
-	if(ix>=0x7f800000) return  one/(x+x*x);
-        if(ix==0) return -one/zero;
-        if(hx<0) return zero/zero;
+	if(ix>=0x7f800000) return  vone/(x+x*x);
+	if(ix==0) return -one/vzero;
+	if(hx<0) return vzero/vzero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
         /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
          * where x0 = x-pi/4
@@ -136,14 +141,14 @@ __ieee754_y0f(float x)
                     if ((s*c)<zero) cc = z/ss;
                     else            ss = z/cc;
                 }
-                if(ix>0x80000000) z = (invsqrtpi*ss)/sqrtf(x);
+                if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
                 else {
                     u = pzerof(x); v = qzerof(x);
                     z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
                 }
                 return z;
 	}
-	if(ix<=0x32000000) {	/* x < 2**-27 */
+	if(ix<=0x39000000) {	/* x < 2**-13 */
 	    return(u00 + tpi*__ieee754_logf(x));
 	}
 	z = x*x;
@@ -224,7 +229,8 @@ static const float pS2[5] = {
   1.4657617569e+01, /* 0x416a859a */
 };
 
-	static float pzerof(float x)
+static __inline float
+pzerof(float x)
 {
 	const float *p,*q;
 	float z,r,s;
@@ -232,9 +238,9 @@ static const float pS2[5] = {
 	GET_FLOAT_WORD(ix,x);
 	ix &= 0x7fffffff;
 	if(ix>=0x41000000)     {p = pR8; q= pS8;}
-	else if(ix>=0x40f71c58){p = pR5; q= pS5;}
-	else if(ix>=0x4036db68){p = pR3; q= pS3;}
-	else if(ix>=0x40000000){p = pR2; q= pS2;}
+	else if(ix>=0x409173eb){p = pR5; q= pS5;}
+	else if(ix>=0x4036d917){p = pR3; q= pS3;}
+	else                   {p = pR2; q= pS2;}	/* ix>=0x40000000 */
 	z = one/(x*x);
 	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@@ -319,7 +325,8 @@ static const float qS2[6] = {
  -5.3109550476e+00, /* 0xc0a9f358 */
 };
 
-	static float qzerof(float x)
+static __inline float
+qzerof(float x)
 {
 	const float *p,*q;
 	float s,r,z;
@@ -327,9 +334,9 @@ static const float qS2[6] = {
 	GET_FLOAT_WORD(ix,x);
 	ix &= 0x7fffffff;
 	if(ix>=0x41000000)     {p = qR8; q= qS8;}
-	else if(ix>=0x40f71c58){p = qR5; q= qS5;}
-	else if(ix>=0x4036db68){p = qR3; q= qS3;}
-	else if(ix>=0x40000000){p = qR2; q= qS2;}
+	else if(ix>=0x409173eb){p = qR5; q= qS5;}
+	else if(ix>=0x4036d917){p = qR3; q= qS3;}
+	else                   {p = qR2; q= qS2;}	/* ix>=0x40000000 */
 	z = one/(x*x);
 	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
diff --git a/lib/msun/src/e_j1.c b/lib/msun/src/e_j1.c
index 63800ad..74a7c69 100644
--- a/lib/msun/src/e_j1.c
+++ b/lib/msun/src/e_j1.c
@@ -62,7 +62,9 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
-static double pone(double), qone(double);
+static __inline double pone(double), qone(double);
+
+static const volatile double vone = 1, vzero = 0;
 
 static const double
 huge    = 1e300,
@@ -147,10 +149,16 @@ __ieee754_y1(double x)
 
 	EXTRACT_WORDS(hx,lx,x);
         ix = 0x7fffffff&hx;
-    /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
-	if(ix>=0x7ff00000) return  one/(x+x*x); 
-        if((ix|lx)==0) return -one/zero;
-        if(hx<0) return zero/zero;
+	/*
+	 * y1(NaN) = NaN.
+	 * y1(Inf) = 0.
+	 * y1(-Inf) = NaN and raise invalid exception.
+	 */
+	if(ix>=0x7ff00000) return  vone/(x+x*x); 
+	/* y1(+-0) = -inf and raise divide-by-zero exception. */
+        if((ix|lx)==0) return -one/vzero;
+	/* y1(x<0) = NaN and raise invalid exception. */
+        if(hx<0) return vzero/vzero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
                 s = sin(x);
                 c = cos(x);
@@ -262,7 +270,8 @@ static const double ps2[5] = {
   8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
 };
 
-	static double pone(double x)
+static __inline double
+pone(double x)
 {
 	const double *p,*q;
 	double z,r,s;
@@ -272,7 +281,7 @@ static const double ps2[5] = {
         if(ix>=0x40200000)     {p = pr8; q= ps8;}
         else if(ix>=0x40122E8B){p = pr5; q= ps5;}
         else if(ix>=0x4006DB6D){p = pr3; q= ps3;}
-        else if(ix>=0x40000000){p = pr2; q= ps2;}
+	else                   {p = pr2; q= ps2;}	/* ix>=0x40000000 */
         z = one/(x*x);
         r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
         s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@@ -358,7 +367,8 @@ static const double qs2[6] = {
  -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
 };
 
-	static double qone(double x)
+static __inline double
+qone(double x)
 {
 	const double *p,*q;
 	double  s,r,z;
@@ -368,7 +378,7 @@ static const double qs2[6] = {
 	if(ix>=0x40200000)     {p = qr8; q= qs8;}
 	else if(ix>=0x40122E8B){p = qr5; q= qs5;}
 	else if(ix>=0x4006DB6D){p = qr3; q= qs3;}
-	else if(ix>=0x40000000){p = qr2; q= qs2;}
+	else                   {p = qr2; q= qs2;}	/* ix>=0x40000000 */
 	z = one/(x*x);
 	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
diff --git a/lib/msun/src/e_j1f.c b/lib/msun/src/e_j1f.c
index 88e2d83..ec7f381 100644
--- a/lib/msun/src/e_j1f.c
+++ b/lib/msun/src/e_j1f.c
@@ -16,10 +16,16 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
+/*
+ * See e_j1.c for complete comments.
+ */
+
 #include "math.h"
 #include "math_private.h"
 
-static float ponef(float), qonef(float);
+static __inline float ponef(float), qonef(float);
+
+static const volatile float vone = 1, vzero = 0;
 
 static const float
 huge    = 1e30,
@@ -63,7 +69,7 @@ __ieee754_j1f(float x)
 	 * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
 	 * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
 	 */
-		if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(y);
+		if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(y); /* |x|>2**49 */
 		else {
 		    u = ponef(y); v = qonef(y);
 		    z = invsqrtpi*(u*cc-v*ss)/sqrtf(y);
@@ -71,7 +77,7 @@ __ieee754_j1f(float x)
 		if(hx<0) return -z;
 		else  	 return  z;
 	}
-	if(ix<0x32000000) {	/* |x|<2**-27 */
+	if(ix<0x39000000) {	/* |x|<2**-13 */
 	    if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */
 	}
 	z = x*x;
@@ -104,10 +110,9 @@ __ieee754_y1f(float x)
 
 	GET_FLOAT_WORD(hx,x);
         ix = 0x7fffffff&hx;
-    /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
-	if(ix>=0x7f800000) return  one/(x+x*x);
-        if(ix==0) return -one/zero;
-        if(hx<0) return zero/zero;
+	if(ix>=0x7f800000) return  vone/(x+x*x);
+	if(ix==0) return -one/vzero;
+	if(hx<0) return vzero/vzero;
         if(ix >= 0x40000000) {  /* |x| >= 2.0 */
                 s = sinf(x);
                 c = cosf(x);
@@ -129,14 +134,14 @@ __ieee754_y1f(float x)
          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
          * to compute the worse one.
          */
-                if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x);
+                if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
                 else {
                     u = ponef(x); v = qonef(x);
                     z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
                 }
                 return z;
         }
-        if(ix<=0x24800000) {    /* x < 2**-54 */
+        if(ix<=0x33000000) {    /* x < 2**-25 */
             return(-tpi/x);
         }
         z = x*x;
@@ -219,7 +224,8 @@ static const float ps2[5] = {
   8.3646392822e+00, /* 0x4105d590 */
 };
 
-	static float ponef(float x)
+static __inline float
+ponef(float x)
 {
 	const float *p,*q;
 	float z,r,s;
@@ -227,9 +233,9 @@ static const float ps2[5] = {
 	GET_FLOAT_WORD(ix,x);
 	ix &= 0x7fffffff;
         if(ix>=0x41000000)     {p = pr8; q= ps8;}
-        else if(ix>=0x40f71c58){p = pr5; q= ps5;}
-        else if(ix>=0x4036db68){p = pr3; q= ps3;}
-        else if(ix>=0x40000000){p = pr2; q= ps2;}
+        else if(ix>=0x409173eb){p = pr5; q= ps5;}
+        else if(ix>=0x4036d917){p = pr3; q= ps3;}
+	else                   {p = pr2; q= ps2;}	/* ix>=0x40000000 */
         z = one/(x*x);
         r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
         s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@@ -315,17 +321,18 @@ static const float qs2[6] = {
  -4.9594988823e+00, /* 0xc09eb437 */
 };
 
-	static float qonef(float x)
+static __inline float
+qonef(float x)
 {
 	const float *p,*q;
 	float  s,r,z;
 	int32_t ix;
 	GET_FLOAT_WORD(ix,x);
 	ix &= 0x7fffffff;
-	if(ix>=0x40200000)     {p = qr8; q= qs8;}
-	else if(ix>=0x40f71c58){p = qr5; q= qs5;}
-	else if(ix>=0x4036db68){p = qr3; q= qs3;}
-	else if(ix>=0x40000000){p = qr2; q= qs2;}
+	if(ix>=0x41000000)     {p = qr8; q= qs8;}
+	else if(ix>=0x409173eb){p = qr5; q= qs5;}
+	else if(ix>=0x4036d917){p = qr3; q= qs3;}
+	else                   {p = qr2; q= qs2;}	/* ix>=0x40000000 */
 	z = one/(x*x);
 	r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
 	s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
diff --git a/lib/msun/src/e_jn.c b/lib/msun/src/e_jn.c
index 8b0bc62..eefd4ff 100644
--- a/lib/msun/src/e_jn.c
+++ b/lib/msun/src/e_jn.c
@@ -43,6 +43,8 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+static const volatile double vone = 1, vzero = 0;
+
 static const double
 invsqrtpi=  5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
 two   =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
@@ -220,10 +222,12 @@ __ieee754_yn(int n, double x)
 
 	EXTRACT_WORDS(hx,lx,x);
 	ix = 0x7fffffff&hx;
-    /* if Y(n,NaN) is NaN */
+	/* yn(n,NaN) = NaN */
 	if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
-	if((ix|lx)==0) return -one/zero;
-	if(hx<0) return zero/zero;
+	/* yn(n,+-0) = -inf and raise divide-by-zero exception. */
+	if((ix|lx)==0) return -one/vzero;
+	/* yn(n,x<0) = NaN and raise invalid exception. */
+	if(hx<0) return vzero/vzero;
 	sign = 1;
 	if(n<0){
 		n = -n;
diff --git a/lib/msun/src/e_jnf.c b/lib/msun/src/e_jnf.c
index f564aec..965feeb 100644
--- a/lib/msun/src/e_jnf.c
+++ b/lib/msun/src/e_jnf.c
@@ -16,9 +16,15 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
+/*
+ * See e_jn.c for complete comments.
+ */
+
 #include "math.h"
 #include "math_private.h"
 
+static const volatile float vone = 1, vzero = 0;
+
 static const float
 two   =  2.0000000000e+00, /* 0x40000000 */
 one   =  1.0000000000e+00; /* 0x3F800000 */
@@ -172,10 +178,9 @@ __ieee754_ynf(int n, float x)
 
 	GET_FLOAT_WORD(hx,x);
 	ix = 0x7fffffff&hx;
-    /* if Y(n,NaN) is NaN */
 	if(ix>0x7f800000) return x+x;
-	if(ix==0) return -one/zero;
-	if(hx<0) return zero/zero;
+	if(ix==0) return -one/vzero;
+	if(hx<0) return vzero/vzero;
 	sign = 1;
 	if(n<0){
 		n = -n;
diff --git a/lib/msun/src/k_exp.c b/lib/msun/src/k_exp.c
index f592f69..ecef54c 100644
--- a/lib/msun/src/k_exp.c
+++ b/lib/msun/src/k_exp.c
@@ -103,6 +103,6 @@ __ldexp_cexp(double complex z, int expt)
 	half_expt = expt - half_expt;
 	INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
 
-	return (cpack(cos(y) * exp_x * scale1 * scale2,
+	return (CMPLX(cos(y) * exp_x * scale1 * scale2,
 	    sin(y) * exp_x * scale1 * scale2));
 }
diff --git a/lib/msun/src/k_expf.c b/lib/msun/src/k_expf.c
index 548a008..f8c254d 100644
--- a/lib/msun/src/k_expf.c
+++ b/lib/msun/src/k_expf.c
@@ -82,6 +82,6 @@ __ldexp_cexpf(float complex z, int expt)
 	half_expt = expt - half_expt;
 	SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23);
 
-	return (cpackf(cosf(y) * exp_x * scale1 * scale2,
+	return (CMPLXF(cosf(y) * exp_x * scale1 * scale2,
 	    sinf(y) * exp_x * scale1 * scale2));
 }
diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h
index 8af2c65..afaf201 100644
--- a/lib/msun/src/math_private.h
+++ b/lib/msun/src/math_private.h
@@ -454,9 +454,15 @@ typedef union {
  * (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
-cpackf(float x, float y)
+CMPLXF(float x, float y)
 {
 	float_complex z;
 
@@ -464,9 +470,11 @@ cpackf(float x, float y)
 	IMAGPART(z) = y;
 	return (z.f);
 }
+#endif
 
+#ifndef CMPLX
 static __inline double complex
-cpack(double x, double y)
+CMPLX(double x, double y)
 {
 	double_complex z;
 
@@ -474,9 +482,11 @@ cpack(double x, double y)
 	IMAGPART(z) = y;
 	return (z.f);
 }
+#endif
 
+#ifndef CMPLXL
 static __inline long double complex
-cpackl(long double x, long double y)
+CMPLXL(long double x, long double y)
 {
 	long_double_complex z;
 
@@ -484,6 +494,8 @@ cpackl(long double x, long double y)
 	IMAGPART(z) = y;
 	return (z.f);
 }
+#endif
+
 #endif /* _COMPLEX_H */
  
 #ifdef __GNUCLIKE_ASM
diff --git a/lib/msun/src/s_ccosh.c b/lib/msun/src/s_ccosh.c
index 9ea962b..fbe15fc 100644
--- a/lib/msun/src/s_ccosh.c
+++ b/lib/msun/src/s_ccosh.c
@@ -62,23 +62,23 @@ ccosh(double complex z)
 	/* Handle the nearly-non-exceptional cases where x and y are finite. */
 	if (ix < 0x7ff00000 && iy < 0x7ff00000) {
 		if ((iy | ly) == 0)
-			return (cpack(cosh(x), x * y));
+			return (CMPLX(cosh(x), x * y));
 		if (ix < 0x40360000)	/* small x: normal case */
-			return (cpack(cosh(x) * cos(y), sinh(x) * sin(y)));
+			return (CMPLX(cosh(x) * cos(y), sinh(x) * sin(y)));
 
 		/* |x| >= 22, so cosh(x) ~= exp(|x|) */
 		if (ix < 0x40862e42) {
 			/* x < 710: exp(|x|) won't overflow */
 			h = exp(fabs(x)) * 0.5;
-			return (cpack(h * cos(y), copysign(h, x) * sin(y)));
+			return (CMPLX(h * cos(y), copysign(h, x) * sin(y)));
 		} else if (ix < 0x4096bbaa) {
 			/* x < 1455: scale to avoid overflow */
-			z = __ldexp_cexp(cpack(fabs(x), y), -1);
-			return (cpack(creal(z), cimag(z) * copysign(1, x)));
+			z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
+			return (CMPLX(creal(z), cimag(z) * copysign(1, x)));
 		} else {
 			/* x >= 1455: the result always overflows */
 			h = huge * x;
-			return (cpack(h * h * cos(y), h * sin(y)));
+			return (CMPLX(h * h * cos(y), h * sin(y)));
 		}
 	}
 
@@ -92,7 +92,7 @@ ccosh(double complex z)
 	 * the same as d(NaN).
 	 */
 	if ((ix | lx) == 0 && iy >= 0x7ff00000)
-		return (cpack(y - y, copysign(0, x * (y - y))));
+		return (CMPLX(y - y, copysign(0, x * (y - y))));
 
 	/*
 	 * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
@@ -102,8 +102,8 @@ ccosh(double complex z)
 	 */
 	if ((iy | ly) == 0 && ix >= 0x7ff00000) {
 		if (((hx & 0xfffff) | lx) == 0)
-			return (cpack(x * x, copysign(0, x) * y));
-		return (cpack(x * x, copysign(0, (x + x) * y)));
+			return (CMPLX(x * x, copysign(0, x) * y));
+		return (CMPLX(x * x, copysign(0, (x + x) * y)));
 	}
 
 	/*
@@ -115,7 +115,7 @@ ccosh(double complex z)
 	 * nonzero x.  Choice = don't raise (except for signaling NaNs).
 	 */
 	if (ix < 0x7ff00000 && iy >= 0x7ff00000)
-		return (cpack(y - y, x * (y - y)));
+		return (CMPLX(y - y, x * (y - y)));
 
 	/*
 	 * cosh(+-Inf + I NaN)  = +Inf + I d(NaN).
@@ -128,8 +128,8 @@ ccosh(double complex z)
 	 */
 	if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
 		if (iy >= 0x7ff00000)
-			return (cpack(x * x, x * (y - y)));
-		return (cpack((x * x) * cos(y), x * sin(y)));
+			return (CMPLX(x * x, x * (y - y)));
+		return (CMPLX((x * x) * cos(y), x * sin(y)));
 	}
 
 	/*
@@ -143,7 +143,7 @@ ccosh(double complex z)
 	 * Optionally raises the invalid floating-point exception for finite
 	 * nonzero y.  Choice = don't raise (except for signaling NaNs).
 	 */
-	return (cpack((x * x) * (y - y), (x + x) * (y - y)));
+	return (CMPLX((x * x) * (y - y), (x + x) * (y - y)));
 }
 
 double complex
@@ -151,5 +151,5 @@ ccos(double complex z)
 {
 
 	/* ccos(z) = ccosh(I * z) */
-	return (ccosh(cpack(-cimag(z), creal(z))));
+	return (ccosh(CMPLX(-cimag(z), creal(z))));
 }
diff --git a/lib/msun/src/s_ccoshf.c b/lib/msun/src/s_ccoshf.c
index 1de9ad4..fe8cf89 100644
--- a/lib/msun/src/s_ccoshf.c
+++ b/lib/msun/src/s_ccoshf.c
@@ -55,50 +55,50 @@ ccoshf(float complex z)
 
 	if (ix < 0x7f800000 && iy < 0x7f800000) {
 		if (iy == 0)
-			return (cpackf(coshf(x), x * y));
+			return (CMPLXF(coshf(x), x * y));
 		if (ix < 0x41100000)	/* small x: normal case */
-			return (cpackf(coshf(x) * cosf(y), sinhf(x) * sinf(y)));
+			return (CMPLXF(coshf(x) * cosf(y), sinhf(x) * sinf(y)));
 
 		/* |x| >= 9, so cosh(x) ~= exp(|x|) */
 		if (ix < 0x42b17218) {
 			/* x < 88.7: expf(|x|) won't overflow */
 			h = expf(fabsf(x)) * 0.5f;
-			return (cpackf(h * cosf(y), copysignf(h, x) * sinf(y)));
+			return (CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y)));
 		} else if (ix < 0x4340b1e7) {
 			/* x < 192.7: scale to avoid overflow */
-			z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
-			return (cpackf(crealf(z), cimagf(z) * copysignf(1, x)));
+			z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
+			return (CMPLXF(crealf(z), cimagf(z) * copysignf(1, x)));
 		} else {
 			/* x >= 192.7: the result always overflows */
 			h = huge * x;
-			return (cpackf(h * h * cosf(y), h * sinf(y)));
+			return (CMPLXF(h * h * cosf(y), h * sinf(y)));
 		}
 	}
 
 	if (ix == 0 && iy >= 0x7f800000)
-		return (cpackf(y - y, copysignf(0, x * (y - y))));
+		return (CMPLXF(y - y, copysignf(0, x * (y - y))));
 
 	if (iy == 0 && ix >= 0x7f800000) {
 		if ((hx & 0x7fffff) == 0)
-			return (cpackf(x * x, copysignf(0, x) * y));
-		return (cpackf(x * x, copysignf(0, (x + x) * y)));
+			return (CMPLXF(x * x, copysignf(0, x) * y));
+		return (CMPLXF(x * x, copysignf(0, (x + x) * y)));
 	}
 
 	if (ix < 0x7f800000 && iy >= 0x7f800000)
-		return (cpackf(y - y, x * (y - y)));
+		return (CMPLXF(y - y, x * (y - y)));
 
 	if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
 		if (iy >= 0x7f800000)
-			return (cpackf(x * x, x * (y - y)));
-		return (cpackf((x * x) * cosf(y), x * sinf(y)));
+			return (CMPLXF(x * x, x * (y - y)));
+		return (CMPLXF((x * x) * cosf(y), x * sinf(y)));
 	}
 
-	return (cpackf((x * x) * (y - y), (x + x) * (y - y)));
+	return (CMPLXF((x * x) * (y - y), (x + x) * (y - y)));
 }
 
 float complex
 ccosf(float complex z)
 {
 
-	return (ccoshf(cpackf(-cimagf(z), crealf(z))));
+	return (ccoshf(CMPLXF(-cimagf(z), crealf(z))));
 }
diff --git a/lib/msun/src/s_cexp.c b/lib/msun/src/s_cexp.c
index abe178f..9e11d51 100644
--- a/lib/msun/src/s_cexp.c
+++ b/lib/msun/src/s_cexp.c
@@ -50,22 +50,22 @@ cexp(double complex z)
 
 	/* cexp(x + I 0) = exp(x) + I 0 */
 	if ((hy | ly) == 0)
-		return (cpack(exp(x), y));
+		return (CMPLX(exp(x), y));
 	EXTRACT_WORDS(hx, lx, x);
 	/* cexp(0 + I y) = cos(y) + I sin(y) */
 	if (((hx & 0x7fffffff) | lx) == 0)
-		return (cpack(cos(y), sin(y)));
+		return (CMPLX(cos(y), sin(y)));
 
 	if (hy >= 0x7ff00000) {
 		if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
 			/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
-			return (cpack(y - y, y - y));
+			return (CMPLX(y - y, y - y));
 		} else if (hx & 0x80000000) {
 			/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
-			return (cpack(0.0, 0.0));
+			return (CMPLX(0.0, 0.0));
 		} else {
 			/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
-			return (cpack(x, y - y));
+			return (CMPLX(x, y - y));
 		}
 	}
 
@@ -84,6 +84,6 @@ cexp(double complex z)
 		 *  -  x = NaN (spurious inexact exception from y)
 		 */
 		exp_x = exp(x);
-		return (cpack(exp_x * cos(y), exp_x * sin(y)));
+		return (CMPLX(exp_x * cos(y), exp_x * sin(y)));
 	}
 }
diff --git a/lib/msun/src/s_cexpf.c b/lib/msun/src/s_cexpf.c
index 0e30d08..0138cd1 100644
--- a/lib/msun/src/s_cexpf.c
+++ b/lib/msun/src/s_cexpf.c
@@ -50,22 +50,22 @@ cexpf(float complex z)
 
 	/* cexp(x + I 0) = exp(x) + I 0 */
 	if (hy == 0)
-		return (cpackf(expf(x), y));
+		return (CMPLXF(expf(x), y));
 	GET_FLOAT_WORD(hx, x);
 	/* cexp(0 + I y) = cos(y) + I sin(y) */
 	if ((hx & 0x7fffffff) == 0)
-		return (cpackf(cosf(y), sinf(y)));
+		return (CMPLXF(cosf(y), sinf(y)));
 
 	if (hy >= 0x7f800000) {
 		if ((hx & 0x7fffffff) != 0x7f800000) {
 			/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
-			return (cpackf(y - y, y - y));
+			return (CMPLXF(y - y, y - y));
 		} else if (hx & 0x80000000) {
 			/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
-			return (cpackf(0.0, 0.0));
+			return (CMPLXF(0.0, 0.0));
 		} else {
 			/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
-			return (cpackf(x, y - y));
+			return (CMPLXF(x, y - y));
 		}
 	}
 
@@ -84,6 +84,6 @@ cexpf(float complex z)
 		 *  -  x = NaN (spurious inexact exception from y)
 		 */
 		exp_x = expf(x);
-		return (cpackf(exp_x * cosf(y), exp_x * sinf(y)));
+		return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y)));
 	}
 }
diff --git a/lib/msun/src/s_conj.c b/lib/msun/src/s_conj.c
index 5770c29..d0dd051 100644
--- a/lib/msun/src/s_conj.c
+++ b/lib/msun/src/s_conj.c
@@ -34,5 +34,5 @@ double complex
 conj(double complex z)
 {
 
-	return (cpack(creal(z), -cimag(z)));
+	return (CMPLX(creal(z), -cimag(z)));
 }
diff --git a/lib/msun/src/s_conjf.c b/lib/msun/src/s_conjf.c
index b090760..c4cf127 100644
--- a/lib/msun/src/s_conjf.c
+++ b/lib/msun/src/s_conjf.c
@@ -34,5 +34,5 @@ float complex
 conjf(float complex z)
 {
 
-	return (cpackf(crealf(z), -cimagf(z)));
+	return (CMPLXF(crealf(z), -cimagf(z)));
 }
diff --git a/lib/msun/src/s_conjl.c b/lib/msun/src/s_conjl.c
index 0e431ef..a4604b6 100644
--- a/lib/msun/src/s_conjl.c
+++ b/lib/msun/src/s_conjl.c
@@ -34,5 +34,5 @@ long double complex
 conjl(long double complex z)
 {
 
-	return (cpackl(creall(z), -cimagl(z)));
+	return (CMPLXL(creall(z), -cimagl(z)));
 }
diff --git a/lib/msun/src/s_cproj.c b/lib/msun/src/s_cproj.c
index 8e9404c..2f0db6e 100644
--- a/lib/msun/src/s_cproj.c
+++ b/lib/msun/src/s_cproj.c
@@ -39,7 +39,7 @@ cproj(double complex z)
 	if (!isinf(creal(z)) && !isinf(cimag(z)))
 		return (z);
 	else
-		return (cpack(INFINITY, copysign(0.0, cimag(z))));
+		return (CMPLX(INFINITY, copysign(0.0, cimag(z))));
 }
 
 #if LDBL_MANT_DIG == 53
diff --git a/lib/msun/src/s_cprojf.c b/lib/msun/src/s_cprojf.c
index 68ea77b..bd27bdc 100644
--- a/lib/msun/src/s_cprojf.c
+++ b/lib/msun/src/s_cprojf.c
@@ -39,5 +39,5 @@ cprojf(float complex z)
 	if (!isinf(crealf(z)) && !isinf(cimagf(z)))
 		return (z);
 	else
-		return (cpackf(INFINITY, copysignf(0.0, cimagf(z))));
+		return (CMPLXF(INFINITY, copysignf(0.0, cimagf(z))));
 }
diff --git a/lib/msun/src/s_cprojl.c b/lib/msun/src/s_cprojl.c
index 07385bc..5cccc91 100644
--- a/lib/msun/src/s_cprojl.c
+++ b/lib/msun/src/s_cprojl.c
@@ -39,5 +39,5 @@ cprojl(long double complex z)
 	if (!isinf(creall(z)) && !isinf(cimagl(z)))
 		return (z);
 	else
-		return (cpackl(INFINITY, copysignl(0.0, cimagl(z))));
+		return (CMPLXL(INFINITY, copysignl(0.0, cimagl(z))));
 }
diff --git a/lib/msun/src/s_csinh.c b/lib/msun/src/s_csinh.c
index c192f30..37f46da 100644
--- a/lib/msun/src/s_csinh.c
+++ b/lib/msun/src/s_csinh.c
@@ -62,23 +62,23 @@ csinh(double complex z)
 	/* Handle the nearly-non-exceptional cases where x and y are finite. */
 	if (ix < 0x7ff00000 && iy < 0x7ff00000) {
 		if ((iy | ly) == 0)
-			return (cpack(sinh(x), y));
+			return (CMPLX(sinh(x), y));
 		if (ix < 0x40360000)	/* small x: normal case */
-			return (cpack(sinh(x) * cos(y), cosh(x) * sin(y)));
+			return (CMPLX(sinh(x) * cos(y), cosh(x) * sin(y)));
 
 		/* |x| >= 22, so cosh(x) ~= exp(|x|) */
 		if (ix < 0x40862e42) {
 			/* x < 710: exp(|x|) won't overflow */
 			h = exp(fabs(x)) * 0.5;
-			return (cpack(copysign(h, x) * cos(y), h * sin(y)));
+			return (CMPLX(copysign(h, x) * cos(y), h * sin(y)));
 		} else if (ix < 0x4096bbaa) {
 			/* x < 1455: scale to avoid overflow */
-			z = __ldexp_cexp(cpack(fabs(x), y), -1);
-			return (cpack(creal(z) * copysign(1, x), cimag(z)));
+			z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
+			return (CMPLX(creal(z) * copysign(1, x), cimag(z)));
 		} else {
 			/* x >= 1455: the result always overflows */
 			h = huge * x;
-			return (cpack(h * cos(y), h * h * sin(y)));
+			return (CMPLX(h * cos(y), h * h * sin(y)));
 		}
 	}
 
@@ -92,7 +92,7 @@ csinh(double complex z)
 	 * the same as d(NaN).
 	 */
 	if ((ix | lx) == 0 && iy >= 0x7ff00000)
-		return (cpack(copysign(0, x * (y - y)), y - y));
+		return (CMPLX(copysign(0, x * (y - y)), y - y));
 
 	/*
 	 * sinh(+-Inf +- I 0) = +-Inf + I +-0.
@@ -101,8 +101,8 @@ csinh(double complex z)
 	 */
 	if ((iy | ly) == 0 && ix >= 0x7ff00000) {
 		if (((hx & 0xfffff) | lx) == 0)
-			return (cpack(x, y));
-		return (cpack(x, copysign(0, y)));
+			return (CMPLX(x, y));
+		return (CMPLX(x, copysign(0, y)));
 	}
 
 	/*
@@ -114,7 +114,7 @@ csinh(double complex z)
 	 * nonzero x.  Choice = don't raise (except for signaling NaNs).
 	 */
 	if (ix < 0x7ff00000 && iy >= 0x7ff00000)
-		return (cpack(y - y, x * (y - y)));
+		return (CMPLX(y - y, x * (y - y)));
 
 	/*
 	 * sinh(+-Inf + I NaN)  = +-Inf + I d(NaN).
@@ -129,8 +129,8 @@ csinh(double complex z)
 	 */
 	if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
 		if (iy >= 0x7ff00000)
-			return (cpack(x * x, x * (y - y)));
-		return (cpack(x * cos(y), INFINITY * sin(y)));
+			return (CMPLX(x * x, x * (y - y)));
+		return (CMPLX(x * cos(y), INFINITY * sin(y)));
 	}
 
 	/*
@@ -144,7 +144,7 @@ csinh(double complex z)
 	 * Optionally raises the invalid floating-point exception for finite
 	 * nonzero y.  Choice = don't raise (except for signaling NaNs).
 	 */
-	return (cpack((x * x) * (y - y), (x + x) * (y - y)));
+	return (CMPLX((x * x) * (y - y), (x + x) * (y - y)));
 }
 
 double complex
@@ -152,6 +152,6 @@ csin(double complex z)
 {
 
 	/* csin(z) = -I * csinh(I * z) */
-	z = csinh(cpack(-cimag(z), creal(z)));
-	return (cpack(cimag(z), -creal(z)));
+	z = csinh(CMPLX(-cimag(z), creal(z)));
+	return (CMPLX(cimag(z), -creal(z)));
 }
diff --git a/lib/msun/src/s_csinhf.c b/lib/msun/src/s_csinhf.c
index c523125..6348272 100644
--- a/lib/msun/src/s_csinhf.c
+++ b/lib/msun/src/s_csinhf.c
@@ -55,51 +55,51 @@ csinhf(float complex z)
 
 	if (ix < 0x7f800000 && iy < 0x7f800000) {
 		if (iy == 0)
-			return (cpackf(sinhf(x), y));
+			return (CMPLXF(sinhf(x), y));
 		if (ix < 0x41100000)	/* small x: normal case */
-			return (cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y)));
+			return (CMPLXF(sinhf(x) * cosf(y), coshf(x) * sinf(y)));
 
 		/* |x| >= 9, so cosh(x) ~= exp(|x|) */
 		if (ix < 0x42b17218) {
 			/* x < 88.7: expf(|x|) won't overflow */
 			h = expf(fabsf(x)) * 0.5f;
-			return (cpackf(copysignf(h, x) * cosf(y), h * sinf(y)));
+			return (CMPLXF(copysignf(h, x) * cosf(y), h * sinf(y)));
 		} else if (ix < 0x4340b1e7) {
 			/* x < 192.7: scale to avoid overflow */
-			z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
-			return (cpackf(crealf(z) * copysignf(1, x), cimagf(z)));
+			z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
+			return (CMPLXF(crealf(z) * copysignf(1, x), cimagf(z)));
 		} else {
 			/* x >= 192.7: the result always overflows */
 			h = huge * x;
-			return (cpackf(h * cosf(y), h * h * sinf(y)));
+			return (CMPLXF(h * cosf(y), h * h * sinf(y)));
 		}
 	}
 
 	if (ix == 0 && iy >= 0x7f800000)
-		return (cpackf(copysignf(0, x * (y - y)), y - y));
+		return (CMPLXF(copysignf(0, x * (y - y)), y - y));
 
 	if (iy == 0 && ix >= 0x7f800000) {
 		if ((hx & 0x7fffff) == 0)
-			return (cpackf(x, y));
-		return (cpackf(x, copysignf(0, y)));
+			return (CMPLXF(x, y));
+		return (CMPLXF(x, copysignf(0, y)));
 	}
 
 	if (ix < 0x7f800000 && iy >= 0x7f800000)
-		return (cpackf(y - y, x * (y - y)));
+		return (CMPLXF(y - y, x * (y - y)));
 
 	if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
 		if (iy >= 0x7f800000)
-			return (cpackf(x * x, x * (y - y)));
-		return (cpackf(x * cosf(y), INFINITY * sinf(y)));
+			return (CMPLXF(x * x, x * (y - y)));
+		return (CMPLXF(x * cosf(y), INFINITY * sinf(y)));
 	}
 
-	return (cpackf((x * x) * (y - y), (x + x) * (y - y)));
+	return (CMPLXF((x * x) * (y - y), (x + x) * (y - y)));
 }
 
 float complex
 csinf(float complex z)
 {
 
-	z = csinhf(cpackf(-cimagf(z), crealf(z)));
-	return (cpackf(cimagf(z), -crealf(z)));
+	z = csinhf(CMPLXF(-cimagf(z), crealf(z)));
+	return (CMPLXF(cimagf(z), -crealf(z)));
 }
diff --git a/lib/msun/src/s_csqrt.c b/lib/msun/src/s_csqrt.c
index 18a7ae3..75e12d8 100644
--- a/lib/msun/src/s_csqrt.c
+++ b/lib/msun/src/s_csqrt.c
@@ -58,12 +58,12 @@ csqrt(double complex z)
 
 	/* Handle special cases. */
 	if (z == 0)
-		return (cpack(0, b));
+		return (CMPLX(0, b));
 	if (isinf(b))
-		return (cpack(INFINITY, b));
+		return (CMPLX(INFINITY, b));
 	if (isnan(a)) {
 		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
-		return (cpack(a, t));	/* return NaN + NaN i */
+		return (CMPLX(a, t));	/* return NaN + NaN i */
 	}
 	if (isinf(a)) {
 		/*
@@ -73,9 +73,9 @@ csqrt(double complex z)
 		 * csqrt(-inf + y i)   = 0   +  inf i
 		 */
 		if (signbit(a))
-			return (cpack(fabs(b - b), copysign(a, b)));
+			return (CMPLX(fabs(b - b), copysign(a, b)));
 		else
-			return (cpack(a, copysign(b - b, b)));
+			return (CMPLX(a, copysign(b - b, b)));
 	}
 	/*
 	 * The remaining special case (b is NaN) is handled just fine by
@@ -94,10 +94,10 @@ csqrt(double complex z)
 	/* Algorithm 312, CACM vol 10, Oct 1967. */
 	if (a >= 0) {
 		t = sqrt((a + hypot(a, b)) * 0.5);
-		result = cpack(t, b / (2 * t));
+		result = CMPLX(t, b / (2 * t));
 	} else {
 		t = sqrt((-a + hypot(a, b)) * 0.5);
-		result = cpack(fabs(b) / (2 * t), copysign(t, b));
+		result = CMPLX(fabs(b) / (2 * t), copysign(t, b));
 	}
 
 	/* Rescale. */
diff --git a/lib/msun/src/s_csqrtf.c b/lib/msun/src/s_csqrtf.c
index da7fe18..7ee513f 100644
--- a/lib/msun/src/s_csqrtf.c
+++ b/lib/msun/src/s_csqrtf.c
@@ -49,12 +49,12 @@ csqrtf(float complex z)
 
 	/* Handle special cases. */
 	if (z == 0)
-		return (cpackf(0, b));
+		return (CMPLXF(0, b));
 	if (isinf(b))
-		return (cpackf(INFINITY, b));
+		return (CMPLXF(INFINITY, b));
 	if (isnan(a)) {
 		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
-		return (cpackf(a, t));	/* return NaN + NaN i */
+		return (CMPLXF(a, t));	/* return NaN + NaN i */
 	}
 	if (isinf(a)) {
 		/*
@@ -64,9 +64,9 @@ csqrtf(float complex z)
 		 * csqrtf(-inf + y i)   = 0   +  inf i
 		 */
 		if (signbit(a))
-			return (cpackf(fabsf(b - b), copysignf(a, b)));
+			return (CMPLXF(fabsf(b - b), copysignf(a, b)));
 		else
-			return (cpackf(a, copysignf(b - b, b)));
+			return (CMPLXF(a, copysignf(b - b, b)));
 	}
 	/*
 	 * The remaining special case (b is NaN) is handled just fine by
@@ -80,9 +80,9 @@ csqrtf(float complex z)
 	 */
 	if (a >= 0) {
 		t = sqrt((a + hypot(a, b)) * 0.5);
-		return (cpackf(t, b / (2.0 * t)));
+		return (CMPLXF(t, b / (2.0 * t)));
 	} else {
 		t = sqrt((-a + hypot(a, b)) * 0.5);
-		return (cpackf(fabsf(b) / (2.0 * t), copysignf(t, b)));
+		return (CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b)));
 	}
 }
diff --git a/lib/msun/src/s_csqrtl.c b/lib/msun/src/s_csqrtl.c
index dd18e1e..ea2ef27 100644
--- a/lib/msun/src/s_csqrtl.c
+++ b/lib/msun/src/s_csqrtl.c
@@ -58,12 +58,12 @@ csqrtl(long double complex z)
 
 	/* Handle special cases. */
 	if (z == 0)
-		return (cpackl(0, b));
+		return (CMPLXL(0, b));
 	if (isinf(b))
-		return (cpackl(INFINITY, b));
+		return (CMPLXL(INFINITY, b));
 	if (isnan(a)) {
 		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
-		return (cpackl(a, t));	/* return NaN + NaN i */
+		return (CMPLXL(a, t));	/* return NaN + NaN i */
 	}
 	if (isinf(a)) {
 		/*
@@ -73,9 +73,9 @@ csqrtl(long double complex z)
 		 * csqrt(-inf + y i)   = 0   +  inf i
 		 */
 		if (signbit(a))
-			return (cpackl(fabsl(b - b), copysignl(a, b)));
+			return (CMPLXL(fabsl(b - b), copysignl(a, b)));
 		else
-			return (cpackl(a, copysignl(b - b, b)));
+			return (CMPLXL(a, copysignl(b - b, b)));
 	}
 	/*
 	 * The remaining special case (b is NaN) is handled just fine by
@@ -94,10 +94,10 @@ csqrtl(long double complex z)
 	/* Algorithm 312, CACM vol 10, Oct 1967. */
 	if (a >= 0) {
 		t = sqrtl((a + hypotl(a, b)) * 0.5);
-		result = cpackl(t, b / (2 * t));
+		result = CMPLXL(t, b / (2 * t));
 	} else {
 		t = sqrtl((-a + hypotl(a, b)) * 0.5);
-		result = cpackl(fabsl(b) / (2 * t), copysignl(t, b));
+		result = CMPLXL(fabsl(b) / (2 * t), copysignl(t, b));
 	}
 
 	/* Rescale. */
diff --git a/lib/msun/src/s_ctanh.c b/lib/msun/src/s_ctanh.c
index d427e28..b15c814 100644
--- a/lib/msun/src/s_ctanh.c
+++ b/lib/msun/src/s_ctanh.c
@@ -102,9 +102,9 @@ ctanh(double complex z)
 	 */
 	if (ix >= 0x7ff00000) {
 		if ((ix & 0xfffff) | lx)	/* x is NaN */
-			return (cpack(x, (y == 0 ? y : x * y)));
+			return (CMPLX(x, (y == 0 ? y : x * y)));
 		SET_HIGH_WORD(x, hx - 0x40000000);	/* x = copysign(1, x) */
-		return (cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y))));
+		return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y))));
 	}
 
 	/*
@@ -112,7 +112,7 @@ ctanh(double complex z)
 	 * ctanh(x +- i Inf) = NaN + i NaN
 	 */
 	if (!isfinite(y))
-		return (cpack(y - y, y - y));
+		return (CMPLX(y - y, y - y));
 
 	/*
 	 * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the
@@ -121,7 +121,7 @@ ctanh(double complex z)
 	 */
 	if (ix >= 0x40360000) {	/* x >= 22 */
 		double exp_mx = exp(-fabs(x));
-		return (cpack(copysign(1, x),
+		return (CMPLX(copysign(1, x),
 		    4 * sin(y) * cos(y) * exp_mx * exp_mx));
 	}
 
@@ -131,7 +131,7 @@ ctanh(double complex z)
 	s = sinh(x);
 	rho = sqrt(1 + s * s);	/* = cosh(x) */
 	denom = 1 + beta * s * s;
-	return (cpack((beta * rho * s) / denom, t / denom));
+	return (CMPLX((beta * rho * s) / denom, t / denom));
 }
 
 double complex
@@ -139,6 +139,6 @@ ctan(double complex z)
 {
 
 	/* ctan(z) = -I * ctanh(I * z) */
-	z = ctanh(cpack(-cimag(z), creal(z)));
-	return (cpack(cimag(z), -creal(z)));
+	z = ctanh(CMPLX(-cimag(z), creal(z)));
+	return (CMPLX(cimag(z), -creal(z)));
 }
diff --git a/lib/msun/src/s_ctanhf.c b/lib/msun/src/s_ctanhf.c
index 4be28d8..9b0cda1 100644
--- a/lib/msun/src/s_ctanhf.c
+++ b/lib/msun/src/s_ctanhf.c
@@ -51,18 +51,18 @@ ctanhf(float complex z)
 
 	if (ix >= 0x7f800000) {
 		if (ix & 0x7fffff)
-			return (cpackf(x, (y == 0 ? y : x * y)));
+			return (CMPLXF(x, (y == 0 ? y : x * y)));
 		SET_FLOAT_WORD(x, hx - 0x40000000);
-		return (cpackf(x,
+		return (CMPLXF(x,
 		    copysignf(0, isinf(y) ? y : sinf(y) * cosf(y))));
 	}
 
 	if (!isfinite(y))
-		return (cpackf(y - y, y - y));
+		return (CMPLXF(y - y, y - y));
 
 	if (ix >= 0x41300000) {	/* x >= 11 */
 		float exp_mx = expf(-fabsf(x));
-		return (cpackf(copysignf(1, x),
+		return (CMPLXF(copysignf(1, x),
 		    4 * sinf(y) * cosf(y) * exp_mx * exp_mx));
 	}
 
@@ -71,14 +71,14 @@ ctanhf(float complex z)
 	s = sinhf(x);
 	rho = sqrtf(1 + s * s);
 	denom = 1 + beta * s * s;
-	return (cpackf((beta * rho * s) / denom, t / denom));
+	return (CMPLXF((beta * rho * s) / denom, t / denom));
 }
 
 float complex
 ctanf(float complex z)
 {
 
-	z = ctanhf(cpackf(-cimagf(z), crealf(z)));
-	return (cpackf(cimagf(z), -crealf(z)));
+	z = ctanhf(CMPLXF(-cimagf(z), crealf(z)));
+	return (CMPLXF(cimagf(z), -crealf(z)));
 }
 
diff --git a/lib/msun/src/s_scalbln.c b/lib/msun/src/s_scalbln.c
index d609d4e..5e4e42e 100644
--- a/lib/msun/src/s_scalbln.c
+++ b/lib/msun/src/s_scalbln.c
@@ -27,50 +27,28 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-#include <limits.h>
 #include <math.h>
 
+#define	NMAX	65536
+#define	NMIN	-65536
+
 double
-scalbln (double x, long n)
+scalbln(double x, long n)
 {
-	int in;
 
-	in = (int)n;
-	if (in != n) {
-		if (n > 0)
-			in = INT_MAX;
-		else
-			in = INT_MIN;
-	}
-	return (scalbn(x, in));
+	return (scalbn(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
 }
 
 float
-scalblnf (float x, long n)
+scalblnf(float x, long n)
 {
-	int in;
 
-	in = (int)n;
-	if (in != n) {
-		if (n > 0)
-			in = INT_MAX;
-		else
-			in = INT_MIN;
-	}
-	return (scalbnf(x, in));
+	return (scalbnf(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
 }
 
 long double
-scalblnl (long double x, long n)
+scalblnl(long double x, long n)
 {
-	int in;
 
-	in = (int)n;
-	if (in != n) {
-		if (n > 0)
-			in = INT_MAX;
-		else
-			in = INT_MIN;
-	}
-	return (scalbnl(x, (int)n));
+	return (scalbnl(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
 }
diff --git a/lib/msun/tests/Makefile b/lib/msun/tests/Makefile
index 4261e48..0479cfb 100644
--- a/lib/msun/tests/Makefile
+++ b/lib/msun/tests/Makefile
@@ -6,13 +6,12 @@ TESTSRC=	${SRCTOP}/contrib/netbsd-tests/lib/libm
 
 TESTSDIR=	${TESTSBASE}/lib/msun
 
-.if ${MACHINE} == "sparc" || ${MACHINE} == "i386" \
-	|| ${MACHINE} == "amd64" || ${MACHINE_CPU} == "arm"	\
-	|| ${MACHINE} == "sparc64"
+# All architectures on FreeBSD have fenv.h
 CFLAGS+=	-DHAVE_FENV_H
-.endif
 
-.if ${MACHINE} == "amd64" || ${MACHINE} == "i386"
+# Not sure why this isn't defined for all architectures, since most
+# have long double.
+.if ${MACHINE_CPUARCH} == "amd64" || ${MACHINE_CPUARCH} == "i386"
 CFLAGS+=	-D__HAVE_LONG_DOUBLE
 .endif
 
@@ -41,8 +40,7 @@ NETBSD_ATF_TESTS_C+=	tanh_test
 
 CSTD=		c99
 
-LDADD+=		-lm
-DPADD+=		${LIBM}
+LIBADD+=	m
 #COPTS+=	-Wfloat-equal
 
 # Copied from lib/msun/Makefile