Add logl, log2l, log10l, and log1pl.

Submitted by:	bde

11 files changed, 1809 insertions, 18 deletions

diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index 642799d..8746420 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -99,8 +99,8 @@ COMMON_SRCS+=	e_acosl.c e_asinl.c e_atan2l.c e_fmodl.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_csqrtl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
-	s_frexpl.c s_logbl.c s_nanl.c s_nextafterl.c s_nexttoward.c \
-	s_remquol.c s_rintl.c s_scalbnl.c \
+	s_frexpl.c s_logbl.c s_logl.c s_nanl.c s_nextafterl.c \
+	s_nexttoward.c s_remquol.c s_rintl.c s_scalbnl.c \
 	s_sinl.c s_tanl.c s_truncl.c w_cabsl.c
 .endif
 
@@ -187,7 +187,9 @@ MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3
 MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3
 MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 lgamma.3 lgammaf.3 \
 	lgamma.3 tgamma.3 lgamma.3 tgammaf.3
-MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log1p.3 log.3 log1pf.3 log.3 logf.3 log.3 log2.3 log.3 log2f.3
+MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log10l.3 log.3 \
+	log1p.3 log.3 log1pf.3 log.3 log1pl.3 log.3 logf.3 log.3 logl.3 \
+	log.3 log2.3 log.3 log2f.3 log.3 log2l.3
 MLINKS+=lrint.3 llrint.3 lrint.3 llrintf.3 lrint.3 llrintl.3 \
 	lrint.3 lrintf.3 lrint.3 lrintl.3
 MLINKS+=lround.3 llround.3 lround.3 llroundf.3 lround.3 llroundl.3 \
diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map
index 38c5941..5e99e76 100644
--- a/lib/msun/Symbol.map
+++ b/lib/msun/Symbol.map
@@ -262,4 +262,8 @@ FBSD_1.3 {
 	ctanh;
 	ctanhf;
 	expl;
+	log10l;
+	log1pl;
+	log2l;
+	logl;
 };
diff --git a/lib/msun/ld128/s_logl.c b/lib/msun/ld128/s_logl.c
new file mode 100644
index 0000000..391d623
--- /dev/null
+++ b/lib/msun/ld128/s_logl.c
@@ -0,0 +1,737 @@
+/*-
+ * 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.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/**
+ * Implementation of the natural logarithm of x for 128-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
+
+#include "fpmath.h"
+#include "math.h"
+#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 ~[-1.1577e-37, 1.1582e-37]:
+ * |log(1 + d)/d - p(d)| < 2**-122.7
+ */
+static const long double
+P2 = -0.5L,
+P3 =  3.33333333333333333333333333333233795e-1L,	/*  0x15555555555555555555555554d42.0p-114L */
+P4 = -2.49999999999999999999999999941139296e-1L,	/* -0x1ffffffffffffffffffffffdab14e.0p-115L */
+P5 =  2.00000000000000000000000085468039943e-1L,	/*  0x19999999999999999999a6d3567f4.0p-115L */
+P6 = -1.66666666666666666666696142372698408e-1L,	/* -0x15555555555555555567267a58e13.0p-115L */
+P7 =  1.42857142857142857119522943477166120e-1L,	/*  0x1249249249249248ed79a0ae434de.0p-115L */
+P8 = -1.24999999999999994863289015033581301e-1L;	/* -0x1fffffffffffffa13e91765e46140.0p-116L */
+/* Double precision gives ~ 53 + log2(P9 * max(|d|)**8) ~= 120 bits. */
+static const double
+P9 =  1.1111111111111401e-1,		/*  0x1c71c71c71c7ed.0p-56 */
+P10 = -1.0000000000040135e-1,		/* -0x199999999a0a92.0p-56 */
+P11 =  9.0909090728136258e-2,		/*  0x1745d173962111.0p-56 */
+P12 = -8.3333318851855284e-2,		/* -0x1555551722c7a3.0p-56 */
+P13 =  7.6928634666404178e-2,		/*  0x13b1985204a4ae.0p-56 */
+P14 = -7.1626810078462499e-2;		/* -0x12562276cdc5d0.0p-56 */
+
+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) */
+	/* The compiler will insert 8 bytes of padding here. */
+	long double F_lo;		/* next 113 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).
+	 *
+	 * The relevant range of dk is [-16445, 16383].  The maximum number
+	 * of bits in F_hi(i) that works is very dependent on i but has
+	 * a minimum of 93.  We only need about 12 bits in F_hi(i) for
+	 * it to provide enough extra precision.
+	 *
+	 * 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, -0x14ee431dae6674afa0c4bfe16e8fd.0p-144L,
+     0xfc0000.0p-24,  0x8102b3.0p-29, -0x1db29ee2d83717be918e1119642ab.0p-144L,
+     0xfa0000.0p-24,  0xc24929.0p-29,  0x1191957d173697cf302cc9476f561.0p-143L,
+     0xf80000.0p-24,  0x820aec.0p-28,  0x13ce8888e02e78eba9b1113bc1c18.0p-142L,
+     0xf60000.0p-24,  0xa33577.0p-28, -0x17a4382ce6eb7bfa509bec8da5f22.0p-142L,
+     0xf48000.0p-24,  0xbc42cb.0p-28, -0x172a21161a107674986dcdca6709c.0p-143L,
+     0xf30000.0p-24,  0xd57797.0p-28, -0x1e09de07cb958897a3ea46e84abb3.0p-142L,
+     0xf10000.0p-24,  0xf7518e.0p-28,  0x1ae1eec1b036c484993c549c4bf40.0p-151L,
+     0xef0000.0p-24,  0x8cb9df.0p-27, -0x1d7355325d560d9e9ab3d6ebab580.0p-141L,
+     0xed8000.0p-24,  0x999ec0.0p-27, -0x1f9f02d256d5037108f4ec21e48cd.0p-142L,
+     0xec0000.0p-24,  0xa6988b.0p-27, -0x16fc0a9d12c17a70f7a684c596b12.0p-143L,
+     0xea0000.0p-24,  0xb80698.0p-27,  0x15d581c1e8da99ded322fb08b8462.0p-141L,
+     0xe80000.0p-24,  0xc99af3.0p-27, -0x1535b3ba8f150ae09996d7bb4653e.0p-143L,
+     0xe70000.0p-24,  0xd273b2.0p-27,  0x163786f5251aefe0ded34c8318f52.0p-145L,
+     0xe50000.0p-24,  0xe442c0.0p-27,  0x1bc4b2368e32d56699c1799a244d4.0p-144L,
+     0xe38000.0p-24,  0xf1b83f.0p-27,  0x1c6090f684e6766abceccab1d7174.0p-141L,
+     0xe20000.0p-24,  0xff448a.0p-27, -0x1890aa69ac9f4215f93936b709efb.0p-142L,
+     0xe08000.0p-24,  0x8673f6.0p-26,  0x1b9985194b6affd511b534b72a28e.0p-140L,
+     0xdf0000.0p-24,  0x8d515c.0p-26, -0x1dc08d61c6ef1d9b2ef7e68680598.0p-143L,
+     0xdd8000.0p-24,  0x943a9e.0p-26, -0x1f72a2dac729b3f46662238a9425a.0p-142L,
+     0xdc0000.0p-24,  0x9b2fe6.0p-26, -0x1fd4dfd3a0afb9691aed4d5e3df94.0p-140L,
+     0xda8000.0p-24,  0xa2315d.0p-26, -0x11b26121629c46c186384993e1c93.0p-142L,
+     0xd90000.0p-24,  0xa93f2f.0p-26,  0x1286d633e8e5697dc6a402a56fce1.0p-141L,
+     0xd78000.0p-24,  0xb05988.0p-26,  0x16128eba9367707ebfa540e45350c.0p-144L,
+     0xd60000.0p-24,  0xb78094.0p-26,  0x16ead577390d31ef0f4c9d43f79b2.0p-140L,
+     0xd50000.0p-24,  0xbc4c6c.0p-26,  0x151131ccf7c7b75e7d900b521c48d.0p-141L,
+     0xd38000.0p-24,  0xc3890a.0p-26, -0x115e2cd714bd06508aeb00d2ae3e9.0p-140L,
+     0xd20000.0p-24,  0xcad2d7.0p-26, -0x1847f406ebd3af80485c2f409633c.0p-142L,
+     0xd10000.0p-24,  0xcfb620.0p-26,  0x1c2259904d686581799fbce0b5f19.0p-141L,
+     0xcf8000.0p-24,  0xd71653.0p-26,  0x1ece57a8d5ae54f550444ecf8b995.0p-140L,
+     0xce0000.0p-24,  0xde843a.0p-26, -0x1f109d4bc4595412b5d2517aaac13.0p-141L,
+     0xcd0000.0p-24,  0xe37fde.0p-26,  0x1bc03dc271a74d3a85b5b43c0e727.0p-141L,
+     0xcb8000.0p-24,  0xeb050c.0p-26, -0x1bf2badc0df841a71b79dd5645b46.0p-145L,
+     0xca0000.0p-24,  0xf29878.0p-26, -0x18efededd89fbe0bcfbe6d6db9f66.0p-147L,
+     0xc90000.0p-24,  0xf7ad6f.0p-26,  0x1373ff977baa6911c7bafcb4d84fb.0p-141L,
+     0xc80000.0p-24,  0xfcc8e3.0p-26,  0x196766f2fb328337cc050c6d83b22.0p-140L,
+     0xc68000.0p-24,  0x823f30.0p-25,  0x19bd076f7c434e5fcf1a212e2a91e.0p-139L,
+     0xc58000.0p-24,  0x84d52c.0p-25, -0x1a327257af0f465e5ecab5f2a6f81.0p-139L,
+     0xc40000.0p-24,  0x88bc74.0p-25,  0x113f23def19c5a0fe396f40f1dda9.0p-141L,
+     0xc30000.0p-24,  0x8b5ae6.0p-25,  0x1759f6e6b37de945a049a962e66c6.0p-139L,
+     0xc20000.0p-24,  0x8dfccb.0p-25,  0x1ad35ca6ed5147bdb6ddcaf59c425.0p-141L,
+     0xc10000.0p-24,  0x90a22b.0p-25,  0x1a1d71a87deba46bae9827221dc98.0p-139L,
+     0xbf8000.0p-24,  0x94a0d8.0p-25, -0x139e5210c2b730e28aba001a9b5e0.0p-140L,
+     0xbe8000.0p-24,  0x974f16.0p-25, -0x18f6ebcff3ed72e23e13431adc4a5.0p-141L,
+     0xbd8000.0p-24,  0x9a00f1.0p-25, -0x1aa268be39aab7148e8d80caa10b7.0p-139L,
+     0xbc8000.0p-24,  0x9cb672.0p-25, -0x14c8815839c5663663d15faed7771.0p-139L,
+     0xbb0000.0p-24,  0xa0cda1.0p-25,  0x1eaf46390dbb2438273918db7df5c.0p-141L,
+     0xba0000.0p-24,  0xa38c6e.0p-25,  0x138e20d831f698298adddd7f32686.0p-141L,
+     0xb90000.0p-24,  0xa64f05.0p-25, -0x1e8d3c41123615b147a5d47bc208f.0p-142L,
+     0xb80000.0p-24,  0xa91570.0p-25,  0x1ce28f5f3840b263acb4351104631.0p-140L,
+     0xb70000.0p-24,  0xabdfbb.0p-25, -0x186e5c0a42423457e22d8c650b355.0p-139L,
+     0xb60000.0p-24,  0xaeadef.0p-25, -0x14d41a0b2a08a465dc513b13f567d.0p-143L,
+     0xb50000.0p-24,  0xb18018.0p-25,  0x16755892770633947ffe651e7352f.0p-139L,
+     0xb40000.0p-24,  0xb45642.0p-25, -0x16395ebe59b15228bfe8798d10ff0.0p-142L,
+     0xb30000.0p-24,  0xb73077.0p-25,  0x1abc65c8595f088b61a335f5b688c.0p-140L,
+     0xb20000.0p-24,  0xba0ec4.0p-25, -0x1273089d3dad88e7d353e9967d548.0p-139L,
+     0xb10000.0p-24,  0xbcf133.0p-25,  0x10f9f67b1f4bbf45de06ecebfaf6d.0p-139L,
+     0xb00000.0p-24,  0xbfd7d2.0p-25, -0x109fab904864092b34edda19a831e.0p-140L,
+     0xaf0000.0p-24,  0xc2c2ac.0p-25, -0x1124680aa43333221d8a9b475a6ba.0p-139L,
+     0xae8000.0p-24,  0xc439b3.0p-25, -0x1f360cc4710fbfe24b633f4e8d84d.0p-140L,
+     0xad8000.0p-24,  0xc72afd.0p-25, -0x132d91f21d89c89c45003fc5d7807.0p-140L,
+     0xac8000.0p-24,  0xca20a2.0p-25, -0x16bf9b4d1f8da8002f2449e174504.0p-139L,
+     0xab8000.0p-24,  0xcd1aae.0p-25,  0x19deb5ce6a6a8717d5626e16acc7d.0p-141L,
+     0xaa8000.0p-24,  0xd0192f.0p-25,  0x1a29fb48f7d3ca87dabf351aa41f4.0p-139L,
+     0xaa0000.0p-24,  0xd19a20.0p-25,  0x1127d3c6457f9d79f51dcc73014c9.0p-141L,
+     0xa90000.0p-24,  0xd49f6a.0p-25, -0x1ba930e486a0ac42d1bf9199188e7.0p-141L,
+     0xa80000.0p-24,  0xd7a94b.0p-25, -0x1b6e645f31549dd1160bcc45c7e2c.0p-139L,
+     0xa70000.0p-24,  0xdab7d0.0p-25,  0x1118a425494b610665377f15625b6.0p-140L,
+     0xa68000.0p-24,  0xdc40d5.0p-25,  0x1966f24d29d3a2d1b2176010478be.0p-140L,
+     0xa58000.0p-24,  0xdf566d.0p-25, -0x1d8e52eb2248f0c95dd83626d7333.0p-142L,
+     0xa48000.0p-24,  0xe270ce.0p-25, -0x1ee370f96e6b67ccb006a5b9890ea.0p-140L,
+     0xa40000.0p-24,  0xe3ffce.0p-25,  0x1d155324911f56db28da4d629d00a.0p-140L,
+     0xa30000.0p-24,  0xe72179.0p-25, -0x1fe6e2f2f867d8f4d60c713346641.0p-140L,
+     0xa20000.0p-24,  0xea4812.0p-25,  0x1b7be9add7f4d3b3d406b6cbf3ce5.0p-140L,
+     0xa18000.0p-24,  0xebdd3d.0p-25,  0x1b3cfb3f7511dd73692609040ccc2.0p-139L,
+     0xa08000.0p-24,  0xef0b5b.0p-25, -0x1220de1f7301901b8ad85c25afd09.0p-139L,
+     0xa00000.0p-24,  0xf0a451.0p-25, -0x176364c9ac81cc8a4dfb804de6867.0p-140L,
+     0x9f0000.0p-24,  0xf3da16.0p-25,  0x1eed6b9aafac8d42f78d3e65d3727.0p-141L,
+     0x9e8000.0p-24,  0xf576e9.0p-25,  0x1d593218675af269647b783d88999.0p-139L,
+     0x9d8000.0p-24,  0xf8b47c.0p-25, -0x13e8eb7da053e063714615f7cc91d.0p-144L,
+     0x9d0000.0p-24,  0xfa553f.0p-25,  0x1c063259bcade02951686d5373aec.0p-139L,
+     0x9c0000.0p-24,  0xfd9ac5.0p-25,  0x1ef491085fa3c1649349630531502.0p-139L,
+     0x9b8000.0p-24,  0xff3f8c.0p-25,  0x1d607a7c2b8c5320619fb9433d841.0p-139L,
+     0x9a8000.0p-24,  0x814697.0p-24, -0x12ad3817004f3f0bdff99f932b273.0p-138L,
+     0x9a0000.0p-24,  0x821b06.0p-24, -0x189fc53117f9e54e78103a2bc1767.0p-141L,
+     0x990000.0p-24,  0x83c5f8.0p-24,  0x14cf15a048907b7d7f47ddb45c5a3.0p-139L,
+     0x988000.0p-24,  0x849c7d.0p-24,  0x1cbb1d35fb82873b04a9af1dd692c.0p-138L,
+     0x978000.0p-24,  0x864ba6.0p-24,  0x1128639b814f9b9770d8cb6573540.0p-138L,
+     0x970000.0p-24,  0x87244c.0p-24,  0x184733853300f002e836dfd47bd41.0p-139L,
+     0x968000.0p-24,  0x87fdaa.0p-24,  0x109d23aef77dd5cd7cc94306fb3ff.0p-140L,
+     0x958000.0p-24,  0x89b293.0p-24, -0x1a81ef367a59de2b41eeebd550702.0p-138L,
+     0x950000.0p-24,  0x8a8e20.0p-24, -0x121ad3dbb2f45275c917a30df4ac9.0p-138L,
+     0x948000.0p-24,  0x8b6a6a.0p-24, -0x1cfb981628af71a89df4e6df2e93b.0p-139L,
+     0x938000.0p-24,  0x8d253a.0p-24, -0x1d21730ea76cfdec367828734cae5.0p-139L,
+     0x930000.0p-24,  0x8e03c2.0p-24,  0x135cc00e566f76b87333891e0dec4.0p-138L,
+     0x928000.0p-24,  0x8ee30d.0p-24, -0x10fcb5df257a263e3bf446c6e3f69.0p-140L,
+     0x918000.0p-24,  0x90a3ee.0p-24, -0x16e171b15433d723a4c7380a448d8.0p-139L,
+     0x910000.0p-24,  0x918587.0p-24, -0x1d050da07f3236f330972da2a7a87.0p-139L,
+     0x908000.0p-24,  0x9267e7.0p-24,  0x1be03669a5268d21148c6002becd3.0p-139L,
+     0x8f8000.0p-24,  0x942f04.0p-24,  0x10b28e0e26c336af90e00533323ba.0p-139L,
+     0x8f0000.0p-24,  0x9513c3.0p-24,  0x1a1d820da57cf2f105a89060046aa.0p-138L,
+     0x8e8000.0p-24,  0x95f950.0p-24, -0x19ef8f13ae3cf162409d8ea99d4c0.0p-139L,
+     0x8e0000.0p-24,  0x96dfab.0p-24, -0x109e417a6e507b9dc10dac743ad7a.0p-138L,
+     0x8d0000.0p-24,  0x98aed2.0p-24,  0x10d01a2c5b0e97c4990b23d9ac1f5.0p-139L,
+     0x8c8000.0p-24,  0x9997a2.0p-24, -0x1d6a50d4b61ea74540bdd2aa99a42.0p-138L,
+     0x8c0000.0p-24,  0x9a8145.0p-24,  0x1b3b190b83f9527e6aba8f2d783c1.0p-138L,
+     0x8b8000.0p-24,  0x9b6bbf.0p-24,  0x13a69fad7e7abe7ba81c664c107e0.0p-138L,
+     0x8b0000.0p-24,  0x9c5711.0p-24, -0x11cd12316f576aad348ae79867223.0p-138L,
+     0x8a8000.0p-24,  0x9d433b.0p-24,  0x1c95c444b807a246726b304ccae56.0p-139L,
+     0x898000.0p-24,  0x9f1e22.0p-24, -0x1b9c224ea698c2f9b47466d6123fe.0p-139L,
+     0x890000.0p-24,  0xa00ce1.0p-24,  0x125ca93186cf0f38b4619a2483399.0p-141L,
+     0x888000.0p-24,  0xa0fc80.0p-24, -0x1ee38a7bc228b3597043be78eaf49.0p-139L,
+     0x880000.0p-24,  0xa1ed00.0p-24, -0x1a0db876613d204147dc69a07a649.0p-138L,
+     0x878000.0p-24,  0xa2de62.0p-24,  0x193224e8516c008d3602a7b41c6e8.0p-139L,
+     0x870000.0p-24,  0xa3d0a9.0p-24,  0x1fa28b4d2541aca7d5844606b2421.0p-139L,
+     0x868000.0p-24,  0xa4c3d6.0p-24,  0x1c1b5760fb4571acbcfb03f16daf4.0p-138L,
+     0x858000.0p-24,  0xa6acea.0p-24,  0x1fed5d0f65949c0a345ad743ae1ae.0p-140L,
+     0x850000.0p-24,  0xa7a2d4.0p-24,  0x1ad270c9d749362382a7688479e24.0p-140L,
+     0x848000.0p-24,  0xa899ab.0p-24,  0x199ff15ce532661ea9643a3a2d378.0p-139L,
+     0x840000.0p-24,  0xa99171.0p-24,  0x1a19e15ccc45d257530a682b80490.0p-139L,
+     0x838000.0p-24,  0xaa8a28.0p-24, -0x121a14ec532b35ba3e1f868fd0b5e.0p-140L,
+     0x830000.0p-24,  0xab83d1.0p-24,  0x1aee319980bff3303dd481779df69.0p-139L,
+     0x828000.0p-24,  0xac7e6f.0p-24, -0x18ffd9e3900345a85d2d86161742e.0p-140L,
+     0x820000.0p-24,  0xad7a03.0p-24, -0x1e4db102ce29f79b026b64b42caa1.0p-140L,
+     0x818000.0p-24,  0xae768f.0p-24,  0x17c35c55a04a82ab19f77652d977a.0p-141L,
+     0x810000.0p-24,  0xaf7415.0p-24,  0x1448324047019b48d7b98c1cf7234.0p-138L,
+     0x808000.0p-24,  0xb07298.0p-24, -0x1750ee3915a197e9c7359dd94152f.0p-138L,
+     0x800000.0p-24,  0xb17218.0p-24, -0x105c610ca86c3898cff81a12a17e2.0p-141L,
+};
+
+#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, val_hi, val_lo;
+	double dd, dk;
+	uint64_t lx, llx;
+	int i, k;
+	uint16_t hx;
+
+	EXTRACT_LDBL128_WORDS(hx, lx, llx, 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 | llx) == 0)
+			RETURN1(rp, -1 / zero);	/* log(+-0) = -Inf */
+		if (hx != 0)
+			/* log(neg or NaN) = qNaN: */
+			RETURN1(rp, (x - x) / zero);
+		x *= 0x1.0p113;		/* subnormal; scale up x */
+		EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+		k = -16383 - 113;
+	} else if (hx >= 0x7fff)
+		RETURN1(rp, x + x);	/* log(Inf or NaN) = Inf or qNaN */
+#ifndef STRUCT_RETURN
+	ENTERI();
+#endif
+	k += hx;
+	dk = k;
+
+	/* Scale x to be in [1, 2). */
+	SET_LDBL_EXPSIGN(x, 0x3fff);
+
+	/* 0 <= i <= INTERVALS: */
+#define	L2I	(49 - LOG2_INTERVALS)
+	i = (lx + (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;
+		double x_lo;
+
+		/*
+		 * 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.
+		 */
+		INSERT_LDBL128_WORDS(x_hi, 0x3fff, lx,
+		    llx & 0xffffffffff000000ULL);
+		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.
+	 */
+	/*
+	 * Use double precision operations wherever possible, since long
+	 * double operations are emulated and are very slow on the only
+	 * known machines that support ld128 (sparc64).  Also, don't try
+	 * to improve parallelism by increasing the number of operations,
+	 * since any parallelism on such machines is needed for the
+	 * emulation.  Horner's method is good for this, and is also good
+	 * for accuracy.  Horner's method doesn't handle the `lo' term
+	 * well, either for efficiency or accuracy.  However, for accuracy
+	 * we evaluate d * d * P2 separately to take advantage of
+	 * by P2 being exact, and this gives a good place to sum the 'lo'
+	 * term too.
+	 */
+	dd = (double)d;
+	val_lo = d * d * d * (P3 +
+	    d * (P4 + d * (P5 + d * (P6 + d * (P7 + d * (P8 +
+	    dd * (P9 + dd * (P10 + dd * (P11 + dd * (P12 + dd * (P13 +
+	    dd * P14))))))))))) + (F_lo(i) + dk * ln2_lo) + d * d * 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, f_lo, val_hi, val_lo;
+	long double f_hi, twopminusk;
+	double d_lo, dd, dk;
+	uint64_t lx, llx;
+	int i, k;
+	int16_t ax, hx;
+
+	DOPRINT_START(&x);
+	EXTRACT_LDBL128_WORDS(hx, lx, llx, 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 | llx) == 0)
+				RETURNP(-1 / zero);	/* log1p(-1) = -Inf */
+			/* log1p(x < 1, or x NaN) = qNaN: */
+			RETURNP((x - x) / (x - x));
+		}
+		if (ax <= 0x3f8d) {	/* |x| < 2**-113 */
+			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 */
+	} else if (hx < 0x40e1) {	/* 1 <= x < 2**226 */
+		f_hi = x;
+		f_lo = 1;
+	} else {			/* 2**226 <= x < +Inf */
+		f_hi = x;
+		f_lo = 0;		/* avoid underflow of the P3 term */
+	}
+	ENTERI();
+	x = f_hi + f_lo;
+	f_lo = (f_hi - x) + f_lo;
+
+	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+	k = -16383;
+
+	k += hx;
+	dk = k;
+
+	SET_LDBL_EXPSIGN(x, 0x3fff);
+	twopminusk = 1;
+	SET_LDBL_EXPSIGN(twopminusk, 0x7ffe - (hx & 0x7fff));
+	f_lo *= twopminusk;
+
+	i = (lx + (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;
+		double x_lo;
+
+		INSERT_LDBL128_WORDS(x_hi, 0x3fff, lx,
+		    llx & 0xffffffffff000000ULL);
+		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;
+
+	dd = (double)d;
+	val_lo = d * d * d * (P3 +
+	    d * (P4 + d * (P5 + d * (P6 + d * (P7 + d * (P8 +
+	    dd * (P9 + dd * (P10 + dd * (P11 + dd * (P12 + dd * (P13 +
+	    dd * P14))))))))))) + (F_lo(i) + dk * ln2_lo + d_lo) + d * d * 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);
+}
+
+/*
+ * 29+113 bit decompositions.  The bits are distributed so that the products
+ * of the hi terms are exact in double precision.  The types are chosen so
+ * that the products of the hi terms are done in at least double precision,
+ * without any explicit conversions.  More natural choices would require a
+ * slow long double precision multiplication.
+ */
+static const double
+invln10_hi =  4.3429448176175356e-1,		/*  0x1bcb7b15000000.0p-54 */
+invln2_hi =  1.4426950402557850e0;		/*  0x17154765000000.0p-52 */
+static const long double
+invln10_lo =  1.41498268538580090791605082294397000e-10L,	/*  0x137287195355baaafad33dc323ee3.0p-145L */
+invln2_lo =  6.33178418956604368501892137426645911e-10L;	/*  0x15c17f0bbbe87fed0691d3e88eb57.0p-143L */
+
+long double
+log10l(long double x)
+{
+	struct ld r;
+	long double lo;
+	float hi;
+
+	ENTERI();
+	DOPRINT_START(&x);
+	k_logl(x, &r);
+	if (!r.lo_set)
+		RETURNPI(r.hi);
+	_2sumF(r.hi, r.lo);
+	hi = 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 lo;
+	float hi;
+
+	ENTERI();
+	DOPRINT_START(&x);
+	k_logl(x, &r);
+	if (!r.lo_set)
+		RETURNPI(r.hi);
+	_2sumF(r.hi, r.lo);
+	hi = 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/msun/ld80/s_logl.c b/lib/msun/ld80/s_logl.c
new file mode 100644
index 0000000..3a35753
--- /dev/null
+++ b/lib/msun/ld80/s_logl.c
@@ -0,0 +1,717 @@
+/*-
+ * 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.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/**
+ * 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/msun/man/log.3 b/lib/msun/man/log.3
index b9fd83c..b08e692 100644
--- a/lib/msun/man/log.3
+++ b/lib/msun/man/log.3
@@ -24,7 +24,7 @@
 .\"
 .\" $FreeBSD$
 .\"
-.Dd December 5, 2010
+.Dd June 3, 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/msun/src/e_log.c b/lib/msun/src/e_log.c
index 19f1ab9..68bc107 100644
--- a/lib/msun/src/e_log.c
+++ b/lib/msun/src/e_log.c
@@ -65,6 +65,8 @@ __FBSDID("$FreeBSD$");
  * to produce the hexadecimal values shown.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -139,3 +141,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/msun/src/e_log10.c b/lib/msun/src/e_log10.c
index a795129..3c89ed2 100644
--- a/lib/msun/src/e_log10.c
+++ b/lib/msun/src/e_log10.c
@@ -22,6 +22,8 @@ __FBSDID("$FreeBSD$");
  * in not-quite-routine extra precision.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 #include "k_log.h"
@@ -86,3 +88,7 @@ __ieee754_log10(double x)
 
 	return val_lo + val_hi;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(log10, log10l);
+#endif
diff --git a/lib/msun/src/e_log2.c b/lib/msun/src/e_log2.c
index 010bc2c..d737f04 100644
--- a/lib/msun/src/e_log2.c
+++ b/lib/msun/src/e_log2.c
@@ -109,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/msun/src/math.h b/lib/msun/src/math.h
index c6cee13..d578e41 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -418,7 +418,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 */
@@ -464,10 +468,6 @@ 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/msun/src/math_private.h b/lib/msun/src/math_private.h
index 5662df0..8ebc7fb 100644
--- a/lib/msun/src/math_private.h
+++ b/lib/msun/src/math_private.h
@@ -188,6 +188,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)					\
@@ -197,6 +224,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)					\
@@ -261,6 +315,110 @@ do {								\
 #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().
  */
 void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
@@ -370,6 +528,140 @@ irintl(long double x)
 
 #endif /* __GNUCLIKE_ASM */
 
+#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/msun/src/s_log1p.c b/lib/msun/src/s_log1p.c
index 55d352c..3cc77bd 100644
--- a/lib/msun/src/s_log1p.c
+++ b/lib/msun/src/s_log1p.c
@@ -174,3 +174,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