From a8a6a9161a123c760a51db32942722249315c1ba Mon Sep 17 00:00:00 2001
From: John Marino <draco@marino.st>
Date: Wed, 12 Jun 2013 00:28:04 +0200
Subject: [PATCH] libm: Add several new functions and symbol versioning

The following long double functions were added to the math library:
  logl
  log2l
  log10l
  log1pl
  expm1l
  acoshl
  asinhl
  atanhl

In addition, the FreeBSD functionality that creates symbol versioning
for libraries was adapted for FreeBSD.  The first version is called
"DFLY36.0".  If it is necessary to create a new version of the 3.5 or
3.6 branch, the number after the decimal will be incremented.  The 3.7
branch will start with "DFLY38.0" if it needs its own version.

libm was baselined with all symbols being the same version: DFLY36.0.
With symbol versioning, it will not be necessary to increment the major
version anymore, so this library shall always be known as libm.so.4 from
this point on.
---
 lib/libc/Versions.def        |  19 +
 lib/libm/Makefile            |  28 +-
 lib/libm/Symbol.map          | 260 +++++++++++++
 lib/libm/i386/Makefile.inc   |   2 +
 lib/libm/i386/Symbol.map     |  14 +
 lib/libm/ld80/s_expl.c       | 233 ++++++++++--
 lib/libm/ld80/s_logl.c       | 716 +++++++++++++++++++++++++++++++++++
 lib/libm/man/acosh.3         |  35 +-
 lib/libm/man/asinh.3         |  35 +-
 lib/libm/man/atanh.3         |  35 +-
 lib/libm/man/exp.3           |  12 +-
 lib/libm/man/log.3           |  43 ++-
 lib/libm/man/sin.3           |   2 +-
 lib/libm/src/catrig.c        |  30 +-
 lib/libm/src/catrigf.c       |  43 ++-
 lib/libm/src/e_acosh.c       |   9 +-
 lib/libm/src/e_acoshl.c      |  87 +++++
 lib/libm/src/e_atanh.c       |   8 +-
 lib/libm/src/e_atanhl.c      |  72 ++++
 lib/libm/src/e_log.c         |   8 +-
 lib/libm/src/e_log10.c       |   8 +-
 lib/libm/src/e_log2.c        |   8 +-
 lib/libm/src/math.h          |  18 +-
 lib/libm/src/math_private.h  | 296 ++++++++++++++-
 lib/libm/src/s_asinh.c       |   8 +-
 lib/libm/src/s_asinhl.c      |  90 +++++
 lib/libm/src/s_expm1.c       |   6 +-
 lib/libm/src/s_log1p.c       |   6 +-
 lib/libm/x86_64/Makefile.inc |   2 +
 lib/libm/x86_64/Symbol.map   |  12 +
 share/mk/bsd.lib.mk          |   9 +
 share/mk/bsd.symver.mk       |  29 ++
 share/mk/version_gen.awk     | 249 ++++++++++++
 33 files changed, 2277 insertions(+), 155 deletions(-)
 create mode 100644 lib/libc/Versions.def
 create mode 100644 lib/libm/Symbol.map
 create mode 100644 lib/libm/i386/Symbol.map
 create mode 100644 lib/libm/ld80/s_logl.c
 create mode 100644 lib/libm/src/e_acoshl.c
 create mode 100644 lib/libm/src/e_atanhl.c
 create mode 100644 lib/libm/src/s_asinhl.c
 create mode 100644 lib/libm/x86_64/Symbol.map
 create mode 100644 share/mk/bsd.symver.mk
 create mode 100644 share/mk/version_gen.awk

diff --git a/lib/libc/Versions.def b/lib/libc/Versions.def
new file mode 100644
index 0000000000..d5041162b2
--- /dev/null
+++ b/lib/libc/Versions.def
@@ -0,0 +1,19 @@
+#
+# Note: Whenever bumping the DFLY version, always make
+#       DFLYprivate_1.0 depend on the new DFLY version.
+#       This will keep it at the end of the dependency chain.
+#
+
+# This is our first version; it depends on no other.
+# This version was first added to DragonFly 3.5 
+DFLY36.0 {
+};
+
+# This is our private namespace.  Any global interfaces that are
+# strictly for use only by other DragonFly applications and libraries
+# are listed here.  We use a separate namespace so we can write
+# simple ABI-checking tools.
+#
+# Please do NOT increment the version of this namespace.
+DFLYprivate_1.0 {
+} DFLY36.0;
diff --git a/lib/libm/Makefile b/lib/libm/Makefile
index 240fef6d4e..a6700dae30 100644
--- a/lib/libm/Makefile
+++ b/lib/libm/Makefile
@@ -62,16 +62,21 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
 # Location of fpmath.h and _fpmath.h
 LIBCDIR=${.CURDIR}/../libc
 CFLAGS+=-I${.CURDIR}/src -I${LIBCDIR}/include -I${LIBCDIR}/${MACHINE_ARCH}
+SYM_MAPS+=	${.CURDIR}/Symbol.map
+
+VERSION_DEF=	${LIBCDIR}/Versions.def
+SYMBOL_MAPS=	${SYM_MAPS}
 
 # C99 long double functions
 COMMON_SRCS+=	s_copysignl.c s_fabsl.c s_modfl.c
 # If long double != double use these; otherwise, we alias the double versions.
-COMMON_SRCS+=	e_acosl.c e_asinl.c e_atan2l.c e_fmodl.c \
-	e_hypotl.c \
+COMMON_SRCS+=	e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \
+	e_fmodl.c e_hypotl.c \
 	invtrig.c k_cosl.c k_sinl.c k_tanl.c \
-	s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \
+	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \
 	s_csqrtl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
-	s_frexpl.c s_nanl.c s_nextafterl.c s_nexttoward.c \
+	s_frexpl.c s_logl.c s_nanl.c s_nextafterl.c \
+	s_nexttoward.c \
 	s_sinl.c s_tanl.c s_truncl.c w_cabsl.c
 
 # C99 complex functions
@@ -107,10 +112,12 @@ MAN=	acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
 
 MLINKS+=acos.3 acosf.3 \
 	acos.3 acosl.3
-MLINKS+=acosh.3 acoshf.3
+MLINKS+=acosh.3 acoshf.3 \
+	acosh.3 acoshl.3
 MLINKS+=asin.3 asinf.3 \
 	asin.3 asinl.3
-MLINKS+=asinh.3 asinhf.3
+MLINKS+=asinh.3 asinhf.3 \
+	asinh.3 asinhl.3
 MLINKS+=atan.3 atanf.3 \
 	atan.3 atanl.3
 MLINKS+=atan2.3 atan2f.3 \
@@ -118,7 +125,8 @@ MLINKS+=atan2.3 atan2f.3 \
 	atan2.3 carg.3 \
 	atan2.3 cargf.3 \
 	atan2.3 cargl.3
-MLINKS+=atanh.3 atanhf.3
+MLINKS+=atanh.3 atanhf.3 \
+	atanh.3 atanhl.3
 MLINKS+=cacos.3 cacosf.3 \
 	cacos.3 cacosh.3 \
 	cacos.3 cacoshf.3 \
@@ -171,6 +179,7 @@ MLINKS+=exp.3 expf.3 \
 	exp.3 exp2l.3 \
 	exp.3 expm1.3 \
 	exp.3 expm1f.3 \
+	exp.3 expm1l.3 \
 	exp.3 pow.3 \
 	exp.3 powf.3
 MLINKS+=fabs.3 fabsf.3 \
@@ -237,10 +246,13 @@ MLINKS+=log.3 logf.3 \
 	log.3 logl.3 \
 	log.3 log10.3 \
 	log.3 log10f.3 \
+	log.3 log10l.3 \
 	log.3 log2.3 \
 	log.3 log2f.3 \
+	log.3 log2l.e \
 	log.3 log1p.3 \
-	log.3 log1pf.3
+	log.3 log1pf.3 \
+	log.3 log1pl.3
 MLINKS+=lrint.3 llrint.3 \
 	lrint.3 llrintf.3 \
 	lrint.3 llrintl.3 \
diff --git a/lib/libm/Symbol.map b/lib/libm/Symbol.map
new file mode 100644
index 0000000000..6129c8e375
--- /dev/null
+++ b/lib/libm/Symbol.map
@@ -0,0 +1,260 @@
+/*
+ * $FreeBSD: head/lib/msun/Symbol.map 251599 2013-06-10 06:04:58Z das $
+ */
+
+DFLY36.0 {
+	__fe_dfl_env;
+	tgamma;
+	acos;
+	acosf;
+	acosh;
+	acoshf;
+	asin;
+	asinf;
+	atan2;
+	atan2f;
+	atanh;
+	atanhf;
+	cosh;
+	coshf;
+	exp;
+	expf;
+	fmod;
+	fmodf;
+	gamma;
+	gamma_r;
+	gammaf;
+	gammaf_r;
+	hypot;
+	hypotf;
+	j0;
+	y0;
+	j0f;
+	y0f;
+	j1;
+	y1;
+	j1f;
+	y1f;
+	jn;
+	yn;
+	jnf;
+	ynf;
+	lgamma;
+	lgamma_r;
+	lgammaf;
+	lgammaf_r;
+	log;
+	log10;
+	log10f;
+	logf;
+	pow;
+	powf;
+	remainder;
+	remainderf;
+	scalb;
+	scalbf;
+	sinh;
+	sinhf;
+	sqrt;
+	sqrtf;
+	asinh;
+	asinhf;
+	atan;
+	atanf;
+	cbrt;
+	cbrtf;
+	ceil;
+	ceilf;
+	ceill;
+	cimag;
+	cimagf;
+	cimagl;
+	conj;
+	conjf;
+	conjl;
+	copysign;
+	copysignf;
+	copysignl;
+	cos;
+	cosf;
+	creal;
+	crealf;
+	creall;
+	erf;
+	erfc;
+	erff;
+	erfcf;
+	exp2;
+	exp2f;
+	expm1;
+	expm1f;
+	fabs;
+	fabsf;
+	fabsl;
+	fdim;
+	fdimf;
+	fdiml;
+	finite;
+	finitef;
+	floor;
+	floorf;
+	floorl;
+	fma;
+	fmaf;
+	fmal;
+	fmax;
+	fmaxf;
+	fmaxl;
+	fmin;
+	fminf;
+	fminl;
+	frexp;
+	frexpf;
+	frexpl;
+	ilogb;
+	ilogbf;
+	ilogbl;
+	__isfinite;
+	__isfinitef;
+	__isfinitel;
+	isnanf;
+	__isnanl;
+	__isnormal;
+	__isnormalf;
+	__isnormall;
+	llrint;
+	llrintf;
+	llround;
+	llroundf;
+	llroundl;
+	log1p;
+	log1pf;
+	logb;
+	logbf;
+	lrint;
+	lrintf;
+	lround;
+	lroundf;
+	lroundl;
+	modff;
+	modfl;
+	nearbyint;
+	nearbyintf;
+	nextafter;
+	nexttoward;
+	nexttowardl;
+	nextafterl;
+	nextafterf;
+	nexttowardf;
+	remquo;
+	remquof;
+	rint;
+	rintf;
+	round;
+	roundf;
+	roundl;
+	scalbln;
+	scalblnf;
+	scalblnl;
+	scalbn;
+	scalbnl;
+	scalbnf;
+	ldexpf;
+	ldexpl;
+	__signbit;
+	__signbitf;
+	__signbitl;
+	signgam;
+	significand;
+	significandf;
+	sin;
+	sinf;
+	tan;
+	tanf;
+	tanh;
+	tanhf;
+	trunc;
+	truncf;
+	truncl;
+	cabs;
+	cabsf;
+	drem;
+	dremf;
+	carg;
+	cargf;
+	csqrt;
+	csqrtf;
+	logbl;
+	nan;
+	nanf;
+	nanl;
+	llrintl;
+	lrintl;
+	nearbyintl;
+	rintl;
+	exp2l;
+	sinl;
+	cosl;
+	tanl;
+	tgammaf;
+	sqrtl;
+	hypotl;
+	cabsl;
+	csqrtl;
+	remquol;
+	remainderl;
+	fmodl;
+	acosl;
+	asinl;
+	atan2l;
+	atanl;
+	cargl;
+	cproj;
+	cprojf;
+	cprojl;
+	__isnanf;
+	cbrtl;
+	cexp;
+	cexpf;
+	log2;
+	log2f;
+	feclearexcept;
+	fegetexceptflag;
+	fetestexcept;
+	fegetround;
+	fesetround;
+	fesetenv;
+	acoshl;
+	asinhl;
+	atanhl;
+	cacos;
+	cacosf;
+	cacosh;
+	cacoshf;
+	casin;
+	casinf;
+	casinh;
+	casinhf;
+	catan;
+	catanf;
+	catanh;
+	catanhf;
+	csin;
+	csinf;
+	csinh;
+	csinhf;
+	ccos;
+	ccosf;
+	ccosh;
+	ccoshf;
+	ctan;
+	ctanf;
+	ctanh;
+	ctanhf;
+	expl;
+	expm1l;
+	log10l;
+	log1pl;
+	log2l;
+	logl;
+};
diff --git a/lib/libm/i386/Makefile.inc b/lib/libm/i386/Makefile.inc
index 6769cee59e..9397a29735 100644
--- a/lib/libm/i386/Makefile.inc
+++ b/lib/libm/i386/Makefile.inc
@@ -16,3 +16,5 @@ ARCH_SRCS+= e_log10f.S e_logf.S e_remainderf.S \
 ARCH_SRCS+= e_remainderl.S e_sqrtl.S s_ceill.S s_copysignl.S \
 	    s_floorl.S s_llrintl.S \
 	    s_logbl.S s_lrintl.S s_remquol.S s_rintl.S s_scalbnl.S s_truncl.S
+
+SYM_MAPS += ${.CURDIR}/i386/Symbol.map
diff --git a/lib/libm/i386/Symbol.map b/lib/libm/i386/Symbol.map
new file mode 100644
index 0000000000..aa0ebc5fad
--- /dev/null
+++ b/lib/libm/i386/Symbol.map
@@ -0,0 +1,14 @@
+/*
+ * $FreeBSD: head/lib/msun/i387/Symbol.map 226218 2011-10-10 15:43:09Z das $
+ */
+DFLY36.0 {
+	__has_sse;
+	__test_sse;
+	fesetexceptflag;
+	feraiseexcept;
+	fegetenv;
+	feholdexcept;
+	feupdateenv;
+	feenableexcept;
+	fedisableexcept;
+};
diff --git a/lib/libm/ld80/s_expl.c b/lib/libm/ld80/s_expl.c
index 48a92ab7e9..3155ed45f2 100644
--- a/lib/libm/ld80/s_expl.c
+++ b/lib/libm/ld80/s_expl.c
@@ -1,5 +1,5 @@
 /*-
- * Copyright (c) 2009-2012 Steven G. Kargl
+ * Copyright (c) 2009-2013 Steven G. Kargl
  * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -25,10 +25,10 @@
  *
  * Optimized by Bruce D. Evans.
  *
- * $FreeBSD: head/lib/msun/ld80/s_expl.c 241516 2012-10-13 19:53:11Z kargl $
+ * $FreeBSD: head/lib/msun/ld80/s_expl.c 251343 2013-06-03 19:51:32Z kargl $
  */
 
-/*-
+/**
  * Compute the exponential of x for Intel 80-bit format.  This is based on:
  *
  *   PTP Tang, "Table-driven implementation of the exponential function
@@ -49,6 +49,7 @@
 #include "math_private.h"
 
 #define	INTERVALS	128
+#define	LOG2_INTERVALS	7
 #define	BIAS	(LDBL_MAX_EXP - 1)
 
 static const long double
@@ -59,9 +60,12 @@ static volatile const long double tiny = 0x1p-10000L;
 
 static const union IEEEl2bits
 /* log(2**16384 - 0.5) rounded towards zero: */
-o_threshold = LD80C(0xb17217f7d1cf79ab, 13,  11356.5234062941439488L),
+/* log(2**16384 - 0.5 + 1) rounded towards zero for expm1l() is the same: */
+o_thresholdu = LD80C(0xb17217f7d1cf79ab, 13,  11356.5234062941439488L),
+#define o_threshold	 (o_thresholdu.e)
 /* log(2**(-16381-64-1)) rounded towards zero: */
-u_threshold = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L);
+u_thresholdu = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L);
+#define u_threshold	 (u_thresholdu.e)
 
 static const double
 /*
@@ -77,11 +81,11 @@ L2 = -3.2819649005320973e-13,		/* -0x1718432a1b0e26.0p-94 */
  * |exp(x) - p(x)| < 2**-77.2
  * (0.002708 is ln2/(2*INTERVALS) rounded up a little).
  */
-P2 =  0.5,
-P3 =  1.6666666666666119e-1,		/*  0x15555555555490.0p-55 */
-P4 =  4.1666666666665887e-2,		/*  0x155555555554e5.0p-57 */
-P5 =  8.3333354987869413e-3,		/*  0x1111115b789919.0p-59 */
-P6 =  1.3888891738560272e-3;		/*  0x16c16c651633ae.0p-62 */
+A2 =  0.5,
+A3 =  1.6666666666666119e-1,		/*  0x15555555555490.0p-55 */
+A4 =  4.1666666666665887e-2,		/*  0x155555555554e5.0p-57 */
+A5 =  8.3333354987869413e-3,		/*  0x1111115b789919.0p-59 */
+A6 =  1.3888891738560272e-3;		/*  0x16c16c651633ae.0p-62 */
 
 /*
  * 2^(i/INTERVALS) for i in [0,INTERVALS] is represented by two values where
@@ -95,8 +99,7 @@ P6 =  1.3888891738560272e-3;		/*  0x16c16c651633ae.0p-62 */
 static const struct {
 	double	hi;
 	double	lo;
-/* XXX should rename 's'. */
-} s[INTERVALS] = {
+} tbl[INTERVALS] = {
 	0x1p+0, 0x0p+0,
 	0x1.0163da9fb3335p+0, 0x1.b61299ab8cdb7p-54,
 	0x1.02c9a3e778060p+0, 0x1.dcdef95949ef4p-53,
@@ -231,7 +234,8 @@ long double
 expl(long double x)
 {
 	union IEEEl2bits u, v;
-	long double fn, q, r, r1, r2, t, t23, t45, twopk, twopkp10000, z;
+	long double fn, q, r, r1, r2, t, twopk, twopkp10000;
+	long double z;
 	int k, n, n2;
 	uint16_t hx, ix;
 
@@ -241,23 +245,21 @@ expl(long double x)
 	ix = hx & 0x7fff;
 	if (ix >= BIAS + 13) {		/* |x| >= 8192 or x is NaN */
 		if (ix == BIAS + LDBL_MAX_EXP) {
-			if (hx & 0x8000 && u.xbits.man == 1ULL << 63)
-				return (0.0L);	/* x is -Inf */
-			return (x + x); /* x is +Inf, NaN or unsupported */
+			if (hx & 0x8000)  /* x is -Inf, -NaN or unsupported */
+				return (-1 / x);
+ 			return (x + x);	/* x is +Inf, +NaN or unsupported */
 		}
-		if (x > o_threshold.e)
+		if (x > o_threshold)
 			return (huge * huge);
-		if (x < u_threshold.e)
+		if (x < u_threshold)
 			return (tiny * tiny);
-	} else if (ix < BIAS - 66) {	/* |x| < 0x1p-66 */
-					/* includes pseudo-denormals */
-		if (huge + x > 1.0L)	/* trigger inexact iff x != 0 */
-			return (1.0L + x);
+	} else if (ix < BIAS - 65) {	/* |x| < 0x1p-65 (includes pseudos) */
+		return (1 + x);		/* 1 with inexact iff x != 0 */
 	}
 
 	ENTERI();
 
-	/* Reduce x to (k*ln2 + midpoint[n2] + r1 + r2). */
+	/* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
 	/* Use a specialized rint() to get fn.  Assume round-to-nearest. */
 	fn = x * INV_L + 0x1.8p63 - 0x1.8p63;
 	r = x - fn * L1 - fn * L2;	/* r = r1 + r2 done independently. */
@@ -269,12 +271,13 @@ expl(long double x)
 	n  = (int)fn;
 #endif
 	n2 = (unsigned)n % INTERVALS;
-	k = (n - n2) / INTERVALS;
+	/* Depend on the sign bit being propagated: */
+	k = n >> LOG2_INTERVALS;
 	r1 = x - fn * L1;
-	r2 = -fn * L2;
+	r2 = fn * -L2;
 
 	/* Prepare scale factors. */
-	v.xbits.man = 1ULL << 63;
+	v.e = 1;
 	if (k >= LDBL_MIN_EXP) {
 		v.xbits.expsign = BIAS + k;
 		twopk = v.e;
@@ -283,21 +286,183 @@ expl(long double x)
 		twopkp10000 = v.e;
 	}
 
-	/* Evaluate expl(midpoint[n2] + r1 + r2) = s[n2] * expl(r1 + r2). */
-	/* Here q = q(r), not q(r1), since r1 is lopped like L1. */
-	t45 = r * P5 + P4;
+	/* Evaluate expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2). */
 	z = r * r;
-	t23 = r * P3 + P2;
-	q = r2 + z * t23 + z * z * t45 + z * z * z * P6;
-	t = (long double)s[n2].lo + s[n2].hi;
-	t = s[n2].lo + t * (q + r1) + s[n2].hi;
+	q = r2 + z * (A2 + r * A3) + z * z * (A4 + r * A5) + z * z * z * A6;
+	t = (long double)tbl[n2].lo + tbl[n2].hi;
+	t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
 
 	/* Scale by 2**k. */
 	if (k >= LDBL_MIN_EXP) {
 		if (k == LDBL_MAX_EXP)
-			RETURNI(t * 2.0L * 0x1p16383L);
+			RETURNI(t * 2 * 0x1p16383L);
 		RETURNI(t * twopk);
 	} else {
 		RETURNI(t * twopkp10000 * twom10000);
 	}
 }
+
+/**
+ * Compute expm1l(x) for Intel 80-bit format.  This is based on:
+ *
+ *   PTP Tang, "Table-driven implementation of the Expm1 function
+ *   in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 18,
+ *   211-222 (1992).
+ */
+
+/*
+ * Our T1 and T2 are chosen to be approximately the points where method
+ * A and method B have the same accuracy.  Tang's T1 and T2 are the
+ * points where method A's accuracy changes by a full bit.  For Tang,
+ * this drop in accuracy makes method A immediately less accurate than
+ * method B, but our larger INTERVALS makes method A 2 bits more
+ * accurate so it remains the most accurate method significantly
+ * closer to the origin despite losing the full bit in our extended
+ * range for it.
+ */
+static const double
+T1 = -0.1659,				/* ~-30.625/128 * log(2) */
+T2 =  0.1659;				/* ~30.625/128 * log(2) */
+
+/*
+ * Domain [-0.1659, 0.1659], range ~[-1.2027e-22, 3.4417e-22]:
+ * |(exp(x)-1-x-x**2/2)/x - p(x)| < 2**-71.2
+ */
+static const union IEEEl2bits
+B3 = LD80C(0xaaaaaaaaaaaaaaab, -3,  1.66666666666666666671e-1L),
+B4 = LD80C(0xaaaaaaaaaaaaaaac, -5,  4.16666666666666666712e-2L);
+
+static const double
+B5  =  8.3333333333333245e-3,		/*  0x1.111111111110cp-7 */
+B6  =  1.3888888888888861e-3,		/*  0x1.6c16c16c16c0ap-10 */
+B7  =  1.9841269841532042e-4,		/*  0x1.a01a01a0319f9p-13 */
+B8  =  2.4801587302069236e-5,		/*  0x1.a01a01a03cbbcp-16 */
+B9  =  2.7557316558468562e-6,		/*  0x1.71de37fd33d67p-19 */
+B10 =  2.7557315829785151e-7,		/*  0x1.27e4f91418144p-22 */
+B11 =  2.5063168199779829e-8,		/*  0x1.ae94fabdc6b27p-26 */
+B12 =  2.0887164654459567e-9;		/*  0x1.1f122d6413fe1p-29 */
+
+long double
+expm1l(long double x)
+{
+	union IEEEl2bits u, v;
+	long double fn, hx2_hi, hx2_lo, q, r, r1, r2, t, twomk, twopk, x_hi;
+	long double x_lo, x2, z;
+	long double x4;
+	int k, n, n2;
+	uint16_t hx, ix;
+
+	/* Filter out exceptional cases. */
+	u.e = x;
+	hx = u.xbits.expsign;
+	ix = hx & 0x7fff;
+	if (ix >= BIAS + 6) {		/* |x| >= 64 or x is NaN */
+		if (ix == BIAS + LDBL_MAX_EXP) {
+			if (hx & 0x8000)  /* x is -Inf, -NaN or unsupported */
+				return (-1 / x - 1);
+			return (x + x);	/* x is +Inf, +NaN or unsupported */
+		}
+		if (x > o_threshold)
+			return (huge * huge);
+		/*
+		 * expm1l() never underflows, but it must avoid
+		 * unrepresentable large negative exponents.  We used a
+		 * much smaller threshold for large |x| above than in
+		 * expl() so as to handle not so large negative exponents
+		 * in the same way as large ones here.
+		 */
+		if (hx & 0x8000)	/* x <= -64 */
+			return (tiny - 1);	/* good for x < -65ln2 - eps */
+	}
+
+	ENTERI();
+
+	if (T1 < x && x < T2) {
+		if (ix < BIAS - 64) {	/* |x| < 0x1p-64 (includes pseudos) */
+			/* x (rounded) with inexact if x != 0: */
+			RETURNI(x == 0 ? x :
+			    (0x1p100 * x + fabsl(x)) * 0x1p-100);
+		}
+
+		x2 = x * x;
+		x4 = x2 * x2;
+		q = x4 * (x2 * (x4 *
+		    /*
+		     * XXX the number of terms is no longer good for
+		     * pairwise grouping of all except B3, and the
+		     * grouping is no longer from highest down.
+		     */
+		    (x2 *            B12  + (x * B11 + B10)) +
+		    (x2 * (x * B9 +  B8) +  (x * B7 +  B6))) +
+			  (x * B5 +  B4.e)) + x2 * x * B3.e;
+
+		x_hi = (float)x;
+		x_lo = x - x_hi;
+		hx2_hi = x_hi * x_hi / 2;
+		hx2_lo = x_lo * (x + x_hi) / 2;
+		if (ix >= BIAS - 7)
+			RETURNI(hx2_lo + x_lo + q + (hx2_hi + x_hi));
+		else
+			RETURNI(hx2_lo + q + hx2_hi + x);
+	}
+
+	/* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
+	/* Use a specialized rint() to get fn.  Assume round-to-nearest. */
+	fn = x * INV_L + 0x1.8p63 - 0x1.8p63;
+#if defined(HAVE_EFFICIENT_IRINTL)
+	n = irintl(fn);
+#elif defined(HAVE_EFFICIENT_IRINT)
+	n = irint(fn);
+#else
+	n = (int)fn;
+#endif
+	n2 = (unsigned)n % INTERVALS;
+	k = n >> LOG2_INTERVALS;
+	r1 = x - fn * L1;
+	r2 = fn * -L2;
+	r = r1 + r2;
+
+	/* Prepare scale factor. */
+	v.e = 1;
+	v.xbits.expsign = BIAS + k;
+	twopk = v.e;
+
+	/*
+	 * Evaluate lower terms of
+	 * expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2).
+	 */
+	z = r * r;
+	q = r2 + z * (A2 + r * A3) + z * z * (A4 + r * A5) + z * z * z * A6;
+
+	t = (long double)tbl[n2].lo + tbl[n2].hi;
+
+	if (k == 0) {
+		t = tbl[n2].lo * (r1 + 1) + t * q + tbl[n2].hi * r1 +
+		    (tbl[n2].hi - 1);
+		RETURNI(t);
+	}
+	if (k == -1) {
+		t = tbl[n2].lo * (r1 + 1) + t * q + tbl[n2].hi * r1 + 
+		    (tbl[n2].hi - 2);
+		RETURNI(t / 2);
+	}
+	if (k < -7) {
+		t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+		RETURNI(t * twopk - 1);
+	}
+	if (k > 2 * LDBL_MANT_DIG - 1) {
+		t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+		if (k == LDBL_MAX_EXP)
+			RETURNI(t * 2 * 0x1p16383L - 1);
+		RETURNI(t * twopk - 1);
+	}
+
+	v.xbits.expsign = BIAS - k;
+	twomk = v.e;
+
+	if (k > LDBL_MANT_DIG - 1)
+		t = tbl[n2].lo - twomk + t * (q + r1) + tbl[n2].hi;
+	else
+		t = tbl[n2].lo + t * (q + r1) + (tbl[n2].hi - twomk);
+	RETURNI(t * twopk);
+}
diff --git a/lib/libm/ld80/s_logl.c b/lib/libm/ld80/s_logl.c
new file mode 100644
index 0000000000..e8b41ca89d
--- /dev/null
+++ b/lib/libm/ld80/s_logl.c
@@ -0,0 +1,716 @@
+/*-
+ * Copyright (c) 2007-2013 Bruce D. Evans
+ * 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.
+ *
+ * $FreeBSD: head/lib/msun/ld80/s_logl.c 251292 2013-06-03 09:14:31Z das $
+ */
+
+/**
+ * Implementation of the natural logarithm of x for Intel 80-bit format.
+ *
+ * First decompose x into its base 2 representation:
+ *
+ *    log(x) = log(X * 2**k), where X is in [1, 2)
+ *           = log(X) + k * log(2).
+ *
+ * Let X = X_i + e, where X_i is the center of one of the intervals
+ * [-1.0/256, 1.0/256), [1.0/256, 3.0/256), .... [2.0-1.0/256, 2.0+1.0/256)
+ * and X is in this interval.  Then
+ *
+ *    log(X) = log(X_i + e)
+ *           = log(X_i * (1 + e / X_i))
+ *           = log(X_i) + log(1 + e / X_i).
+ *
+ * The values log(X_i) are tabulated below.  Let d = e / X_i and use
+ *
+ *    log(1 + d) = p(d)
+ *
+ * where p(d) = d - 0.5*d*d + ... is a special minimax polynomial of
+ * suitably high degree.
+ *
+ * To get sufficiently small roundoff errors, k * log(2), log(X_i), and
+ * sometimes (if |k| is not large) the first term in p(d) must be evaluated
+ * and added up in extra precision.  Extra precision is not needed for the
+ * rest of p(d).  In the worst case when k = 0 and log(X_i) is 0, the final
+ * error is controlled mainly by the error in the second term in p(d).  The
+ * error in this term itself is at most 0.5 ulps from the d*d operation in
+ * it.  The error in this term relative to the first term is thus at most
+ * 0.5 * |-0.5| * |d| < 1.0/1024 ulps.  We aim for an accumulated error of
+ * at most twice this at the point of the final rounding step.  Thus the
+ * final error should be at most 0.5 + 1.0/512 = 0.5020 ulps.  Exhaustive
+ * testing of a float variant of this function showed a maximum final error
+ * of 0.5008 ulps.  Non-exhaustive testing of a double variant of this
+ * function showed a maximum final error of 0.5078 ulps (near 1+1.0/256).
+ *
+ * We made the maximum of |d| (and thus the total relative error and the
+ * degree of p(d)) small by using a large number of intervals.  Using
+ * centers of intervals instead of endpoints reduces this maximum by a
+ * factor of 2 for a given number of intervals.  p(d) is special only
+ * in beginning with the Taylor coefficients 0 + 1*d, which tends to happen
+ * naturally.  The most accurate minimax polynomial of a given degree might
+ * be different, but then we wouldn't want it since we would have to do
+ * extra work to avoid roundoff error (especially for P0*d instead of d).
+ */
+
+#ifdef DEBUG
+#include <assert.h>
+#include <fenv.h>
+#endif
+
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#define	i386_SSE_GOOD
+#ifndef NO_STRUCT_RETURN
+#define	STRUCT_RETURN
+#endif
+#include "math_private.h"
+
+#if !defined(NO_UTAB) && !defined(NO_UTABL)
+#define	USE_UTAB
+#endif
+
+/*
+ * Domain [-0.005280, 0.004838], range ~[-5.1736e-22, 5.1738e-22]:
+ * |log(1 + d)/d - p(d)| < 2**-70.7
+ */
+static const double
+P2 = -0.5,
+P3 =  3.3333333333333359e-1,		/*  0x1555555555555a.0p-54 */
+P4 = -2.5000000000004424e-1,		/* -0x1000000000031d.0p-54 */
+P5 =  1.9999999992970016e-1,		/*  0x1999999972f3c7.0p-55 */
+P6 = -1.6666666072191585e-1,		/* -0x15555548912c09.0p-55 */
+P7 =  1.4286227413310518e-1,		/*  0x12494f9d9def91.0p-55 */
+P8 = -1.2518388626763144e-1;		/* -0x1006068cc0b97c.0p-55 */
+
+static volatile const double zero = 0;
+
+#define	INTERVALS	128
+#define	LOG2_INTERVALS	7
+#define	TSIZE		(INTERVALS + 1)
+#define	G(i)		(T[(i)].G)
+#define	F_hi(i)		(T[(i)].F_hi)
+#define	F_lo(i)		(T[(i)].F_lo)
+#define	ln2_hi		F_hi(TSIZE - 1)
+#define	ln2_lo		F_lo(TSIZE - 1)
+#define	E(i)		(U[(i)].E)
+#define	H(i)		(U[(i)].H)
+
+static const struct {
+	float	G;			/* 1/(1 + i/128) rounded to 8/9 bits */
+	float	F_hi;			/* log(1 / G_i) rounded (see below) */
+	double	F_lo;			/* next 53 bits for log(1 / G_i) */
+} T[TSIZE] = {
+	/*
+	 * ln2_hi and each F_hi(i) are rounded to a number of bits that
+	 * makes F_hi(i) + dk*ln2_hi exact for all i and all dk.
+	 *
+	 * The last entry (for X just below 2) is used to define ln2_hi
+	 * and ln2_lo, to ensure that F_hi(i) and F_lo(i) cancel exactly
+	 * with dk*ln2_hi and dk*ln2_lo, respectively, when dk = -1.
+	 * This is needed for accuracy when x is just below 1.  (To avoid
+	 * special cases, such x are "reduced" strangely to X just below
+	 * 2 and dk = -1, and then the exact cancellation is needed
+	 * because any the error from any non-exactness would be too
+	 * large).
+	 *
+	 * We want to share this table between double precision and ld80,
+	 * so the relevant range of dk is the larger one of ld80
+	 * ([-16445, 16383]) and the relevant exactness requirement is
+	 * the stricter one of double precision.  The maximum number of
+	 * bits in F_hi(i) that works is very dependent on i but has
+	 * a minimum of 33.  We only need about 12 bits in F_hi(i) for
+	 * it to provide enough extra precision in double precision (11
+	 * more than that are required for ld80).
+	 *
+	 * We round F_hi(i) to 24 bits so that it can have type float,
+	 * mainly to minimize the size of the table.  Using all 24 bits
+	 * in a float for it automatically satisfies the above constraints.
+	 */
+	 0x800000.0p-23,  0,               0,
+	 0xfe0000.0p-24,  0x8080ac.0p-30, -0x14ee431dae6675.0p-84,
+	 0xfc0000.0p-24,  0x8102b3.0p-29, -0x1db29ee2d83718.0p-84,
+	 0xfa0000.0p-24,  0xc24929.0p-29,  0x1191957d173698.0p-83,
+	 0xf80000.0p-24,  0x820aec.0p-28,  0x13ce8888e02e79.0p-82,
+	 0xf60000.0p-24,  0xa33577.0p-28, -0x17a4382ce6eb7c.0p-82,
+	 0xf48000.0p-24,  0xbc42cb.0p-28, -0x172a21161a1076.0p-83,
+	 0xf30000.0p-24,  0xd57797.0p-28, -0x1e09de07cb9589.0p-82,
+	 0xf10000.0p-24,  0xf7518e.0p-28,  0x1ae1eec1b036c5.0p-91,
+	 0xef0000.0p-24,  0x8cb9df.0p-27, -0x1d7355325d560e.0p-81,
+	 0xed8000.0p-24,  0x999ec0.0p-27, -0x1f9f02d256d503.0p-82,
+	 0xec0000.0p-24,  0xa6988b.0p-27, -0x16fc0a9d12c17a.0p-83,
+	 0xea0000.0p-24,  0xb80698.0p-27,  0x15d581c1e8da9a.0p-81,
+	 0xe80000.0p-24,  0xc99af3.0p-27, -0x1535b3ba8f150b.0p-83,
+	 0xe70000.0p-24,  0xd273b2.0p-27,  0x163786f5251af0.0p-85,
+	 0xe50000.0p-24,  0xe442c0.0p-27,  0x1bc4b2368e32d5.0p-84,
+	 0xe38000.0p-24,  0xf1b83f.0p-27,  0x1c6090f684e676.0p-81,
+	 0xe20000.0p-24,  0xff448a.0p-27, -0x1890aa69ac9f42.0p-82,
+	 0xe08000.0p-24,  0x8673f6.0p-26,  0x1b9985194b6b00.0p-80,
+	 0xdf0000.0p-24,  0x8d515c.0p-26, -0x1dc08d61c6ef1e.0p-83,
+	 0xdd8000.0p-24,  0x943a9e.0p-26, -0x1f72a2dac729b4.0p-82,
+	 0xdc0000.0p-24,  0x9b2fe6.0p-26, -0x1fd4dfd3a0afb9.0p-80,
+	 0xda8000.0p-24,  0xa2315d.0p-26, -0x11b26121629c47.0p-82,
+	 0xd90000.0p-24,  0xa93f2f.0p-26,  0x1286d633e8e569.0p-81,
+	 0xd78000.0p-24,  0xb05988.0p-26,  0x16128eba936770.0p-84,
+	 0xd60000.0p-24,  0xb78094.0p-26,  0x16ead577390d32.0p-80,
+	 0xd50000.0p-24,  0xbc4c6c.0p-26,  0x151131ccf7c7b7.0p-81,
+	 0xd38000.0p-24,  0xc3890a.0p-26, -0x115e2cd714bd06.0p-80,
+	 0xd20000.0p-24,  0xcad2d7.0p-26, -0x1847f406ebd3b0.0p-82,
+	 0xd10000.0p-24,  0xcfb620.0p-26,  0x1c2259904d6866.0p-81,
+	 0xcf8000.0p-24,  0xd71653.0p-26,  0x1ece57a8d5ae55.0p-80,
+	 0xce0000.0p-24,  0xde843a.0p-26, -0x1f109d4bc45954.0p-81,
+	 0xcd0000.0p-24,  0xe37fde.0p-26,  0x1bc03dc271a74d.0p-81,
+	 0xcb8000.0p-24,  0xeb050c.0p-26, -0x1bf2badc0df842.0p-85,
+	 0xca0000.0p-24,  0xf29878.0p-26, -0x18efededd89fbe.0p-87,
+	 0xc90000.0p-24,  0xf7ad6f.0p-26,  0x1373ff977baa69.0p-81,
+	 0xc80000.0p-24,  0xfcc8e3.0p-26,  0x196766f2fb3283.0p-80,
+	 0xc68000.0p-24,  0x823f30.0p-25,  0x19bd076f7c434e.0p-79,
+	 0xc58000.0p-24,  0x84d52c.0p-25, -0x1a327257af0f46.0p-79,
+	 0xc40000.0p-24,  0x88bc74.0p-25,  0x113f23def19c5a.0p-81,
+	 0xc30000.0p-24,  0x8b5ae6.0p-25,  0x1759f6e6b37de9.0p-79,
+	 0xc20000.0p-24,  0x8dfccb.0p-25,  0x1ad35ca6ed5148.0p-81,
+	 0xc10000.0p-24,  0x90a22b.0p-25,  0x1a1d71a87deba4.0p-79,
+	 0xbf8000.0p-24,  0x94a0d8.0p-25, -0x139e5210c2b731.0p-80,
+	 0xbe8000.0p-24,  0x974f16.0p-25, -0x18f6ebcff3ed73.0p-81,
+	 0xbd8000.0p-24,  0x9a00f1.0p-25, -0x1aa268be39aab7.0p-79,
+	 0xbc8000.0p-24,  0x9cb672.0p-25, -0x14c8815839c566.0p-79,
+	 0xbb0000.0p-24,  0xa0cda1.0p-25,  0x1eaf46390dbb24.0p-81,
+	 0xba0000.0p-24,  0xa38c6e.0p-25,  0x138e20d831f698.0p-81,
+	 0xb90000.0p-24,  0xa64f05.0p-25, -0x1e8d3c41123616.0p-82,
+	 0xb80000.0p-24,  0xa91570.0p-25,  0x1ce28f5f3840b2.0p-80,
+	 0xb70000.0p-24,  0xabdfbb.0p-25, -0x186e5c0a424234.0p-79,
+	 0xb60000.0p-24,  0xaeadef.0p-25, -0x14d41a0b2a08a4.0p-83,
+	 0xb50000.0p-24,  0xb18018.0p-25,  0x16755892770634.0p-79,
+	 0xb40000.0p-24,  0xb45642.0p-25, -0x16395ebe59b152.0p-82,
+	 0xb30000.0p-24,  0xb73077.0p-25,  0x1abc65c8595f09.0p-80,
+	 0xb20000.0p-24,  0xba0ec4.0p-25, -0x1273089d3dad89.0p-79,
+	 0xb10000.0p-24,  0xbcf133.0p-25,  0x10f9f67b1f4bbf.0p-79,
+	 0xb00000.0p-24,  0xbfd7d2.0p-25, -0x109fab90486409.0p-80,
+	 0xaf0000.0p-24,  0xc2c2ac.0p-25, -0x1124680aa43333.0p-79,
+	 0xae8000.0p-24,  0xc439b3.0p-25, -0x1f360cc4710fc0.0p-80,
+	 0xad8000.0p-24,  0xc72afd.0p-25, -0x132d91f21d89c9.0p-80,
+	 0xac8000.0p-24,  0xca20a2.0p-25, -0x16bf9b4d1f8da8.0p-79,
+	 0xab8000.0p-24,  0xcd1aae.0p-25,  0x19deb5ce6a6a87.0p-81,
+	 0xaa8000.0p-24,  0xd0192f.0p-25,  0x1a29fb48f7d3cb.0p-79,
+	 0xaa0000.0p-24,  0xd19a20.0p-25,  0x1127d3c6457f9d.0p-81,
+	 0xa90000.0p-24,  0xd49f6a.0p-25, -0x1ba930e486a0ac.0p-81,
+	 0xa80000.0p-24,  0xd7a94b.0p-25, -0x1b6e645f31549e.0p-79,
+	 0xa70000.0p-24,  0xdab7d0.0p-25,  0x1118a425494b61.0p-80,
+	 0xa68000.0p-24,  0xdc40d5.0p-25,  0x1966f24d29d3a3.0p-80,
+	 0xa58000.0p-24,  0xdf566d.0p-25, -0x1d8e52eb2248f1.0p-82,
+	 0xa48000.0p-24,  0xe270ce.0p-25, -0x1ee370f96e6b68.0p-80,
+	 0xa40000.0p-24,  0xe3ffce.0p-25,  0x1d155324911f57.0p-80,
+	 0xa30000.0p-24,  0xe72179.0p-25, -0x1fe6e2f2f867d9.0p-80,
+	 0xa20000.0p-24,  0xea4812.0p-25,  0x1b7be9add7f4d4.0p-80,
+	 0xa18000.0p-24,  0xebdd3d.0p-25,  0x1b3cfb3f7511dd.0p-79,
+	 0xa08000.0p-24,  0xef0b5b.0p-25, -0x1220de1f730190.0p-79,
+	 0xa00000.0p-24,  0xf0a451.0p-25, -0x176364c9ac81cd.0p-80,
+	 0x9f0000.0p-24,  0xf3da16.0p-25,  0x1eed6b9aafac8d.0p-81,
+	 0x9e8000.0p-24,  0xf576e9.0p-25,  0x1d593218675af2.0p-79,
+	 0x9d8000.0p-24,  0xf8b47c.0p-25, -0x13e8eb7da053e0.0p-84,
+	 0x9d0000.0p-24,  0xfa553f.0p-25,  0x1c063259bcade0.0p-79,
+	 0x9c0000.0p-24,  0xfd9ac5.0p-25,  0x1ef491085fa3c1.0p-79,
+	 0x9b8000.0p-24,  0xff3f8c.0p-25,  0x1d607a7c2b8c53.0p-79,
+	 0x9a8000.0p-24,  0x814697.0p-24, -0x12ad3817004f3f.0p-78,
+	 0x9a0000.0p-24,  0x821b06.0p-24, -0x189fc53117f9e5.0p-81,
+	 0x990000.0p-24,  0x83c5f8.0p-24,  0x14cf15a048907b.0p-79,
+	 0x988000.0p-24,  0x849c7d.0p-24,  0x1cbb1d35fb8287.0p-78,
+	 0x978000.0p-24,  0x864ba6.0p-24,  0x1128639b814f9c.0p-78,
+	 0x970000.0p-24,  0x87244c.0p-24,  0x184733853300f0.0p-79,
+	 0x968000.0p-24,  0x87fdaa.0p-24,  0x109d23aef77dd6.0p-80,
+	 0x958000.0p-24,  0x89b293.0p-24, -0x1a81ef367a59de.0p-78,
+	 0x950000.0p-24,  0x8a8e20.0p-24, -0x121ad3dbb2f452.0p-78,
+	 0x948000.0p-24,  0x8b6a6a.0p-24, -0x1cfb981628af72.0p-79,
+	 0x938000.0p-24,  0x8d253a.0p-24, -0x1d21730ea76cfe.0p-79,
+	 0x930000.0p-24,  0x8e03c2.0p-24,  0x135cc00e566f77.0p-78,
+	 0x928000.0p-24,  0x8ee30d.0p-24, -0x10fcb5df257a26.0p-80,
+	 0x918000.0p-24,  0x90a3ee.0p-24, -0x16e171b15433d7.0p-79,
+	 0x910000.0p-24,  0x918587.0p-24, -0x1d050da07f3237.0p-79,
+	 0x908000.0p-24,  0x9267e7.0p-24,  0x1be03669a5268d.0p-79,
+	 0x8f8000.0p-24,  0x942f04.0p-24,  0x10b28e0e26c337.0p-79,
+	 0x8f0000.0p-24,  0x9513c3.0p-24,  0x1a1d820da57cf3.0p-78,
+	 0x8e8000.0p-24,  0x95f950.0p-24, -0x19ef8f13ae3cf1.0p-79,
+	 0x8e0000.0p-24,  0x96dfab.0p-24, -0x109e417a6e507c.0p-78,
+	 0x8d0000.0p-24,  0x98aed2.0p-24,  0x10d01a2c5b0e98.0p-79,
+	 0x8c8000.0p-24,  0x9997a2.0p-24, -0x1d6a50d4b61ea7.0p-78,
+	 0x8c0000.0p-24,  0x9a8145.0p-24,  0x1b3b190b83f952.0p-78,
+	 0x8b8000.0p-24,  0x9b6bbf.0p-24,  0x13a69fad7e7abe.0p-78,
+	 0x8b0000.0p-24,  0x9c5711.0p-24, -0x11cd12316f576b.0p-78,
+	 0x8a8000.0p-24,  0x9d433b.0p-24,  0x1c95c444b807a2.0p-79,
+	 0x898000.0p-24,  0x9f1e22.0p-24, -0x1b9c224ea698c3.0p-79,
+	 0x890000.0p-24,  0xa00ce1.0p-24,  0x125ca93186cf0f.0p-81,
+	 0x888000.0p-24,  0xa0fc80.0p-24, -0x1ee38a7bc228b3.0p-79,
+	 0x880000.0p-24,  0xa1ed00.0p-24, -0x1a0db876613d20.0p-78,
+	 0x878000.0p-24,  0xa2de62.0p-24,  0x193224e8516c01.0p-79,
+	 0x870000.0p-24,  0xa3d0a9.0p-24,  0x1fa28b4d2541ad.0p-79,
+	 0x868000.0p-24,  0xa4c3d6.0p-24,  0x1c1b5760fb4572.0p-78,
+	 0x858000.0p-24,  0xa6acea.0p-24,  0x1fed5d0f65949c.0p-80,
+	 0x850000.0p-24,  0xa7a2d4.0p-24,  0x1ad270c9d74936.0p-80,
+	 0x848000.0p-24,  0xa899ab.0p-24,  0x199ff15ce53266.0p-79,
+	 0x840000.0p-24,  0xa99171.0p-24,  0x1a19e15ccc45d2.0p-79,
+	 0x838000.0p-24,  0xaa8a28.0p-24, -0x121a14ec532b36.0p-80,
+	 0x830000.0p-24,  0xab83d1.0p-24,  0x1aee319980bff3.0p-79,
+	 0x828000.0p-24,  0xac7e6f.0p-24, -0x18ffd9e3900346.0p-80,
+	 0x820000.0p-24,  0xad7a03.0p-24, -0x1e4db102ce29f8.0p-80,
+	 0x818000.0p-24,  0xae768f.0p-24,  0x17c35c55a04a83.0p-81,
+	 0x810000.0p-24,  0xaf7415.0p-24,  0x1448324047019b.0p-78,
+	 0x808000.0p-24,  0xb07298.0p-24, -0x1750ee3915a198.0p-78,
+	 0x800000.0p-24,  0xb17218.0p-24, -0x105c610ca86c39.0p-81,
+};
+
+#ifdef USE_UTAB
+static const struct {
+	float	H;			/* 1 + i/INTERVALS (exact) */
+	float	E;			/* H(i) * G(i) - 1 (exact) */
+} U[TSIZE] = {
+	 0x800000.0p-23,  0,
+	 0x810000.0p-23, -0x800000.0p-37,
+	 0x820000.0p-23, -0x800000.0p-35,
+	 0x830000.0p-23, -0x900000.0p-34,
+	 0x840000.0p-23, -0x800000.0p-33,
+	 0x850000.0p-23, -0xc80000.0p-33,
+	 0x860000.0p-23, -0xa00000.0p-36,
+	 0x870000.0p-23,  0x940000.0p-33,
+	 0x880000.0p-23,  0x800000.0p-35,
+	 0x890000.0p-23, -0xc80000.0p-34,
+	 0x8a0000.0p-23,  0xe00000.0p-36,
+	 0x8b0000.0p-23,  0x900000.0p-33,
+	 0x8c0000.0p-23, -0x800000.0p-35,
+	 0x8d0000.0p-23, -0xe00000.0p-33,
+	 0x8e0000.0p-23,  0x880000.0p-33,
+	 0x8f0000.0p-23, -0xa80000.0p-34,
+	 0x900000.0p-23, -0x800000.0p-35,
+	 0x910000.0p-23,  0x800000.0p-37,
+	 0x920000.0p-23,  0x900000.0p-35,
+	 0x930000.0p-23,  0xd00000.0p-35,
+	 0x940000.0p-23,  0xe00000.0p-35,
+	 0x950000.0p-23,  0xc00000.0p-35,
+	 0x960000.0p-23,  0xe00000.0p-36,
+	 0x970000.0p-23, -0x800000.0p-38,
+	 0x980000.0p-23, -0xc00000.0p-35,
+	 0x990000.0p-23, -0xd00000.0p-34,
+	 0x9a0000.0p-23,  0x880000.0p-33,
+	 0x9b0000.0p-23,  0xe80000.0p-35,
+	 0x9c0000.0p-23, -0x800000.0p-35,
+	 0x9d0000.0p-23,  0xb40000.0p-33,
+	 0x9e0000.0p-23,  0x880000.0p-34,
+	 0x9f0000.0p-23, -0xe00000.0p-35,
+	 0xa00000.0p-23,  0x800000.0p-33,
+	 0xa10000.0p-23, -0x900000.0p-36,
+	 0xa20000.0p-23, -0xb00000.0p-33,
+	 0xa30000.0p-23, -0xa00000.0p-36,
+	 0xa40000.0p-23,  0x800000.0p-33,
+	 0xa50000.0p-23, -0xf80000.0p-35,
+	 0xa60000.0p-23,  0x880000.0p-34,
+	 0xa70000.0p-23, -0x900000.0p-33,
+	 0xa80000.0p-23, -0x800000.0p-35,
+	 0xa90000.0p-23,  0x900000.0p-34,
+	 0xaa0000.0p-23,  0xa80000.0p-33,
+	 0xab0000.0p-23, -0xac0000.0p-34,
+	 0xac0000.0p-23, -0x800000.0p-37,
+	 0xad0000.0p-23,  0xf80000.0p-35,
+	 0xae0000.0p-23,  0xf80000.0p-34,
+	 0xaf0000.0p-23, -0xac0000.0p-33,
+	 0xb00000.0p-23, -0x800000.0p-33,
+	 0xb10000.0p-23, -0xb80000.0p-34,
+	 0xb20000.0p-23, -0x800000.0p-34,
+	 0xb30000.0p-23, -0xb00000.0p-35,
+	 0xb40000.0p-23, -0x800000.0p-35,
+	 0xb50000.0p-23, -0xe00000.0p-36,
+	 0xb60000.0p-23, -0x800000.0p-35,
+	 0xb70000.0p-23, -0xb00000.0p-35,
+	 0xb80000.0p-23, -0x800000.0p-34,
+	 0xb90000.0p-23, -0xb80000.0p-34,
+	 0xba0000.0p-23, -0x800000.0p-33,
+	 0xbb0000.0p-23, -0xac0000.0p-33,
+	 0xbc0000.0p-23,  0x980000.0p-33,
+	 0xbd0000.0p-23,  0xbc0000.0p-34,
+	 0xbe0000.0p-23,  0xe00000.0p-36,
+	 0xbf0000.0p-23, -0xb80000.0p-35,
+	 0xc00000.0p-23, -0x800000.0p-33,
+	 0xc10000.0p-23,  0xa80000.0p-33,
+	 0xc20000.0p-23,  0x900000.0p-34,
+	 0xc30000.0p-23, -0x800000.0p-35,
+	 0xc40000.0p-23, -0x900000.0p-33,
+	 0xc50000.0p-23,  0x820000.0p-33,
+	 0xc60000.0p-23,  0x800000.0p-38,
+	 0xc70000.0p-23, -0x820000.0p-33,
+	 0xc80000.0p-23,  0x800000.0p-33,
+	 0xc90000.0p-23, -0xa00000.0p-36,
+	 0xca0000.0p-23, -0xb00000.0p-33,
+	 0xcb0000.0p-23,  0x840000.0p-34,
+	 0xcc0000.0p-23, -0xd00000.0p-34,
+	 0xcd0000.0p-23,  0x800000.0p-33,
+	 0xce0000.0p-23, -0xe00000.0p-35,
+	 0xcf0000.0p-23,  0xa60000.0p-33,
+	 0xd00000.0p-23, -0x800000.0p-35,
+	 0xd10000.0p-23,  0xb40000.0p-33,
+	 0xd20000.0p-23, -0x800000.0p-35,
+	 0xd30000.0p-23,  0xaa0000.0p-33,
+	 0xd40000.0p-23, -0xe00000.0p-35,
+	 0xd50000.0p-23,  0x880000.0p-33,
+	 0xd60000.0p-23, -0xd00000.0p-34,
+	 0xd70000.0p-23,  0x9c0000.0p-34,
+	 0xd80000.0p-23, -0xb00000.0p-33,
+	 0xd90000.0p-23, -0x800000.0p-38,
+	 0xda0000.0p-23,  0xa40000.0p-33,
+	 0xdb0000.0p-23, -0xdc0000.0p-34,
+	 0xdc0000.0p-23,  0xc00000.0p-35,
+	 0xdd0000.0p-23,  0xca0000.0p-33,
+	 0xde0000.0p-23, -0xb80000.0p-34,
+	 0xdf0000.0p-23,  0xd00000.0p-35,
+	 0xe00000.0p-23,  0xc00000.0p-33,
+	 0xe10000.0p-23, -0xf40000.0p-34,
+	 0xe20000.0p-23,  0x800000.0p-37,
+	 0xe30000.0p-23,  0x860000.0p-33,
+	 0xe40000.0p-23, -0xc80000.0p-33,
+	 0xe50000.0p-23, -0xa80000.0p-34,
+	 0xe60000.0p-23,  0xe00000.0p-36,
+	 0xe70000.0p-23,  0x880000.0p-33,
+	 0xe80000.0p-23, -0xe00000.0p-33,
+	 0xe90000.0p-23, -0xfc0000.0p-34,
+	 0xea0000.0p-23, -0x800000.0p-35,
+	 0xeb0000.0p-23,  0xe80000.0p-35,
+	 0xec0000.0p-23,  0x900000.0p-33,
+	 0xed0000.0p-23,  0xe20000.0p-33,
+	 0xee0000.0p-23, -0xac0000.0p-33,
+	 0xef0000.0p-23, -0xc80000.0p-34,
+	 0xf00000.0p-23, -0x800000.0p-35,
+	 0xf10000.0p-23,  0x800000.0p-35,
+	 0xf20000.0p-23,  0xb80000.0p-34,
+	 0xf30000.0p-23,  0x940000.0p-33,
+	 0xf40000.0p-23,  0xc80000.0p-33,
+	 0xf50000.0p-23, -0xf20000.0p-33,
+	 0xf60000.0p-23, -0xc80000.0p-33,
+	 0xf70000.0p-23, -0xa20000.0p-33,
+	 0xf80000.0p-23, -0x800000.0p-33,
+	 0xf90000.0p-23, -0xc40000.0p-34,
+	 0xfa0000.0p-23, -0x900000.0p-34,
+	 0xfb0000.0p-23, -0xc80000.0p-35,
+	 0xfc0000.0p-23, -0x800000.0p-35,
+	 0xfd0000.0p-23, -0x900000.0p-36,
+	 0xfe0000.0p-23, -0x800000.0p-37,
+	 0xff0000.0p-23, -0x800000.0p-39,
+	 0x800000.0p-22,  0,
+};
+#endif /* USE_UTAB */
+
+#ifdef STRUCT_RETURN
+#define	RETURN1(rp, v) do {	\
+	(rp)->hi = (v);		\
+	(rp)->lo_set = 0;	\
+	return;			\
+} while (0)
+
+#define	RETURN2(rp, h, l) do {	\
+	(rp)->hi = (h);		\
+	(rp)->lo = (l);		\
+	(rp)->lo_set = 1;	\
+	return;			\
+} while (0)
+
+struct ld {
+	long double hi;
+	long double lo;
+	int	lo_set;
+};
+#else
+#define	RETURN1(rp, v)	RETURNF(v)
+#define	RETURN2(rp, h, l)	RETURNI((h) + (l))
+#endif
+
+#ifdef STRUCT_RETURN
+static inline __always_inline void
+k_logl(long double x, struct ld *rp)
+#else
+long double
+logl(long double x)
+#endif
+{
+	long double d, dk, val_hi, val_lo, z;
+	uint64_t ix, lx;
+	int i, k;
+	uint16_t hx;
+
+	EXTRACT_LDBL80_WORDS(hx, lx, x);
+	k = -16383;
+#if 0 /* Hard to do efficiently.  Don't do it until we support all modes. */
+	if (x == 1)
+		RETURN1(rp, 0);		/* log(1) = +0 in all rounding modes */
+#endif
+	if (hx == 0 || hx >= 0x8000) {	/* zero, negative or subnormal? */
+		if (((hx & 0x7fff) | lx) == 0)
+			RETURN1(rp, -1 / zero);	/* log(+-0) = -Inf */
+		if (hx != 0)
+			/* log(neg or [pseudo-]NaN) = qNaN: */
+			RETURN1(rp, (x - x) / zero);
+		x *= 0x1.0p65;		/* subnormal; scale up x */
+					/* including pseudo-subnormals */
+		EXTRACT_LDBL80_WORDS(hx, lx, x);
+		k = -16383 - 65;
+	} else if (hx >= 0x7fff || (lx & 0x8000000000000000ULL) == 0)
+		RETURN1(rp, x + x);	/* log(Inf or NaN) = Inf or qNaN */
+					/* log(pseudo-Inf) = qNaN */
+					/* log(pseudo-NaN) = qNaN */
+					/* log(unnormal) = qNaN */
+#ifndef STRUCT_RETURN
+	ENTERI();
+#endif
+	k += hx;
+	ix = lx & 0x7fffffffffffffffULL;
+	dk = k;
+
+	/* Scale x to be in [1, 2). */
+	SET_LDBL_EXPSIGN(x, 0x3fff);
+
+	/* 0 <= i <= INTERVALS: */
+#define	L2I	(64 - LOG2_INTERVALS)
+	i = (ix + (1LL << (L2I - 2))) >> (L2I - 1);
+
+	/*
+	 * -0.005280 < d < 0.004838.  In particular, the infinite-
+	 * precision |d| is <= 2**-7.  Rounding of G(i) to 8 bits
+	 * ensures that d is representable without extra precision for
+	 * this bound on |d| (since when this calculation is expressed
+	 * as x*G(i)-1, the multiplication needs as many extra bits as
+	 * G(i) has and the subtraction cancels 8 bits).  But for
+	 * most i (107 cases out of 129), the infinite-precision |d|
+	 * is <= 2**-8.  G(i) is rounded to 9 bits for such i to give
+	 * better accuracy (this works by improving the bound on |d|,
+	 * which in turn allows rounding to 9 bits in more cases).
+	 * This is only important when the original x is near 1 -- it
+	 * lets us avoid using a special method to give the desired
+	 * accuracy for such x.
+	 */
+	if (0)
+		d = x * G(i) - 1;
+	else {
+#ifdef USE_UTAB
+		d = (x - H(i)) * G(i) + E(i);
+#else
+		long double x_hi, x_lo;
+		float fx_hi;
+
+		/*
+		 * Split x into x_hi + x_lo to calculate x*G(i)-1 exactly.
+		 * G(i) has at most 9 bits, so the splitting point is not
+		 * critical.
+		 */
+		SET_FLOAT_WORD(fx_hi, (lx >> 40) | 0x3f800000);
+		x_hi = fx_hi;
+		x_lo = x - x_hi;
+		d = x_hi * G(i) - 1 + x_lo * G(i);
+#endif
+	}
+
+	/*
+	 * Our algorithm depends on exact cancellation of F_lo(i) and
+	 * F_hi(i) with dk*ln_2_lo and dk*ln2_hi when k is -1 and i is
+	 * at the end of the table.  This and other technical complications
+	 * make it difficult to avoid the double scaling in (dk*ln2) *
+	 * log(base) for base != e without losing more accuracy and/or
+	 * efficiency than is gained.
+	 */
+	z = d * d;
+	val_lo = z * d * z * (z * (d * P8 + P7) + (d * P6 + P5)) +
+	    (F_lo(i) + dk * ln2_lo + z * d * (d * P4 + P3)) + z * P2;
+	val_hi = d;
+#ifdef DEBUG
+	if (fetestexcept(FE_UNDERFLOW))
+		breakpoint();
+#endif
+
+	_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
+	RETURN2(rp, val_hi, val_lo);
+}
+
+long double
+log1pl(long double x)
+{
+	long double d, d_hi, d_lo, dk, f_lo, val_hi, val_lo, z;
+	long double f_hi, twopminusk;
+	uint64_t ix, lx;
+	int i, k;
+	int16_t ax, hx;
+
+	DOPRINT_START(&x);
+	EXTRACT_LDBL80_WORDS(hx, lx, x);
+	if (hx < 0x3fff) {		/* x < 1, or x neg NaN */
+		ax = hx & 0x7fff;
+		if (ax >= 0x3fff) {	/* x <= -1, or x neg NaN */
+			if (ax == 0x3fff && lx == 0x8000000000000000ULL)
+				RETURNP(-1 / zero);	/* log1p(-1) = -Inf */
+			/* log1p(x < 1, or x [pseudo-]NaN) = qNaN: */
+			RETURNP((x - x) / (x - x));
+		}
+		if (ax <= 0x3fbe) {	/* |x| < 2**-64 */
+			if ((int)x == 0)
+				RETURNP(x);	/* x with inexact if x != 0 */
+		}
+		f_hi = 1;
+		f_lo = x;
+	} else if (hx >= 0x7fff) {	/* x +Inf or non-neg NaN */
+		RETURNP(x + x);		/* log1p(Inf or NaN) = Inf or qNaN */
+					/* log1p(pseudo-Inf) = qNaN */
+					/* log1p(pseudo-NaN) = qNaN */
+					/* log1p(unnormal) = qNaN */
+	} else if (hx < 0x407f) {	/* 1 <= x < 2**128 */
+		f_hi = x;
+		f_lo = 1;
+	} else {			/* 2**128 <= x < +Inf */
+		f_hi = x;
+		f_lo = 0;		/* avoid underflow of the P5 term */
+	}
+	ENTERI();
+	x = f_hi + f_lo;
+	f_lo = (f_hi - x) + f_lo;
+
+	EXTRACT_LDBL80_WORDS(hx, lx, x);
+	k = -16383;
+
+	k += hx;
+	ix = lx & 0x7fffffffffffffffULL;
+	dk = k;
+
+	SET_LDBL_EXPSIGN(x, 0x3fff);
+	twopminusk = 1;
+	SET_LDBL_EXPSIGN(twopminusk, 0x7ffe - (hx & 0x7fff));
+	f_lo *= twopminusk;
+
+	i = (ix + (1LL << (L2I - 2))) >> (L2I - 1);
+
+	/*
+	 * x*G(i)-1 (with a reduced x) can be represented exactly, as
+	 * above, but now we need to evaluate the polynomial on d =
+	 * (x+f_lo)*G(i)-1 and extra precision is needed for that.
+	 * Since x+x_lo is a hi+lo decomposition and subtracting 1
+	 * doesn't lose too many bits, an inexact calculation for
+	 * f_lo*G(i) is good enough.
+	 */
+	if (0)
+		d_hi = x * G(i) - 1;
+	else {
+#ifdef USE_UTAB
+		d_hi = (x - H(i)) * G(i) + E(i);
+#else
+		long double x_hi, x_lo;
+		float fx_hi;
+
+		SET_FLOAT_WORD(fx_hi, (lx >> 40) | 0x3f800000);
+		x_hi = fx_hi;
+		x_lo = x - x_hi;
+		d_hi = x_hi * G(i) - 1 + x_lo * G(i);
+#endif
+	}
+	d_lo = f_lo * G(i);
+
+	/*
+	 * This is _2sumF(d_hi, d_lo) inlined.  The condition
+	 * (d_hi == 0 || |d_hi| >= |d_lo|) for using _2sumF() is not
+	 * always satisifed, so it is not clear that this works, but
+	 * it works in practice.  It works even if it gives a wrong
+	 * normalized d_lo, since |d_lo| > |d_hi| implies that i is
+	 * nonzero and d is tiny, so the F(i) term dominates d_lo.
+	 * In float precision:
+	 * (By exhaustive testing, the worst case is d_hi = 0x1.bp-25.
+	 * And if d is only a little tinier than that, we would have
+	 * another underflow problem for the P3 term; this is also ruled
+	 * out by exhaustive testing.)
+	 */
+	d = d_hi + d_lo;
+	d_lo = d_hi - d + d_lo;
+	d_hi = d;
+
+	z = d * d;
+	val_lo = z * d * z * (z * (d * P8 + P7) + (d * P6 + P5)) +
+	    (F_lo(i) + dk * ln2_lo + d_lo + z * d * (d * P4 + P3)) + z * P2;
+	val_hi = d_hi;
+#ifdef DEBUG
+	if (fetestexcept(FE_UNDERFLOW))
+		breakpoint();
+#endif
+
+	_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
+	RETURN2PI(val_hi, val_lo);
+}
+
+#ifdef STRUCT_RETURN
+
+long double
+logl(long double x)
+{
+	struct ld r;
+
+	ENTERI();
+	DOPRINT_START(&x);
+	k_logl(x, &r);
+	RETURNSPI(&r);
+}
+
+static const double
+invln10_hi =  4.3429448190317999e-1,		/*  0x1bcb7b1526e000.0p-54 */
+invln10_lo =  7.1842412889749798e-14,		/*  0x1438ca9aadd558.0p-96 */
+invln2_hi =  1.4426950408887933e0,		/*  0x171547652b8000.0p-52 */
+invln2_lo =  1.7010652264631490e-13;		/*  0x17f0bbbe87fed0.0p-95 */
+
+long double
+log10l(long double x)
+{
+	struct ld r;
+	long double hi, lo;
+
+	ENTERI();
+	DOPRINT_START(&x);
+	k_logl(x, &r);
+	if (!r.lo_set)
+		RETURNPI(r.hi);
+	_2sumF(r.hi, r.lo);
+	hi = (float)r.hi;
+	lo = r.lo + (r.hi - hi);
+	RETURN2PI(invln10_hi * hi,
+	    (invln10_lo + invln10_hi) * lo + invln10_lo * hi);
+}
+
+long double
+log2l(long double x)
+{
+	struct ld r;
+	long double hi, lo;
+
+	ENTERI();
+	DOPRINT_START(&x);
+	k_logl(x, &r);
+	if (!r.lo_set)
+		RETURNPI(r.hi);
+	_2sumF(r.hi, r.lo);
+	hi = (float)r.hi;
+	lo = r.lo + (r.hi - hi);
+	RETURN2PI(invln2_hi * hi,
+	    (invln2_lo + invln2_hi) * lo + invln2_lo * hi);
+}
+
+#endif /* STRUCT_RETURN */
diff --git a/lib/libm/man/acosh.3 b/lib/libm/man/acosh.3
index 8db27ea782..2bd787dbc9 100644
--- a/lib/libm/man/acosh.3
+++ b/lib/libm/man/acosh.3
@@ -26,14 +26,15 @@
 .\" SUCH DAMAGE.
 .\"
 .\"     from: @(#)acosh.3	5.2 (Berkeley) 5/6/91
-.\" $FreeBSD: head/lib/msun/man/acosh.3 165906 2007-01-09 01:02:06Z imp $
+.\" $FreeBSD: head/lib/msun/man/acosh.3 251599 2013-06-10 06:04:58Z das $
 .\"
-.Dd January 14, 2005
+.Dd June 11, 2013
 .Dt ACOSH 3
 .Os
 .Sh NAME
 .Nm acosh ,
-.Nm acoshf
+.Nm acoshf ,
+.Nm acoshl
 .Nd inverse hyperbolic cosine functions
 .Sh LIBRARY
 .Lb libm
@@ -43,11 +44,14 @@
 .Fn acosh "double x"
 .Ft float
 .Fn acoshf "float x"
+.Ft long double
+.Fn acoshl "long double x"
 .Sh DESCRIPTION
 The
-.Fn acosh
-and the
-.Fn acoshf
+.Fn acosh ,
+.Fn acoshf ,
+and
+.Fn acoshl
 functions compute the inverse hyperbolic cosine
 of the real
 argument
@@ -55,11 +59,7 @@ argument
 For a discussion of error due to roundoff, see
 .Xr math 3 .
 .Sh RETURN VALUES
-The
-.Fn acosh
-and the
-.Fn acoshf
-functions
+These functions
 return the inverse hyperbolic cosine of
 .Ar x .
 If the argument is less than 1,
@@ -73,6 +73,13 @@ raises an invalid exception and returns an \*(Na.
 .Xr math 3
 .Sh HISTORY
 The
-.Fn acosh
-function appeared in
-.Bx 4.3 .
+.Fn acosh ,
+.Fn acoshf ,
+and
+.Fn acoshl
+functions appeared in
+.Bx 4.3 ,
+.Fx 2.0 ,
+and
+.Dx 3.5 ,
+respectively.
diff --git a/lib/libm/man/asinh.3 b/lib/libm/man/asinh.3
index 77ac19e7ec..c08ba30aec 100644
--- a/lib/libm/man/asinh.3
+++ b/lib/libm/man/asinh.3
@@ -26,14 +26,15 @@
 .\" SUCH DAMAGE.
 .\"
 .\"     from: @(#)asinh.3	6.4 (Berkeley) 5/6/91
-.\" $FreeBSD: head/lib/msun/man/asinh.3 165906 2007-01-09 01:02:06Z imp $
+.\" $FreeBSD: head/lib/msun/man/asinh.3 251599 2013-06-10 06:04:58Z das $
 .\"
-.Dd May 6, 1991
+.Dd June 11, 2013
 .Dt ASINH 3
 .Os
 .Sh NAME
 .Nm asinh ,
-.Nm asinhf
+.Nm asinhf ,
+.Nm asinhl
 .Nd inverse hyperbolic sine functions
 .Sh LIBRARY
 .Lb libm
@@ -43,11 +44,14 @@
 .Fn asinh "double x"
 .Ft float
 .Fn asinhf "float x"
+.Ft long double
+.Fn asinhl "long double x"
 .Sh DESCRIPTION
 The
-.Fn asinh
-and the
-.Fn asinhf
+.Fn asinh ,
+.Fn asinhf ,
+and
+.Fn asinhl
 functions compute the inverse hyperbolic sine
 of the real
 argument
@@ -55,11 +59,7 @@ argument
 For a discussion of error due to roundoff, see
 .Xr math 3 .
 .Sh RETURN VALUES
-The
-.Fn asinh
-and the
-.Fn asinhf
-functions
+These functions
 return the inverse hyperbolic sine of
 .Ar x .
 .Sh SEE ALSO
@@ -69,6 +69,13 @@ return the inverse hyperbolic sine of
 .Xr math 3
 .Sh HISTORY
 The
-.Fn asinh
-function appeared in
-.Bx 4.3 .
+.Fn asinh ,
+.Fn asinhf ,
+and
+.Fn asinhl
+functions appeared in
+.Bx 4.3 ,
+.Fx 2.0 ,
+and
+.Dx 3.5 ,
+respectively.
diff --git a/lib/libm/man/atanh.3 b/lib/libm/man/atanh.3
index 9c1480a24f..bb6d967344 100644
--- a/lib/libm/man/atanh.3
+++ b/lib/libm/man/atanh.3
@@ -26,14 +26,15 @@
 .\" SUCH DAMAGE.
 .\"
 .\"     from: @(#)atanh.3	5.2 (Berkeley) 5/6/91
-.\" $FreeBSD: head/lib/msun/man/atanh.3 165906 2007-01-09 01:02:06Z imp $
+.\" $FreeBSD: head/lib/msun/man/atanh.3 251599 2013-06-10 06:04:58Z das $
 .\"
-.Dd January 14, 2005
+.Dd June 11, 2013
 .Dt ATANH 3
 .Os
 .Sh NAME
 .Nm atanh ,
-.Nm atanhf
+.Nm atanhf ,
+.Nm atanhl
 .Nd inverse hyperbolic tangent functions
 .Sh LIBRARY
 .Lb libm
@@ -43,11 +44,14 @@
 .Fn atanh "double x"
 .Ft float
 .Fn atanhf "float x"
+.Ft long double
+.Fn atanhl "long double x"
 .Sh DESCRIPTION
 The
-.Fn atanh
-and the
-.Fn atanhf
+.Fn atanh ,
+.Fn atanhf ,
+and
+.Fn atanhl
 functions compute the inverse hyperbolic tangent
 of the real
 argument
@@ -55,11 +59,7 @@ argument
 For a discussion of error due to roundoff, see
 .Xr math 3 .
 .Sh RETURN VALUES
-The
-.Fn atanh
-and the
-.Fn atanhf
-functions
+These functions
 return the inverse hyperbolic tangent of
 .Ar x
 if successful.
@@ -76,6 +76,13 @@ If
 .Xr math 3
 .Sh HISTORY
 The
-.Fn atanh
-function appeared in
-.Bx 4.3 .
+.Fn atanh ,
+.Fn atanhf ,
+and
+.Fn atanhl
+functions appeared in
+.Bx 4.3 ,
+.Fx 2.0 ,
+and
+.Dx 3.5 ,
+respectively.
diff --git a/lib/libm/man/exp.3 b/lib/libm/man/exp.3
index 143aae2a80..29af562010 100644
--- a/lib/libm/man/exp.3
+++ b/lib/libm/man/exp.3
@@ -26,9 +26,9 @@
 .\" SUCH DAMAGE.
 .\"
 .\"     from: @(#)exp.3	6.12 (Berkeley) 7/31/91
-.\" $FreeBSD: head/lib/msun/man/exp.3 238722 2012-07-23 19:13:55Z kargl $
+.\" $FreeBSD: head/lib/msun/man/exp.3 251343 2013-06-03 19:51:32Z kargl $
 .\"
-.Dd July 10, 2012
+.Dd June 11, 2013
 .Dt EXP 3
 .Os
 .Sh NAME
@@ -41,6 +41,7 @@
 .Nm exp2l ,
 .Nm expm1 ,
 .Nm expm1f ,
+.Nm expm1l ,
 .Nm pow ,
 .Nm powf
 .Nd exponential and power functions
@@ -64,6 +65,8 @@
 .Fn expm1 "double x"
 .Ft float
 .Fn expm1f "float x"
+.Ft long double
+.Fn expm1l "long double x"
 .Ft double
 .Fn pow "double x" "double y"
 .Ft float
@@ -88,9 +91,10 @@ functions compute the base 2 exponential of the given argument
 .Fa x .
 .Pp
 The
-.Fn expm1
+.Fn expm1 ,
+.Fn expm1f ,
 and the
-.Fn expm1f
+.Fn expm1l
 functions compute the value exp(x)\-1 accurately even for tiny argument
 .Fa x .
 .Pp
diff --git a/lib/libm/man/log.3 b/lib/libm/man/log.3
index c068b0a1e1..a18e73243c 100644
--- a/lib/libm/man/log.3
+++ b/lib/libm/man/log.3
@@ -22,9 +22,9 @@
 .\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 .\" SUCH DAMAGE.
 .\"
-.\" $FreeBSD: head/lib/msun/man/log.3 216211 2010-12-05 22:11:22Z das $
+.\" $FreeBSD: head/lib/msun/man/log.3 251292 2013-06-03 09:14:31Z das $
 .\"
-.Dd December 21, 2011
+.Dd June 11, 2013
 .Dt LOG 3
 .Os
 .Sh NAME
@@ -33,10 +33,13 @@
 .Nm logl ,
 .Nm log10 ,
 .Nm log10f ,
+.Nm log10l ,
 .Nm log2 ,
 .Nm log2f ,
+.Nm log2l ,
 .Nm log1p ,
-.Nm log1pf
+.Nm log1pf ,
+.Nm log1pl
 .Nd logarithm functions
 .Sh LIBRARY
 .Lb libm
@@ -46,43 +49,55 @@
 .Fn log "double x"
 .Ft float
 .Fn logf "float x"
+.Ft long double
+.Fn logl "long double x"
 .Ft double
 .Fn log10 "double x"
 .Ft float
 .Fn log10f "float x"
+.Ft long double
+.Fn log10l "long double x"
 .Ft double
 .Fn log2 "double x"
 .Ft float
 .Fn log2f "float x"
+.Ft long double
+.Fn log2l "long double x"
 .Ft double
 .Fn log1p "double x"
 .Ft float
 .Fn log1pf "float x"
+.Ft long double
+.Fn log1pl "long double x"
 .Sh DESCRIPTION
 The
-.Fn log
+.Fn log ,
+.Fn logf ,
 and
-.Fn logf
+.Fn logl
 functions compute the natural logarithm of
 .Fa x .
 .Pp
 The
-.Fn log10
+.Fn log10 ,
+.Fn log10f ,
 and
-.Fn log10f
+.Fn log10l
 functions compute the logarithm base 10 of
 .Fa x ,
 while
-.Fn log2
+.Fn log2 ,
+.Fn log2f ,
 and
-.Fn log2f
+.Fn log2l
 compute the logarithm base 2 of
 .Fa x .
 .Pp
 The
-.Fn log1p
+.Fn log1p ,
+.Fn log1pf ,
 and
-.Fn log1pf
+.Fn log1pl
 functions compute the natural logarithm of
 .No "1 + x" .
 Computing the natural logarithm as
@@ -107,12 +122,16 @@ results in an invalid exception and a return value of \*(Na.
 The
 .Fn log ,
 .Fn logf ,
+.Fn logl ,
 .Fn log10 ,
 .Fn log10f ,
+.Fn log10l ,
 .Fn log2 ,
 .Fn log2f ,
+.Fn log2l ,
 .Fn log1p ,
+.Fn log1pf ,
 and
-.Fn log1pf
+.Fn log1pl
 functions conform to
 .St -isoC-99 .
diff --git a/lib/libm/man/sin.3 b/lib/libm/man/sin.3
index d5497ac608..049cca53d7 100644
--- a/lib/libm/man/sin.3
+++ b/lib/libm/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/libm/src/catrig.c b/lib/libm/src/catrig.c
index fd1312d2af..f407f10608 100644
--- a/lib/libm/src/catrig.c
+++ b/lib/libm/src/catrig.c
@@ -23,7 +23,7 @@
  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  * SUCH DAMAGE.
  *
- * $FreeBSD: head/lib/msun/src/catrig.c 251121 2013-05-30 04:49:26Z das $
+ * $FreeBSD: head/lib/msun/src/catrig.c 251404 2013-06-05 05:33:01Z das $
  */
 
 #include <complex.h>
@@ -173,7 +173,7 @@ do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
 		 * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y).
 		 * rx = log1p(Am1 + sqrt(Am1*(A+1)))
 		 */
-		if (y == 1 && x < DBL_EPSILON*DBL_EPSILON / 128) {
+		if (y == 1 && x < DBL_EPSILON * DBL_EPSILON / 128) {
 			/*
 			 * fp is of order x^2, and fm = x/2.
 			 * A = 1 (inexactly).
@@ -192,7 +192,7 @@ do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
 			 * A = 1 (inexactly).
 			 */
 			*rx = x / sqrt((1 - y) * (1 + y));
-		} else /* if (y > 1) */ {
+		} else {		/* if (y > 1) */
 			/*
 			 * A-1 = y-1 (inexactly).
 			 */
@@ -252,7 +252,7 @@ do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
 			*sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y /
 				sqrt((y + 1) * (y - 1));
 			*new_y = y * (4 / DBL_EPSILON / DBL_EPSILON);
-		} else /* if (y < 1) */ {
+		} else {		/* if (y < 1) */
 			/*
 			 * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and
 			 * A = 1 (inexactly).
@@ -297,7 +297,6 @@ casinh(double complex z)
 		 * C99 leaves it optional whether to raise invalid if one of
 		 * the arguments is not NaN, so we opt not to raise it.
 		 */
-		/* Bruce Evans tells me this is the way to do this: */
 		return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
 	}
 
@@ -336,6 +335,7 @@ double complex
 casin(double complex z)
 {
 	double complex w = casinh(cpack(cimag(z), creal(z)));
+
 	return (cpack(cimag(w), creal(w)));
 }
 
@@ -401,17 +401,17 @@ cacos(double complex z)
 	/* All remaining cases are inexact. */
 	raise_inexact();
 
-	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON/4)
+	if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
 		return (cpack(pio2_hi - (x - pio2_lo), -y));
 
 	do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
 	if (B_is_usable) {
-		if (sx==0)
+		if (sx == 0)
 			rx = acos(B);
 		else
 			rx = acos(-B);
 	} else {
-		if (sx==0)
+		if (sx == 0)
 			rx = atan2(sqrt_A2mx2, new_x);
 		else
 			rx = atan2(sqrt_A2mx2, -new_x);
@@ -486,10 +486,6 @@ clog_for_large_values(double complex z)
 	return (cpack(log(ax * ax + ay * ay) / 2, atan2(y, x)));
 }
 
-/*
- *=============================================================================
- */
-
 /*
  *				=================
  *				| catanh, catan |
@@ -510,6 +506,7 @@ sum_squares(double x, double y)
 	/* Avoid underflow when y is small. */
 	if (y < SQRT_MIN)
 		return (x * x);
+
 	return (x * x + y * y);
 }
 
@@ -588,10 +585,9 @@ catanh(double complex z)
 		if (isinf(x))
 			return (cpack(copysign(0, x), y + y));
 		/* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
-		if (isinf(y)) {
+		if (isinf(y))
 			return (cpack(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
@@ -600,10 +596,9 @@ catanh(double complex z)
 		return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
 	}
 
-	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
 		return (cpack(real_part_reciprocal(x, y),
 			      copysign(pio2_hi + pio2_lo, y)));
-	}
 
 	if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
 		/*
@@ -616,7 +611,7 @@ catanh(double complex z)
 	}
 
 	if (ax == 1 && ay < DBL_EPSILON)
-		rx = (log(ay) - m_ln2) / -2;
+		rx = (m_ln2 - log(ay)) / 2;
 	else
 		rx = log1p(4 * ax / sum_squares(ax - 1, ay)) / 4;
 
@@ -638,5 +633,6 @@ double complex
 catan(double complex z)
 {
 	double complex w = catanh(cpack(cimag(z), creal(z)));
+
 	return (cpack(cimag(w), creal(w)));
 }
diff --git a/lib/libm/src/catrigf.c b/lib/libm/src/catrigf.c
index f44a319bad..5c787c91ef 100644
--- a/lib/libm/src/catrigf.c
+++ b/lib/libm/src/catrigf.c
@@ -23,7 +23,7 @@
  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  * SUCH DAMAGE.
  *
- * $FreeBSD: head/lib/msun/src/catrigf.c 251121 2013-05-30 04:49:26Z das $
+ * $FreeBSD: head/lib/msun/src/catrigf.c 251404 2013-06-05 05:33:01Z das $
  */
 
 /*
@@ -33,7 +33,11 @@
  * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
  * http://dl.acm.org/citation.cfm?id=275324.
  *
- * The code for catrig.c contains complete comments.
+ * See catrig.c for complete comments.
+ *
+ * XXX comments were removed automatically, and even short ones on the right
+ * of statements were removed (all of them), contrary to normal style.  Only
+ * a few comments on the right of declarations remain.
  */
 
 #include <complex.h>
@@ -100,7 +104,7 @@ do_hard_work(float x, float y, float *rx, int *B_is_usable, float *B,
 			Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
 			*rx = log1pf(Am1 + sqrtf(Am1 * (A + 1)));
 		} else if (y < 1) {
-			*rx = x / sqrtf((1 - y)*(1 + y));
+			*rx = x / sqrtf((1 - y) * (1 + y));
 		} else {
 			*rx = log1pf((y - 1) + sqrtf((y - 1) * (y + 1)));
 		}
@@ -188,6 +192,7 @@ float complex
 casinf(float complex z)
 {
 	float complex w = casinhf(cpackf(cimagf(z), crealf(z)));
+
 	return (cpackf(cimagf(w), crealf(w)));
 }
 
@@ -211,7 +216,8 @@ cacosf(float complex z)
 			return (cpackf(y + y, -INFINITY));
 		if (isinf(y))
 			return (cpackf(x + x, -y));
-		if (x == 0) return (cpackf(pio2_hi + pio2_lo, y + y));
+		if (x == 0)
+			return (cpackf(pio2_hi + pio2_lo, y + y));
 		return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
 	}
 
@@ -234,17 +240,17 @@ cacosf(float complex z)
 
 	do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
 	if (B_is_usable) {
-		if (sx==0)
+		if (sx == 0)
 			rx = acosf(B);
 		else
 			rx = acosf(-B);
 	} else {
-		if (sx==0)
+		if (sx == 0)
 			rx = atan2f(sqrt_A2mx2, new_x);
 		else
 			rx = atan2f(sqrt_A2mx2, -new_x);
 	}
-	if (sy==0)
+	if (sy == 0)
 		ry = -ry;
 	return (cpackf(rx, ry));
 }
@@ -283,10 +289,9 @@ clog_for_large_values(float complex z)
 		ay = t;
 	}
 
-	if (ax > FLT_MAX / 2) {
+	if (ax > FLT_MAX / 2)
 		return (cpackf(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)));
@@ -299,8 +304,9 @@ sum_squares(float x, float y)
 {
 
 	if (y < SQRT_MIN)
-		return (x*x);
-	return (x*x + y*y);
+		return (x * x);
+
+	return (x * x + y * y);
 }
 
 static inline float
@@ -317,9 +323,9 @@ real_part_reciprocal(float x, float y)
 #define	BIAS	(FLT_MAX_EXP - 1)
 #define	CUTOFF	(FLT_MANT_DIG / 2 + 1)
 	if (ix - iy >= CUTOFF << 23 || isinf(x))
-		return (1/x);
+		return (1 / x);
 	if (iy - ix >= CUTOFF << 23)
-		return (x/y/y);
+		return (x / y / y);
 	if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23)
 		return (x / (x * x + y * y));
 	SET_FLOAT_WORD(scale, 0x7f800000 - ix);
@@ -346,18 +352,16 @@ catanhf(float complex z)
 
 	if (isnan(x) || isnan(y)) {
 		if (isinf(x))
-			return (cpackf(copysignf(0, x), y+y));
-		if (isinf(y)) {
+			return (cpackf(copysignf(0, x), y + y));
+		if (isinf(y))
 			return (cpackf(copysignf(0, x),
 				       copysignf(pio2_hi + pio2_lo, y)));
-		}
 		return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
 	}
 
-	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
+	if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
 		return (cpackf(real_part_reciprocal(x, y),
 			       copysignf(pio2_hi + pio2_lo, y)));
-	}
 
 	if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
 		raise_inexact();
@@ -365,7 +369,7 @@ catanhf(float complex z)
 	}
 
 	if (ax == 1 && ay < FLT_EPSILON)
-		rx = (logf(ay) - m_ln2) / -2;
+		rx = (m_ln2 - logf(ay)) / 2;
 	else
 		rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4;
 
@@ -383,5 +387,6 @@ float complex
 catanf(float complex z)
 {
 	float complex w = catanhf(cpackf(cimagf(z), crealf(z)));
+
 	return (cpackf(cimagf(w), crealf(w)));
 }
diff --git a/lib/libm/src/e_acosh.c b/lib/libm/src/e_acosh.c
index a14e09cb03..8de9689ea3 100644
--- a/lib/libm/src/e_acosh.c
+++ b/lib/libm/src/e_acosh.c
@@ -1,6 +1,6 @@
 
 /* @(#)e_acosh.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das $ */
+/* $FreeBSD: head/lib/msun/src/e_acosh.c 251599 2013-06-10 06:04:58Z das $
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -10,6 +10,7 @@
  * software is freely granted, provided that this notice 
  * is preserved.
  * ====================================================
+ *
  */
 
 /* __ieee754_acosh(x)
@@ -26,6 +27,8 @@
  *	acosh(NaN) is NaN without signal.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -57,3 +60,7 @@ __ieee754_acosh(double x)
 	    return log1p(t+sqrt(2.0*t+t*t));
 	}
 }
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(acosh, acoshl);
+#endif
diff --git a/lib/libm/src/e_acoshl.c b/lib/libm/src/e_acoshl.c
new file mode 100644
index 0000000000..e25cd9c2be
--- /dev/null
+++ b/lib/libm/src/e_acoshl.c
@@ -0,0 +1,87 @@
+/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)e_acosh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ * $FreeBSD: head/lib/msun/src/e_acoshl.c 251599 2013-06-10 06:04:58Z das $
+ */
+
+/*
+ * See e_acosh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_LARGE is the threshold above which we use acosh(x) ~= log(2x). */
+#if LDBL_MANT_DIG == 64
+#define	EXP_LARGE	34
+#elif LDBL_MANT_DIG == 113
+#define	EXP_LARGE	58
+#else
+#error "Unsupported long double format"
+#endif
+
+#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 double
+one	= 1.0;
+
+#if LDBL_MANT_DIG == 64
+static const union IEEEl2bits
+u_ln2 =  LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
+#define	ln2	u_ln2.e
+#elif LDBL_MANT_DIG == 113
+static const long double
+ln2 =  6.93147180559945309417232121458176568e-1L;	/* 0x162e42fefa39ef35793c7673007e6.0p-113 */
+#else
+#error "Unsupported long double format"
+#endif
+
+long double
+acoshl(long double x)
+{
+	long double t;
+	int16_t hx;
+
+	ENTERI();
+	GET_LDBL_EXPSIGN(hx, x);
+	if (hx < 0x3fff) {		/* x < 1, or misnormal */
+	    RETURNI((x-x)/(x-x));
+	} else if (hx >= BIAS + EXP_LARGE) { /* x >= LARGE */
+	    if (hx >= 0x7fff) {		/* x is inf, NaN or misnormal */
+	        RETURNI(x+x);
+	    } else 
+		RETURNI(logl(x)+ln2);	/* acosh(huge)=log(2x), or misnormal */
+	} else if (hx == 0x3fff && x == 1) {
+	    RETURNI(0.0);		/* acosh(1) = 0 */
+	} else if (hx >= 0x4000) {	/* LARGE > x >= 2, or misnormal */
+	    t=x*x;
+	    RETURNI(logl(2.0*x-one/(x+sqrtl(t-one))));
+	} else {			/* 1<x<2 */
+	    t = x-one;
+	    RETURNI(log1pl(t+sqrtl(2.0*t+t*t)));
+	}
+}
diff --git a/lib/libm/src/e_atanh.c b/lib/libm/src/e_atanh.c
index c91f1add14..12a6f87675 100644
--- a/lib/libm/src/e_atanh.c
+++ b/lib/libm/src/e_atanh.c
@@ -1,6 +1,6 @@
 
 /* @(#)e_atanh.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das $ */
+/* $FreeBSD: head/lib/msun/src/e_atanh.c 251599 2013-06-10 06:04:58Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -31,6 +31,8 @@
  *
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -58,3 +60,7 @@ __ieee754_atanh(double x)
 	    t = 0.5*log1p((x+x)/(one-x));
 	if(hx>=0) return t; else return -t;
 }
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(atanh, atanhl);
+#endif
diff --git a/lib/libm/src/e_atanhl.c b/lib/libm/src/e_atanhl.c
new file mode 100644
index 0000000000..be26e0c186
--- /dev/null
+++ b/lib/libm/src/e_atanhl.c
@@ -0,0 +1,72 @@
+/* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)e_atanh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ * $FreeBSD: head/lib/msun/src/e_atanhl.c 251599 2013-06-10 06:04:58Z das $
+ */
+
+/*
+ * See e_atanh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_TINY is the threshold below which we use atanh(x) ~= x. */
+#if LDBL_MANT_DIG == 64
+#define	EXP_TINY	-34
+#elif LDBL_MANT_DIG == 113
+#define	EXP_TINY	-58
+#else
+#error "Unsupported long double format"
+#endif
+
+#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 double one = 1.0, huge = 1e300;
+static const double zero = 0.0;
+
+long double
+atanhl(long double x)
+{
+	long double t;
+	uint16_t hx, ix;
+
+	ENTERI();
+	GET_LDBL_EXPSIGN(hx, x);
+	ix = hx & 0x7fff;
+	if (ix >= 0x3fff)		/* |x| >= 1, or NaN or misnormal */
+	    RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x));
+	if (ix < BIAS + EXP_TINY && (huge + x) > zero)
+	    RETURNI(x);			/* x is tiny */
+	SET_LDBL_EXPSIGN(x, ix);
+	if (ix < 0x3ffe) {		/* |x| < 0.5, or misnormal */
+	    t = x+x;
+	    t = 0.5*log1pl(t+t*x/(one-x));
+	} else 
+	    t = 0.5*log1pl((x+x)/(one-x));
+	RETURNI((hx & 0x8000) == 0 ? t : -t);
+}
diff --git a/lib/libm/src/e_log.c b/lib/libm/src/e_log.c
index dd875c1065..2fc3f0695d 100644
--- a/lib/libm/src/e_log.c
+++ b/lib/libm/src/e_log.c
@@ -1,6 +1,6 @@
 
 /* @(#)e_log.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_log.c 251024 2013-05-27 08:50:10Z das $ */
+/* $FreeBSD: head/lib/msun/src/e_log.c 251292 2013-06-03 09:14:31Z das $
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -63,6 +63,8 @@
  * to produce the hexadecimal values shown.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -137,3 +139,7 @@ __ieee754_log(double x)
 		     return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
 	}
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log, logl);
+#endif
diff --git a/lib/libm/src/e_log10.c b/lib/libm/src/e_log10.c
index d744fb5b53..ccb6373a78 100644
--- a/lib/libm/src/e_log10.c
+++ b/lib/libm/src/e_log10.c
@@ -1,6 +1,6 @@
 
 /* @(#)e_log10.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_log10.c 251024 2013-05-27 08:50:10Z das $ */
+/* $FreeBSD: head/lib/msun/src/e_log10.c 251292 2013-06-03 09:14:31Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -20,6 +20,8 @@
  * in not-quite-routine extra precision.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 #include "k_log.h"
@@ -84,3 +86,7 @@ __ieee754_log10(double x)
 
 	return val_lo + val_hi;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log10, log10l);
+#endif
diff --git a/lib/libm/src/e_log2.c b/lib/libm/src/e_log2.c
index a7388e5ce7..080e46bbba 100644
--- a/lib/libm/src/e_log2.c
+++ b/lib/libm/src/e_log2.c
@@ -1,6 +1,6 @@
 
 /* @(#)e_log10.c 1.3 95/01/18 */
-/* $FreeBSD: head/lib/msun/src/e_log2.c 251024 2013-05-27 08:50:10Z das $ */
+/* $FreeBSD: head/lib/msun/src/e_log2.c 251404 2013-06-05 05:33:01Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -22,6 +22,8 @@
  * in not-quite-routine extra precision.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 #include "k_log.h"
@@ -107,3 +109,7 @@ __ieee754_log2(double x)
 
 	return val_lo + val_hi;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log2, log2l);
+#endif
diff --git a/lib/libm/src/math.h b/lib/libm/src/math.h
index 3768ac35b5..f9befc7f28 100644
--- a/lib/libm/src/math.h
+++ b/lib/libm/src/math.h
@@ -11,7 +11,7 @@
 
 /*
  * from: @(#)fdlibm.h 5.1 93/09/24
- * $FreeBSD: head/lib/msun/src/math.h 238722 2012-07-23 19:13:55Z kargl $
+ * $FreeBSD: head/lib/msun/src/math.h 251599 2013-06-10 06:04:58Z das $
  */
 
 #ifndef _MATH_H_
@@ -396,9 +396,12 @@ float	significandf(float);
  * long double versions of ISO/POSIX math functions
  */
 #if __ISO_C_VISIBLE >= 1999
+long double	acoshl(long double);
 long double	acosl(long double);
+long double	asinhl(long double);
 long double	asinl(long double);
 long double	atan2l(long double, long double);
+long double	atanhl(long double);
 long double	atanl(long double);
 long double	cbrtl(long double);
 long double	ceill(long double);
@@ -406,6 +409,7 @@ long double	copysignl(long double, long double) __pure2;
 long double	cosl(long double);
 long double	exp2l(long double);
 long double	expl(long double);
+long double	expm1l(long double);
 long double	fabsl(long double) __pure2;
 long double	fdiml(long double, long double);
 long double	floorl(long double);
@@ -419,7 +423,11 @@ int		ilogbl(long double) __pure2;
 long double	ldexpl(long double, int);
 long long	llrintl(long double);
 long long	llroundl(long double);
+long double	log10l(long double);
+long double	log1pl(long double);
+long double	log2l(long double);
 long double	logbl(long double);
+long double	logl(long double);
 long		lrintl(long double);
 long		lroundl(long double);
 long double	modfl(long double, long double *); /* fundamentally !__pure2 */
@@ -457,18 +465,10 @@ __END_DECLS
  */
 __BEGIN_DECLS
 
-long double	acoshl(long double);
-long double	asinhl(long double);
-long double	atanhl(long double);
 long double	coshl(long double);
 long double	erfcl(long double);
 long double	erfl(long double);
-long double	expm1l(long double);
 long double	lgammal(long double);
-long double	log10l(long double);
-long double	log1pl(long double);
-long double	log2l(long double);
-long double	logl(long double);
 long double	powl(long double, long double);
 long double	sinhl(long double);
 long double	tanhl(long double);
diff --git a/lib/libm/src/math_private.h b/lib/libm/src/math_private.h
index 34e8f5967a..d4623aff07 100644
--- a/lib/libm/src/math_private.h
+++ b/lib/libm/src/math_private.h
@@ -11,7 +11,7 @@
 
 /*
  * from: @(#)fdlibm.h 5.1 93/09/24
- * $FreeBSD: head/lib/msun/src/math_private.h 241051 2012-09-29 16:40:12Z kargl $
+ * $FreeBSD: head/lib/msun/src/math_private.h 251292 2013-06-03 09:14:31Z das $
  */
 
 #ifndef _MATH_PRIVATE_H_
@@ -180,6 +180,33 @@ do {								\
   (d) = sf_u.value;						\
 } while (0)
 
+/*
+ * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long
+ * double.
+ */
+
+#define	EXTRACT_LDBL80_WORDS(ix0,ix1,d)				\
+do {								\
+  union IEEEl2bits ew_u;					\
+  ew_u.e = (d);							\
+  (ix0) = ew_u.xbits.expsign;					\
+  (ix1) = ew_u.xbits.man;					\
+} while (0)
+
+/*
+ * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit
+ * long double.
+ */
+
+#define	EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d)			\
+do {								\
+  union IEEEl2bits ew_u;					\
+  ew_u.e = (d);							\
+  (ix0) = ew_u.xbits.expsign;					\
+  (ix1) = ew_u.xbits.manh;					\
+  (ix2) = ew_u.xbits.manl;					\
+} while (0)
+
 /* Get expsign as a 16 bit int from a long double.  */
 
 #define	GET_LDBL_EXPSIGN(i,d)					\
@@ -189,6 +216,33 @@ do {								\
   (i) = ge_u.xbits.expsign;					\
 } while (0)
 
+/*
+ * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int
+ * mantissa.
+ */
+
+#define	INSERT_LDBL80_WORDS(d,ix0,ix1)				\
+do {								\
+  union IEEEl2bits iw_u;					\
+  iw_u.xbits.expsign = (ix0);					\
+  iw_u.xbits.man = (ix1);					\
+  (d) = iw_u.e;							\
+} while (0)
+
+/*
+ * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints
+ * comprising the mantissa.
+ */
+
+#define	INSERT_LDBL128_WORDS(d,ix0,ix1,ix2)			\
+do {								\
+  union IEEEl2bits iw_u;					\
+  iw_u.xbits.expsign = (ix0);					\
+  iw_u.xbits.manh = (ix1);					\
+  iw_u.xbits.manl = (ix2);					\
+  (d) = iw_u.e;							\
+} while (0)
+
 /* Set expsign of a long double from a 16 bit int.  */
 
 #define	SET_LDBL_EXPSIGN(d,v)					\
@@ -252,6 +306,110 @@ do {								\
 /* Default return statement if hack*_t() is not used. */
 #define      RETURNF(v)      return (v)
 
+/*
+ * 2sum gives the same result as 2sumF without requiring |a| >= |b| or
+ * a == 0, but is slower.
+ */
+#define	_2sum(a, b) do {	\
+	__typeof(a) __s, __w;	\
+				\
+	__w = (a) + (b);	\
+	__s = __w - (a);	\
+	(b) = ((a) - (__w - __s)) + ((b) - __s); \
+	(a) = __w;		\
+} while (0)
+
+/*
+ * 2sumF algorithm.
+ *
+ * "Normalize" the terms in the infinite-precision expression a + b for
+ * the sum of 2 floating point values so that b is as small as possible
+ * relative to 'a'.  (The resulting 'a' is the value of the expression in
+ * the same precision as 'a' and the resulting b is the rounding error.)
+ * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and
+ * exponent overflow or underflow must not occur.  This uses a Theorem of
+ * Dekker (1971).  See Knuth (1981) 4.2.2 Theorem C.  The name "TwoSum"
+ * is apparently due to Skewchuk (1997).
+ *
+ * For this to always work, assignment of a + b to 'a' must not retain any
+ * extra precision in a + b.  This is required by C standards but broken
+ * in many compilers.  The brokenness cannot be worked around using
+ * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this
+ * algorithm would be destroyed by non-null strict assignments.  (The
+ * compilers are correct to be broken -- the efficiency of all floating
+ * point code calculations would be destroyed similarly if they forced the
+ * conversions.)
+ *
+ * Fortunately, a case that works well can usually be arranged by building
+ * any extra precision into the type of 'a' -- 'a' should have type float_t,
+ * double_t or long double.  b's type should be no larger than 'a's type.
+ * Callers should use these types with scopes as large as possible, to
+ * reduce their own extra-precision and efficiciency problems.  In
+ * particular, they shouldn't convert back and forth just to call here.
+ */
+#ifdef DEBUG
+#define	_2sumF(a, b) do {				\
+	__typeof(a) __w;				\
+	volatile __typeof(a) __ia, __ib, __r, __vw;	\
+							\
+	__ia = (a);					\
+	__ib = (b);					\
+	assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib));	\
+							\
+	__w = (a) + (b);				\
+	(b) = ((a) - __w) + (b);			\
+	(a) = __w;					\
+							\
+	/* The next 2 assertions are weak if (a) is already long double. */ \
+	assert((long double)__ia + __ib == (long double)(a) + (b));	\
+	__vw = __ia + __ib;				\
+	__r = __ia - __vw;				\
+	__r += __ib;					\
+	assert(__vw == (a) && __r == (b));		\
+} while (0)
+#else /* !DEBUG */
+#define	_2sumF(a, b) do {	\
+	__typeof(a) __w;	\
+				\
+	__w = (a) + (b);	\
+	(b) = ((a) - __w) + (b); \
+	(a) = __w;		\
+} while (0)
+#endif /* DEBUG */
+
+/*
+ * Set x += c, where x is represented in extra precision as a + b.
+ * x must be sufficiently normalized and sufficiently larger than c,
+ * and the result is then sufficiently normalized.
+ *
+ * The details of ordering are that |a| must be >= |c| (so that (a, c)
+ * can be normalized without extra work to swap 'a' with c).  The details of
+ * the normalization are that b must be small relative to the normalized 'a'.
+ * Normalization of (a, c) makes the normalized c tiny relative to the
+ * normalized a, so b remains small relative to 'a' in the result.  However,
+ * b need not ever be tiny relative to 'a'.  For example, b might be about
+ * 2**20 times smaller than 'a' to give about 20 extra bits of precision.
+ * That is usually enough, and adding c (which by normalization is about
+ * 2**53 times smaller than a) cannot change b significantly.  However,
+ * cancellation of 'a' with c in normalization of (a, c) may reduce 'a'
+ * significantly relative to b.  The caller must ensure that significant
+ * cancellation doesn't occur, either by having c of the same sign as 'a',
+ * or by having |c| a few percent smaller than |a|.  Pre-normalization of
+ * (a, b) may help.
+ *
+ * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2
+ * exercise 19).  We gain considerable efficiency by requiring the terms to
+ * be sufficiently normalized and sufficiently increasing.
+ */
+#define	_3sumF(a, b, c) do {	\
+	__typeof(a) __tmp;	\
+				\
+	__tmp = (c);		\
+	_2sumF(__tmp, (a));	\
+	(b) += (a);		\
+	(a) = __tmp;		\
+} while (0)
+
 /*
  * Common routine to process the arguments to nan(), nanf(), and nanl().
  */
@@ -348,6 +506,7 @@ irint(double x)
 #define	HAVE_EFFICIENT_IRINT
 #endif
 
+#if defined(__x86_64__) || defined(__i386__)
 static __inline int
 irintl(long double x)
 {
@@ -357,9 +516,144 @@ irintl(long double x)
 	return (n);
 }
 #define	HAVE_EFFICIENT_IRINTL
+#endif
 
 #endif /* __GNUC__ */
 
+#ifdef DEBUG
+#if defined(__amd64__) || defined(__i386__)
+#define	breakpoint()	asm("int $3")
+#else
+#include <signal.h>
+
+#define	breakpoint()	raise(SIGTRAP)
+#endif
+#endif
+
+/* Write a pari script to test things externally. */
+#ifdef DOPRINT
+#include <stdio.h>
+
+#ifndef DOPRINT_SWIZZLE
+#define	DOPRINT_SWIZZLE		0
+#endif
+
+#ifdef DOPRINT_LD80
+
+#define	DOPRINT_START(xp) do {						\
+	uint64_t __lx;							\
+	uint16_t __hx;							\
+									\
+	/* Hack to give more-problematic args. */			\
+	EXTRACT_LDBL80_WORDS(__hx, __lx, *xp);				\
+	__lx ^= DOPRINT_SWIZZLE;					\
+	INSERT_LDBL80_WORDS(*xp, __hx, __lx);				\
+	printf("x = %.21Lg; ", (long double)*xp);			\
+} while (0)
+#define	DOPRINT_END1(v)							\
+	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define	DOPRINT_END2(hi, lo)						\
+	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
+	    (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_D64)
+
+#define	DOPRINT_START(xp) do {						\
+	uint32_t __hx, __lx;						\
+									\
+	EXTRACT_WORDS(__hx, __lx, *xp);					\
+	__lx ^= DOPRINT_SWIZZLE;					\
+	INSERT_WORDS(*xp, __hx, __lx);					\
+	printf("x = %.21Lg; ", (long double)*xp);			\
+} while (0)
+#define	DOPRINT_END1(v)							\
+	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define	DOPRINT_END2(hi, lo)						\
+	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
+	    (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_F32)
+
+#define	DOPRINT_START(xp) do {						\
+	uint32_t __hx;							\
+									\
+	GET_FLOAT_WORD(__hx, *xp);					\
+	__hx ^= DOPRINT_SWIZZLE;					\
+	SET_FLOAT_WORD(*xp, __hx);					\
+	printf("x = %.21Lg; ", (long double)*xp);			\
+} while (0)
+#define	DOPRINT_END1(v)							\
+	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define	DOPRINT_END2(hi, lo)						\
+	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
+	    (long double)(hi), (long double)(lo))
+
+#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
+
+#ifndef DOPRINT_SWIZZLE_HIGH
+#define	DOPRINT_SWIZZLE_HIGH	0
+#endif
+
+#define	DOPRINT_START(xp) do {						\
+	uint64_t __lx, __llx;						\
+	uint16_t __hx;							\
+									\
+	EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp);			\
+	__llx ^= DOPRINT_SWIZZLE;					\
+	__lx ^= DOPRINT_SWIZZLE_HIGH;					\
+	INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx);			\
+	printf("x = %.36Lg; ", (long double)*xp);					\
+} while (0)
+#define	DOPRINT_END1(v)							\
+	printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define	DOPRINT_END2(hi, lo)						\
+	printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n",		\
+	    (long double)(hi), (long double)(lo))
+
+#endif /* DOPRINT_LD80 */
+
+#else /* !DOPRINT */
+#define	DOPRINT_START(xp)
+#define	DOPRINT_END1(v)
+#define	DOPRINT_END2(hi, lo)
+#endif /* DOPRINT */
+
+#define	RETURNP(x) do {			\
+	DOPRINT_END1(x);		\
+	RETURNF(x);			\
+} while (0)
+#define	RETURNPI(x) do {		\
+	DOPRINT_END1(x);		\
+	RETURNI(x);			\
+} while (0)
+#define	RETURN2P(x, y) do {		\
+	DOPRINT_END2((x), (y));		\
+	RETURNF((x) + (y));		\
+} while (0)
+#define	RETURN2PI(x, y) do {		\
+	DOPRINT_END2((x), (y));		\
+	RETURNI((x) + (y));		\
+} while (0)
+#ifdef STRUCT_RETURN
+#define	RETURNSP(rp) do {		\
+	if (!(rp)->lo_set)		\
+		RETURNP((rp)->hi);	\
+	RETURN2P((rp)->hi, (rp)->lo);	\
+} while (0)
+#define	RETURNSPI(rp) do {		\
+	if (!(rp)->lo_set)		\
+		RETURNPI((rp)->hi);	\
+	RETURN2PI((rp)->hi, (rp)->lo);	\
+} while (0)
+#endif
+#define	SUM2P(x, y) ({			\
+	const __typeof (x) __x = (x);	\
+	const __typeof (y) __y = (y);	\
+					\
+	DOPRINT_END2(__x, __y);		\
+	__x + __y;			\
+})
+
 /*
  * ieee style elementary functions
  *
diff --git a/lib/libm/src/s_asinh.c b/lib/libm/src/s_asinh.c
index 50720bc6f3..38ac00e342 100644
--- a/lib/libm/src/s_asinh.c
+++ b/lib/libm/src/s_asinh.c
@@ -1,5 +1,5 @@
 /* @(#)s_asinh.c 5.1 93/09/24 */
-/* $FreeBSD: head/lib/msun/src/s_asinh.c 176451 2008-02-22 02:30:36Z das $ */
+/* $FreeBSD: head/lib/msun/src/s_asinh.c 251599 2013-06-10 06:04:58Z das $
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -22,6 +22,8 @@
  *		 := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -52,3 +54,7 @@ asinh(double x)
 	}
 	if(hx>0) return w; else return -w;
 }
+
+#if LDBL_MANT_DIG == 53
+__weak_reference(asinh, asinhl);
+#endif
diff --git a/lib/libm/src/s_asinhl.c b/lib/libm/src/s_asinhl.c
new file mode 100644
index 0000000000..d6577b85e2
--- /dev/null
+++ b/lib/libm/src/s_asinhl.c
@@ -0,0 +1,90 @@
+/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */
+
+/* @(#)s_asinh.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.
+ * ====================================================
+ *
+ * $FreeBSD: head/lib/msun/src/s_asinhl.c 251599 2013-06-10 06:04:58Z das $
+ */
+
+/*
+ * See s_asinh.c for complete comments.
+ *
+ * Converted to long double by David Schultz <das@FreeBSD.ORG> and
+ * Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* EXP_LARGE is the threshold above which we use asinh(x) ~= log(2x). */
+/* EXP_TINY is the threshold below which we use asinh(x) ~= x. */
+#if LDBL_MANT_DIG == 64
+#define	EXP_LARGE	34
+#define	EXP_TINY	-34
+#elif LDBL_MANT_DIG == 113
+#define	EXP_LARGE	58
+#define	EXP_TINY	-58
+#else
+#error "Unsupported long double format"
+#endif
+
+#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 double
+one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+huge=  1.00000000000000000000e+300;
+
+#if LDBL_MANT_DIG == 64
+static const union IEEEl2bits
+u_ln2 =  LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
+#define	ln2	u_ln2.e
+#elif LDBL_MANT_DIG == 113
+static const long double
+ln2 =  6.93147180559945309417232121458176568e-1L;	/* 0x162e42fefa39ef35793c7673007e6.0p-113 */
+#else
+#error "Unsupported long double format"
+#endif
+
+long double
+asinhl(long double x)
+{
+	long double t, w;
+	uint16_t hx, ix;
+
+	ENTERI();
+	GET_LDBL_EXPSIGN(hx, x);
+	ix = hx & 0x7fff;
+	if (ix >= 0x7fff) RETURNI(x+x);	/* x is inf, NaN or misnormal */
+	if (ix < BIAS + EXP_TINY) {	/* |x| < TINY, or misnormal */
+	    if (huge + x > one) RETURNI(x);	/* return x inexact except 0 */
+	}
+	if (ix >= BIAS + EXP_LARGE) {	/* |x| >= LARGE, or misnormal */
+	    w = logl(fabsl(x))+ln2;
+	} else if (ix >= 0x4000) {	/* LARGE > |x| >= 2.0, or misnormal */
+	    t = fabsl(x);
+	    w = logl(2.0*t+one/(sqrtl(x*x+one)+t));
+	} else {		/* 2.0 > |x| >= TINY, or misnormal */
+	    t = x*x;
+	    w =log1pl(fabsl(x)+t/(one+sqrtl(one+t)));
+	}
+	RETURNI((hx & 0x8000) == 0 ? w : -w);
+}
diff --git a/lib/libm/src/s_expm1.c b/lib/libm/src/s_expm1.c
index dfae2af85b..89ac6742c9 100644
--- a/lib/libm/src/s_expm1.c
+++ b/lib/libm/src/s_expm1.c
@@ -1,5 +1,5 @@
 /* @(#)s_expm1.c 5.1 93/09/24 */
-/* $FreeBSD: head/lib/msun/src/s_expm1.c 251024 2013-05-27 08:50:10Z das $ */
+/* $FreeBSD: head/lib/msun/src/s_expm1.c 251343 2013-06-03 19:51:32Z kargl $
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -214,3 +214,7 @@ expm1(double x)
 	}
 	return y;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(expm1, expm1l);
+#endif
diff --git a/lib/libm/src/s_log1p.c b/lib/libm/src/s_log1p.c
index 7d68280be5..7cf2a406b1 100644
--- a/lib/libm/src/s_log1p.c
+++ b/lib/libm/src/s_log1p.c
@@ -1,5 +1,5 @@
 /* @(#)s_log1p.c 5.1 93/09/24 */
-/* $FreeBSD: head/lib/msun/src/s_log1p.c 251024 2013-05-27 08:50:10Z das $ */
+/* $FreeBSD: head/lib/msun/src/s_log1p.c 251292 2013-06-03 09:14:31Z das $ */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -172,3 +172,7 @@ log1p(double x)
 	if(k==0) return f-(hfsq-s*(hfsq+R)); else
 		 return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log1p, log1pl);
+#endif
diff --git a/lib/libm/x86_64/Makefile.inc b/lib/libm/x86_64/Makefile.inc
index abd2e850fa..b69d369e81 100644
--- a/lib/libm/x86_64/Makefile.inc
+++ b/lib/libm/x86_64/Makefile.inc
@@ -6,3 +6,5 @@ ARCH_SRCS = e_remainder.S e_remainderf.S e_remainderl.S \
 	    s_logbl.S s_lrint.S s_lrintf.S s_lrintl.S \
 	    s_remquo.S s_remquof.S s_remquol.S \
 	    s_rintl.S s_scalbn.S s_scalbnf.S s_scalbnl.S
+
+SYM_MAPS += ${.CURDIR}/x86_64/Symbol.map
diff --git a/lib/libm/x86_64/Symbol.map b/lib/libm/x86_64/Symbol.map
new file mode 100644
index 0000000000..e889f36652
--- /dev/null
+++ b/lib/libm/x86_64/Symbol.map
@@ -0,0 +1,12 @@
+/*
+ * $FreeBSD: head/lib/msun/amd64/Symbol.map 226218 2011-10-10 15:43:09Z das $
+ */
+DFLY36.0 {
+	fesetexceptflag;
+	feraiseexcept;
+	fegetenv;
+	feholdexcept;
+	feupdateenv;
+	feenableexcept;
+	fedisableexcept;
+};
diff --git a/share/mk/bsd.lib.mk b/share/mk/bsd.lib.mk
index 336a4704e5..380007148c 100644
--- a/share/mk/bsd.lib.mk
+++ b/share/mk/bsd.lib.mk
@@ -109,6 +109,15 @@ PO_CXXFLAGS=${CXXFLAGS:N-ffunction-sections}
 
 all: objwarn
 
+.include <bsd.symver.mk>
+
+# Allow libraries to specify their own version map or have it
+# automatically generated (see bsd.symver.mk above).
+.if !defined(NO_SYMVER) && !empty(VERSION_MAP)
+${SHLIB_NAME}:	${VERSION_MAP}
+LDFLAGS+=	-Wl,--version-script=${VERSION_MAP}
+.endif
+
 .include <bsd.patch.mk>
 
 .if defined(LIB) && !empty(LIB) || defined(SHLIB_NAME)
diff --git a/share/mk/bsd.symver.mk b/share/mk/bsd.symver.mk
new file mode 100644
index 0000000000..29d347185b
--- /dev/null
+++ b/share/mk/bsd.symver.mk
@@ -0,0 +1,29 @@
+# $FreeBSD$
+
+.if !target(__<bsd.symver.mk>__)
+__<bsd.symver.mk>__:
+
+.include <bsd.init.mk>
+
+# Generate the version map given the version definitions
+# and symbol maps.
+.if !defined(NO_SYMVER) && !empty(VERSION_DEF) && !empty(SYMBOL_MAPS)
+# Find the awk script that generates the version map.
+VERSION_GEN?=	version_gen.awk
+VERSION_MAP?=	Version.map
+
+CLEANFILES+=	${VERSION_MAP}
+
+.if exists(${.PARSEDIR}/${VERSION_GEN})
+_vgen:=	${.PARSEDIR}/${VERSION_GEN}
+.else
+.error	${VERSION_GEN} not found in ${.PARSEDIR}
+.endif
+
+# Run the symbol maps through the C preprocessor before passing
+# them to the symbol version generator.
+${VERSION_MAP}: ${VERSION_DEF} ${_vgen} ${SYMBOL_MAPS}
+	cat ${SYMBOL_MAPS} | ${CPP} - - \
+	    | awk -v vfile=${VERSION_DEF} -f ${_vgen} > ${.TARGET}
+.endif	# !empty(VERSION_DEF) && !empty(SYMBOL_MAPS)
+.endif  # !target(__<bsd.symver.mk>__)
diff --git a/share/mk/version_gen.awk b/share/mk/version_gen.awk
new file mode 100644
index 0000000000..93fbc4fc91
--- /dev/null
+++ b/share/mk/version_gen.awk
@@ -0,0 +1,249 @@
+#
+# Copyright (C) 2006 Daniel M. Eischen.  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 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 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$
+#
+
+#
+# Make a list of all the library versions listed in the master file.
+#
+#   versions[] - array indexed by version name, contains number
+#                of symbols (+ 1) found for each version.
+#   successors[] - array index by version name, contains successor
+#                  version name.
+#   symbols[][] - array index by [version name, symbol index], contains
+#                 names of symbols defined for each version.
+#   names[] - array index is symbol name and value is its first version seen,
+#	      used to check for duplicate symbols and warn about them.
+#
+BEGIN {
+	brackets = 0;
+	errors = warns = 0;
+	version_count = 0;
+	current_version = "";
+	stderr = "/dev/stderr";
+	while (getline < vfile) {
+		# Strip comments.
+		sub("#.*$", "", $0);
+
+		# Strip leading and trailing whitespace.
+		sub("^[ \t]+", "", $0);
+		sub("[ \t]+$", "", $0);
+
+		if (/^[a-zA-Z0-9._]+[ \t]*{$/) {
+			# Strip brace.
+			sub("{", "", $1);
+			brackets++;
+			symver = $1;
+			versions[symver] = 1;
+			successors[symver] = "";
+			generated[symver] = 0;
+			version_count++;
+		}
+		else if (/^}[ \t]*[a-zA-Z0-9._]+[ \t]*;$/) {
+			v = $1 != "}" ? $1 : $2;
+			# Strip brace.
+			sub("}", "", v);
+			# Strip semicolon.
+			sub(";", "", v);
+			if (symver == "") {
+				printf("File %s: Unmatched bracket.\n",
+				vfile) > stderr;
+				errors++;
+			}
+			else if (versions[v] != 1) {
+				printf("File %s: `%s' has unknown " \
+				    "successor `%s'.\n",
+				    vfile, symver, v) > stderr;
+				errors++;
+			}
+			else
+				successors[symver] = v;
+			brackets--;
+		}
+		else if (/^}[ \t]*;$/) {
+			if (symver == "") {
+				printf("File %s: Unmatched bracket.\n",
+				    vfile) > stderr;
+				errors++;
+			}
+			# No successor
+			brackets--;
+		}
+		else if (/^}$/) {
+			printf("File %s: Missing final semicolon.\n",
+			    vfile) > stderr;
+			errors++;
+		}
+		else if (/^$/)
+			;  # Ignore blank lines.
+		else {
+			printf("File %s: Unknown directive: `%s'.\n",
+			    vfile, $0) > stderr;
+			errors++;
+		}
+	}
+	brackets = 0;
+}
+
+{
+	# Set meaningful filename for diagnostics.
+	filename = FILENAME != "" ? FILENAME : "<stdin>";
+
+	# Delete comments, preceding and trailing whitespace, then
+	# consume blank lines.
+	sub("#.*$", "", $0);
+	sub("^[ \t]+", "", $0);
+	sub("[ \t]+$", "", $0);
+	if ($0 == "")
+		next;
+}
+
+/^[a-zA-Z0-9._]+[ \t]*{$/ {
+	# Strip bracket from version name.
+	sub("{", "", $1);
+	if (current_version != "") {
+		printf("File %s, line %d: Illegal nesting detected.\n",
+		    filename, FNR) > stderr;
+		errors++;
+	}
+	else if (versions[$1] == 0) {
+		printf("File %s, line %d: Undefined " \
+		    "library version `%s'.\n", filename, FNR, $1) > stderr;
+		errors++;
+		# Remove this entry from the versions.
+		delete versions[$1];
+	}
+	else
+		current_version = $1;
+	brackets++;
+	next;
+}
+
+/^[a-zA-Z0-9._]+[ \t]*;$/ {
+	# Strip semicolon.
+	sub(";", "", $1);
+	if (current_version != "") {
+		count = versions[current_version];
+		versions[current_version]++;
+		symbols[current_version, count] = $1;
+		if ($1 in names && names[$1] != current_version) {
+			#
+			# A graver case when a dup symbol appears under
+			# different versions in the map.  That can result
+			# in subtle problems with the library later.
+			#
+			printf("File %s, line %d: Duplicated symbol `%s' " \
+			    "in version `%s', first seen in `%s'. " \
+			    "Did you forget to move it to ObsoleteVersions?\n",
+			    filename, FNR, $1,
+			    current_version, names[$1]) > stderr;
+			errors++;
+		}
+		else if (names[$1] == current_version) {
+			#
+			# A harmless case: a dup symbol with the same version.
+			#
+			printf("File %s, line %d: warning: " \
+			    "Duplicated symbol `%s' in version `%s'.\n",
+			    filename, FNR, $1, current_version) > stderr;
+			warns++;
+		}
+		else
+			names[$1] = current_version;
+	}
+	else {
+		printf("File %s, line %d: Symbol `%s' outside version scope.\n",
+		    filename, FNR, $1) > stderr;
+		errors++;
+	}
+	next;
+}
+
+/^}[ \t]*;$/ {
+	brackets--;
+	if (brackets < 0) {
+		printf("File %s, line %d: Unmatched bracket.\n",
+		    filename, FNR, $1) > stderr;
+		errors++;
+		brackets = 0;	# Reset
+	}
+	current_version = "";
+	next;
+}
+
+
+{
+	printf("File %s, line %d: Unknown directive: `%s'.\n",
+	    filename, FNR, $0) > stderr;
+	errors++;
+}
+
+function print_version(v)
+{
+	# This function is recursive, so return if this version
+	# has already been printed.  Otherwise, if there is an
+	# ancestral version, recursively print its symbols before
+	# printing the symbols for this version.
+	#
+	if (generated[v] == 1)
+		return;
+	if (successors[v] != "")
+		print_version(successors[v]);
+
+	printf("%s {\n", v);
+
+	# The version count is always one more that actual,
+	# so the loop ranges from 1 to n-1.
+	#
+	for (i = 1; i < versions[v]; i++) {
+		if (i == 1)
+			printf("global:\n");
+		printf("\t%s;\n", symbols[v, i]);
+	}
+
+	version_count--;
+	if (version_count == 0) {
+		printf("local:\n");
+		printf("\t*;\n");
+	}
+	if (successors[v] == "")
+		printf("};\n");
+	else
+		printf("} %s;\n", successors[v]);
+	printf("\n");
+
+	generated[v] = 1;
+    }
+
+END {
+	if (errors) {
+		printf("%d error(s) total.\n", errors) > stderr;
+		exit(1);
+	}
+	# OK, no errors.
+	for (v in versions) {
+		print_version(v);
+	}
+}
-- 
2.41.0