diff --git a/libm/upstream-freebsd/lib/msun/src/e_hypotl.c b/libm/upstream-freebsd/lib/msun/src/e_hypotl.c
index 9189b1f..fc43538 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_hypotl.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_hypotl.c
@@ -82,7 +82,7 @@
 	        man_t manh, manl;
 		GET_LDBL_MAN(manh,manl,b);
 		if((manh|manl)==0) return a;
-		t1=0;
+		t1=1;
 		SET_HIGH_WORD(t1,ESW(MAX_EXP-2));	/* t1=2^(MAX_EXP-2) */
 		b *= t1;
 		a *= t1;
diff --git a/libm/upstream-freebsd/lib/msun/src/e_powf.c b/libm/upstream-freebsd/lib/msun/src/e_powf.c
index 53f1d37..122da45 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_powf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_powf.c
@@ -136,7 +136,7 @@
     /* |y| is huge */
 	if(iy>0x4d000000) { /* if |y| > 2**27 */
 	/* over/underflow if x is not close to one */
-	    if(ix<0x3f7ffff8) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
+	    if(ix<0x3f7ffff6) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
 	    if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny;
 	/* now |1-x| is tiny <= 2**-20, suffice to compute
 	   log(x) by x-x^2/2+x^3/3-x^4/4 */
diff --git a/libm/upstream-freebsd/lib/msun/src/e_sqrt.c b/libm/upstream-freebsd/lib/msun/src/e_sqrt.c
index 12fb56e..37351a4 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_sqrt.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_sqrt.c
@@ -14,6 +14,18 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef USE_BUILTIN_SQRT
+double
+__ieee754_sqrt(double x)
+{
+	return (__builtin_sqrt(x));
+}
+#else
 /* __ieee754_sqrt(x)
  * Return correctly rounded sqrt.
  *           ------------------------------------------
@@ -84,11 +96,6 @@
  *---------------
  */
 
-#include <float.h>
-
-#include "math.h"
-#include "math_private.h"
-
 static	const double	one	= 1.0, tiny=1.0e-300;
 
 double
@@ -187,6 +194,7 @@
 	INSERT_WORDS(z,ix0,ix1);
 	return z;
 }
+#endif
 
 #if (LDBL_MANT_DIG == 53)
 __weak_reference(sqrt, sqrtl);
diff --git a/libm/upstream-freebsd/lib/msun/src/e_sqrtf.c b/libm/upstream-freebsd/lib/msun/src/e_sqrtf.c
index 7eba4d0..06e5d62 100644
--- a/libm/upstream-freebsd/lib/msun/src/e_sqrtf.c
+++ b/libm/upstream-freebsd/lib/msun/src/e_sqrtf.c
@@ -20,6 +20,13 @@
 #include "math.h"
 #include "math_private.h"
 
+#ifdef USE_BUILTIN_SQRTF
+float
+__ieee754_sqrtf(float x)
+{
+	return (__builtin_sqrtf(x));
+}
+#else
 static	const float	one	= 1.0, tiny=1.0e-30;
 
 float
@@ -87,3 +94,4 @@
 	SET_FLOAT_WORD(z,ix);
 	return z;
 }
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/k_cospi.h b/libm/upstream-freebsd/lib/msun/src/k_cospi.h
new file mode 100644
index 0000000..32985be
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/k_cospi.h
@@ -0,0 +1,44 @@
+/*-
+ * Copyright (c) 2017 Steven G. Kargl
+ * 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 unmodified, 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 ``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 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 basic kernel for x in [0,0.25].  To use the kernel for cos(x), the
+ * argument to __kernel_cospi() must be multiplied by pi.
+ */
+
+static inline double
+__kernel_cospi(double x)
+{
+	double_t hi, lo;
+
+	hi = (float)x;
+	lo = x - hi;
+	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
+	hi *= pi_hi;
+	_2sumF(hi, lo);
+	return (__kernel_cos(hi, lo));
+}
+
diff --git a/libm/upstream-freebsd/lib/msun/src/k_sincosl.h b/libm/upstream-freebsd/lib/msun/src/k_sincosl.h
index 4d4dc69..6425f14 100644
--- a/libm/upstream-freebsd/lib/msun/src/k_sincosl.h
+++ b/libm/upstream-freebsd/lib/msun/src/k_sincosl.h
@@ -76,13 +76,6 @@
 #elif LDBL_MANT_DIG == 113	/* ld128 version of k_sincosl.c. */
 
 static const long double
-C1 =  0.04166666666666666666666666666666658424671L,
-C2 = -0.001388888888888888888888888888863490893732L,
-C3 =  0.00002480158730158730158730158600795304914210L,
-C4 = -0.2755731922398589065255474947078934284324e-6L,
-C5 =  0.2087675698786809897659225313136400793948e-8L,
-C6 = -0.1147074559772972315817149986812031204775e-10L,
-C7 =  0.4779477332386808976875457937252120293400e-13L,
 S1 = -0.16666666666666666666666666666666666606732416116558L,
 S2 =  0.0083333333333333333333333333333331135404851288270047L,
 S3 = -0.00019841269841269841269841269839935785325638310428717L,
@@ -93,15 +86,25 @@
 S8 =  0.28114572543451292625024967174638477283187397621303e-14L;
 
 static const double
-C8  = -0.1561920696721507929516718307820958119868e-15,
-C9  =  0.4110317413744594971475941557607804508039e-18,
-C10 = -0.8896592467191938803288521958313920156409e-21,
-C11 =  0.1601061435794535138244346256065192782581e-23,
 S9  = -0.82206352458348947812512122163446202498005154296863e-17,
 S10 =  0.19572940011906109418080609928334380560135358385256e-19,
 S11 = -0.38680813379701966970673724299207480965452616911420e-22,
 S12 =  0.64038150078671872796678569586315881020659912139412e-25;
 
+static const long double
+C1 =  4.16666666666666666666666666666666667e-02L,
+C2 = -1.38888888888888888888888888888888834e-03L,
+C3 =  2.48015873015873015873015873015446795e-05L,
+C4 = -2.75573192239858906525573190949988493e-07L,
+C5 =  2.08767569878680989792098886701451072e-09L,
+C6 = -1.14707455977297247136657111139971865e-11L,
+C7 =  4.77947733238738518870113294139830239e-14L,
+C8 = -1.56192069685858079920640872925306403e-16L,
+C9 =  4.11031762320473354032038893429515732e-19L,
+C10= -8.89679121027589608738005163931958096e-22L,
+C11=  1.61171797801314301767074036661901531e-24L,
+C12= -2.46748624357670948912574279501044295e-27L;
+
 static inline void
 __kernel_sincosl(long double x, long double y, int iy, long double *sn, 
     long double *cs)
@@ -120,12 +123,12 @@
 	if (iy == 0)
 		*sn = x + v * (S1 + z * r);
 	else
-		*cs = x - ((z * (y / 2 - v * r) - y) - v * S1);
+		*sn = x - ((z * (y / 2 - v * r) - y) - v * S1);
 
 	hz = z / 2;
 	w = 1 - hz;
 	r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 + 
-	    z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * C11))))))))));
+	    z * (C7 + z * (C8 + z * (C9 + z * (C10 + z * (C11+z*C12)))))))))));
 
 	*cs =  w + (((1 - w) - hz) + (z * r - x * y));
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/k_sinpi.h b/libm/upstream-freebsd/lib/msun/src/k_sinpi.h
new file mode 100644
index 0000000..66152ba
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/k_sinpi.h
@@ -0,0 +1,43 @@
+/*-
+ * Copyright (c) 2017 Steven G. Kargl
+ * 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 unmodified, 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 ``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 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 basic kernel for x in [0,0.25].  To use the kernel for sin(x), the
+ * argument to __kernel_sinpi() must be multiplied by pi.
+ */
+
+static inline double
+__kernel_sinpi(double x)
+{
+	double_t hi, lo;
+
+	hi = (float)x;
+	lo = x - hi;
+	lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
+	hi *= pi_hi;
+	_2sumF(hi, lo);
+	return (__kernel_sin(hi, lo, 1));
+}
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cexp.c b/libm/upstream-freebsd/lib/msun/src/s_cexp.c
index 2ef8ba1..a1f853e 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_cexp.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_cexp.c
@@ -30,6 +30,7 @@
 __FBSDID("$FreeBSD$");
 
 #include <complex.h>
+#include <float.h>
 #include <math.h>
 
 #include "math_private.h"
@@ -41,7 +42,7 @@
 double complex
 cexp(double complex z)
 {
-	double x, y, exp_x;
+	double c, exp_x, s, x, y;
 	uint32_t hx, hy, lx, ly;
 
 	x = creal(z);
@@ -55,8 +56,10 @@
 		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 (CMPLX(cos(y), sin(y)));
+	if (((hx & 0x7fffffff) | lx) == 0) {
+		sincos(y, &s, &c);
+		return (CMPLX(c, s));
+	}
 
 	if (hy >= 0x7ff00000) {
 		if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
@@ -86,6 +89,11 @@
 		 *  -  x = NaN (spurious inexact exception from y)
 		 */
 		exp_x = exp(x);
-		return (CMPLX(exp_x * cos(y), exp_x * sin(y)));
+		sincos(y, &s, &c);
+		return (CMPLX(exp_x * c, exp_x * s));
 	}
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(cexp, cexpl);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cexpf.c b/libm/upstream-freebsd/lib/msun/src/s_cexpf.c
index b815c99..d905b74 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_cexpf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_cexpf.c
@@ -41,7 +41,7 @@
 float complex
 cexpf(float complex z)
 {
-	float x, y, exp_x;
+	float c, exp_x, s, x, y;
 	uint32_t hx, hy;
 
 	x = crealf(z);
@@ -55,8 +55,10 @@
 		return (CMPLXF(expf(x), y));
 	GET_FLOAT_WORD(hx, x);
 	/* cexp(0 + I y) = cos(y) + I sin(y) */
-	if ((hx & 0x7fffffff) == 0)
-		return (CMPLXF(cosf(y), sinf(y)));
+	if ((hx & 0x7fffffff) == 0) {
+		sincosf(y, &s, &c);
+		return (CMPLXF(c, s));
+	}
 
 	if (hy >= 0x7f800000) {
 		if ((hx & 0x7fffffff) != 0x7f800000) {
@@ -86,6 +88,7 @@
 		 *  -  x = NaN (spurious inexact exception from y)
 		 */
 		exp_x = expf(x);
-		return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y)));
+		sincosf(y, &s, &c);
+		return (CMPLXF(exp_x * c, exp_x * s));
 	}
 }
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cosl.c b/libm/upstream-freebsd/lib/msun/src/s_cosl.c
index 46a2e86..3d06648 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_cosl.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_cosl.c
@@ -39,12 +39,17 @@
 #include <ieeefp.h>
 #endif
 
+#include "fpmath.h"
 #include "math.h"
 #include "math_private.h"
 #if LDBL_MANT_DIG == 64
 #include "../ld80/e_rem_pio2l.h"
+static const union IEEEl2bits
+pio4u = LD80C(0xc90fdaa22168c235, -00001,  7.85398163397448309628e-01L);
+#define	pio4	(pio4u.e)
 #elif LDBL_MANT_DIG == 113
 #include "../ld128/e_rem_pio2l.h"
+long double pio4 =  7.85398163397448309615660845819875721e-1L;
 #else
 #error "Unsupported long double format"
 #endif
@@ -71,7 +76,7 @@
 	ENTERI();
 
 	/* Optimize the case where x is already within range. */
-	if (z.e < M_PI_4)
+	if (z.e < pio4)
 		RETURNI(__kernel_cosl(z.e, 0));
 
 	e0 = __ieee754_rem_pio2l(x, y);
diff --git a/libm/upstream-freebsd/lib/msun/src/s_cospi.c b/libm/upstream-freebsd/lib/msun/src/s_cospi.c
new file mode 100644
index 0000000..860219e
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/s_cospi.c
@@ -0,0 +1,152 @@
+/*-
+ * Copyright (c) 2017 Steven G. Kargl
+ * 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 unmodified, 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 ``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 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.
+ */
+
+/**
+ * cospi(x) computes cos(pi*x) without multiplication by pi (almost).  First,
+ * note that cospi(-x) = cospi(x), so the algorithm considers only |x|.  The
+ * method used depends on the magnitude of x.
+ *
+ * 1. For small |x|, cospi(x) = 1 with FE_INEXACT raised where a sloppy
+ *    threshold is used.  The threshold is |x| < 0x1pN with N = -(P/2+M).
+ *    P is the precision of the floating-point type and M = 2 to 4.
+ *
+ * 2. For |x| < 1, argument reduction is not required and sinpi(x) is 
+ *    computed by calling a kernel that leverages the kernels for sin(x)
+ *    ans cos(x).  See k_sinpi.c and k_cospi.c for details.
+ *
+ * 3. For 1 <= |x| < 0x1p(P-1), argument reduction is required where
+ *    |x| = j0 + r with j0 an integer and the remainder r satisfies
+ *    0 <= r < 1.  With the given domain, a simplified inline floor(x)
+ *    is used.  Also, note the following identity
+ *
+ *    cospi(x) = cos(pi*(j0+r))
+ *             = cos(pi*j0) * cos(pi*r) - sin(pi*j0) * sin(pi*r)
+ *             = cos(pi*j0) * cos(pi*r)
+ *             = +-cospi(r)
+ *
+ *    If j0 is even, then cos(pi*j0) = 1. If j0 is odd, then cos(pi*j0) = -1.
+ *    cospi(r) is then computed via an appropriate kernel.
+ *
+ * 4. For |x| >= 0x1p(P-1), |x| is integral and cospi(x) = 1.
+ *
+ * 5. Special cases:
+ *
+ *    cospi(+-0) = 1.
+ *    cospi(n.5) = 0 for n an integer.
+ *    cospi(+-inf) = nan.  Raises the "invalid" floating-point exception.
+ *    cospi(nan) = nan.  Raises the "invalid" floating-point exception.
+ */
+
+#include <float.h>
+#include "math.h"
+#include "math_private.h"
+
+static const double
+pi_hi = 3.1415926814079285e+00,	/* 0x400921fb 0x58000000 */
+pi_lo =-2.7818135228334233e-08;	/* 0xbe5dde97 0x3dcb3b3a */
+
+#include "k_cospi.h"
+#include "k_sinpi.h"
+
+volatile static const double vzero = 0;
+
+double
+cospi(double x)
+{
+	double ax, c;
+	uint32_t hx, ix, j0, lx;
+
+	EXTRACT_WORDS(hx, lx, x);
+	ix = hx & 0x7fffffff;
+	INSERT_WORDS(ax, ix, lx);
+
+	if (ix < 0x3ff00000) {			/* |x| < 1 */
+		if (ix < 0x3fd00000) {		/* |x| < 0.25 */
+			if (ix < 0x3e200000) {	/* |x| < 0x1p-29 */
+				if ((int)ax == 0)
+					return (1);
+			}
+			return (__kernel_cospi(ax));
+		}
+
+		if (ix < 0x3fe00000)		/* |x| < 0.5 */
+			c = __kernel_sinpi(0.5 - ax);
+		else if (ix < 0x3fe80000){	/* |x| < 0.75 */
+			if (ax == 0.5)
+				return (0);
+			c = -__kernel_sinpi(ax - 0.5);
+		} else
+			c = -__kernel_cospi(1 - ax);
+		return (c);
+	}
+
+	if (ix < 0x43300000) {		/* 1 <= |x| < 0x1p52 */
+		/* Determine integer part of ax. */
+		j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
+		if (j0 < 20) {
+			ix &= ~(0x000fffff >> j0);
+			lx = 0;
+		} else {
+			lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
+		}
+		INSERT_WORDS(x, ix, lx);
+
+		ax -= x;
+		EXTRACT_WORDS(ix, lx, ax);
+
+
+		if (ix < 0x3fe00000) {		/* |x| < 0.5 */
+			if (ix < 0x3fd00000)	/* |x| < 0.25 */
+				c = ix == 0 ? 1 : __kernel_cospi(ax);
+			else 
+				c = __kernel_sinpi(0.5 - ax);
+		} else {
+			if (ix < 0x3fe80000) {	/* |x| < 0.75 */
+				if (ax == 0.5)
+					return (0);
+				c = -__kernel_sinpi(ax - 0.5);
+			} else
+				c = -__kernel_cospi(1 - ax);
+		}
+
+		if (j0 > 30)
+			x -= 0x1p30;
+		j0 = (uint32_t)x;
+		return (j0 & 1 ? -c : c);
+	}
+
+	if (ix >= 0x7f800000)
+		return (vzero / vzero);
+
+	/*
+	 * |x| >= 0x1p52 is always an even integer, so return 1.
+	 */
+	return (1);
+}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(cospi, cospil);
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_ctanh.c b/libm/upstream-freebsd/lib/msun/src/s_ctanh.c
index 88afeb5..e5840a1 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_ctanh.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_ctanh.c
@@ -111,11 +111,13 @@
 	}
 
 	/*
-	 * ctanh(x + I NaN) = d(NaN) + I d(NaN)
-	 * ctanh(x +- I Inf) = dNaN + I dNaN
+	 * ctanh(+-0 + i NAN) = +-0 + i NaN
+	 * ctanh(+-0 +- i Inf) = +-0 + i NaN
+	 * ctanh(x + i NAN) = NaN + i NaN
+	 * ctanh(x +- i Inf) = NaN + i NaN
 	 */
 	if (!isfinite(y))
-		return (CMPLX(y - y, y - y));
+		return (CMPLX(x ? y - y : x, y - y));
 
 	/*
 	 * ctanh(+-huge +- I y) ~= +-1 +- I 2sin(2y)/exp(2x), using the
diff --git a/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c b/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c
index d2bd0b6..c46f86d 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_ctanhf.c
@@ -61,7 +61,7 @@
 	}
 
 	if (!isfinite(y))
-		return (CMPLXF(y - y, y - y));
+		return (CMPLXF(ix ? y - y : x, y - y));
 
 	if (ix >= 0x41300000) {	/* |x| >= 11 */
 		float exp_mx = expf(-fabsf(x));
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fma.c b/libm/upstream-freebsd/lib/msun/src/s_fma.c
index 41a6424..95cffd0 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fma.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fma.c
@@ -35,6 +35,13 @@
 
 #include "math_private.h"
 
+#ifdef USE_BUILTIN_FMA
+double
+fma(double x, double y, double z)
+{
+	return (__builtin_fma(x, y, z));
+}
+#else
 /*
  * A struct dd represents a floating-point number with twice the precision
  * of a double.  We maintain the invariant that "hi" stores the 53 high-order
@@ -284,6 +291,7 @@
 	else
 		return (add_and_denormalize(r.hi, adj, spread));
 }
+#endif /* !USE_BUILTIN_FMA */
 
 #if (LDBL_MANT_DIG == 53)
 __weak_reference(fma, fmal);
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fmaf.c b/libm/upstream-freebsd/lib/msun/src/s_fmaf.c
index 389cf1b..4591cc2 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fmaf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fmaf.c
@@ -34,6 +34,13 @@
 #include "math.h"
 #include "math_private.h"
 
+#ifdef USE_BUILTIN_FMAF
+float
+fmaf(float x, float y, float z)
+{
+	return (__builtin_fmaf(x, y, z));
+}
+#else
 /*
  * Fused multiply-add: Compute x * y + z with a single rounding error.
  *
@@ -69,3 +76,4 @@
 		SET_LOW_WORD(adjusted_result, lr + 1);
 	return (adjusted_result);
 }
+#endif /* !USE_BUILTIN_FMAF */
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fmax.c b/libm/upstream-freebsd/lib/msun/src/s_fmax.c
index 0c234bc..b53b1e6 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fmax.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fmax.c
@@ -34,6 +34,13 @@
 
 #include "fpmath.h"
 
+#ifdef USE_BUILTIN_FMAX
+double
+fmax(double x, double y)
+{
+	return (__builtin_fmax(x, y));
+}
+#else
 double
 fmax(double x, double y)
 {
@@ -54,6 +61,7 @@
 
 	return (x > y ? x : y);
 }
+#endif
 
 #if (LDBL_MANT_DIG == 53)
 __weak_reference(fmax, fmaxl);
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fmaxf.c b/libm/upstream-freebsd/lib/msun/src/s_fmaxf.c
index 8e9d1ba..8d3d14f 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fmaxf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fmaxf.c
@@ -33,6 +33,13 @@
 
 #include "fpmath.h"
 
+#ifdef USE_BUILTIN_FMAXF
+float
+fmaxf(float x, float y)
+{
+	return (__builtin_fmaxf(x, y));
+}
+#else
 float
 fmaxf(float x, float y)
 {
@@ -53,3 +60,4 @@
 
 	return (x > y ? x : y);
 }
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fmin.c b/libm/upstream-freebsd/lib/msun/src/s_fmin.c
index d7f24c1..53f36c1 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fmin.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fmin.c
@@ -34,6 +34,13 @@
 
 #include "fpmath.h"
 
+#ifdef USE_BUILTIN_FMIN
+double
+fmin(double x, double y)
+{
+	return (__builtin_fmin(x, y));
+}
+#else
 double
 fmin(double x, double y)
 {
@@ -54,6 +61,7 @@
 
 	return (x < y ? x : y);
 }
+#endif
 
 #if (LDBL_MANT_DIG == 53)
 __weak_reference(fmin, fminl);
diff --git a/libm/upstream-freebsd/lib/msun/src/s_fminf.c b/libm/upstream-freebsd/lib/msun/src/s_fminf.c
index 2583167..58b6a48 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_fminf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_fminf.c
@@ -33,6 +33,13 @@
 
 #include "fpmath.h"
 
+#ifdef USE_BUILTIN_FMINF
+float
+fminf(float x, float y)
+{
+	return (__builtin_fminf(x, y));
+}
+#else
 float
 fminf(float x, float y)
 {
@@ -53,3 +60,4 @@
 
 	return (x < y ? x : y);
 }
+#endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_lround.c b/libm/upstream-freebsd/lib/msun/src/s_lround.c
index 66d9183..1dd8e69 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_lround.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_lround.c
@@ -49,9 +49,11 @@
  * that everything is in range.  At compile time, INRANGE(x) should reduce to
  * two floating-point comparisons in the former case, or TRUE otherwise.
  */
+static const type type_min = (type)DTYPE_MIN;
+static const type type_max = (type)DTYPE_MAX;
 static const type dtype_min = (type)DTYPE_MIN - 0.5;
 static const type dtype_max = (type)DTYPE_MAX + 0.5;
-#define	INRANGE(x)	(dtype_max - (type)DTYPE_MAX != 0.5 || \
+#define	INRANGE(x)	(dtype_max - type_max != 0.5 || \
 			 ((x) > dtype_min && (x) < dtype_max))
 
 dtype
diff --git a/libm/upstream-freebsd/lib/msun/src/s_scalbn.c b/libm/upstream-freebsd/lib/msun/src/s_scalbn.c
index b048b05..2d4f7a3 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_scalbn.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_scalbn.c
@@ -1,66 +1,47 @@
-/* @(#)s_scalbn.c 5.1 93/09/24 */
 /*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2005-2020 Rich Felker, et al.
  *
- * 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.
- * ====================================================
+ * SPDX-License-Identifier: MIT
+ *
+ * Please see https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+ * for all contributors to musl.
  */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
-
-/*
- * scalbn (double x, int n)
- * scalbn(x,n) returns x* 2**n  computed by  exponent
- * manipulation rather than by actually performing an
- * exponentiation or a multiplication.
- */
-
 #include <float.h>
+#include <math.h>
+#include <stdint.h>
 
-#include "math.h"
-#include "math_private.h"
-
-static const double
-two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
-twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
-huge   = 1.0e+300,
-tiny   = 1.0e-300;
-
-double
-scalbn (double x, int n)
+double scalbn(double x, int n)
 {
-	int32_t k,hx,lx;
-	EXTRACT_WORDS(hx,lx,x);
-        k = (hx&0x7ff00000)>>20;		/* extract exponent */
-        if (k==0) {				/* 0 or subnormal x */
-            if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
-	    x *= two54;
-	    GET_HIGH_WORD(hx,x);
-	    k = ((hx&0x7ff00000)>>20) - 54;
-            if (n< -50000) return tiny*x; 	/*underflow*/
-	    }
-        if (k==0x7ff) return x+x;		/* NaN or Inf */
-        k = k+n;
-        if (k >  0x7fe) return huge*copysign(huge,x); /* overflow  */
-        if (k > 0) 				/* normal result */
-	    {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
-        if (k <= -54) {
-            if (n > 50000) 	/* in case integer overflow in n+k */
-		return huge*copysign(huge,x);	/*overflow*/
-	    else
-		return tiny*copysign(tiny,x); 	/*underflow*/
+	union {double f; uint64_t i;} u;
+	double_t y = x;
+
+	if (n > 1023) {
+		y *= 0x1p1023;
+		n -= 1023;
+		if (n > 1023) {
+			y *= 0x1p1023;
+			n -= 1023;
+			if (n > 1023)
+				n = 1023;
+		}
+	} else if (n < -1022) {
+		/* make sure final n < -53 to avoid double
+		   rounding in the subnormal range */
+		y *= 0x1p-1022 * 0x1p53;
+		n += 1022 - 53;
+		if (n < -1022) {
+			y *= 0x1p-1022 * 0x1p53;
+			n += 1022 - 53;
+			if (n < -1022)
+				n = -1022;
+		}
 	}
-        k += 54;				/* subnormal result */
-	SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
-        return x*twom54;
+	u.i = (uint64_t)(0x3ff+n)<<52;
+	x = y * u.f;
+	return x;
 }
 
-#if (LDBL_MANT_DIG == 53)
+#if (LDBL_MANT_DIG == 53) && !defined(scalbn)
 __weak_reference(scalbn, ldexpl);
 __weak_reference(scalbn, scalbnl);
 #endif
diff --git a/libm/upstream-freebsd/lib/msun/src/s_scalbnf.c b/libm/upstream-freebsd/lib/msun/src/s_scalbnf.c
index 21d001c..8cf1e01 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_scalbnf.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_scalbnf.c
@@ -1,57 +1,41 @@
-/* s_scalbnf.c -- float version of s_scalbn.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
 /*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2005-2020 Rich Felker, et al.
  *
- * 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.
- * ====================================================
+ * SPDX-License-Identifier: MIT
+ *
+ * Please see https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+ * for all contributors to musl.
  */
+#include <math.h>
+#include <stdint.h>
 
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
-
-#include "math.h"
-#include "math_private.h"
-
-static const float
-two25   =  3.355443200e+07,	/* 0x4c000000 */
-twom25  =  2.9802322388e-08,	/* 0x33000000 */
-huge   = 1.0e+30,
-tiny   = 1.0e-30;
-
-float
-scalbnf (float x, int n)
+float scalbnf(float x, int n)
 {
-	int32_t k,ix;
-	GET_FLOAT_WORD(ix,x);
-        k = (ix&0x7f800000)>>23;		/* extract exponent */
-        if (k==0) {				/* 0 or subnormal x */
-            if ((ix&0x7fffffff)==0) return x; /* +-0 */
-	    x *= two25;
-	    GET_FLOAT_WORD(ix,x);
-	    k = ((ix&0x7f800000)>>23) - 25;
-            if (n< -50000) return tiny*x; 	/*underflow*/
-	    }
-        if (k==0xff) return x+x;		/* NaN or Inf */
-        k = k+n;
-        if (k >  0xfe) return huge*copysignf(huge,x); /* overflow  */
-        if (k > 0) 				/* normal result */
-	    {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
-        if (k <= -25) {
-            if (n > 50000) 	/* in case integer overflow in n+k */
-		return huge*copysignf(huge,x);	/*overflow*/
-	    else
-		return tiny*copysignf(tiny,x);	/*underflow*/
+	union {float f; uint32_t i;} u;
+	float_t y = x;
+
+	if (n > 127) {
+		y *= 0x1p127f;
+		n -= 127;
+		if (n > 127) {
+			y *= 0x1p127f;
+			n -= 127;
+			if (n > 127)
+				n = 127;
+		}
+	} else if (n < -126) {
+		y *= 0x1p-126f * 0x1p24f;
+		n += 126 - 24;
+		if (n < -126) {
+			y *= 0x1p-126f * 0x1p24f;
+			n += 126 - 24;
+			if (n < -126)
+				n = -126;
+		}
 	}
-        k += 25;				/* subnormal result */
-	SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
-        return x*twom25;
+	u.i = (uint32_t)(0x7f+n)<<23;
+	x = y * u.f;
+	return x;
 }
 
 __strong_reference(scalbnf, ldexpf);
diff --git a/libm/upstream-freebsd/lib/msun/src/s_scalbnl.c b/libm/upstream-freebsd/lib/msun/src/s_scalbnl.c
index 28b0cf9..6044c1b 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_scalbnl.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_scalbnl.c
@@ -1,71 +1,49 @@
-/* @(#)s_scalbn.c 5.1 93/09/24 */
 /*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2005-2020 Rich Felker, et al.
  *
- * 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.
- * ====================================================
+ * SPDX-License-Identifier: MIT
+ *
+ * Please see https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+ * for all contributors to musl.
  */
-
-#include <sys/cdefs.h>
-__FBSDID("$FreeBSD$");
-
+#include <math.h>
+#include <float.h>
+#include "math_private.h"
+#include "fpmath.h"
 /*
  * scalbnl (long double x, int n)
  * scalbnl(x,n) returns x* 2**n  computed by  exponent
  * manipulation rather than by actually performing an
  * exponentiation or a multiplication.
  */
-
-/*
- * We assume that a long double has a 15-bit exponent.  On systems
- * where long double is the same as double, scalbnl() is an alias
- * for scalbn(), so we don't use this routine.
- */
-
-#include <float.h>
-#include <math.h>
-
-#include "fpmath.h"
-
-#if LDBL_MAX_EXP != 0x4000
-#error "Unsupported long double format"
-#endif
-
-static const long double
-huge = 0x1p16000L,
-tiny = 0x1p-16000L;
-
-long double
-scalbnl (long double x, int n)
+#if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
+long double scalbnl(long double x, int n)
 {
 	union IEEEl2bits u;
-	int k;
-	u.e = x;
-        k = u.bits.exp;				/* extract exponent */
-        if (k==0) {				/* 0 or subnormal x */
-            if ((u.bits.manh|u.bits.manl)==0) return x;	/* +-0 */
-	    u.e *= 0x1p+128;
-	    k = u.bits.exp - 128;
-            if (n< -50000) return tiny*x; 	/*underflow*/
-	    }
-        if (k==0x7fff) return x+x;		/* NaN or Inf */
-        k = k+n;
-        if (k >= 0x7fff) return huge*copysignl(huge,x); /* overflow  */
-        if (k > 0) 				/* normal result */
-	    {u.bits.exp = k; return u.e;}
-        if (k <= -128) {
-            if (n > 50000) 	/* in case integer overflow in n+k */
-		return huge*copysign(huge,x);	/*overflow*/
-	    else
-		return tiny*copysign(tiny,x); 	/*underflow*/
-	}
-        k += 128;				/* subnormal result */
-	u.bits.exp = k;
-        return u.e*0x1p-128;
-}
 
+	if (n > 16383) {
+		x *= 0x1p16383L;
+		n -= 16383;
+		if (n > 16383) {
+			x *= 0x1p16383L;
+			n -= 16383;
+			if (n > 16383)
+				n = 16383;
+		}
+	} else if (n < -16382) {
+		x *= 0x1p-16382L * 0x1p113L;
+		n += 16382 - 113;
+		if (n < -16382) {
+			x *= 0x1p-16382L * 0x1p113L;
+			n += 16382 - 113;
+			if (n < -16382)
+				n = -16382;
+		}
+	}
+	u.e = 1.0;
+	u.xbits.expsign = 0x3fff + n;
+	return x * u.e;
+}
 __strong_reference(scalbnl, ldexpl);
+#endif
+
diff --git a/libm/upstream-freebsd/lib/msun/src/s_sincosl.c b/libm/upstream-freebsd/lib/msun/src/s_sincosl.c
index aef36c2..3dd3457 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_sincosl.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_sincosl.c
@@ -50,11 +50,10 @@
 sincosl(long double x, long double *sn, long double *cs)
 {
 	union IEEEl2bits z;
-	int e0, sgn;
+	int e0;
 	long double y[2];
 
 	z.e = x;
-	sgn = z.bits.sign;
 	z.bits.sign = 0;
 
 	ENTERV();
diff --git a/libm/upstream-freebsd/lib/msun/src/s_sinpi.c b/libm/upstream-freebsd/lib/msun/src/s_sinpi.c
new file mode 100644
index 0000000..858459a
--- /dev/null
+++ b/libm/upstream-freebsd/lib/msun/src/s_sinpi.c
@@ -0,0 +1,169 @@
+/*-
+ * Copyright (c) 2017 Steven G. Kargl
+ * 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 unmodified, 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 ``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 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.
+ */
+
+/**
+ * sinpi(x) computes sin(pi*x) without multiplication by pi (almost).  First,
+ * note that sinpi(-x) = -sinpi(x), so the algorithm considers only |x| and
+ * includes reflection symmetry by considering the sign of x on output.  The
+ * method used depends on the magnitude of x.
+ *
+ * 1. For small |x|, sinpi(x) = pi * x where a sloppy threshold is used.  The
+ *    threshold is |x| < 0x1pN with N = -(P/2+M).  P is the precision of the
+ *    floating-point type and M = 2 to 4.  To achieve high accuracy, pi is 
+ *    decomposed into high and low parts with the high part containing a
+ *    number of trailing zero bits.  x is also split into high and low parts.
+ *
+ * 2. For |x| < 1, argument reduction is not required and sinpi(x) is 
+ *    computed by calling a kernel that leverages the kernels for sin(x)
+ *    ans cos(x).  See k_sinpi.c and k_cospi.c for details.
+ *
+ * 3. For 1 <= |x| < 0x1p(P-1), argument reduction is required where
+ *    |x| = j0 + r with j0 an integer and the remainder r satisfies
+ *    0 <= r < 1.  With the given domain, a simplified inline floor(x)
+ *    is used.  Also, note the following identity
+ *
+ *    sinpi(x) = sin(pi*(j0+r))
+ *             = sin(pi*j0) * cos(pi*r) + cos(pi*j0) * sin(pi*r)
+ *             = cos(pi*j0) * sin(pi*r)
+ *             = +-sinpi(r)
+ *
+ *    If j0 is even, then cos(pi*j0) = 1. If j0 is odd, then cos(pi*j0) = -1.
+ *    sinpi(r) is then computed via an appropriate kernel.
+ *
+ * 4. For |x| >= 0x1p(P-1), |x| is integral and sinpi(x) = copysign(0,x).
+ *
+ * 5. Special cases:
+ *
+ *    sinpi(+-0) = +-0
+ *    sinpi(+-n) = +-0, for positive integers n.
+ *    sinpi(+-inf) = nan.  Raises the "invalid" floating-point exception.
+ *    sinpi(nan) = nan.  Raises the "invalid" floating-point exception.
+ */
+
+#include <float.h>
+#include "math.h"
+#include "math_private.h"
+
+static const double
+pi_hi = 3.1415926814079285e+00,	/* 0x400921fb 0x58000000 */
+pi_lo =-2.7818135228334233e-08;	/* 0xbe5dde97 0x3dcb3b3a */
+
+#include "k_cospi.h"
+#include "k_sinpi.h"
+
+volatile static const double vzero = 0;
+
+double
+sinpi(double x)
+{
+	double ax, hi, lo, s;
+	uint32_t hx, ix, j0, lx;
+
+	EXTRACT_WORDS(hx, lx, x);
+	ix = hx & 0x7fffffff;
+	INSERT_WORDS(ax, ix, lx);
+
+	if (ix < 0x3ff00000) {			/* |x| < 1 */
+		if (ix < 0x3fd00000) {		/* |x| < 0.25 */
+			if (ix < 0x3e200000) {	/* |x| < 0x1p-29 */
+				if (x == 0)
+					return (x);
+				/*
+				 * To avoid issues with subnormal values,
+				 * scale the computation and rescale on 
+				 * return.
+				 */
+				INSERT_WORDS(hi, hx, 0);
+				hi *= 0x1p53;
+				lo = x * 0x1p53 - hi;
+				s = (pi_lo + pi_hi) * lo + pi_lo * hi +
+				    pi_hi * hi;
+				return (s * 0x1p-53);
+			}
+
+			s = __kernel_sinpi(ax);
+			return ((hx & 0x80000000) ? -s : s);
+		}
+
+		if (ix < 0x3fe00000)		/* |x| < 0.5 */
+			s = __kernel_cospi(0.5 - ax);
+		else if (ix < 0x3fe80000)	/* |x| < 0.75 */
+			s = __kernel_cospi(ax - 0.5);
+		else
+			s = __kernel_sinpi(1 - ax);
+		return ((hx & 0x80000000) ? -s : s);
+	}
+
+	if (ix < 0x43300000) {			/* 1 <= |x| < 0x1p52 */
+		/* Determine integer part of ax. */
+		j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
+		if (j0 < 20) {
+			ix &= ~(0x000fffff >> j0);
+			lx = 0;
+		} else {
+			lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
+		}
+		INSERT_WORDS(x, ix, lx);
+
+		ax -= x;
+		EXTRACT_WORDS(ix, lx, ax);
+
+		if (ix == 0)
+			s = 0;
+		else {
+			if (ix < 0x3fe00000) {		/* |x| < 0.5 */
+				if (ix < 0x3fd00000)	/* |x| < 0.25 */
+					s = __kernel_sinpi(ax);
+				else 
+					s = __kernel_cospi(0.5 - ax);
+			} else {
+				if (ix < 0x3fe80000)	/* |x| < 0.75 */
+					s = __kernel_cospi(ax - 0.5);
+				else
+					s = __kernel_sinpi(1 - ax);
+			}
+
+			if (j0 > 30)
+				x -= 0x1p30;
+			j0 = (uint32_t)x;
+			if (j0 & 1) s = -s;
+		}
+
+		return ((hx & 0x80000000) ? -s : s);
+	}
+
+	if (ix >= 0x7f800000)
+		return (vzero / vzero);
+
+	/*
+	 * |x| >= 0x1p52 is always an integer, so return +-0.
+	 */
+	return (copysign(0, x));
+}
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(sinpi, sinpil);
+#endif