MFC r257770 r257818 r257823 r260066 r260067 r260089 r260145 r268587 r268588

    r268589 r268590 r268593 r268597 r269758 r270845 r270847 r270893 r270932
    r270947 r271147

Merge libm work by kargl, bde and das from the past few months.
Besides optimisations and small bug fixes this includes new implementations
for C99 functions expl, coshl, sinhl, tanhl, erfl and erfcl.

Approved by:	re (kib)

17 files changed, 669 insertions, 218 deletions

diff --git a/lib/msun/src/e_cosh.c b/lib/msun/src/e_cosh.c
index a363695..246b5fb 100644
--- a/lib/msun/src/e_cosh.c
+++ b/lib/msun/src/e_cosh.c
@@ -35,6 +35,8 @@ __FBSDID("$FreeBSD$");
  *	only cosh(0)=1 is exact for finite x.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -77,3 +79,7 @@ __ieee754_cosh(double x)
     /* |x| > overflowthresold, cosh(x) overflow */
 	return huge*huge;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(cosh, coshl);
+#endif
diff --git a/lib/msun/src/e_coshl.c b/lib/msun/src/e_coshl.c
new file mode 100644
index 0000000..0a21277
--- /dev/null
+++ b/lib/msun/src/e_coshl.c
@@ -0,0 +1,130 @@
+/* from: FreeBSD: head/lib/msun/src/e_coshl.c XXX */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_cosh.c for complete comments.
+ *
+ * Converted to long double by Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+#include "k_expl.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+static const volatile long double huge = 0x1p10000L, tiny = 0x1p-10000L;
+#if LDBL_MANT_DIG == 64
+/*
+ * Domain [-1, 1], range ~[-1.8211e-21, 1.8211e-21]:
+ * |cosh(x) - c(x)| < 2**-68.8
+ */
+static const union IEEEl2bits
+C4u = LD80C(0xaaaaaaaaaaaaac78, -5,  4.16666666666666682297e-2L);
+#define	C4	C4u.e
+static const double
+C2  =  0.5,
+C6  =  1.3888888888888616e-3,		/*  0x16c16c16c16b99.0p-62 */
+C8  =  2.4801587301767953e-5,		/*  0x1a01a01a027061.0p-68 */
+C10 =  2.7557319163300398e-7,		/*  0x127e4fb6c9b55f.0p-74 */
+C12 =  2.0876768371393075e-9,		/*  0x11eed99406a3f4.0p-81 */
+C14 =  1.1469537039374480e-11,		/*  0x1938c67cd18c48.0p-89 */
+C16 =  4.8473490896852041e-14;		/*  0x1b49c429701e45.0p-97 */
+#elif LDBL_MANT_DIG == 113
+/*
+ * Domain [-1, 1], range ~[-2.3194e-37, 2.3194e-37]:
+ * |cosh(x) - c(x)| < 2**-121.69
+ */
+static const long double
+C4  =  4.16666666666666666666666666666666225e-2L,	/*  0x1555555555555555555555555554e.0p-117L */
+C6  =  1.38888888888888888888888888889434831e-3L,	/*  0x16c16c16c16c16c16c16c16c1dd7a.0p-122L */
+C8  =  2.48015873015873015873015871870962089e-5L,	/*  0x1a01a01a01a01a01a01a017af2756.0p-128L */
+C10 =  2.75573192239858906525574318600800201e-7L,	/*  0x127e4fb7789f5c72ef01c8a040640.0p-134L */
+C12 =  2.08767569878680989791444691755468269e-9L,	/*  0x11eed8eff8d897b543d0679607399.0p-141L */
+C14=  1.14707455977297247387801189650495351e-11L,	/*  0x193974a8c07c9d24ae169a7fa9b54.0p-149L */
+C16 =  4.77947733238737883626416876486279985e-14L;	/*  0x1ae7f3e733b814d4e1b90f5727fe4.0p-157L */
+static const double
+C2  =  0.5,
+C18 =  1.5619206968597871e-16,		/*  0x16827863b9900b.0p-105 */
+C20 =  4.1103176218528049e-19,		/*  0x1e542ba3d3c269.0p-114 */
+C22 =  8.8967926401641701e-22,		/*  0x10ce399542a014.0p-122 */
+C24 =  1.6116681626523904e-24,		/*  0x1f2c981d1f0cb7.0p-132 */
+C26 =  2.5022374732804632e-27;		/*  0x18c7ecf8b2c4a0.0p-141 */
+#else
+#error "Unsupported long double format"
+#endif /* LDBL_MANT_DIG == 64 */
+
+/* log(2**16385 - 0.5) rounded up: */
+static const float
+o_threshold =  1.13572168e4;		/*  0xb174de.0p-10 */
+
+long double
+coshl(long double x)
+{
+	long double hi,lo,x2,x4;
+	double dx2;
+	uint16_t ix;
+
+	GET_LDBL_EXPSIGN(ix,x);
+	ix &= 0x7fff;
+
+    /* x is INF or NaN */
+	if(ix>=0x7fff) return x*x;
+
+	ENTERI();
+
+    /* |x| < 1, return 1 or c(x) */
+	if(ix<0x3fff) {
+	    if (ix<BIAS-(LDBL_MANT_DIG+1)/2) 	/* |x| < TINY */
+		RETURNI(1+tiny);	/* cosh(tiny) = 1(+) with inexact */
+	    x2 = x*x;
+#if LDBL_MANT_DIG == 64
+	    x4 = x2*x2;
+	    RETURNI(((C16*x2 + C14)*x4 + (C12*x2 + C10))*(x4*x4*x2) +
+		((C8*x2 + C6)*x2 + C4)*x4 + C2*x2 + 1);
+#elif LDBL_MANT_DIG == 113
+	    dx2 = x2;
+	    RETURNI((((((((((((C26*dx2 + C24)*dx2 + C22)*dx2 +
+		C20)*x2 + C18)*x2 +
+		C16)*x2 + C14)*x2 + C12)*x2 + C10)*x2 + C8)*x2 + C6)*x2 +
+		C4)*(x2*x2) + C2*x2 + 1);
+#endif
+	}
+
+    /* |x| in [1, 64), return accurate exp(|x|)/2+1/exp(|x|)/2 */
+	if (ix < 0x4005) {
+	    k_hexpl(fabsl(x), &hi, &lo);
+	    RETURNI(lo + 0.25/(hi + lo) + hi);
+	}
+
+    /* |x| in [64, o_threshold], return correctly-overflowing exp(|x|)/2 */
+	if (fabsl(x) <= o_threshold)
+	    RETURNI(hexpl(fabsl(x)));
+
+    /* |x| > o_threshold, cosh(x) overflow */
+	RETURNI(huge*huge);
+}
diff --git a/lib/msun/src/e_lgamma_r.c b/lib/msun/src/e_lgamma_r.c
index 1cff592..7a95ea4 100644
--- a/lib/msun/src/e_lgamma_r.c
+++ b/lib/msun/src/e_lgamma_r.c
@@ -86,8 +86,10 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+static const volatile double vzero = 0;
+
 static const double
-two52=  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
+zero=  0.00000000000000000000e+00,
 half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 pi  =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
@@ -154,39 +156,35 @@ w4  = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
 w5  =  8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
 w6  = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
 
-static const double zero=  0.00000000000000000000e+00;
-
-	static double sin_pi(double x)
+/*
+ * Compute sin(pi*x) without actually doing the pi*x multiplication.
+ * sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is
+ * the precision of x.
+ */
+static double
+sin_pi(double x)
 {
+	volatile double vz;
 	double y,z;
-	int n,ix;
+	int n;
+
+	y = -x;
 
-	GET_HIGH_WORD(ix,x);
-	ix &= 0x7fffffff;
+	vz = y+0x1p52;			/* depend on 0 <= y < 0x1p52 */
+	z = vz-0x1p52;			/* rint(y) for the above range */
+	if (z == y)
+	    return zero;
 
-	if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0);
-	y = -x;		/* x is assume negative */
+	vz = y+0x1p50;
+	GET_LOW_WORD(n,vz);		/* bits for rounded y (units 0.25) */
+	z = vz-0x1p50;			/* y rounded to a multiple of 0.25 */
+	if (z > y) {
+	    z -= 0.25;			/* adjust to round down */
+	    n--;
+	}
+	n &= 7;				/* octant of y mod 2 */
+	y = y - z + n * 0.25;		/* y mod 2 */
 
-    /*
-     * argument reduction, make sure inexact flag not raised if input
-     * is an integer
-     */
-	z = floor(y);
-	if(z!=y) {				/* inexact anyway */
-	    y  *= 0.5;
-	    y   = 2.0*(y - floor(y));		/* y = |x| mod 2.0 */
-	    n   = (int) (y*4.0);
-	} else {
-            if(ix>=0x43400000) {
-                y = zero; n = 0;                 /* y must be even */
-            } else {
-                if(ix<0x43300000) z = y+two52;	/* exact */
-		GET_LOW_WORD(n,z);
-		n &= 1;
-                y  = n;
-                n<<= 2;
-            }
-        }
 	switch (n) {
 	    case 0:   y =  __kernel_sin(pi*y,zero,0); break;
 	    case 1:   
@@ -206,7 +204,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
 {
 	double t,y,z,nadj,p,p1,p2,p3,q,r,w;
 	int32_t hx;
-	int i,lx,ix;
+	int i,ix,lx;
 
 	EXTRACT_WORDS(hx,lx,x);
 
@@ -214,7 +212,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
 	if(ix>=0x7ff00000) return x*x;
-	if((ix|lx)==0) return one/zero;
+	if((ix|lx)==0) return one/vzero;
 	if(ix<0x3b900000) {	/* |x|<2**-70, return -log(|x|) */
 	    if(hx<0) {
 	        *signgamp = -1;
@@ -223,9 +221,9 @@ __ieee754_lgamma_r(double x, int *signgamp)
 	}
 	if(hx<0) {
 	    if(ix>=0x43300000) 	/* |x|>=2**52, must be -integer */
-		return one/zero;
+		return one/vzero;
 	    t = sin_pi(x);
-	    if(t==zero) return one/zero; /* -integer */
+	    if(t==zero) return one/vzero; /* -integer */
 	    nadj = __ieee754_log(pi/fabs(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
diff --git a/lib/msun/src/e_lgammaf_r.c b/lib/msun/src/e_lgammaf_r.c
index e2d90ef..9a7ab39 100644
--- a/lib/msun/src/e_lgammaf_r.c
+++ b/lib/msun/src/e_lgammaf_r.c
@@ -19,8 +19,10 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+static const volatile float vzero = 0;
+
 static const float
-two23=  8.3886080000e+06, /* 0x4b000000 */
+zero=  0.0000000000e+00,
 half=  5.0000000000e-01, /* 0x3f000000 */
 one =  1.0000000000e+00, /* 0x3f800000 */
 pi  =  3.1415927410e+00, /* 0x40490fdb */
@@ -87,39 +89,30 @@ w4  = -5.9518753551e-04, /* 0xba1c065c */
 w5  =  8.3633989561e-04, /* 0x3a5b3dd2 */
 w6  = -1.6309292987e-03; /* 0xbad5c4e8 */
 
-static const float zero=  0.0000000000e+00;
-
-	static float sin_pif(float x)
+static float
+sin_pif(float x)
 {
+	volatile float vz;
 	float y,z;
-	int n,ix;
+	int n;
 
-	GET_FLOAT_WORD(ix,x);
-	ix &= 0x7fffffff;
+	y = -x;
 
-	if(ix<0x3e800000) return __kernel_sindf(pi*x);
-	y = -x;		/* x is assume negative */
+	vz = y+0x1p23F;			/* depend on 0 <= y < 0x1p23 */
+	z = vz-0x1p23F;			/* rintf(y) for the above range */
+	if (z == y)
+	    return zero;
+
+	vz = y+0x1p21F;
+	GET_FLOAT_WORD(n,vz);		/* bits for rounded y (units 0.25) */
+	z = vz-0x1p21F;			/* y rounded to a multiple of 0.25 */
+	if (z > y) {
+	    z -= 0.25F;			/* adjust to round down */
+	    n--;
+	}
+	n &= 7;				/* octant of y mod 2 */
+	y = y - z + n * 0.25F;		/* y mod 2 */
 
-    /*
-     * argument reduction, make sure inexact flag not raised if input
-     * is an integer
-     */
-	z = floorf(y);
-	if(z!=y) {				/* inexact anyway */
-	    y  *= (float)0.5;
-	    y   = (float)2.0*(y - floorf(y));	/* y = |x| mod 2.0 */
-	    n   = (int) (y*(float)4.0);
-	} else {
-            if(ix>=0x4b800000) {
-                y = zero; n = 0;                 /* y must be even */
-            } else {
-                if(ix<0x4b000000) z = y+two23;	/* exact */
-		GET_FLOAT_WORD(n,z);
-		n &= 1;
-                y  = n;
-                n<<= 2;
-            }
-        }
 	switch (n) {
 	    case 0:   y =  __kernel_sindf(pi*y); break;
 	    case 1:
@@ -147,7 +140,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
 	if(ix>=0x7f800000) return x*x;
-	if(ix==0) return one/zero;
+	if(ix==0) return one/vzero;
 	if(ix<0x35000000) {	/* |x|<2**-21, return -log(|x|) */
 	    if(hx<0) {
 	        *signgamp = -1;
@@ -156,9 +149,9 @@ __ieee754_lgammaf_r(float x, int *signgamp)
 	}
 	if(hx<0) {
 	    if(ix>=0x4b000000) 	/* |x|>=2**23, must be -integer */
-		return one/zero;
+		return one/vzero;
 	    t = sin_pif(x);
-	    if(t==zero) return one/zero; /* -integer */
+	    if(t==zero) return one/vzero; /* -integer */
 	    nadj = __ieee754_logf(pi/fabsf(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
diff --git a/lib/msun/src/e_pow.c b/lib/msun/src/e_pow.c
index 7607a4a..d54af9d 100644
--- a/lib/msun/src/e_pow.c
+++ b/lib/msun/src/e_pow.c
@@ -19,20 +19,20 @@ __FBSDID("$FreeBSD$");
  *	1. Compute and return log2(x) in two pieces:
  *		log2(x) = w1 + w2,
  *	   where w1 has 53-24 = 29 bit trailing zeros.
- *	2. Perform y*log2(x) = n+y' by simulating muti-precision 
+ *	2. Perform y*log2(x) = n+y' by simulating multi-precision 
  *	   arithmetic, where |y'|<=0.5.
  *	3. Return x**y = 2**n*exp(y'*log2)
  *
  * Special cases:
  *	1.  (anything) ** 0  is 1
  *	2.  (anything) ** 1  is itself
- *	3.  (anything) ** NAN is NAN
+ *	3.  (anything) ** NAN is NAN except 1 ** NAN = 1
  *	4.  NAN ** (anything except 0) is NAN
  *	5.  +-(|x| > 1) **  +INF is +INF
  *	6.  +-(|x| > 1) **  -INF is +0
  *	7.  +-(|x| < 1) **  +INF is +0
  *	8.  +-(|x| < 1) **  -INF is +INF
- *	9.  +-1         ** +-INF is NAN
+ *	9.  +-1         ** +-INF is 1
  *	10. +0 ** (+anything except 0, NAN)               is +0
  *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
  *	12. +0 ** (-anything except 0, NAN)               is +INF
@@ -141,7 +141,7 @@ __ieee754_pow(double x, double y)
 	if(ly==0) { 	
 	    if (iy==0x7ff00000) {	/* y is +-inf */
 	        if(((ix-0x3ff00000)|lx)==0)
-		    return  one;	/* (-1)**+-inf is NaN */
+		    return  one;	/* (-1)**+-inf is 1 */
 	        else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
 		    return (hy>=0)? y: zero;
 	        else			/* (|x|<1)**-,+inf = inf,0 */
diff --git a/lib/msun/src/e_sinh.c b/lib/msun/src/e_sinh.c
index 17442d0..6c01f4a 100644
--- a/lib/msun/src/e_sinh.c
+++ b/lib/msun/src/e_sinh.c
@@ -32,6 +32,8 @@ __FBSDID("$FreeBSD$");
  *	only sinh(0)=0 is exact for finite x.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -71,3 +73,7 @@ __ieee754_sinh(double x)
     /* |x| > overflowthresold, sinh(x) overflow */
 	return x*shuge;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(sinh, sinhl);
+#endif
diff --git a/lib/msun/src/e_sinhl.c b/lib/msun/src/e_sinhl.c
new file mode 100644
index 0000000..ce7e333
--- /dev/null
+++ b/lib/msun/src/e_sinhl.c
@@ -0,0 +1,131 @@
+/* from: FreeBSD: head/lib/msun/src/e_sinhl.c XXX */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_sinh.c for complete comments.
+ *
+ * Converted to long double by Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+#include "k_expl.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+static const long double shuge = 0x1p16383L;
+#if LDBL_MANT_DIG == 64
+/*
+ * Domain [-1, 1], range ~[-6.6749e-22, 6.6749e-22]:
+ * |sinh(x)/x - s(x)| < 2**-70.3
+ */
+static const union IEEEl2bits
+S3u = LD80C(0xaaaaaaaaaaaaaaaa, -3,  1.66666666666666666658e-1L);
+#define	S3	S3u.e
+static const double
+S5  =  8.3333333333333332e-3,		/*  0x11111111111111.0p-59 */
+S7  =  1.9841269841270074e-4,		/*  0x1a01a01a01a070.0p-65 */
+S9  =  2.7557319223873889e-6,		/*  0x171de3a5565fe6.0p-71 */
+S11 =  2.5052108406704084e-8,		/*  0x1ae6456857530f.0p-78 */
+S13 =  1.6059042748655297e-10,		/*  0x161245fa910697.0p-85 */
+S15 =  7.6470006914396920e-13,		/*  0x1ae7ce4eff2792.0p-93 */
+S17 =  2.8346142308424267e-15;		/*  0x19882ce789ffc6.0p-101 */
+#elif LDBL_MANT_DIG == 113
+/*
+ * Domain [-1, 1], range ~[-2.9673e-36, 2.9673e-36]:
+ * |sinh(x)/x - s(x)| < 2**-118.0
+ */
+static const long double
+S3  =  1.66666666666666666666666666666666033e-1L,	/*  0x1555555555555555555555555553b.0p-115L */
+S5  =  8.33333333333333333333333333337643193e-3L,	/*  0x111111111111111111111111180f5.0p-119L */
+S7  =  1.98412698412698412698412697391263199e-4L,	/*  0x1a01a01a01a01a01a01a0176aad11.0p-125L */
+S9  =  2.75573192239858906525574406205464218e-6L,	/*  0x171de3a556c7338faac243aaa9592.0p-131L */
+S11 =  2.50521083854417187749675637460977997e-8L,	/*  0x1ae64567f544e38fe59b3380d7413.0p-138L */
+S13 =  1.60590438368216146368737762431552702e-10L,	/*  0x16124613a86d098059c7620850fc2.0p-145L */
+S15 =  7.64716373181980539786802470969096440e-13L,	/*  0x1ae7f3e733b814193af09ce723043.0p-153L */
+S17 =  2.81145725434775409870584280722701574e-15L;	/*  0x1952c77030c36898c3fd0b6dfc562.0p-161L */
+static const double
+S19=  8.2206352435411005e-18,		/*  0x12f49b4662b86d.0p-109 */
+S21=  1.9572943931418891e-20,		/*  0x171b8f2fab9628.0p-118 */
+S23 =  3.8679983530666939e-23,		/*  0x17617002b73afc.0p-127 */
+S25 =  6.5067867911512749e-26;		/*  0x1423352626048a.0p-136 */
+#else
+#error "Unsupported long double format"
+#endif /* LDBL_MANT_DIG == 64 */
+
+/* log(2**16385 - 0.5) rounded up: */
+static const float
+o_threshold =  1.13572168e4;		/*  0xb174de.0p-10 */
+
+long double
+sinhl(long double x)
+{
+	long double hi,lo,x2,x4;
+	double dx2,s;
+	int16_t ix,jx;
+
+	GET_LDBL_EXPSIGN(jx,x);
+	ix = jx&0x7fff;
+
+    /* x is INF or NaN */
+	if(ix>=0x7fff) return x+x;
+
+	ENTERI();
+
+	s = 1;
+	if (jx<0) s = -1;
+
+    /* |x| < 64, return x, s(x), or accurate s*(exp(|x|)/2-1/exp(|x|)/2) */
+	if (ix<0x4005) {		/* |x|<64 */
+	    if (ix<BIAS-(LDBL_MANT_DIG+1)/2) 	/* |x|<TINY */
+		if(shuge+x>1) RETURNI(x);  /* sinh(tiny) = tiny with inexact */
+	    if (ix<0x3fff) {		/* |x|<1 */
+		x2 = x*x;
+#if LDBL_MANT_DIG == 64
+		x4 = x2*x2;
+		RETURNI(((S17*x2 + S15)*x4 + (S13*x2 + S11))*(x2*x*x4*x4) +
+		    ((S9*x2 + S7)*x2 + S5)*(x2*x*x2) + S3*(x2*x) + x);
+#elif LDBL_MANT_DIG == 113
+		dx2 = x2;
+		RETURNI(((((((((((S25*dx2 + S23)*dx2 +
+		    S21)*x2 + S19)*x2 +
+		    S17)*x2 + S15)*x2 + S13)*x2 + S11)*x2 + S9)*x2 + S7)*x2 +
+		    S5)* (x2*x*x2) +
+		    S3*(x2*x) + x);
+#endif
+	    }
+	    k_hexpl(fabsl(x), &hi, &lo);
+	    RETURNI(s*(lo - 0.25/(hi + lo) + hi));
+	}
+
+    /* |x| in [64, o_threshold], return correctly-overflowing s*exp(|x|)/2 */
+	if (fabsl(x) <= o_threshold)
+	    RETURNI(s*hexpl(fabsl(x)));
+
+    /* |x| > o_threshold, sinh(x) overflow */
+	return x*shuge;
+}
diff --git a/lib/msun/src/imprecise.c b/lib/msun/src/imprecise.c
index a7503bf..92fb2d0 100644
--- a/lib/msun/src/imprecise.c
+++ b/lib/msun/src/imprecise.c
@@ -60,10 +60,5 @@ DECLARE_WEAK(powl);
 	long double imprecise_ ## f ## l(long double v) { return f(v); }\
 	DECLARE_WEAK(f ## l)
 
-DECLARE_IMPRECISE(cosh);
-DECLARE_IMPRECISE(erfc);
-DECLARE_IMPRECISE(erf);
 DECLARE_IMPRECISE(lgamma);
-DECLARE_IMPRECISE(sinh);
-DECLARE_IMPRECISE(tanh);
 DECLARE_IMPRECISE(tgamma);
diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h
index 1bd931c..3ab76f8 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -451,7 +451,10 @@ long double	atanl(long double);
 long double	cbrtl(long double);
 long double	ceill(long double);
 long double	copysignl(long double, long double) __pure2;
+long double	coshl(long double);
 long double	cosl(long double);
+long double	erfcl(long double);
+long double	erfl(long double);
 long double	exp2l(long double);
 long double	expl(long double);
 long double	expm1l(long double);
@@ -466,6 +469,7 @@ long double	frexpl(long double value, int *); /* fundamentally !__pure2 */
 long double	hypotl(long double, long double);
 int		ilogbl(long double) __pure2;
 long double	ldexpl(long double, int);
+long double	lgammal(long double);
 long long	llrintl(long double);
 long long	llroundl(long double);
 long double	log10l(long double);
@@ -482,45 +486,22 @@ long double	nextafterl(long double, long double);
 double		nexttoward(double, long double);
 float		nexttowardf(float, long double);
 long double	nexttowardl(long double, long double);
+long double	powl(long double, long double);
 long double	remainderl(long double, long double);
 long double	remquol(long double, long double, int *);
 long double	rintl(long double);
 long double	roundl(long double);
 long double	scalblnl(long double, long);
 long double	scalbnl(long double, int);
+long double	sinhl(long double);
 long double	sinl(long double);
 long double	sqrtl(long double);
+long double	tanhl(long double);
 long double	tanl(long double);
+long double	tgammal(long double);
 long double	truncl(long double);
 
 #endif /* __ISO_C_VISIBLE >= 1999 */
 __END_DECLS
 
 #endif /* !_MATH_H_ */
-
-/* separate header for cmath */
-#ifndef _MATH_EXTRA_H_
-#if __ISO_C_VISIBLE >= 1999
-#if _DECLARE_C99_LDBL_MATH
-
-#define _MATH_EXTRA_H_
-
-/*
- * extra long double versions of math functions for C99 and cmath
- */
-__BEGIN_DECLS
-
-long double	coshl(long double);
-long double	erfcl(long double);
-long double	erfl(long double);
-long double	lgammal(long double);
-long double	powl(long double, long double);
-long double	sinhl(long double);
-long double	tanhl(long double);
-long double	tgammal(long double);
-
-__END_DECLS
-
-#endif /* !_DECLARE_C99_LDBL_MATH */
-#endif /* __ISO_C_VISIBLE >= 1999 */
-#endif /* !_MATH_EXTRA_H_ */
diff --git a/lib/msun/src/s_erf.c b/lib/msun/src/s_erf.c
index 854767b..e1d63bc 100644
--- a/lib/msun/src/s_erf.c
+++ b/lib/msun/src/s_erf.c
@@ -111,18 +111,25 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile double tiny= 1e-300;
+
 static const double
-tiny	    = 1e-300,
-half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
-one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
-two =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
-	/* c = (float)0.84506291151 */
+half= 0.5,
+one = 1,
+two = 2,
+/* c = (float)0.84506291151 */
 erx =  8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
 /*
- * Coefficients for approximation to  erf on [0,0.84375]
+ * In the domain [0, 2**-28], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-84.
  */
 efx =  1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
 efx8=  1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
+/*
+ * Coefficients for approximation to erf on [0,0.84375]
+ */
 pp0  =  1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
 pp1  = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
 pp2  = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
@@ -134,7 +141,7 @@ qq3  =  5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
 qq4  =  1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
 qq5  = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
 /*
- * Coefficients for approximation to  erf  in [0.84375,1.25]
+ * Coefficients for approximation to erf in [0.84375,1.25]
  */
 pa0  = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
 pa1  =  4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
@@ -150,7 +157,7 @@ qa4  =  1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
 qa5  =  1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
 qa6  =  1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
 /*
- * Coefficients for approximation to  erfc in [1.25,1/0.35]
+ * Coefficients for approximation to erfc in [1.25,1/0.35]
  */
 ra0  = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
 ra1  = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
@@ -169,7 +176,7 @@ sa6  =  1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
 sa7  =  6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
 sa8  = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
 /*
- * Coefficients for approximation to  erfc in [1/.35,28]
+ * Coefficients for approximation to erfc in [1/.35,28]
  */
 rb0  = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
 rb1  = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
@@ -222,15 +229,12 @@ erf(double x)
 	x = fabs(x);
  	s = one/(x*x);
 	if(ix< 0x4006DB6E) {	/* |x| < 1/0.35 */
-	    R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
-				ra5+s*(ra6+s*ra7))))));
-	    S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
-				sa5+s*(sa6+s*(sa7+s*sa8)))))));
+	    R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
+	    S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		s*sa8)))))));
 	} else {	/* |x| >= 1/0.35 */
-	    R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
-				rb5+s*rb6)))));
-	    S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
-				sb5+s*(sb6+s*sb7))))));
+	    R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
+	    S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
 	}
 	z  = x;
 	SET_LOW_WORD(z,0);
@@ -238,6 +242,10 @@ erf(double x)
 	if(hx>=0) return one-r/x; else return  r/x-one;
 }
 
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(erf, erfl);
+#endif
+
 double
 erfc(double x)
 {
@@ -279,23 +287,23 @@ erfc(double x)
 	    x = fabs(x);
  	    s = one/(x*x);
 	    if(ix< 0x4006DB6D) {	/* |x| < 1/.35 ~ 2.857143*/
-	        R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
-				ra5+s*(ra6+s*ra7))))));
-	        S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
-				sa5+s*(sa6+s*(sa7+s*sa8)))))));
+		R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
+		S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		    s*sa8)))))));
 	    } else {			/* |x| >= 1/.35 ~ 2.857143 */
 		if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
-	        R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
-				rb5+s*rb6)))));
-	        S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
-				sb5+s*(sb6+s*sb7))))));
+		R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
+		S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
 	    }
 	    z  = x;
 	    SET_LOW_WORD(z,0);
-	    r  =  __ieee754_exp(-z*z-0.5625)*
-			__ieee754_exp((z-x)*(z+x)+R/S);
+	    r  =  __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
 	    if(hx>0) return r/x; else return two-r/x;
 	} else {
 	    if(hx>0) return tiny*tiny; else return two-tiny;
 	}
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(erfc, erfcl);
+#endif
diff --git a/lib/msun/src/s_erff.c b/lib/msun/src/s_erff.c
index b97ca1d..d6cfbd2 100644
--- a/lib/msun/src/s_erff.c
+++ b/lib/msun/src/s_erff.c
@@ -19,64 +19,63 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile float tiny = 1e-30;
+
 static const float
-tiny	    = 1e-30,
-half=  5.0000000000e-01, /* 0x3F000000 */
-one =  1.0000000000e+00, /* 0x3F800000 */
-two =  2.0000000000e+00, /* 0x40000000 */
+half= 0.5,
+one = 1,
+two = 2,
+erx = 8.42697144e-01,			/* 0x3f57bb00 */
 /*
- * Coefficients for approximation to erf on [0,0.84375]
+ * In the domain [0, 2**-14], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-42.
  */
-efx =  1.2837916613e-01, /* 0x3e0375d4 */
-efx8=  1.0270333290e+00, /* 0x3f8375d4 */
+efx = 1.28379166e-01, /* 0x3e0375d4 */
+efx8= 1.02703333e+00, /* 0x3f8375d4 */
 /*
- *  Domain [0, 0.84375], range ~[-5.4446e-10,5.5197e-10]:
- *  |(erf(x) - x)/x - p(x)/q(x)| < 2**-31.
+ * Domain [0, 0.84375], range ~[-5.4419e-10, 5.5179e-10]:
+ * |(erf(x) - x)/x - pp(x)/qq(x)| < 2**-31
  */
-pp0  =  1.28379166e-01F, /*  0x1.06eba8p-3 */
-pp1  = -3.36030394e-01F, /* -0x1.58185ap-2 */
-pp2  = -1.86260219e-03F, /* -0x1.e8451ep-10 */
-qq1  =  3.12324286e-01F, /*  0x1.3fd1f0p-2 */
-qq2  =  2.16070302e-02F, /*  0x1.620274p-6 */
-qq3  = -1.98859419e-03F, /* -0x1.04a626p-9 */
+pp0  =  1.28379166e-01, /* 0x3e0375d4 */
+pp1  = -3.36030394e-01, /* 0xbeac0c2d */
+pp2  = -1.86261395e-03, /* 0xbaf422f4 */
+qq1  =  3.12324315e-01, /* 0x3e9fe8f9 */
+qq2  =  2.16070414e-02, /* 0x3cb10140 */
+qq3  = -1.98859372e-03, /* 0xbb025311 */
 /*
- * Domain [0.84375, 1.25], range ~[-1.953e-11,1.940e-11]:
- * |(erf(x) - erx) - p(x)/q(x)| < 2**-36.
+ * Domain [0.84375, 1.25], range ~[-1.023e-9, 1.023e-9]:
+ * |(erf(x) - erx) - pa(x)/qa(x)| < 2**-31
  */
-erx  =  8.42697144e-01F, /*  0x1.af7600p-1.  erf(1) rounded to 16 bits. */
-pa0  =  3.64939137e-06F, /*  0x1.e9d022p-19 */
-pa1  =  4.15109694e-01F, /*  0x1.a91284p-2 */
-pa2  = -1.65179938e-01F, /* -0x1.5249dcp-3 */
-pa3  =  1.10914491e-01F, /*  0x1.c64e46p-4 */
-qa1  =  6.02074385e-01F, /*  0x1.344318p-1 */
-qa2  =  5.35934687e-01F, /*  0x1.126608p-1 */
-qa3  =  1.68576106e-01F, /*  0x1.593e6ep-3 */
-qa4  =  5.62181212e-02F, /*  0x1.cc89f2p-5 */
+pa0  =  3.65041046e-06, /* 0x3674f993 */
+pa1  =  4.15109307e-01, /* 0x3ed48935 */
+pa2  = -2.09395722e-01, /* 0xbe566bd5 */
+pa3  =  8.67677554e-02, /* 0x3db1b34b */
+qa1  =  4.95560974e-01, /* 0x3efdba2b */
+qa2  =  3.71248513e-01, /* 0x3ebe1449 */
+qa3  =  3.92478965e-02, /* 0x3d20c267 */
 /*
- * Domain [1.25,1/0.35], range ~[-7.043e-10,7.457e-10]:
- * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-30
+ * Domain [1.25,1/0.35], range ~[-4.821e-9, 4.927e-9]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - ra(x)/sa(x)| < 2**-28
  */
-ra0  = -9.87132732e-03F, /* -0x1.4376b2p-7 */
-ra1  = -5.53605914e-01F, /* -0x1.1b723cp-1 */
-ra2  = -2.17589188e+00F, /* -0x1.1683a0p+1 */
-ra3  = -1.43268085e+00F, /* -0x1.6ec42cp+0 */
-sa1  =  5.45995426e+00F, /*  0x1.5d6fe4p+2 */
-sa2  =  6.69798088e+00F, /*  0x1.acabb8p+2 */
-sa3  =  1.43113089e+00F, /*  0x1.6e5e98p+0 */
-sa4  = -5.77397496e-02F, /* -0x1.d90108p-5 */
+ra0  = -9.88156721e-03, /* 0xbc21e64c */
+ra1  = -5.43658376e-01, /* 0xbf0b2d32 */
+ra2  = -1.66828310e+00, /* 0xbfd58a4d */
+ra3  = -6.91554189e-01, /* 0xbf3109b2 */
+sa1  =  4.48581553e+00, /* 0x408f8bcd */
+sa2  =  4.10799170e+00, /* 0x408374ab */
+sa3  =  5.53855181e-01, /* 0x3f0dc974 */
 /*
- * Domain [1/0.35, 11], range ~[-2.264e-13,2.336e-13]:
- * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-42
+ * Domain [2.85715, 11], range ~[-1.484e-9, 1.505e-9]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rb(x)/sb(x)| < 2**-30
  */
-rb0  = -9.86494310e-03F, /* -0x1.434124p-7 */
-rb1  = -6.25171244e-01F, /* -0x1.401672p-1 */
-rb2  = -6.16498327e+00F, /* -0x1.8a8f16p+2 */
-rb3  = -1.66696873e+01F, /* -0x1.0ab70ap+4 */
-rb4  = -9.53764343e+00F, /* -0x1.313460p+3 */
-sb1  =  1.26884899e+01F, /*  0x1.96081cp+3 */
-sb2  =  4.51839523e+01F, /*  0x1.6978bcp+5 */
-sb3  =  4.72810211e+01F, /*  0x1.7a3f88p+5 */
-sb4  =  8.93033314e+00F; /*  0x1.1dc54ap+3 */
+rb0  = -9.86496918e-03, /* 0xbc21a0ae */
+rb1  = -5.48049808e-01, /* 0xbf0c4cfe */
+rb2  = -1.84115684e+00, /* 0xbfebab07 */
+sb1  =  4.87132740e+00, /* 0x409be1ea */
+sb2  =  3.04982710e+00, /* 0x4043305e */
+sb3  = -7.61900663e-01; /* 0xbf430bec */
 
 float
 erff(float x)
@@ -85,9 +84,9 @@ erff(float x)
 	float R,S,P,Q,s,y,z,r;
 	GET_FLOAT_WORD(hx,x);
 	ix = hx&0x7fffffff;
-	if(ix>=0x7f800000) {		/* erf(nan)=nan */
+	if(ix>=0x7f800000) {		/* erff(nan)=nan */
 	    i = ((u_int32_t)hx>>31)<<1;
-	    return (float)(1-i)+one/x;	/* erf(+-inf)=+-1 */
+	    return (float)(1-i)+one/x;	/* erff(+-inf)=+-1 */
 	}
 
 	if(ix < 0x3f580000) {		/* |x|<0.84375 */
@@ -105,7 +104,7 @@ erff(float x)
 	if(ix < 0x3fa00000) {		/* 0.84375 <= |x| < 1.25 */
 	    s = fabsf(x)-one;
 	    P = pa0+s*(pa1+s*(pa2+s*pa3));
-	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*qa4)));
+	    Q = one+s*(qa1+s*(qa2+s*qa3));
 	    if(hx>=0) return erx + P/Q; else return -erx - P/Q;
 	}
 	if (ix >= 0x40800000) {		/* inf>|x|>=4 */
@@ -113,12 +112,12 @@ erff(float x)
 	}
 	x = fabsf(x);
  	s = one/(x*x);
-	if(ix< 0x4036DB6E) {	/* |x| < 1/0.35 */
+	if(ix< 0x4036db8c) {	/* |x| < 2.85715 ~ 1/0.35 */
 	    R=ra0+s*(ra1+s*(ra2+s*ra3));
-	    S=one+s*(sa1+s*(sa2+s*(sa3+s*sa4)));
-	} else {	/* |x| >= 1/0.35 */
-	    R=rb0+s*(rb1+s*(rb2+s*(rb3+s*rb4)));
-	    S=one+s*(sb1+s*(sb2+s*(sb3+s*sb4)));
+	    S=one+s*(sa1+s*(sa2+s*sa3));
+	} else {	/* |x| >= 2.85715 ~ 1/0.35 */
+	    R=rb0+s*(rb1+s*rb2);
+	    S=one+s*(sb1+s*(sb2+s*sb3));
 	}
 	SET_FLOAT_WORD(z,hx&0xffffe000);
 	r  = expf(-z*z-0.5625F)*expf((z-x)*(z+x)+R/S);
@@ -132,8 +131,8 @@ erfcf(float x)
 	float R,S,P,Q,s,y,z,r;
 	GET_FLOAT_WORD(hx,x);
 	ix = hx&0x7fffffff;
-	if(ix>=0x7f800000) {			/* erfc(nan)=nan */
-						/* erfc(+-inf)=0,2 */
+	if(ix>=0x7f800000) {			/* erfcf(nan)=nan */
+						/* erfcf(+-inf)=0,2 */
 	    return (float)(((u_int32_t)hx>>31)<<1)+one/x;
 	}
 
@@ -155,7 +154,7 @@ erfcf(float x)
 	if(ix < 0x3fa00000) {		/* 0.84375 <= |x| < 1.25 */
 	    s = fabsf(x)-one;
 	    P = pa0+s*(pa1+s*(pa2+s*pa3));
-	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*qa4)));
+	    Q = one+s*(qa1+s*(qa2+s*qa3));
 	    if(hx>=0) {
 	        z  = one-erx; return z - P/Q;
 	    } else {
@@ -165,13 +164,13 @@ erfcf(float x)
 	if (ix < 0x41300000) {		/* |x|<11 */
 	    x = fabsf(x);
  	    s = one/(x*x);
-	    if(ix< 0x4036DB6D) {	/* |x| < 1/.35 ~ 2.857143*/
-	        R=ra0+s*(ra1+s*(ra2+s*ra3));
-	        S=one+s*(sa1+s*(sa2+s*(sa3+s*sa4)));
-	    } else {			/* |x| >= 1/.35 ~ 2.857143 */
+	    if(ix< 0x4036db8c) {	/* |x| < 2.85715 ~ 1/.35 */
+		R=ra0+s*(ra1+s*(ra2+s*ra3));
+		S=one+s*(sa1+s*(sa2+s*sa3));
+	    } else {			/* |x| >= 2.85715 ~ 1/.35 */
 		if(hx<0&&ix>=0x40a00000) return two-tiny;/* x < -5 */
-	        R=rb0+s*(rb1+s*(rb2+s*(rb3+s*rb4)));
-		S=one+s*(sb1+s*(sb2+s*(sb3+s*sb4)));
+		R=rb0+s*(rb1+s*rb2);
+		S=one+s*(sb1+s*(sb2+s*sb3));
 	    }
 	    SET_FLOAT_WORD(z,hx&0xffffe000);
 	    r  = expf(-z*z-0.5625F)*expf((z-x)*(z+x)+R/S);
diff --git a/lib/msun/src/s_round.c b/lib/msun/src/s_round.c
index 65de31b..fab3019 100644
--- a/lib/msun/src/s_round.c
+++ b/lib/msun/src/s_round.c
@@ -27,25 +27,34 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-#include <math.h>
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
 
 double
 round(double x)
 {
 	double t;
+	uint32_t hx;
 
-	if (!isfinite(x))
-		return (x);
+	GET_HIGH_WORD(hx, x);
+	if ((hx & 0x7fffffff) == 0x7ff00000)
+		return (x + x);
 
-	if (x >= 0.0) {
+	if (!(hx & 0x80000000)) {
 		t = floor(x);
 		if (t - x <= -0.5)
-			t += 1.0;
+			t += 1;
 		return (t);
 	} else {
 		t = floor(-x);
 		if (t + x <= -0.5)
-			t += 1.0;
+			t += 1;
 		return (-t);
 	}
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(round, roundl);
+#endif
diff --git a/lib/msun/src/s_roundf.c b/lib/msun/src/s_roundf.c
index 952e8e7..e7e2eb9 100644
--- a/lib/msun/src/s_roundf.c
+++ b/lib/msun/src/s_roundf.c
@@ -27,25 +27,28 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-#include <math.h>
+#include "math.h"
+#include "math_private.h"
 
 float
 roundf(float x)
 {
 	float t;
+	uint32_t hx;
 
-	if (!isfinite(x))
-		return (x);
+	GET_FLOAT_WORD(hx, x);
+	if ((hx & 0x7fffffff) == 0x7f800000)
+		return (x + x);
 
-	if (x >= 0.0) {
+	if (!(hx & 0x80000000)) {
 		t = floorf(x);
-		if (t - x <= -0.5)
-			t += 1.0;
+		if (t - x <= -0.5F)
+			t += 1;
 		return (t);
 	} else {
 		t = floorf(-x);
-		if (t + x <= -0.5)
-			t += 1.0;
+		if (t + x <= -0.5F)
+			t += 1;
 		return (-t);
 	}
 }
diff --git a/lib/msun/src/s_roundl.c b/lib/msun/src/s_roundl.c
index a70b617..2d15e13 100644
--- a/lib/msun/src/s_roundl.c
+++ b/lib/msun/src/s_roundl.c
@@ -27,25 +27,36 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-#include <math.h>
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
 
 long double
 roundl(long double x)
 {
 	long double t;
+	uint16_t hx;
+
+	GET_LDBL_EXPSIGN(hx, x);
+	if ((hx & 0x7fff) == 0x7fff)
+		return (x + x);
 
-	if (!isfinite(x))
-		return (x);
+	ENTERI();
 
-	if (x >= 0.0) {
+	if (!(hx & 0x8000)) {
 		t = floorl(x);
-		if (t - x <= -0.5)
-			t += 1.0;
-		return (t);
+		if (t - x <= -0.5L)
+			t += 1;
+		RETURNI(t);
 	} else {
 		t = floorl(-x);
-		if (t + x <= -0.5)
-			t += 1.0;
-		return (-t);
+		if (t + x <= -0.5L)
+			t += 1;
+		RETURNI(-t);
 	}
 }
diff --git a/lib/msun/src/s_tanh.c b/lib/msun/src/s_tanh.c
index 96e3565..6d26c69 100644
--- a/lib/msun/src/s_tanh.c
+++ b/lib/msun/src/s_tanh.c
@@ -37,10 +37,13 @@ __FBSDID("$FreeBSD$");
  *	only tanh(0)=0 is exact for finite argument.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
-static const double one = 1.0, two = 2.0, tiny = 1.0e-300, huge = 1.0e300;
+static const volatile double tiny = 1.0e-300;
+static const double one = 1.0, two = 2.0, huge = 1.0e300;
 
 double
 tanh(double x)
@@ -75,3 +78,7 @@ tanh(double x)
 	}
 	return (jx>=0)? z: -z;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(tanh, tanhl);
+#endif
diff --git a/lib/msun/src/s_tanhf.c b/lib/msun/src/s_tanhf.c
index 04f09c6..f537be4 100644
--- a/lib/msun/src/s_tanhf.c
+++ b/lib/msun/src/s_tanhf.c
@@ -19,7 +19,9 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
-static const float one=1.0, two=2.0, tiny = 1.0e-30, huge = 1.0e30;
+static const volatile float tiny = 1.0e-30;
+static const float one=1.0, two=2.0, huge = 1.0e30;
+
 float
 tanhf(float x)
 {
diff --git a/lib/msun/src/s_tanhl.c b/lib/msun/src/s_tanhl.c
new file mode 100644
index 0000000..886158b
--- /dev/null
+++ b/lib/msun/src/s_tanhl.c
@@ -0,0 +1,172 @@
+/* from: FreeBSD: head/lib/msun/src/s_tanhl.c XXX */
+
+/* @(#)s_tanh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See s_tanh.c for complete comments.
+ *
+ * Converted to long double by Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "math.h"
+#include "math_private.h"
+#include "fpmath.h"
+#include "k_expl.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+static const volatile double tiny = 1.0e-300;
+static const double one = 1.0;
+#if LDBL_MANT_DIG == 64
+/*
+ * Domain [-0.25, 0.25], range ~[-1.6304e-22, 1.6304e-22]:
+ * |tanh(x)/x - t(x)| < 2**-72.3
+ */
+static const union IEEEl2bits
+T3u = LD80C(0xaaaaaaaaaaaaaa9f, -2, -3.33333333333333333017e-1L);
+#define	T3	T3u.e
+static const double
+T5  =  1.3333333333333314e-1,		/*  0x1111111111110a.0p-55 */
+T7  = -5.3968253968210485e-2,		/* -0x1ba1ba1ba1a1a1.0p-57 */
+T9  =  2.1869488531393817e-2,		/*  0x1664f488172022.0p-58 */
+T11 = -8.8632352345964591e-3,		/* -0x1226e34bc138d5.0p-59 */
+T13 =  3.5921169709993771e-3,		/*  0x1d6d371d3e400f.0p-61 */
+T15 = -1.4555786415756001e-3,		/* -0x17d923aa63814d.0p-62 */
+T17 =  5.8645267876296793e-4,		/*  0x13378589b85aa7.0p-63 */
+T19 = -2.1121033571392224e-4;		/* -0x1baf0af80c4090.0p-65 */
+#elif LDBL_MANT_DIG == 113
+/*
+ * Domain [-0.25, 0.25], range ~[-2.4211e-37, 2.4211e-37]:
+ * |tanh(x)/x - t(x)| < 2**121.6
+ */
+static const long double
+T3 = -3.33333333333333333333333333333332980e-1L,	/* -0x1555555555555555555555555554e.0p-114L */
+T5  =  1.33333333333333333333333333332707260e-1L,	/*  0x1111111111111111111111110ab7b.0p-115L */
+T7  = -5.39682539682539682539682535723482314e-2L,	/* -0x1ba1ba1ba1ba1ba1ba1ba17b5fc98.0p-117L */
+T9  =  2.18694885361552028218693591149061717e-2L,	/*  0x1664f4882c10f9f32d6b1a12a25e5.0p-118L */
+T11 = -8.86323552990219656883762347736381851e-3L,	/* -0x1226e355e6c23c8f5a5a0f386cb4d.0p-119L */
+T13 =  3.59212803657248101358314398220822722e-3L,	/*  0x1d6d3d0e157ddfb403ad3637442c6.0p-121L */
+T15 = -1.45583438705131796512568010348874662e-3L;	/* -0x17da36452b75e150c44cc34253b34.0p-122L */
+static const double
+T17 =  5.9002744094556621e-4,		/*  0x1355824803668e.0p-63 */
+T19 = -2.3912911424260516e-4,		/* -0x1f57d7734c8dde.0p-65 */
+T21 =  9.6915379535512898e-5,		/*  0x1967e18ad6a6ca.0p-66 */
+T23 = -3.9278322983156353e-5,		/* -0x1497d8e6b75729.0p-67 */
+T25 =  1.5918887220143869e-5,		/*  0x10b1319998cafa.0p-68 */
+T27 = -6.4514295231630956e-6,		/* -0x1b0f2b71b218eb.0p-70 */
+T29 =  2.6120754043964365e-6,		/*  0x15e963a3cf3a39.0p-71 */
+T31 = -1.0407567231003314e-6,		/* -0x1176041e656869.0p-72 */
+T33 =  3.4744117554063574e-7;		/*  0x1750fe732cab9c.0p-74 */
+#endif /* LDBL_MANT_DIG == 64 */
+
+static inline long double
+divl(long double a, long double b, long double c, long double d,
+    long double e, long double f)
+{
+	long double inv, r;
+	float fr, fw;
+
+	_2sumF(a, c);
+	b = b + c;
+	_2sumF(d, f);
+	e = e + f;
+
+	inv = 1 / (d + e);
+
+	r = (a + b) * inv;
+	fr = r;
+	r = fr;
+
+	fw = d + e;
+	e = d - fw + e;
+	d = fw;
+
+	r = r + (a - d * r + b - e * r) * inv;
+
+	return r;
+}
+
+long double
+tanhl(long double x)
+{
+	long double hi,lo,s,x2,x4,z;
+	double dx2;
+	int16_t jx,ix;
+
+	GET_LDBL_EXPSIGN(jx,x);
+	ix = jx&0x7fff;
+
+    /* x is INF or NaN */
+	if(ix>=0x7fff) {
+	    if (jx>=0) return one/x+one;    /* tanh(+-inf)=+-1 */
+	    else       return one/x-one;    /* tanh(NaN) = NaN */
+	}
+
+	ENTERI();
+
+    /* |x| < 40 */
+	if (ix < 0x4004 || fabsl(x) < 40) {	/* |x|<40 */
+	    if (__predict_false(ix<BIAS-(LDBL_MANT_DIG+1)/2)) {	/* |x|<TINY */
+		/* tanh(+-0) = +0; tanh(tiny) = tiny(-+) with inexact: */
+		return (x == 0 ? x : (0x1p200 * x - x) * 0x1p-200);
+	    }
+	    if (ix<0x3ffd) {		/* |x|<0.25 */
+		x2 = x*x;
+#if LDBL_MANT_DIG == 64
+		x4 = x2*x2;
+		RETURNI(((T19*x2 + T17)*x4 + (T15*x2 + T13))*(x2*x*x2*x4*x4) +
+		    ((T11*x2 + T9)*x4 + (T7*x2 + T5))*(x2*x*x2) +
+		    T3*(x2*x) + x);
+#elif LDBL_MANT_DIG == 113
+		dx2 = x2;
+#if 0
+		RETURNI(((((((((((((((T33*dx2 + T31)*dx2 + T29)*dx2 + T27)*dx2 +
+		    T25)*x2 + T23)*x2 + T21)*x2 + T19)*x2 + T17)*x2 +
+		    T15)*x2 + T13)*x2 + T11)*x2 + T9)*x2 + T7)*x2 + T5)*
+		    (x2*x*x2) +
+		    T3*(x2*x) + x);
+#else
+		long double q = ((((((((((((((T33*dx2 + T31)*dx2 + T29)*dx2 + T27)*dx2 +
+		    T25)*x2 + T23)*x2 + T21)*x2 + T19)*x2 + T17)*x2 +
+		    T15)*x2 + T13)*x2 + T11)*x2 + T9)*x2 + T7)*x2 + T5)*
+		    (x2*x*x2);
+		RETURNI(q + T3*(x2*x) + x);
+#endif
+#endif
+	    }
+	    k_hexpl(2*fabsl(x), &hi, &lo);
+	    if (ix<0x4001 && fabsl(x) < 1.5)	/* |x|<1.5 */
+		z = divl(hi, lo, -0.5, hi, lo, 0.5);
+	    else
+		z = one - one/(lo+0.5+hi);
+    /* |x| >= 40, return +-1 */
+	} else {
+	    z = one - tiny;		/* raise inexact flag */
+	}
+	s = 1;
+	if (jx<0) s = -1;
+	RETURNI(s*z);
+}