MFC r257770 r257818 r257823 r260066 r260067 r260089 r260145 r268587 r268588

    r268589 r268590 r268593 r268597 r269758 r270845 r270847 r270893 r270932
    r270947 r271147

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

Approved by:	re (kib)

29 files changed, 2140 insertions, 723 deletions

diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index 4bca55b..654d0f5 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -72,7 +72,7 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
 	s_lround.c s_lroundf.c s_lroundl.c s_modff.c \
 	s_nan.c s_nearbyint.c s_nextafter.c s_nextafterf.c \
 	s_nexttowardf.c s_remquo.c s_remquof.c \
-	s_rint.c s_rintf.c s_round.c s_roundf.c s_roundl.c \
+	s_rint.c s_rintf.c s_round.c s_roundf.c \
 	s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
 	s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \
 	s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \
@@ -97,13 +97,14 @@ COMMON_SRCS+=	s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c
 .if ${LDBL_PREC} != 53
 # If long double != double use these; otherwise, we alias the double versions.
 COMMON_SRCS+=	e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \
-	e_fmodl.c e_hypotl.c e_remainderl.c e_sqrtl.c \
+	e_coshl.c e_fmodl.c e_hypotl.c \
+	e_remainderl.c e_sinhl.c e_sqrtl.c \
 	invtrig.c k_cosl.c k_sinl.c k_tanl.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_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.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
+	s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c s_scalbnl.c \
+	s_sinl.c s_tanhl.c s_tanl.c s_truncl.c w_cabsl.c
 .endif
 
 # C99 complex functions
@@ -161,9 +162,9 @@ MLINKS+=cimag.3 cimagf.3 cimag.3 cimagl.3 \
 	cimag.3 creal.3 cimag.3 crealf.3 cimag.3 creall.3
 MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3
 MLINKS+=cos.3 cosf.3 cos.3 cosl.3
-MLINKS+=cosh.3 coshf.3
+MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3
 MLINKS+=csqrt.3 csqrtf.3 csqrt.3 csqrtl.3
-MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3
+MLINKS+=erf.3 erfc.3 erf.3 erff.3 erf.3 erfcf.3 erf.3 erfl.3 erf.3 erfcl.3
 MLINKS+=exp.3 expm1.3 exp.3 expm1f.3 exp.3 expm1l.3 exp.3 pow.3 exp.3 powf.3 \
 	exp.3 exp2.3 exp.3 exp2f.3 exp.3 exp2l.3 exp.3 expf.3 exp.3 expl.3
 MLINKS+=fabs.3 fabsf.3 fabs.3 fabsl.3
@@ -209,11 +210,11 @@ MLINKS+=round.3 roundf.3 round.3 roundl.3
 MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3
 MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
 MLINKS+=sin.3 sinf.3 sin.3 sinl.3
-MLINKS+=sinh.3 sinhf.3
+MLINKS+=sinh.3 sinhf.3 sinh.3 sinhl.3
 MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
 	sqrt.3 sqrtl.3
 MLINKS+=tan.3 tanf.3 tan.3 tanl.3
-MLINKS+=tanh.3 tanhf.3
+MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3
 MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3
 
 .include <bsd.lib.mk>
diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map
index 037659d..e53ca07 100644
--- a/lib/msun/Symbol.map
+++ b/lib/msun/Symbol.map
@@ -260,23 +260,23 @@ FBSD_1.3 {
 	ccosf;
 	ccosh;
 	ccoshf;
+	coshl;
 	ctan;
 	ctanf;
 	ctanh;
 	ctanhf;
+	erfcl;
+	erfl;
 	expl;
 	expm1l;
 	log10l;
 	log1pl;
 	log2l;
 	logl;
+	sinhl;
+	tanhl;
 	/* Implemented as weak aliases for imprecise versions */
-	coshl;
-	erfcl;
-	erfl;
 	lgammal;
 	powl;
-	sinhl;
-	tanhl;
 	tgammal;
 };
diff --git a/lib/msun/ld128/k_expl.h b/lib/msun/ld128/k_expl.h
new file mode 100644
index 0000000..a5668fd
--- /dev/null
+++ b/lib/msun/ld128/k_expl.h
@@ -0,0 +1,328 @@
+/* from: FreeBSD: head/lib/msun/ld128/s_expl.c 251345 2013-06-03 20:09:22Z kargl */
+
+/*-
+ * Copyright (c) 2009-2013 Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice unmodified, this list of conditions, and the following
+ *    disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * ld128 version of k_expl.h.  See ../ld80/s_expl.c for most comments.
+ *
+ * See ../src/e_exp.c and ../src/k_exp.h for precision-independent comments
+ * about the secondary kernels.
+ */
+
+#define	INTERVALS	128
+#define	LOG2_INTERVALS	7
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+static const double
+/*
+ * ln2/INTERVALS = L1+L2 (hi+lo decomposition for multiplication).  L1 must
+ * have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest
+ * bits zero so that multiplication of it by n is exact.
+ */
+INV_L = 1.8466496523378731e+2,		/*  0x171547652b82fe.0p-45 */
+L2 = -1.0253670638894731e-29;		/* -0x1.9ff0342542fc3p-97 */
+static const long double
+/* 0x1.62e42fefa39ef35793c768000000p-8 */
+L1 =  5.41521234812457272982212595914567508e-3L;
+
+/*
+ * XXX values in hex in comments have been lost (or were never present)
+ * from here.
+ */
+static const long double
+/*
+ * Domain [-0.002708, 0.002708], range ~[-2.4021e-38, 2.4234e-38]:
+ * |exp(x) - p(x)| < 2**-124.9
+ * (0.002708 is ln2/(2*INTERVALS) rounded up a little).
+ *
+ * XXX the coeffs aren't very carefully rounded, and I get 3.6 more bits.
+ */
+A2  =  0.5,
+A3  =  1.66666666666666666666666666651085500e-1L,
+A4  =  4.16666666666666666666666666425885320e-2L,
+A5  =  8.33333333333333333334522877160175842e-3L,
+A6  =  1.38888888888888888889971139751596836e-3L;
+
+static const double
+A7  =  1.9841269841269470e-4,		/*  0x1.a01a01a019f91p-13 */
+A8  =  2.4801587301585286e-5,		/*  0x1.71de3ec75a967p-19 */
+A9  =  2.7557324277411235e-6,		/*  0x1.71de3ec75a967p-19 */
+A10 =  2.7557333722375069e-7;		/*  0x1.27e505ab56259p-22 */
+
+static const struct {
+	/*
+	 * hi must be rounded to at most 106 bits so that multiplication
+	 * by r1 in expm1l() is exact, but it is rounded to 88 bits due to
+	 * historical accidents.
+	 *
+	 * XXX it is wasteful to use long double for both hi and lo.  ld128
+	 * exp2l() uses only float for lo (in a very differently organized
+	 * table; ld80 exp2l() is different again.  It uses 2 doubles in a
+	 * table organized like this one.  1 double and 1 float would
+	 * suffice).  There are different packing/locality/alignment/caching
+	 * problems with these methods.
+	 *
+	 * XXX C's bad %a format makes the bits unreadable.  They happen
+	 * to all line up for the hi values 1 before the point and 88
+	 * in 22 nybbles, but for the low values the nybbles are shifted
+	 * randomly.
+	 */
+	long double	hi;
+	long double	lo;
+} tbl[INTERVALS] = {
+	0x1p0L, 0x0p0L,
+	0x1.0163da9fb33356d84a66aep0L, 0x3.36dcdfa4003ec04c360be2404078p-92L,
+	0x1.02c9a3e778060ee6f7cacap0L, 0x4.f7a29bde93d70a2cabc5cb89ba10p-92L,
+	0x1.04315e86e7f84bd738f9a2p0L, 0xd.a47e6ed040bb4bfc05af6455e9b8p-96L,
+	0x1.059b0d31585743ae7c548ep0L, 0xb.68ca417fe53e3495f7df4baf84a0p-92L,
+	0x1.0706b29ddf6ddc6dc403a8p0L, 0x1.d87b27ed07cb8b092ac75e311753p-88L,
+	0x1.0874518759bc808c35f25cp0L, 0x1.9427fa2b041b2d6829d8993a0d01p-88L,
+	0x1.09e3ecac6f3834521e060cp0L, 0x5.84d6b74ba2e023da730e7fccb758p-92L,
+	0x1.0b5586cf9890f6298b92b6p0L, 0x1.1842a98364291408b3ceb0a2a2bbp-88L,
+	0x1.0cc922b7247f7407b705b8p0L, 0x9.3dc5e8aac564e6fe2ef1d431fd98p-92L,
+	0x1.0e3ec32d3d1a2020742e4ep0L, 0x1.8af6a552ac4b358b1129e9f966a4p-88L,
+	0x1.0fb66affed31af232091dcp0L, 0x1.8a1426514e0b627bda694a400a27p-88L,
+	0x1.11301d0125b50a4ebbf1aep0L, 0xd.9318ceac5cc47ab166ee57427178p-92L,
+	0x1.12abdc06c31cbfb92bad32p0L, 0x4.d68e2f7270bdf7cedf94eb1cb818p-92L,
+	0x1.1429aaea92ddfb34101942p0L, 0x1.b2586d01844b389bea7aedd221d4p-88L,
+	0x1.15a98c8a58e512480d573cp0L, 0x1.d5613bf92a2b618ee31b376c2689p-88L,
+	0x1.172b83c7d517adcdf7c8c4p0L, 0x1.0eb14a792035509ff7d758693f24p-88L,
+	0x1.18af9388c8de9bbbf70b9ap0L, 0x3.c2505c97c0102e5f1211941d2840p-92L,
+	0x1.1a35beb6fcb753cb698f68p0L, 0x1.2d1c835a6c30724d5cfae31b84e5p-88L,
+	0x1.1bbe084045cd39ab1e72b4p0L, 0x4.27e35f9acb57e473915519a1b448p-92L,
+	0x1.1d4873168b9aa7805b8028p0L, 0x9.90f07a98b42206e46166cf051d70p-92L,
+	0x1.1ed5022fcd91cb8819ff60p0L, 0x1.121d1e504d36c47474c9b7de6067p-88L,
+	0x1.2063b88628cd63b8eeb028p0L, 0x1.50929d0fc487d21c2b84004264dep-88L,
+	0x1.21f49917ddc962552fd292p0L, 0x9.4bdb4b61ea62477caa1dce823ba0p-92L,
+	0x1.2387a6e75623866c1fadb0p0L, 0x1.c15cb593b0328566902df69e4de2p-88L,
+	0x1.251ce4fb2a63f3582ab7dep0L, 0x9.e94811a9c8afdcf796934bc652d0p-92L,
+	0x1.26b4565e27cdd257a67328p0L, 0x1.d3b249dce4e9186ddd5ff44e6b08p-92L,
+	0x1.284dfe1f5638096cf15cf0p0L, 0x3.ca0967fdaa2e52d7c8106f2e262cp-92L,
+	0x1.29e9df51fdee12c25d15f4p0L, 0x1.a24aa3bca890ac08d203fed80a07p-88L,
+	0x1.2b87fd0dad98ffddea4652p0L, 0x1.8fcab88442fdc3cb6de4519165edp-88L,
+	0x1.2d285a6e4030b40091d536p0L, 0xd.075384589c1cd1b3e4018a6b1348p-92L,
+	0x1.2ecafa93e2f5611ca0f45cp0L, 0x1.523833af611bdcda253c554cf278p-88L,
+	0x1.306fe0a31b7152de8d5a46p0L, 0x3.05c85edecbc27343629f502f1af2p-92L,
+	0x1.32170fc4cd8313539cf1c2p0L, 0x1.008f86dde3220ae17a005b6412bep-88L,
+	0x1.33c08b26416ff4c9c8610cp0L, 0x1.96696bf95d1593039539d94d662bp-88L,
+	0x1.356c55f929ff0c94623476p0L, 0x3.73af38d6d8d6f9506c9bbc93cbc0p-92L,
+	0x1.371a7373aa9caa7145502ep0L, 0x1.4547987e3e12516bf9c699be432fp-88L,
+	0x1.38cae6d05d86585a9cb0d8p0L, 0x1.bed0c853bd30a02790931eb2e8f0p-88L,
+	0x1.3a7db34e59ff6ea1bc9298p0L, 0x1.e0a1d336163fe2f852ceeb134067p-88L,
+	0x1.3c32dc313a8e484001f228p0L, 0xb.58f3775e06ab66353001fae9fca0p-92L,
+	0x1.3dea64c12342235b41223ep0L, 0x1.3d773fba2cb82b8244267c54443fp-92L,
+	0x1.3fa4504ac801ba0bf701aap0L, 0x4.1832fb8c1c8dbdff2c49909e6c60p-92L,
+	0x1.4160a21f72e29f84325b8ep0L, 0x1.3db61fb352f0540e6ba05634413ep-88L,
+	0x1.431f5d950a896dc7044394p0L, 0x1.0ccec81e24b0caff7581ef4127f7p-92L,
+	0x1.44e086061892d03136f408p0L, 0x1.df019fbd4f3b48709b78591d5cb5p-88L,
+	0x1.46a41ed1d005772512f458p0L, 0x1.229d97df404ff21f39c1b594d3a8p-88L,
+	0x1.486a2b5c13cd013c1a3b68p0L, 0x1.062f03c3dd75ce8757f780e6ec99p-88L,
+	0x1.4a32af0d7d3de672d8bcf4p0L, 0x6.f9586461db1d878b1d148bd3ccb8p-92L,
+	0x1.4bfdad5362a271d4397afep0L, 0xc.42e20e0363ba2e159c579f82e4b0p-92L,
+	0x1.4dcb299fddd0d63b36ef1ap0L, 0x9.e0cc484b25a5566d0bd5f58ad238p-92L,
+	0x1.4f9b2769d2ca6ad33d8b68p0L, 0x1.aa073ee55e028497a329a7333dbap-88L,
+	0x1.516daa2cf6641c112f52c8p0L, 0x4.d822190e718226177d7608d20038p-92L,
+	0x1.5342b569d4f81df0a83c48p0L, 0x1.d86a63f4e672a3e429805b049465p-88L,
+	0x1.551a4ca5d920ec52ec6202p0L, 0x4.34ca672645dc6c124d6619a87574p-92L,
+	0x1.56f4736b527da66ecb0046p0L, 0x1.64eb3c00f2f5ab3d801d7cc7272dp-88L,
+	0x1.58d12d497c7fd252bc2b72p0L, 0x1.43bcf2ec936a970d9cc266f0072fp-88L,
+	0x1.5ab07dd48542958c930150p0L, 0x1.91eb345d88d7c81280e069fbdb63p-88L,
+	0x1.5c9268a5946b701c4b1b80p0L, 0x1.6986a203d84e6a4a92f179e71889p-88L,
+	0x1.5e76f15ad21486e9be4c20p0L, 0x3.99766a06548a05829e853bdb2b52p-92L,
+	0x1.605e1b976dc08b076f592ap0L, 0x4.86e3b34ead1b4769df867b9c89ccp-92L,
+	0x1.6247eb03a5584b1f0fa06ep0L, 0x1.d2da42bb1ceaf9f732275b8aef30p-88L,
+	0x1.6434634ccc31fc76f8714cp0L, 0x4.ed9a4e41000307103a18cf7a6e08p-92L,
+	0x1.66238825522249127d9e28p0L, 0x1.b8f314a337f4dc0a3adf1787ff74p-88L,
+	0x1.68155d44ca973081c57226p0L, 0x1.b9f32706bfe4e627d809a85dcc66p-88L,
+	0x1.6a09e667f3bcc908b2fb12p0L, 0x1.66ea957d3e3adec17512775099dap-88L,
+	0x1.6c012750bdabeed76a9980p0L, 0xf.4f33fdeb8b0ecd831106f57b3d00p-96L,
+	0x1.6dfb23c651a2ef220e2cbep0L, 0x1.bbaa834b3f11577ceefbe6c1c411p-92L,
+	0x1.6ff7df9519483cf87e1b4ep0L, 0x1.3e213bff9b702d5aa477c12523cep-88L,
+	0x1.71f75e8ec5f73dd2370f2ep0L, 0xf.0acd6cb434b562d9e8a20adda648p-92L,
+	0x1.73f9a48a58173bd5c9a4e6p0L, 0x8.ab1182ae217f3a7681759553e840p-92L,
+	0x1.75feb564267c8bf6e9aa32p0L, 0x1.a48b27071805e61a17b954a2dad8p-88L,
+	0x1.780694fde5d3f619ae0280p0L, 0x8.58b2bb2bdcf86cd08e35fb04c0f0p-92L,
+	0x1.7a11473eb0186d7d51023ep0L, 0x1.6cda1f5ef42b66977960531e821bp-88L,
+	0x1.7c1ed0130c1327c4933444p0L, 0x1.937562b2dc933d44fc828efd4c9cp-88L,
+	0x1.7e2f336cf4e62105d02ba0p0L, 0x1.5797e170a1427f8fcdf5f3906108p-88L,
+	0x1.80427543e1a11b60de6764p0L, 0x9.a354ea706b8e4d8b718a672bf7c8p-92L,
+	0x1.82589994cce128acf88afap0L, 0xb.34a010f6ad65cbbac0f532d39be0p-92L,
+	0x1.8471a4623c7acce52f6b96p0L, 0x1.c64095370f51f48817914dd78665p-88L,
+	0x1.868d99b4492ec80e41d90ap0L, 0xc.251707484d73f136fb5779656b70p-92L,
+	0x1.88ac7d98a669966530bcdep0L, 0x1.2d4e9d61283ef385de170ab20f96p-88L,
+	0x1.8ace5422aa0db5ba7c55a0p0L, 0x1.92c9bb3e6ed61f2733304a346d8fp-88L,
+	0x1.8cf3216b5448bef2aa1cd0p0L, 0x1.61c55d84a9848f8c453b3ca8c946p-88L,
+	0x1.8f1ae991577362b982745cp0L, 0x7.2ed804efc9b4ae1458ae946099d4p-92L,
+	0x1.9145b0b91ffc588a61b468p0L, 0x1.f6b70e01c2a90229a4c4309ea719p-88L,
+	0x1.93737b0cdc5e4f4501c3f2p0L, 0x5.40a22d2fc4af581b63e8326efe9cp-92L,
+	0x1.95a44cbc8520ee9b483694p0L, 0x1.a0fc6f7c7d61b2b3a22a0eab2cadp-88L,
+	0x1.97d829fde4e4f8b9e920f8p0L, 0x1.1e8bd7edb9d7144b6f6818084cc7p-88L,
+	0x1.9a0f170ca07b9ba3109b8cp0L, 0x4.6737beb19e1eada6825d3c557428p-92L,
+	0x1.9c49182a3f0901c7c46b06p0L, 0x1.1f2be58ddade50c217186c90b457p-88L,
+	0x1.9e86319e323231824ca78ep0L, 0x6.4c6e010f92c082bbadfaf605cfd4p-92L,
+	0x1.a0c667b5de564b29ada8b8p0L, 0xc.ab349aa0422a8da7d4512edac548p-92L,
+	0x1.a309bec4a2d3358c171f76p0L, 0x1.0daad547fa22c26d168ea762d854p-88L,
+	0x1.a5503b23e255c8b424491cp0L, 0xa.f87bc8050a405381703ef7caff50p-92L,
+	0x1.a799e1330b3586f2dfb2b0p0L, 0x1.58f1a98796ce8908ae852236ca94p-88L,
+	0x1.a9e6b5579fdbf43eb243bcp0L, 0x1.ff4c4c58b571cf465caf07b4b9f5p-88L,
+	0x1.ac36bbfd3f379c0db966a2p0L, 0x1.1265fc73e480712d20f8597a8e7bp-88L,
+	0x1.ae89f995ad3ad5e8734d16p0L, 0x1.73205a7fbc3ae675ea440b162d6cp-88L,
+	0x1.b0e07298db66590842acdep0L, 0x1.c6f6ca0e5dcae2aafffa7a0554cbp-88L,
+	0x1.b33a2b84f15faf6bfd0e7ap0L, 0x1.d947c2575781dbb49b1237c87b6ep-88L,
+	0x1.b59728de559398e3881110p0L, 0x1.64873c7171fefc410416be0a6525p-88L,
+	0x1.b7f76f2fb5e46eaa7b081ap0L, 0xb.53c5354c8903c356e4b625aacc28p-92L,
+	0x1.ba5b030a10649840cb3c6ap0L, 0xf.5b47f297203757e1cc6eadc8bad0p-92L,
+	0x1.bcc1e904bc1d2247ba0f44p0L, 0x1.b3d08cd0b20287092bd59be4ad98p-88L,
+	0x1.bf2c25bd71e088408d7024p0L, 0x1.18e3449fa073b356766dfb568ff4p-88L,
+	0x1.c199bdd85529c2220cb12ap0L, 0x9.1ba6679444964a36661240043970p-96L,
+	0x1.c40ab5fffd07a6d14df820p0L, 0xf.1828a5366fd387a7bdd54cdf7300p-92L,
+	0x1.c67f12e57d14b4a2137fd2p0L, 0xf.2b301dd9e6b151a6d1f9d5d5f520p-96L,
+	0x1.c8f6d9406e7b511acbc488p0L, 0x5.c442ddb55820171f319d9e5076a8p-96L,
+	0x1.cb720dcef90691503cbd1ep0L, 0x9.49db761d9559ac0cb6dd3ed599e0p-92L,
+	0x1.cdf0b555dc3f9c44f8958ep0L, 0x1.ac51be515f8c58bdfb6f5740a3a4p-88L,
+	0x1.d072d4a07897b8d0f22f20p0L, 0x1.a158e18fbbfc625f09f4cca40874p-88L,
+	0x1.d2f87080d89f18ade12398p0L, 0x9.ea2025b4c56553f5cdee4c924728p-92L,
+	0x1.d5818dcfba48725da05aeap0L, 0x1.66e0dca9f589f559c0876ff23830p-88L,
+	0x1.d80e316c98397bb84f9d04p0L, 0x8.805f84bec614de269900ddf98d28p-92L,
+	0x1.da9e603db3285708c01a5ap0L, 0x1.6d4c97f6246f0ec614ec95c99392p-88L,
+	0x1.dd321f301b4604b695de3cp0L, 0x6.30a393215299e30d4fb73503c348p-96L,
+	0x1.dfc97337b9b5eb968cac38p0L, 0x1.ed291b7225a944efd5bb5524b927p-88L,
+	0x1.e264614f5a128a12761fa0p0L, 0x1.7ada6467e77f73bf65e04c95e29dp-88L,
+	0x1.e502ee78b3ff6273d13014p0L, 0x1.3991e8f49659e1693be17ae1d2f9p-88L,
+	0x1.e7a51fbc74c834b548b282p0L, 0x1.23786758a84f4956354634a416cep-88L,
+	0x1.ea4afa2a490d9858f73a18p0L, 0xf.5db301f86dea20610ceee13eb7b8p-92L,
+	0x1.ecf482d8e67f08db0312fap0L, 0x1.949cef462010bb4bc4ce72a900dfp-88L,
+	0x1.efa1bee615a27771fd21a8p0L, 0x1.2dac1f6dd5d229ff68e46f27e3dfp-88L,
+	0x1.f252b376bba974e8696fc2p0L, 0x1.6390d4c6ad5476b5162f40e1d9a9p-88L,
+	0x1.f50765b6e4540674f84b76p0L, 0x2.862baff99000dfc4352ba29b8908p-92L,
+	0x1.f7bfdad9cbe138913b4bfep0L, 0x7.2bd95c5ce7280fa4d2344a3f5618p-92L,
+	0x1.fa7c1819e90d82e90a7e74p0L, 0xb.263c1dc060c36f7650b4c0f233a8p-92L,
+	0x1.fd3c22b8f71f10975ba4b2p0L, 0x1.2bcf3a5e12d269d8ad7c1a4a8875p-88L
+};
+
+/*
+ * Kernel for expl(x).  x must be finite and not tiny or huge.
+ * "tiny" is anything that would make us underflow (|A6*x^6| < ~LDBL_MIN).
+ * "huge" is anything that would make fn*L1 inexact (|x| > ~2**17*ln2).
+ */
+static inline void
+__k_expl(long double x, long double *hip, long double *lop, int *kp)
+{
+	long double q, r, r1, t;
+	double dr, fn, r2;
+	int n, n2;
+
+	/* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
+	/* Use a specialized rint() to get fn.  Assume round-to-nearest. */
+	/* XXX assume no extra precision for the additions, as for trig fns. */
+	/* XXX this set of comments is now quadruplicated. */
+	/* XXX but see ../src/e_exp.c for a fix using double_t. */
+	fn = (double)x * INV_L + 0x1.8p52 - 0x1.8p52;
+#if defined(HAVE_EFFICIENT_IRINT)
+	n = irint(fn);
+#else
+	n = (int)fn;
+#endif
+	n2 = (unsigned)n % INTERVALS;
+	/* Depend on the sign bit being propagated: */
+	*kp = n >> LOG2_INTERVALS;
+	r1 = x - fn * L1;
+	r2 = fn * -L2;
+	r = r1 + r2;
+
+	/* Evaluate expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2). */
+	dr = r;
+	q = r2 + r * r * (A2 + r * (A3 + r * (A4 + r * (A5 + r * (A6 +
+	    dr * (A7 + dr * (A8 + dr * (A9 + dr * A10))))))));
+	t = tbl[n2].lo + tbl[n2].hi;
+	*hip = tbl[n2].hi;
+	*lop = tbl[n2].lo + t * (q + r1);
+}
+
+/*
+ * XXX: the rest of the functions are identical for ld80 and ld128.
+ * However, we should use scalbnl() for ld128, since long double
+ * multiplication is very slow on the only supported ld128 arch (sparc64).
+ */
+
+static inline void
+k_hexpl(long double x, long double *hip, long double *lop)
+{
+	float twopkm1;
+	int k;
+
+	__k_expl(x, hip, lop, &k);
+	SET_FLOAT_WORD(twopkm1, 0x3f800000 + ((k - 1) << 23));
+	*hip *= twopkm1;
+	*lop *= twopkm1;
+}
+
+static inline long double
+hexpl(long double x)
+{
+	long double hi, lo, twopkm2;
+	int k;
+
+	twopkm2 = 1;
+	__k_expl(x, &hi, &lo, &k);
+	SET_LDBL_EXPSIGN(twopkm2, BIAS + k - 2);
+	return (lo + hi) * 2 * twopkm2;
+}
+
+#ifdef _COMPLEX_H
+/*
+ * See ../src/k_exp.c for details.
+ */
+static inline long double complex
+__ldexp_cexpl(long double complex z, int expt)
+{
+	long double exp_x, hi, lo;
+	long double x, y, scale1, scale2;
+	int half_expt, k;
+
+	x = creall(z);
+	y = cimagl(z);
+	__k_expl(x, &hi, &lo, &k);
+
+	exp_x = (lo + hi) * 0x1p16382;
+	expt += k - 16382;
+
+	scale1 = 1;
+	half_expt = expt / 2;
+	SET_LDBL_EXPSIGN(scale1, BIAS + half_expt);
+	scale2 = 1;
+	SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
+
+	return (cpackl(cos(y) * exp_x * scale1 * scale2,
+	    sinl(y) * exp_x * scale1 * scale2));
+}
+#endif /* _COMPLEX_H */
diff --git a/lib/msun/ld128/s_erfl.c b/lib/msun/ld128/s_erfl.c
new file mode 100644
index 0000000..e29c969
--- /dev/null
+++ b/lib/msun/ld128/s_erfl.c
@@ -0,0 +1,329 @@
+/* @(#)s_erf.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See s_erf.c for complete comments.
+ *
+ * Converted to long double by Steven G. Kargl.
+ */
+#include <float.h>
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* XXX Prevent compilers from erroneously constant folding these: */
+static const volatile long double tiny = 0x1p-10000L;
+
+static const double
+half= 0.5,
+one = 1,
+two = 2;
+/*
+ * In the domain [0, 2**-40], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-120.
+ */
+static const long double
+efx  =  1.28379167095512573896158903121545167e-01L,	/* 0xecbff6a7, 0x481dd788, 0xb64d21a8, 0xeb06fc3f */
+efx8 =  1.02703333676410059116927122497236133e+00L,	/* 0xecbff6a7, 0x481dd788, 0xb64d21a8, 0xeb06ff3f */
+/*
+ * Domain [0, 0.84375], range ~[-1.919e-38, 1.919e-38]:
+ * |(erf(x) - x)/x - pp(x)/qq(x)| < 2**-125.29
+ */
+pp0  =  1.28379167095512573896158903121545167e-01L,	/* 0x3ffc06eb, 0xa8214db6, 0x88d71d48, 0xa7f6bfec */
+pp1  = -3.14931554396568573802046931159683404e-01L,	/* 0xbffd427d, 0x6ada7263, 0x547eb096, 0x95f37463 */
+pp2  = -5.27514920282183487103576956956725309e-02L,	/* 0xbffab023, 0xe5a271e3, 0xb0e79b01, 0x2f7ac962 */
+pp3  = -1.13202828509005281355609495523452713e-02L,	/* 0xbff872f1, 0x6a5023a1, 0xe08b3884, 0x326af20f */
+pp4  = -9.18626155872522453865998391206048506e-04L,	/* 0xbff4e19f, 0xea5fb024, 0x43247a37, 0xe430b06c */
+pp5  = -7.87518862406176274922506447157284230e-05L,	/* 0xbff14a4f, 0x31a85fe0, 0x7fff2204, 0x09c49b37 */
+pp6  = -3.42357944472240436548115331090560881e-06L,	/* 0xbfeccb81, 0x4b43c336, 0xcd2eb6c2, 0x903f2d87 */
+pp7  = -1.37317432573890412634717890726745428e-07L,	/* 0xbfe826e3, 0x0e915eb6, 0x42aee414, 0xf7e36805 */
+pp8  = -2.71115170113861755855049008732113726e-09L,	/* 0xbfe2749e, 0x2b94fd00, 0xecb4d166, 0x0efb91f8 */
+pp9  = -3.37925756196555959454018189718117864e-11L,	/* 0xbfdc293e, 0x1d9060cb, 0xd043204a, 0x314cd7f0 */
+qq1  =  4.76672625471551170489978555182449450e-01L,	/* 0x3ffde81c, 0xde6531f0, 0x76803bee, 0x526e29e9 */
+qq2  =  1.06713144672281502058807525850732240e-01L,	/* 0x3ffbb518, 0xd7a6bb74, 0xcd9bdd33, 0x7601eee5 */
+qq3  =  1.47747613127513761102189201923147490e-02L,	/* 0x3ff8e423, 0xae527e18, 0xf12cb447, 0x723b4749 */
+qq4  =  1.39939377672028671891148770908874816e-03L,	/* 0x3ff56ed7, 0xba055d84, 0xc21b45c4, 0x388d1812 */
+qq5  =  9.44302939359455241271983309378738276e-05L,	/* 0x3ff18c11, 0xc18c99a4, 0x86d0fe09, 0x46387b4c */
+qq6  =  4.56199342312522842161301671745365650e-06L,	/* 0x3fed3226, 0x73421d05, 0x08875300, 0x32fa1432 */
+qq7  =  1.53019260483764773845294600092361197e-07L,	/* 0x3fe8489b, 0x3a63f627, 0x2b9ad2ce, 0x26516e57 */
+qq8  =  3.25542691121324805094777901250005508e-09L,	/* 0x3fe2bf6c, 0x26d93a29, 0x9142be7c, 0x9f1dd043 */
+qq9  =  3.37405581964478060434410167262684979e-11L;	/* 0x3fdc28c8, 0xfb8fa1be, 0x10e57eec, 0xaa19e49f */
+
+static const long double
+erx  =  8.42700792949714894142232424201210961e-01L,	/* 0x3ffeaf76, 0x7a741088, 0xb0000000, 0x00000000 */
+/*
+ * Domain [0.84375, 1.25], range ~[-2.521e-36, 2.523e-36]:
+ * |(erf(x) - erx) - pa(x)/qa(x)| < 2**-120.15
+ */
+pa0  = -2.48010117891186017024438233323795897e-17L,	/* 0xbfc7c97f, 0x77812279, 0x6c877f22, 0xef4bfb2e */
+pa1  =  4.15107497420594680894327969504526489e-01L,	/* 0x3ffda911, 0xf096fbc2, 0x55662005, 0x2337fa64 */
+pa2  = -3.94180628087084846724448515851892609e-02L,	/* 0xbffa42e9, 0xab54528c, 0xad529da1, 0x6efc2af3 */
+pa3  =  4.48897599625192107295954790681677462e-02L,	/* 0x3ffa6fbc, 0xa65edba1, 0x0e4cbcea, 0x73ef9a31 */
+pa4  =  8.02069252143016600110972019232995528e-02L,	/* 0x3ffb4887, 0x0e8b548e, 0x3230b417, 0x11b553b3 */
+pa5  = -1.02729816533435279443621120242391295e-02L,	/* 0xbff850a0, 0x041de3ee, 0xd5bca6c9, 0x4ef5f9f2 */
+pa6  =  5.70777694530755634864821094419982095e-03L,	/* 0x3ff77610, 0x9b501e10, 0x4c978382, 0x742df68f */
+pa7  =  1.22635150233075521018231779267077071e-03L,	/* 0x3ff5417b, 0x0e623682, 0x60327da0, 0x96b9219e */
+pa8  =  5.36100234820204569428412542856666503e-04L,	/* 0x3ff41912, 0x27ceb4c1, 0x1d3298ec, 0x84ced627 */
+pa9  = -1.97753571846365167177187858667583165e-04L,	/* 0xbff29eb8, 0x23f5bcf3, 0x15c83c46, 0xe4fda98b */
+pa10 =  6.19333039900846970674794789568415105e-05L,	/* 0x3ff103c4, 0x60f88e46, 0xc0c9fb02, 0x13cc7fc1 */
+pa11 = -5.40531400436645861492290270311751349e-06L,	/* 0xbfed6abe, 0x9665f8a8, 0xdd0ad3ba, 0xe5dc0ee3 */
+qa1  =  9.05041313265490487793231810291907851e-01L,	/* 0x3ffecf61, 0x93340222, 0xe9930620, 0xc4e61168 */
+qa2  =  6.79848064708886864767240880834868092e-01L,	/* 0x3ffe5c15, 0x0ba858dc, 0xf7900ae9, 0xfea1e09a */
+qa3  =  4.04720609926471677581066689316516445e-01L,	/* 0x3ffd9e6f, 0x145e9b00, 0x6d8c1749, 0xd2928623 */
+qa4  =  1.69183273898369996364661075664302225e-01L,	/* 0x3ffc5a7c, 0xc2a363c1, 0xd6c19097, 0xef9b4063 */
+qa5  =  7.44476185988067992342479750486764248e-02L,	/* 0x3ffb30ef, 0xfc7259ef, 0x1bcbb089, 0x686dd62d */
+qa6  =  2.02981172725892407200420389604788573e-02L,	/* 0x3ff94c90, 0x7976cb0e, 0x21e1d36b, 0x0f09ca2b */
+qa7  =  6.94281866271607668268269403102277234e-03L,	/* 0x3ff7c701, 0x2b193250, 0xc5d46ecc, 0x374843d8 */
+qa8  =  1.12952275469171559611651594706820034e-03L,	/* 0x3ff52818, 0xfd2a7c06, 0xd13e38fd, 0xda4b34f5 */
+qa9  =  3.13736683241992737197226578597710179e-04L,	/* 0x3ff348fa, 0x0cb48d18, 0x051f849b, 0x135ccf74 */
+qa10 =  1.17037675204033225470121134087771410e-05L,	/* 0x3fee88b6, 0x98f47704, 0xa5d8f8f2, 0xc6422e11 */
+qa11 =  4.61312518293853991439362806880973592e-06L,	/* 0x3fed3594, 0xe31db94f, 0x3592b693, 0xed4386b4 */
+qa12 = -1.02158572037456893687737553657431771e-06L;	/* 0xbfeb123a, 0xd60d9b1e, 0x1f6fdeb9, 0x7dc8410a */
+/*
+ * Domain [1.25,2.85715], range ~[-2.922e-37,2.922e-37]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - ra(x)/sa(x)| < 2**-121.36
+ */
+static const long double
+ra0  = -9.86494292470069009555706994426014461e-03L,	/* 0xbff84341, 0x239e8709, 0xe941b06a, 0xcb4b6ec5 */
+ra1  = -1.13580436992565640457579040117568870e+00L,	/* 0xbfff22c4, 0x133f7c0d, 0x72d5e231, 0x2eb1ee3f */
+ra2  = -4.89744330295291950661185707066921755e+01L,	/* 0xc00487cb, 0xa38b4fc2, 0xc136695b, 0xc1df8047 */
+ra3  = -1.10766149300215937173768072715352140e+03L,	/* 0xc00914ea, 0x55e6beb3, 0xabc50e07, 0xb6e5664d */
+ra4  = -1.49991031232170934967642795601952100e+04L,	/* 0xc00cd4b8, 0xd33243e6, 0xffbf6545, 0x3c57ef6e */
+ra5  = -1.29805749738318462882524181556996692e+05L,	/* 0xc00ffb0d, 0xbfeed9b6, 0x5b2a3ff4, 0xe245bd3c */
+ra6  = -7.42828497044940065828871976644647850e+05L,	/* 0xc0126ab5, 0x8fe7caca, 0x473352d9, 0xcd4e0c90 */
+ra7  = -2.85637299581890734287995171242421106e+06L,	/* 0xc0145cad, 0xa7f76fe7, 0x3e358051, 0x1799f927 */
+ra8  = -7.40674797129824999383748865571026084e+06L,	/* 0xc015c412, 0x6fe29c02, 0x298ad158, 0x7d24e45c */
+ra9  = -1.28653420911930973914078724204151759e+07L,	/* 0xc016889e, 0x7c2eb0dc, 0x95d5863b, 0x0aa34dc3 */
+ra10 = -1.47198163599330179552932489109452638e+07L,	/* 0xc016c136, 0x90b84923, 0xf9bcb497, 0x19bbd0f5 */
+ra11 = -1.07812992258382800318665248311522624e+07L,	/* 0xc0164904, 0xe673a113, 0x35d7f079, 0xe13701f3 */
+ra12 = -4.83545565681708642630419905537756076e+06L,	/* 0xc0152721, 0xfea094a8, 0x869eb39d, 0x413d6f13 */
+ra13 = -1.23956521201673964822976917356685286e+06L,	/* 0xc0132ea0, 0xd3646baa, 0x2fe62b0d, 0xbae5ce85 */
+ra14 = -1.62289333553652417591275333240371812e+05L,	/* 0xc0103cf8, 0xaab1e2d6, 0x4c25e014, 0x248d76ab */
+ra15 = -8.82890392601176969729168894389833110e+03L,	/* 0xc00c13e7, 0x3b3d8f94, 0x6fbda6f6, 0xe7049a82 */
+ra16 = -1.22591866337261720023681535568334619e+02L,	/* 0xc005ea5e, 0x12358891, 0xcfa712c5, 0x77f050d4 */
+sa1  =  6.44508918884710829371852723353794047e+01L,	/* 0x400501cd, 0xb69a6c0f, 0x5716de14, 0x47161af6 */
+sa2  =  1.76118475473171481523704824327358534e+03L,	/* 0x4009b84b, 0xd305829f, 0xc4c771b0, 0xbf1f7f9b */
+sa3  =  2.69448346969488374857087646131950188e+04L,	/* 0x400da503, 0x56bacc05, 0x4fdba68d, 0x2cca27e6 */
+sa4  =  2.56826633369941456778326497384543763e+05L,	/* 0x4010f59d, 0x51124428, 0x69c41de6, 0xbd0d5753 */
+sa5  =  1.60647413092257206847700054645905859e+06L,	/* 0x40138834, 0xa2184244, 0x557a1bed, 0x68c9d556 */
+sa6  =  6.76963075165099718574753447122393797e+06L,	/* 0x40159d2f, 0x7b01b0cc, 0x8bac9e95, 0x5d35d56e */
+sa7  =  1.94295690905361884290986932493647741e+07L,	/* 0x40172878, 0xc1172d61, 0x3068501e, 0x2f3c71da */
+sa8  =  3.79774781017759149060839255547073541e+07L,	/* 0x401821be, 0xc30d06fe, 0x410563d7, 0x032111fd */
+sa9  =  5.00659831846029484248302236457727397e+07L,	/* 0x40187df9, 0x1f97a111, 0xc51d6ac2, 0x4b389793 */
+sa10 =  4.36486287620506484276130525941972541e+07L,	/* 0x40184d03, 0x3a618ae0, 0x2a723357, 0xfa45c60a */
+sa11 =  2.43779678791333894255510508253951934e+07L,	/* 0x401773fa, 0x6fe10ee2, 0xc467850d, 0xc6b7ff30 */
+sa12 =  8.30732360384443202039372372212966542e+06L,	/* 0x4015fb09, 0xee6a5631, 0xdd98de7e, 0x8b00461a */
+sa13 =  1.60160846942050515734192397495105693e+06L,	/* 0x40138704, 0x8782bf13, 0x5b8fb315, 0xa898abe5 */
+sa14 =  1.54255505242533291014555153757001825e+05L,	/* 0x40102d47, 0xc0abc98e, 0x843c9490, 0xb4352440 */
+sa15 =  5.87949220002375547561467275493888824e+03L,	/* 0x400b6f77, 0xe00d21d1, 0xec4d41e8, 0x2f8e1673 */
+sa16 =  4.97272976346793193860385983372237710e+01L;	/* 0x40048dd1, 0x816c1b3f, 0x24f540a6, 0x4cfe03cc */
+/*
+ * Domain [2.85715,9], range ~[-7.886e-37,7.918e-37]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rb(x)/sb(x)| < 2**-120
+ */
+static const long double
+rb0  = -9.86494292470008707171371994479162369e-3L, /* 0xbff84341, 0x239e86f4, 0x2f57e561, 0xf4469360 */
+rb1  = -1.57047326624110727986326503729442830L,    /* 0xbfff920a, 0x8935bf73, 0x8803b894, 0x4656482d */
+rb2  = -1.03228196364885474342132255440317065e2L,  /* 0xc0059ce9, 0xac4ed0ff, 0x2cff0ff7, 0x5e70d1ab */
+rb3  = -3.74000570653418227179358710865224376e3L,  /* 0xc00ad380, 0x2ebf7835, 0xf6b07ed2, 0x861242f7 */
+rb4  = -8.35435477739098044190860390632813956e4L,  /* 0xc00f4657, 0x8c3ae934, 0x3647d7b3, 0x80e76fb7 */
+rb5  = -1.21398672055223642118716640216747152e6L,  /* 0xc0132862, 0x2b8761c8, 0x27d18c0f, 0x137c9463 */
+rb6  = -1.17669175877248796101665344873273970e7L,  /* 0xc0166719, 0x0b2cea46, 0x81f14174, 0x11602ea5 */
+rb7  = -7.66108006086998253606773064264599615e7L,  /* 0xc019243f, 0x3c26f4f0, 0x1cc05241, 0x3b953728 */
+rb8  = -3.32547117558141845968704725353130804e8L,  /* 0xc01b3d24, 0x42d8ee26, 0x24ef6f3b, 0x604a8c65 */
+rb9  = -9.41561252426350696802167711221739746e8L,  /* 0xc01cc0f8, 0xad23692a, 0x8ddb2310, 0xe9937145 */
+rb10 = -1.67157110805390944549427329626281063e9L,  /* 0xc01d8e88, 0x9a903734, 0x09a55fa3, 0xd205c903 */
+rb11 = -1.74339631004410841337645931421427373e9L,  /* 0xc01d9fa8, 0x77582d2a, 0xc183b8ab, 0x7e00cb05 */
+rb12 = -9.57655233596934915727573141357471703e8L,  /* 0xc01cc8a5, 0x460cc685, 0xd0271fa0, 0x6a70e3da */
+rb13 = -2.26320062731339353035254704082495066e8L,  /* 0xc01aafab, 0xd7d76721, 0xc9720e11, 0x6a8bd489 */
+rb14 = -1.42777302996263256686002973851837039e7L,  /* 0xc016b3b8, 0xc499689f, 0x2b88d965, 0xc32414f9 */
+sb1  =  1.08512869705594540211033733976348506e2L,  /* 0x4005b20d, 0x2db7528d, 0x00d20dcb, 0x858f6191 */
+sb2  =  5.02757713761390460534494530537572834e3L,  /* 0x400b3a39, 0x3bf4a690, 0x3025d28d, 0xfd40a891 */
+sb3  =  1.31019107205412870059331647078328430e5L,  /* 0x400fffcb, 0x1b71d05e, 0x3b28361d, 0x2a3c3690 */
+sb4  =  2.13021555152296846166736757455018030e6L,  /* 0x40140409, 0x3c6984df, 0xc4491d7c, 0xb04aa08d */
+sb5  =  2.26649105281820861953868568619768286e7L,  /* 0x401759d6, 0xce8736f0, 0xf28ad037, 0x2a901e0c */
+sb6  =  1.61071939490875921812318684143076081e8L,  /* 0x401a3338, 0x686fb541, 0x6bd27d06, 0x4f95c9ac */
+sb7  =  7.66895673844301852676056750497991966e8L,  /* 0x401c6daf, 0x31cec121, 0x54699126, 0x4bd9bf9e */
+sb8  =  2.41884450436101936436023058196042526e9L,  /* 0x401e2059, 0x46b0b8d7, 0x87b64cbf, 0x78bc296d */
+sb9  =  4.92403055884071695093305291535107666e9L,  /* 0x401f257e, 0xbe5ed739, 0x39e17346, 0xcadd2e55 */
+sb10 =  6.18627786365587486459633615573786416e9L,  /* 0x401f70bb, 0x1be7a7e7, 0x6a45b5ae, 0x607c70f0 */
+sb11 =  4.45898013426501378097430226324743199e9L,  /* 0x401f09c6, 0xa32643d7, 0xf1724620, 0x9ea46c32 */
+sb12 =  1.63006115763329848117160344854224975e9L,  /* 0x401d84a3, 0x0996887f, 0x65a4f43b, 0x978c1d74 */
+sb13 =  2.39216717012421697446304015847567721e8L,  /* 0x401ac845, 0x09a065c2, 0x30095da7, 0x9d72d6ae */
+sb14 =  7.84837329009278694937250358810225609e6L;  /* 0x4015df06, 0xd5290e15, 0x63031fac, 0x4d9c894c */
+/*
+ * Domain [9,108], range ~[-5.324e-38,5.340e-38]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-124
+ */
+static const long double
+rc0  = -9.86494292470008707171367567652935673e-3L, /* 0xbff84341, 0x239e86f4, 0x2f57e55b, 0x1aa10fd3 */
+rc1  = -1.26229447747315096406518846411562266L,    /* 0xbfff4325, 0xbb1aab28, 0xda395cd9, 0xfb861c15 */
+rc2  = -6.13742634438922591780742637728666162e1L,  /* 0xc004eafe, 0x7dd51cd8, 0x3c7c5928, 0x751e50cf */
+rc3  = -1.50455835478908280402912854338421517e3L,  /* 0xc0097823, 0xbc15b9ab, 0x3d60745c, 0x523e80a5 */
+rc4  = -2.04415631865861549920184039902945685e4L,  /* 0xc00d3f66, 0x40b3fc04, 0x5388f2ec, 0xb009e1f0 */
+rc5  = -1.57625662981714582753490610560037638e5L,  /* 0xc01033dc, 0xd4dc95b6, 0xfd4da93b, 0xf355b4a9 */
+rc6  = -6.73473451616752528402917538033283794e5L,  /* 0xc01248d8, 0x2e73a4f9, 0xcded49c5, 0xfa3bfeb7 */
+rc7  = -1.47433165421387483167186683764364857e6L,  /* 0xc01367f1, 0xba77a8f7, 0xcfdd0dbb, 0x25d554b3 */
+rc8  = -1.38811981807868828563794929997744139e6L,  /* 0xc01352e5, 0x7d16d9ad, 0xbbdcbf38, 0x38fbc5ea */
+rc9  = -3.59659700530831825640766479698155060e5L,  /* 0xc0115f3a, 0xecd57f45, 0x21f8ad6c, 0x910a5958 */
+sc1  =  7.72730753022908298637508998072635696e1L,  /* 0x40053517, 0xa10d52bc, 0xdabb55b6, 0xbd0328cd */
+sc2  =  2.36825757341694050500333261769082182e3L,  /* 0x400a2808, 0x3e0a9b42, 0x82977842, 0x9c5de29e */
+sc3  =  3.72210540173034735352888847134073099e4L,  /* 0x400e22ca, 0x1ba827ef, 0xac8390d7, 0x1fc39a41 */
+sc4  =  3.24136032646418336712461033591393412e5L,  /* 0x40113c8a, 0x0216e100, 0xc59d1e44, 0xf0e68d9d */
+sc5  =  1.57836135851134393802505823370009175e6L,  /* 0x40138157, 0x95bc7664, 0x17575961, 0xdbe58eeb */
+sc6  =  4.12881981392063738026679089714182355e6L,  /* 0x4014f801, 0x9e82e8d2, 0xb8b3a70e, 0xfd84185d */
+sc7  =  5.24438427289213488410596395361544142e6L,  /* 0x40154017, 0x81177109, 0x2aa6c3b0, 0x1f106625 */
+sc8  =  2.59909544563616121735963429710382149e6L,  /* 0x40143d45, 0xbb90a9b1, 0x12bf9390, 0xa827a700 */
+sc9  =  2.80930665169282501639651995082335693e5L;  /* 0x40111258, 0xaa92222e, 0xa97e3216, 0xa237fa6c */
+
+long double
+erfl(long double x)
+{
+	long double ax,R,S,P,Q,s,y,z,r;
+	uint64_t lx, llx;
+	int32_t i;
+	uint16_t hx;
+
+	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+
+	if((hx & 0x7fff) == 0x7fff) {	/* erfl(nan)=nan */
+		i = (hx>>15)<<1;
+		return (1-i)+one/x;	/* erfl(+-inf)=+-1 */
+	}
+
+	ax = fabsl(x);
+	if(ax < 0.84375) {
+	    if(ax < 0x1p-40L) {
+	        if(ax < 0x1p-16373L)	
+		    return (8*x+efx8*x)/8;	/* avoid spurious underflow */
+		return x + efx*x;
+	    }
+	    z = x*x;
+	    r = pp0+z*(pp1+z*(pp2+z*(pp3+z*(pp4+z*(pp5+z*(pp6+z*(pp7+
+		z*(pp8+z*pp9))))))));
+	    s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*(qq5+z*(qq6+z*(qq7+
+		z*(qq8+z*qq9))))))));
+	    y = r/s;
+	    return x + x*y;
+	}
+	if(ax < 1.25) {
+	    s = ax-one;
+	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*(pa6+s*(pa7+
+		s*(pa8+s*(pa9+s*(pa10+s*pa11))))))))));
+	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*(qa6+s*(qa7+
+		s*(qa8+s*(qa9+s*(qa10+s*(qa11+s*qa12)))))))))));
+	    if(x>=0) return (erx + P/Q); else return (-erx - P/Q);
+	}
+	if (ax >= 9) {			/* inf>|x|>= 9 */
+	    if(x>=0) return (one-tiny); else return (tiny-one);
+	}
+	s = one/(ax*ax);
+	if(ax < 2.85715) {	/* |x| < 2.85715 */
+	    R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*(ra7+
+		s*(ra8+s*(ra9+s*(ra10+s*(ra11+s*(ra12+s*(ra13+s*(ra14+
+		s*(ra15+s*ra16)))))))))))))));
+	    S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		s*(sa8+s*(sa9+s*(sa10+s*(sa11+s*(sa12+s*(sa13+s*(sa14+
+		s*(sa15+s*sa16)))))))))))))));
+	} else {	/* |x| >= 2.85715 */
+	    R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*(rb6+s*(rb7+
+		s*(rb8+s*(rb9+s*(rb10+s*(rb11+s*(rb12+s*(rb13+
+		s*rb14)))))))))))));
+	    S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*(sb7+
+		s*(sb8+s*(sb9+s*(sb10+s*(sb11+s*(sb12+s*(sb13+
+		s*sb14)))))))))))));
+	}
+	z = (float)ax;
+	r = expl(-z*z-0.5625)*expl((z-ax)*(z+ax)+R/S);
+	if(x>=0) return (one-r/ax); else return (r/ax-one);
+}
+
+long double
+erfcl(long double x)
+{
+	long double ax,R,S,P,Q,s,y,z,r;
+	uint64_t lx, llx;
+	uint16_t hx;
+
+	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+
+	if((hx & 0x7fff) == 0x7fff) {	/* erfcl(nan)=nan */
+					/* erfcl(+-inf)=0,2 */
+	    return ((hx>>15)<<1)+one/x;
+	}
+
+	ax = fabsl(x);
+	if(ax < 0.84375L) {
+	    if(ax < 0x1p-34L)
+		return one-x;
+	    z = x*x;
+	    r = pp0+z*(pp1+z*(pp2+z*(pp3+z*(pp4+z*(pp5+z*(pp6+z*(pp7+
+		z*(pp8+z*pp9))))))));
+	    s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*(qq5+z*(qq6+z*(qq7+
+		z*(qq8+z*qq9))))))));
+	    y = r/s;
+	    if(ax < 0.25L) {  	/* x<1/4 */
+		return one-(x+x*y);
+	    } else {
+		r = x*y;
+		r += (x-half);
+	       return half - r;
+	    }
+	}
+	if(ax < 1.25L) {
+	    s = ax-one;
+	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*(pa6+s*(pa7+
+		    s*(pa8+s*(pa9+s*(pa10+s*pa11))))))))));
+	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*(qa6+s*(qa7+
+		    s*(qa8+s*(qa9+s*(qa10+s*(qa11+s*qa12)))))))))));
+	    if(x>=0) {
+	        z  = one-erx; return z - P/Q;
+	    } else {
+		z = erx+P/Q; return one+z;
+	    }
+	}
+
+	if(ax < 108) {			/* |x| < 108 */
+ 	    s = one/(ax*ax);
+	    if(ax < 2.85715) {		/* |x| < 2.85715 */
+	        R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*(ra7+
+		    s*(ra8+s*(ra9+s*(ra10+s*(ra11+s*(ra12+s*(ra13+s*(ra14+
+		    s*(ra15+s*ra16)))))))))))))));
+	        S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		    s*(sa8+s*(sa9+s*(sa10+s*(sa11+s*(sa12+s*(sa13+s*(sa14+
+		    s*(sa15+s*sa16)))))))))))))));
+	    } else if(ax < 9) {
+		R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*(rb6+s*(rb7+
+		    s*(rb8+s*(rb9+s*(rb10+s*(rb11+s*(rb12+s*(rb13+
+		    s*rb14)))))))))))));
+		S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*(sb7+
+		    s*(sb8+s*(sb9+s*(sb10+s*(sb11+s*(sb12+s*(sb13+
+		    s*sb14)))))))))))));
+	    } else {
+		if(x < -9) return two-tiny;	/* x < -9 */
+		R=rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*(rc5+s*(rc6+s*(rc7+
+		    s*(rc8+s*rc9))))))));
+		S=one+s*(sc1+s*(sc2+s*(sc3+s*(sc4+s*(sc5+s*(sc6+s*(sc7+
+		    s*(sc8+s*sc9))))))));
+	    }
+	    z = (float)ax;
+	    r = expl(-z*z-0.5625)*expl((z-ax)*(z+ax)+R/S);
+	    if(x>0) return r/ax; else return two-r/ax;
+	} else {
+	    if(x>0) return tiny*tiny; else return two-tiny;
+	}
+}
diff --git a/lib/msun/ld128/s_expl.c b/lib/msun/ld128/s_expl.c
index 176c932..1357e70 100644
--- a/lib/msun/ld128/s_expl.c
+++ b/lib/msun/ld128/s_expl.c
@@ -38,16 +38,15 @@ __FBSDID("$FreeBSD$");
 #include "fpmath.h"
 #include "math.h"
 #include "math_private.h"
+#include "k_expl.h"
 
-#define	INTERVALS	128
-#define	LOG2_INTERVALS	7
-#define	BIAS	(LDBL_MAX_EXP - 1)
+/* XXX Prevent compilers from erroneously constant folding these: */
+static const volatile long double
+huge = 0x1p10000L,
+tiny = 0x1p-10000L;
 
 static const long double
-huge = 0x1p10000L,
 twom10000 = 0x1p-10000L;
-/* XXX Prevent gcc from erroneously constant folding this: */
-static volatile const long double tiny = 0x1p-10000L;
 
 static const long double
 /* log(2**16384 - 0.5) rounded towards zero: */
@@ -56,184 +55,16 @@ o_threshold =  11356.523406294143949491931077970763428L,
 /* log(2**(-16381-64-1)) rounded towards zero: */
 u_threshold = -11433.462743336297878837243843452621503L;
 
-static const double
-/*
- * ln2/INTERVALS = L1+L2 (hi+lo decomposition for multiplication).  L1 must
- * have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest
- * bits zero so that multiplication of it by n is exact.
- */
-INV_L = 1.8466496523378731e+2,		/*  0x171547652b82fe.0p-45 */
-L2 = -1.0253670638894731e-29;		/* -0x1.9ff0342542fc3p-97 */
-static const long double
-/* 0x1.62e42fefa39ef35793c768000000p-8 */
-L1 =  5.41521234812457272982212595914567508e-3L;
-
-static const long double
-/*
- * Domain [-0.002708, 0.002708], range ~[-2.4021e-38, 2.4234e-38]:
- * |exp(x) - p(x)| < 2**-124.9
- * (0.002708 is ln2/(2*INTERVALS) rounded up a little).
- */
-A2  =  0.5,
-A3  =  1.66666666666666666666666666651085500e-1L,
-A4  =  4.16666666666666666666666666425885320e-2L,
-A5  =  8.33333333333333333334522877160175842e-3L,
-A6  =  1.38888888888888888889971139751596836e-3L;
-
-static const double
-A7  =  1.9841269841269471e-4,
-A8  =  2.4801587301585284e-5,
-A9  =  2.7557324277411234e-6,
-A10 =  2.7557333722375072e-7;
-
-static const struct {
-	/*
-	 * hi must be rounded to at most 106 bits so that multiplication
-	 * by r1 in expm1l() is exact, but it is rounded to 88 bits due to
-	 * historical accidents.
-	 */
-	long double	hi;
-	long double	lo;
-} tbl[INTERVALS] = {
-	0x1p0L, 0x0p0L,
-	0x1.0163da9fb33356d84a66aep0L, 0x3.36dcdfa4003ec04c360be2404078p-92L,
-	0x1.02c9a3e778060ee6f7cacap0L, 0x4.f7a29bde93d70a2cabc5cb89ba10p-92L,
-	0x1.04315e86e7f84bd738f9a2p0L, 0xd.a47e6ed040bb4bfc05af6455e9b8p-96L,
-	0x1.059b0d31585743ae7c548ep0L, 0xb.68ca417fe53e3495f7df4baf84a0p-92L,
-	0x1.0706b29ddf6ddc6dc403a8p0L, 0x1.d87b27ed07cb8b092ac75e311753p-88L,
-	0x1.0874518759bc808c35f25cp0L, 0x1.9427fa2b041b2d6829d8993a0d01p-88L,
-	0x1.09e3ecac6f3834521e060cp0L, 0x5.84d6b74ba2e023da730e7fccb758p-92L,
-	0x1.0b5586cf9890f6298b92b6p0L, 0x1.1842a98364291408b3ceb0a2a2bbp-88L,
-	0x1.0cc922b7247f7407b705b8p0L, 0x9.3dc5e8aac564e6fe2ef1d431fd98p-92L,
-	0x1.0e3ec32d3d1a2020742e4ep0L, 0x1.8af6a552ac4b358b1129e9f966a4p-88L,
-	0x1.0fb66affed31af232091dcp0L, 0x1.8a1426514e0b627bda694a400a27p-88L,
-	0x1.11301d0125b50a4ebbf1aep0L, 0xd.9318ceac5cc47ab166ee57427178p-92L,
-	0x1.12abdc06c31cbfb92bad32p0L, 0x4.d68e2f7270bdf7cedf94eb1cb818p-92L,
-	0x1.1429aaea92ddfb34101942p0L, 0x1.b2586d01844b389bea7aedd221d4p-88L,
-	0x1.15a98c8a58e512480d573cp0L, 0x1.d5613bf92a2b618ee31b376c2689p-88L,
-	0x1.172b83c7d517adcdf7c8c4p0L, 0x1.0eb14a792035509ff7d758693f24p-88L,
-	0x1.18af9388c8de9bbbf70b9ap0L, 0x3.c2505c97c0102e5f1211941d2840p-92L,
-	0x1.1a35beb6fcb753cb698f68p0L, 0x1.2d1c835a6c30724d5cfae31b84e5p-88L,
-	0x1.1bbe084045cd39ab1e72b4p0L, 0x4.27e35f9acb57e473915519a1b448p-92L,
-	0x1.1d4873168b9aa7805b8028p0L, 0x9.90f07a98b42206e46166cf051d70p-92L,
-	0x1.1ed5022fcd91cb8819ff60p0L, 0x1.121d1e504d36c47474c9b7de6067p-88L,
-	0x1.2063b88628cd63b8eeb028p0L, 0x1.50929d0fc487d21c2b84004264dep-88L,
-	0x1.21f49917ddc962552fd292p0L, 0x9.4bdb4b61ea62477caa1dce823ba0p-92L,
-	0x1.2387a6e75623866c1fadb0p0L, 0x1.c15cb593b0328566902df69e4de2p-88L,
-	0x1.251ce4fb2a63f3582ab7dep0L, 0x9.e94811a9c8afdcf796934bc652d0p-92L,
-	0x1.26b4565e27cdd257a67328p0L, 0x1.d3b249dce4e9186ddd5ff44e6b08p-92L,
-	0x1.284dfe1f5638096cf15cf0p0L, 0x3.ca0967fdaa2e52d7c8106f2e262cp-92L,
-	0x1.29e9df51fdee12c25d15f4p0L, 0x1.a24aa3bca890ac08d203fed80a07p-88L,
-	0x1.2b87fd0dad98ffddea4652p0L, 0x1.8fcab88442fdc3cb6de4519165edp-88L,
-	0x1.2d285a6e4030b40091d536p0L, 0xd.075384589c1cd1b3e4018a6b1348p-92L,
-	0x1.2ecafa93e2f5611ca0f45cp0L, 0x1.523833af611bdcda253c554cf278p-88L,
-	0x1.306fe0a31b7152de8d5a46p0L, 0x3.05c85edecbc27343629f502f1af2p-92L,
-	0x1.32170fc4cd8313539cf1c2p0L, 0x1.008f86dde3220ae17a005b6412bep-88L,
-	0x1.33c08b26416ff4c9c8610cp0L, 0x1.96696bf95d1593039539d94d662bp-88L,
-	0x1.356c55f929ff0c94623476p0L, 0x3.73af38d6d8d6f9506c9bbc93cbc0p-92L,
-	0x1.371a7373aa9caa7145502ep0L, 0x1.4547987e3e12516bf9c699be432fp-88L,
-	0x1.38cae6d05d86585a9cb0d8p0L, 0x1.bed0c853bd30a02790931eb2e8f0p-88L,
-	0x1.3a7db34e59ff6ea1bc9298p0L, 0x1.e0a1d336163fe2f852ceeb134067p-88L,
-	0x1.3c32dc313a8e484001f228p0L, 0xb.58f3775e06ab66353001fae9fca0p-92L,
-	0x1.3dea64c12342235b41223ep0L, 0x1.3d773fba2cb82b8244267c54443fp-92L,
-	0x1.3fa4504ac801ba0bf701aap0L, 0x4.1832fb8c1c8dbdff2c49909e6c60p-92L,
-	0x1.4160a21f72e29f84325b8ep0L, 0x1.3db61fb352f0540e6ba05634413ep-88L,
-	0x1.431f5d950a896dc7044394p0L, 0x1.0ccec81e24b0caff7581ef4127f7p-92L,
-	0x1.44e086061892d03136f408p0L, 0x1.df019fbd4f3b48709b78591d5cb5p-88L,
-	0x1.46a41ed1d005772512f458p0L, 0x1.229d97df404ff21f39c1b594d3a8p-88L,
-	0x1.486a2b5c13cd013c1a3b68p0L, 0x1.062f03c3dd75ce8757f780e6ec99p-88L,
-	0x1.4a32af0d7d3de672d8bcf4p0L, 0x6.f9586461db1d878b1d148bd3ccb8p-92L,
-	0x1.4bfdad5362a271d4397afep0L, 0xc.42e20e0363ba2e159c579f82e4b0p-92L,
-	0x1.4dcb299fddd0d63b36ef1ap0L, 0x9.e0cc484b25a5566d0bd5f58ad238p-92L,
-	0x1.4f9b2769d2ca6ad33d8b68p0L, 0x1.aa073ee55e028497a329a7333dbap-88L,
-	0x1.516daa2cf6641c112f52c8p0L, 0x4.d822190e718226177d7608d20038p-92L,
-	0x1.5342b569d4f81df0a83c48p0L, 0x1.d86a63f4e672a3e429805b049465p-88L,
-	0x1.551a4ca5d920ec52ec6202p0L, 0x4.34ca672645dc6c124d6619a87574p-92L,
-	0x1.56f4736b527da66ecb0046p0L, 0x1.64eb3c00f2f5ab3d801d7cc7272dp-88L,
-	0x1.58d12d497c7fd252bc2b72p0L, 0x1.43bcf2ec936a970d9cc266f0072fp-88L,
-	0x1.5ab07dd48542958c930150p0L, 0x1.91eb345d88d7c81280e069fbdb63p-88L,
-	0x1.5c9268a5946b701c4b1b80p0L, 0x1.6986a203d84e6a4a92f179e71889p-88L,
-	0x1.5e76f15ad21486e9be4c20p0L, 0x3.99766a06548a05829e853bdb2b52p-92L,
-	0x1.605e1b976dc08b076f592ap0L, 0x4.86e3b34ead1b4769df867b9c89ccp-92L,
-	0x1.6247eb03a5584b1f0fa06ep0L, 0x1.d2da42bb1ceaf9f732275b8aef30p-88L,
-	0x1.6434634ccc31fc76f8714cp0L, 0x4.ed9a4e41000307103a18cf7a6e08p-92L,
-	0x1.66238825522249127d9e28p0L, 0x1.b8f314a337f4dc0a3adf1787ff74p-88L,
-	0x1.68155d44ca973081c57226p0L, 0x1.b9f32706bfe4e627d809a85dcc66p-88L,
-	0x1.6a09e667f3bcc908b2fb12p0L, 0x1.66ea957d3e3adec17512775099dap-88L,
-	0x1.6c012750bdabeed76a9980p0L, 0xf.4f33fdeb8b0ecd831106f57b3d00p-96L,
-	0x1.6dfb23c651a2ef220e2cbep0L, 0x1.bbaa834b3f11577ceefbe6c1c411p-92L,
-	0x1.6ff7df9519483cf87e1b4ep0L, 0x1.3e213bff9b702d5aa477c12523cep-88L,
-	0x1.71f75e8ec5f73dd2370f2ep0L, 0xf.0acd6cb434b562d9e8a20adda648p-92L,
-	0x1.73f9a48a58173bd5c9a4e6p0L, 0x8.ab1182ae217f3a7681759553e840p-92L,
-	0x1.75feb564267c8bf6e9aa32p0L, 0x1.a48b27071805e61a17b954a2dad8p-88L,
-	0x1.780694fde5d3f619ae0280p0L, 0x8.58b2bb2bdcf86cd08e35fb04c0f0p-92L,
-	0x1.7a11473eb0186d7d51023ep0L, 0x1.6cda1f5ef42b66977960531e821bp-88L,
-	0x1.7c1ed0130c1327c4933444p0L, 0x1.937562b2dc933d44fc828efd4c9cp-88L,
-	0x1.7e2f336cf4e62105d02ba0p0L, 0x1.5797e170a1427f8fcdf5f3906108p-88L,
-	0x1.80427543e1a11b60de6764p0L, 0x9.a354ea706b8e4d8b718a672bf7c8p-92L,
-	0x1.82589994cce128acf88afap0L, 0xb.34a010f6ad65cbbac0f532d39be0p-92L,
-	0x1.8471a4623c7acce52f6b96p0L, 0x1.c64095370f51f48817914dd78665p-88L,
-	0x1.868d99b4492ec80e41d90ap0L, 0xc.251707484d73f136fb5779656b70p-92L,
-	0x1.88ac7d98a669966530bcdep0L, 0x1.2d4e9d61283ef385de170ab20f96p-88L,
-	0x1.8ace5422aa0db5ba7c55a0p0L, 0x1.92c9bb3e6ed61f2733304a346d8fp-88L,
-	0x1.8cf3216b5448bef2aa1cd0p0L, 0x1.61c55d84a9848f8c453b3ca8c946p-88L,
-	0x1.8f1ae991577362b982745cp0L, 0x7.2ed804efc9b4ae1458ae946099d4p-92L,
-	0x1.9145b0b91ffc588a61b468p0L, 0x1.f6b70e01c2a90229a4c4309ea719p-88L,
-	0x1.93737b0cdc5e4f4501c3f2p0L, 0x5.40a22d2fc4af581b63e8326efe9cp-92L,
-	0x1.95a44cbc8520ee9b483694p0L, 0x1.a0fc6f7c7d61b2b3a22a0eab2cadp-88L,
-	0x1.97d829fde4e4f8b9e920f8p0L, 0x1.1e8bd7edb9d7144b6f6818084cc7p-88L,
-	0x1.9a0f170ca07b9ba3109b8cp0L, 0x4.6737beb19e1eada6825d3c557428p-92L,
-	0x1.9c49182a3f0901c7c46b06p0L, 0x1.1f2be58ddade50c217186c90b457p-88L,
-	0x1.9e86319e323231824ca78ep0L, 0x6.4c6e010f92c082bbadfaf605cfd4p-92L,
-	0x1.a0c667b5de564b29ada8b8p0L, 0xc.ab349aa0422a8da7d4512edac548p-92L,
-	0x1.a309bec4a2d3358c171f76p0L, 0x1.0daad547fa22c26d168ea762d854p-88L,
-	0x1.a5503b23e255c8b424491cp0L, 0xa.f87bc8050a405381703ef7caff50p-92L,
-	0x1.a799e1330b3586f2dfb2b0p0L, 0x1.58f1a98796ce8908ae852236ca94p-88L,
-	0x1.a9e6b5579fdbf43eb243bcp0L, 0x1.ff4c4c58b571cf465caf07b4b9f5p-88L,
-	0x1.ac36bbfd3f379c0db966a2p0L, 0x1.1265fc73e480712d20f8597a8e7bp-88L,
-	0x1.ae89f995ad3ad5e8734d16p0L, 0x1.73205a7fbc3ae675ea440b162d6cp-88L,
-	0x1.b0e07298db66590842acdep0L, 0x1.c6f6ca0e5dcae2aafffa7a0554cbp-88L,
-	0x1.b33a2b84f15faf6bfd0e7ap0L, 0x1.d947c2575781dbb49b1237c87b6ep-88L,
-	0x1.b59728de559398e3881110p0L, 0x1.64873c7171fefc410416be0a6525p-88L,
-	0x1.b7f76f2fb5e46eaa7b081ap0L, 0xb.53c5354c8903c356e4b625aacc28p-92L,
-	0x1.ba5b030a10649840cb3c6ap0L, 0xf.5b47f297203757e1cc6eadc8bad0p-92L,
-	0x1.bcc1e904bc1d2247ba0f44p0L, 0x1.b3d08cd0b20287092bd59be4ad98p-88L,
-	0x1.bf2c25bd71e088408d7024p0L, 0x1.18e3449fa073b356766dfb568ff4p-88L,
-	0x1.c199bdd85529c2220cb12ap0L, 0x9.1ba6679444964a36661240043970p-96L,
-	0x1.c40ab5fffd07a6d14df820p0L, 0xf.1828a5366fd387a7bdd54cdf7300p-92L,
-	0x1.c67f12e57d14b4a2137fd2p0L, 0xf.2b301dd9e6b151a6d1f9d5d5f520p-96L,
-	0x1.c8f6d9406e7b511acbc488p0L, 0x5.c442ddb55820171f319d9e5076a8p-96L,
-	0x1.cb720dcef90691503cbd1ep0L, 0x9.49db761d9559ac0cb6dd3ed599e0p-92L,
-	0x1.cdf0b555dc3f9c44f8958ep0L, 0x1.ac51be515f8c58bdfb6f5740a3a4p-88L,
-	0x1.d072d4a07897b8d0f22f20p0L, 0x1.a158e18fbbfc625f09f4cca40874p-88L,
-	0x1.d2f87080d89f18ade12398p0L, 0x9.ea2025b4c56553f5cdee4c924728p-92L,
-	0x1.d5818dcfba48725da05aeap0L, 0x1.66e0dca9f589f559c0876ff23830p-88L,
-	0x1.d80e316c98397bb84f9d04p0L, 0x8.805f84bec614de269900ddf98d28p-92L,
-	0x1.da9e603db3285708c01a5ap0L, 0x1.6d4c97f6246f0ec614ec95c99392p-88L,
-	0x1.dd321f301b4604b695de3cp0L, 0x6.30a393215299e30d4fb73503c348p-96L,
-	0x1.dfc97337b9b5eb968cac38p0L, 0x1.ed291b7225a944efd5bb5524b927p-88L,
-	0x1.e264614f5a128a12761fa0p0L, 0x1.7ada6467e77f73bf65e04c95e29dp-88L,
-	0x1.e502ee78b3ff6273d13014p0L, 0x1.3991e8f49659e1693be17ae1d2f9p-88L,
-	0x1.e7a51fbc74c834b548b282p0L, 0x1.23786758a84f4956354634a416cep-88L,
-	0x1.ea4afa2a490d9858f73a18p0L, 0xf.5db301f86dea20610ceee13eb7b8p-92L,
-	0x1.ecf482d8e67f08db0312fap0L, 0x1.949cef462010bb4bc4ce72a900dfp-88L,
-	0x1.efa1bee615a27771fd21a8p0L, 0x1.2dac1f6dd5d229ff68e46f27e3dfp-88L,
-	0x1.f252b376bba974e8696fc2p0L, 0x1.6390d4c6ad5476b5162f40e1d9a9p-88L,
-	0x1.f50765b6e4540674f84b76p0L, 0x2.862baff99000dfc4352ba29b8908p-92L,
-	0x1.f7bfdad9cbe138913b4bfep0L, 0x7.2bd95c5ce7280fa4d2344a3f5618p-92L,
-	0x1.fa7c1819e90d82e90a7e74p0L, 0xb.263c1dc060c36f7650b4c0f233a8p-92L,
-	0x1.fd3c22b8f71f10975ba4b2p0L, 0x1.2bcf3a5e12d269d8ad7c1a4a8875p-88L
-};
-
 long double
 expl(long double x)
 {
-	union IEEEl2bits u, v;
-	long double q, r, r1, t, twopk, twopkp10000;
-	double dr, fn, r2;
-	int k, n, n2;
+	union IEEEl2bits u;
+	long double hi, lo, t, twopk;
+	int k;
 	uint16_t hx, ix;
 
+	DOPRINT_START(&x);
+
 	/* Filter out exceptional cases. */
 	u.e = x;
 	hx = u.xbits.expsign;
@@ -241,60 +72,33 @@ expl(long double x)
 	if (ix >= BIAS + 13) {		/* |x| >= 8192 or x is NaN */
 		if (ix == BIAS + LDBL_MAX_EXP) {
 			if (hx & 0x8000)  /* x is -Inf or -NaN */
-				return (-1 / x);
-			return (x + x);	/* x is +Inf or +NaN */
+				RETURNP(-1 / x);
+			RETURNP(x + x);	/* x is +Inf or +NaN */
 		}
 		if (x > o_threshold)
-			return (huge * huge);
+			RETURNP(huge * huge);
 		if (x < u_threshold)
-			return (tiny * tiny);
+			RETURNP(tiny * tiny);
 	} else if (ix < BIAS - 114) {	/* |x| < 0x1p-114 */
-		return (1 + x);		/* 1 with inexact iff x != 0 */
+		RETURN2P(1, x);		/* 1 with inexact iff x != 0 */
 	}
 
 	ENTERI();
 
-	/* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
-	/* Use a specialized rint() to get fn.  Assume round-to-nearest. */
-	/* XXX assume no extra precision for the additions, as for trig fns. */
-	/* XXX this set of comments is now quadruplicated. */
-	fn = (double)x * INV_L + 0x1.8p52 - 0x1.8p52;
-#if 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 factors. */
-	/* XXX sparc64 multiplication is so slow that scalbnl() is faster. */
-	v.e = 1;
-	if (k >= LDBL_MIN_EXP) {
-		v.xbits.expsign = BIAS + k;
-		twopk = v.e;
-	} else {
-		v.xbits.expsign = BIAS + k + 10000;
-		twopkp10000 = v.e;
-	}
-
-	/* Evaluate expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2). */
-	dr = r;
-	q = r2 + r * r * (A2 + r * (A3 + r * (A4 + r * (A5 + r * (A6 +
-	    dr * (A7 + dr * (A8 + dr * (A9 + dr * A10))))))));
-	t = tbl[n2].lo + tbl[n2].hi;
-	t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+	twopk = 1;
+	__k_expl(x, &hi, &lo, &k);
+	t = SUM2P(hi, lo);
 
 	/* Scale by 2**k. */
+	/* XXX sparc64 multiplication is so slow that scalbnl() is faster. */
 	if (k >= LDBL_MIN_EXP) {
 		if (k == LDBL_MAX_EXP)
 			RETURNI(t * 2 * 0x1p16383L);
+		SET_LDBL_EXPSIGN(twopk, BIAS + k);
 		RETURNI(t * twopk);
 	} else {
-		RETURNI(t * twopkp10000 * twom10000);
+		SET_LDBL_EXPSIGN(twopk, BIAS + k + 10000);
+		RETURNI(t * twopk * twom10000);
 	}
 }
 
@@ -312,6 +116,12 @@ expl(long double x)
  * Setting T3 to 0 would require the |x| < 0x1p-113 condition to appear
  * in both subintervals, so set T3 = 2**-5, which places the condition
  * into the [T1, T3] interval.
+ *
+ * XXX we now do this more to (partially) balance the number of terms
+ * in the C and D polys than to avoid checking the condition in both
+ * intervals.
+ *
+ * XXX these micro-optimizations are excessive.
  */
 static const double
 T1 = -0.1659,				/* ~-30.625/128 * log(2) */
@@ -321,6 +131,12 @@ T3 =  0.03125;
 /*
  * Domain [-0.1659, 0.03125], range ~[2.9134e-44, 1.8404e-37]:
  * |(exp(x)-1-x-x**2/2)/x - p(x)| < 2**-122.03
+/*
+ * XXX none of the long double C or D coeffs except C10 is correctly printed.
+ * If you re-print their values in %.35Le format, the result is always
+ * different.  For example, the last 2 digits in C3 should be 59, not 67.
+ * 67 is apparently from rounding an extra-precision value to 36 decimal
+ * places.
  */
 static const long double
 C3  =  1.66666666666666666666666666666666667e-1L,
@@ -335,6 +151,13 @@ C11 =  2.50521083854417203619031960151253944e-8L,
 C12 =  2.08767569878679576457272282566520649e-9L,
 C13 =  1.60590438367252471783548748824255707e-10L;
 
+/*
+ * XXX this has 1 more coeff than needed.
+ * XXX can start the double coeffs but not the double mults at C10.
+ * With my coeffs (C10-C17 double; s = best_s):
+ * Domain [-0.1659, 0.03125], range ~[-1.1976e-37, 1.1976e-37]:
+ * |(exp(x)-1-x-x**2/2)/x - p(x)| ~< 2**-122.65
+ */
 static const double
 C14 =  1.1470745580491932e-11,		/*  0x1.93974a81dae30p-37 */
 C15 =  7.6471620181090468e-13,		/*  0x1.ae7f3820adab1p-41 */
@@ -359,6 +182,13 @@ D11 =  2.50521083855084570046480450935267433e-8L,
 D12 =  2.08767569819738524488686318024854942e-9L,
 D13 =  1.60590442297008495301927448122499313e-10L;
 
+/*
+ * XXX this has 1 more coeff than needed.
+ * XXX can start the double coeffs but not the double mults at D11.
+ * With my coeffs (D11-D16 double):
+ * Domain [0.03125, 0.1659], range ~[-1.1980e-37, 1.1980e-37]:
+ * |(exp(x)-1-x-x**2/2)/x - p(x)| ~< 2**-122.65
+ */
 static const double
 D14 =  1.1470726176204336e-11,		/*  0x1.93971dc395d9ep-37 */
 D15 =  7.6478532249581686e-13,		/*  0x1.ae892e3D16fcep-41 */
@@ -375,6 +205,8 @@ expm1l(long double x)
 	int k, n, n2;
 	uint16_t hx, ix;
 
+	DOPRINT_START(&x);
+
 	/* Filter out exceptional cases. */
 	u.e = x;
 	hx = u.xbits.expsign;
@@ -382,11 +214,11 @@ expm1l(long double x)
 	if (ix >= BIAS + 7) {		/* |x| >= 128 or x is NaN */
 		if (ix == BIAS + LDBL_MAX_EXP) {
 			if (hx & 0x8000)  /* x is -Inf or -NaN */
-				return (-1 / x - 1);
-			return (x + x);	/* x is +Inf or +NaN */
+				RETURNP(-1 / x - 1);
+			RETURNP(x + x);	/* x is +Inf or +NaN */
 		}
 		if (x > o_threshold)
-			return (huge * huge);
+			RETURNP(huge * huge);
 		/*
 		 * expm1l() never underflows, but it must avoid
 		 * unrepresentable large negative exponents.  We used a
@@ -395,7 +227,7 @@ expm1l(long double x)
 		 * in the same way as large ones here.
 		 */
 		if (hx & 0x8000)	/* x <= -128 */
-			return (tiny - 1);	/* good for x < -114ln2 - eps */
+			RETURN2P(tiny, -1);	/* good for x < -114ln2 - eps */
 	}
 
 	ENTERI();
@@ -407,7 +239,7 @@ expm1l(long double x)
 		if (x < T3) {
 			if (ix < BIAS - 113) {	/* |x| < 0x1p-113 */
 				/* x (rounded) with inexact if x != 0: */
-				RETURNI(x == 0 ? x :
+				RETURNPI(x == 0 ? x :
 				    (0x1p200 * x + fabsl(x)) * 0x1p-200);
 			}
 			q = x * x2 * C3 + x2 * x2 * (C4 + x * (C5 + x * (C6 +
@@ -428,9 +260,9 @@ expm1l(long double x)
 		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));
+			RETURN2PI(hx2_hi + x_hi, hx2_lo + x_lo + q);
 		else
-			RETURNI(hx2_lo + q + hx2_hi + x);
+			RETURN2PI(x, hx2_lo + q + hx2_hi);
 	}
 
 	/* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
@@ -463,21 +295,21 @@ expm1l(long double x)
 	t = 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);
+		t = SUM2P(tbl[n2].hi - 1, tbl[n2].lo * (r1 + 1) + t * q +
+		    tbl[n2].hi * r1);
 		RETURNI(t);
 	}
 	if (k == -1) {
-		t = tbl[n2].lo * (r1 + 1) + t * q + tbl[n2].hi * r1 + 
-		    (tbl[n2].hi - 2);
+		t = SUM2P(tbl[n2].hi - 2, tbl[n2].lo * (r1 + 1) + t * q +
+		    tbl[n2].hi * r1);
 		RETURNI(t / 2);
 	}
 	if (k < -7) {
-		t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+		t = SUM2P(tbl[n2].hi, tbl[n2].lo + t * (q + r1));
 		RETURNI(t * twopk - 1);
 	}
 	if (k > 2 * LDBL_MANT_DIG - 1) {
-		t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+		t = SUM2P(tbl[n2].hi, tbl[n2].lo + t * (q + r1));
 		if (k == LDBL_MAX_EXP)
 			RETURNI(t * 2 * 0x1p16383L - 1);
 		RETURNI(t * twopk - 1);
@@ -487,8 +319,8 @@ expm1l(long double x)
 	twomk = v.e;
 
 	if (k > LDBL_MANT_DIG - 1)
-		t = tbl[n2].lo - twomk + t * (q + r1) + tbl[n2].hi;
+		t = SUM2P(tbl[n2].hi, tbl[n2].lo - twomk + t * (q + r1));
 	else
-		t = tbl[n2].lo + t * (q + r1) + (tbl[n2].hi - twomk);
+		t = SUM2P(tbl[n2].hi - twomk, tbl[n2].lo + t * (q + r1));
 	RETURNI(t * twopk);
 }
diff --git a/lib/msun/ld80/k_expl.h b/lib/msun/ld80/k_expl.h
new file mode 100644
index 0000000..ebfb9a8
--- /dev/null
+++ b/lib/msun/ld80/k_expl.h
@@ -0,0 +1,305 @@
+/* from: FreeBSD: head/lib/msun/ld80/s_expl.c 251343 2013-06-03 19:51:32Z kargl */
+
+/*-
+ * Copyright (c) 2009-2013 Steven G. Kargl
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice unmodified, this list of conditions, and the following
+ *    disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See s_expl.c for more comments about __k_expl().
+ *
+ * See ../src/e_exp.c and ../src/k_exp.h for precision-independent comments
+ * about the secondary kernels.
+ */
+
+#define	INTERVALS	128
+#define	LOG2_INTERVALS	7
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+static const double
+/*
+ * ln2/INTERVALS = L1+L2 (hi+lo decomposition for multiplication).  L1 must
+ * have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest
+ * bits zero so that multiplication of it by n is exact.
+ */
+INV_L = 1.8466496523378731e+2,		/*  0x171547652b82fe.0p-45 */
+L1 =  5.4152123484527692e-3,		/*  0x162e42ff000000.0p-60 */
+L2 = -3.2819649005320973e-13,		/* -0x1718432a1b0e26.0p-94 */
+/*
+ * Domain [-0.002708, 0.002708], range ~[-5.7136e-24, 5.7110e-24]:
+ * |exp(x) - p(x)| < 2**-77.2
+ * (0.002708 is ln2/(2*INTERVALS) rounded up a little).
+ */
+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
+ * the first 53 bits of the significand are stored in hi and the next 53
+ * bits are in lo.  Tang's paper states that the trailing 6 bits of hi must
+ * be zero for his algorithm in both single and double precision, because
+ * the table is re-used in the implementation of expm1() where a floating
+ * point addition involving hi must be exact.  Here hi is double, so
+ * converting it to long double gives 11 trailing zero bits.
+ */
+static const struct {
+	double	hi;
+	double	lo;
+} tbl[INTERVALS] = {
+	0x1p+0, 0x0p+0,
+	/*
+	 * XXX hi is rounded down, and the formatting is not quite normal.
+	 * But I rather like both.  The 0x1.*p format is good for 4N+1
+	 * mantissa bits.  Rounding down makes the lo terms positive,
+	 * so that the columnar formatting can be simpler.
+	 */
+	0x1.0163da9fb3335p+0, 0x1.b61299ab8cdb7p-54,
+	0x1.02c9a3e778060p+0, 0x1.dcdef95949ef4p-53,
+	0x1.04315e86e7f84p+0, 0x1.7ae71f3441b49p-53,
+	0x1.059b0d3158574p+0, 0x1.d73e2a475b465p-55,
+	0x1.0706b29ddf6ddp+0, 0x1.8db880753b0f6p-53,
+	0x1.0874518759bc8p+0, 0x1.186be4bb284ffp-57,
+	0x1.09e3ecac6f383p+0, 0x1.1487818316136p-54,
+	0x1.0b5586cf9890fp+0, 0x1.8a62e4adc610bp-54,
+	0x1.0cc922b7247f7p+0, 0x1.01edc16e24f71p-54,
+	0x1.0e3ec32d3d1a2p+0, 0x1.03a1727c57b53p-59,
+	0x1.0fb66affed31ap+0, 0x1.e464123bb1428p-53,
+	0x1.11301d0125b50p+0, 0x1.49d77e35db263p-53,
+	0x1.12abdc06c31cbp+0, 0x1.f72575a649ad2p-53,
+	0x1.1429aaea92ddfp+0, 0x1.66820328764b1p-53,
+	0x1.15a98c8a58e51p+0, 0x1.2406ab9eeab0ap-55,
+	0x1.172b83c7d517ap+0, 0x1.b9bef918a1d63p-53,
+	0x1.18af9388c8de9p+0, 0x1.777ee1734784ap-53,
+	0x1.1a35beb6fcb75p+0, 0x1.e5b4c7b4968e4p-55,
+	0x1.1bbe084045cd3p+0, 0x1.3563ce56884fcp-53,
+	0x1.1d4873168b9aap+0, 0x1.e016e00a2643cp-54,
+	0x1.1ed5022fcd91cp+0, 0x1.71033fec2243ap-53,
+	0x1.2063b88628cd6p+0, 0x1.dc775814a8495p-55,
+	0x1.21f49917ddc96p+0, 0x1.2a97e9494a5eep-55,
+	0x1.2387a6e756238p+0, 0x1.9b07eb6c70573p-54,
+	0x1.251ce4fb2a63fp+0, 0x1.ac155bef4f4a4p-55,
+	0x1.26b4565e27cddp+0, 0x1.2bd339940e9d9p-55,
+	0x1.284dfe1f56380p+0, 0x1.2d9e2b9e07941p-53,
+	0x1.29e9df51fdee1p+0, 0x1.612e8afad1255p-55,
+	0x1.2b87fd0dad98fp+0, 0x1.fbbd48ca71f95p-53,
+	0x1.2d285a6e4030bp+0, 0x1.0024754db41d5p-54,
+	0x1.2ecafa93e2f56p+0, 0x1.1ca0f45d52383p-56,
+	0x1.306fe0a31b715p+0, 0x1.6f46ad23182e4p-55,
+	0x1.32170fc4cd831p+0, 0x1.a9ce78e18047cp-55,
+	0x1.33c08b26416ffp+0, 0x1.32721843659a6p-54,
+	0x1.356c55f929ff0p+0, 0x1.928c468ec6e76p-53,
+	0x1.371a7373aa9cap+0, 0x1.4e28aa05e8a8fp-53,
+	0x1.38cae6d05d865p+0, 0x1.0b53961b37da2p-53,
+	0x1.3a7db34e59ff6p+0, 0x1.d43792533c144p-53,
+	0x1.3c32dc313a8e4p+0, 0x1.08003e4516b1ep-53,
+	0x1.3dea64c123422p+0, 0x1.ada0911f09ebcp-55,
+	0x1.3fa4504ac801bp+0, 0x1.417ee03548306p-53,
+	0x1.4160a21f72e29p+0, 0x1.f0864b71e7b6cp-53,
+	0x1.431f5d950a896p+0, 0x1.b8e088728219ap-53,
+	0x1.44e086061892dp+0, 0x1.89b7a04ef80d0p-59,
+	0x1.46a41ed1d0057p+0, 0x1.c944bd1648a76p-54,
+	0x1.486a2b5c13cd0p+0, 0x1.3c1a3b69062f0p-56,
+	0x1.4a32af0d7d3dep+0, 0x1.9cb62f3d1be56p-54,
+	0x1.4bfdad5362a27p+0, 0x1.d4397afec42e2p-56,
+	0x1.4dcb299fddd0dp+0, 0x1.8ecdbbc6a7833p-54,
+	0x1.4f9b2769d2ca6p+0, 0x1.5a67b16d3540ep-53,
+	0x1.516daa2cf6641p+0, 0x1.8225ea5909b04p-53,
+	0x1.5342b569d4f81p+0, 0x1.be1507893b0d5p-53,
+	0x1.551a4ca5d920ep+0, 0x1.8a5d8c4048699p-53,
+	0x1.56f4736b527dap+0, 0x1.9bb2c011d93adp-54,
+	0x1.58d12d497c7fdp+0, 0x1.295e15b9a1de8p-55,
+	0x1.5ab07dd485429p+0, 0x1.6324c054647adp-54,
+	0x1.5c9268a5946b7p+0, 0x1.c4b1b816986a2p-60,
+	0x1.5e76f15ad2148p+0, 0x1.ba6f93080e65ep-54,
+	0x1.605e1b976dc08p+0, 0x1.60edeb25490dcp-53,
+	0x1.6247eb03a5584p+0, 0x1.63e1f40dfa5b5p-53,
+	0x1.6434634ccc31fp+0, 0x1.8edf0e2989db3p-53,
+	0x1.6623882552224p+0, 0x1.224fb3c5371e6p-53,
+	0x1.68155d44ca973p+0, 0x1.038ae44f73e65p-57,
+	0x1.6a09e667f3bccp+0, 0x1.21165f626cdd5p-53,
+	0x1.6c012750bdabep+0, 0x1.daed533001e9ep-53,
+	0x1.6dfb23c651a2ep+0, 0x1.e441c597c3775p-53,
+	0x1.6ff7df9519483p+0, 0x1.9f0fc369e7c42p-53,
+	0x1.71f75e8ec5f73p+0, 0x1.ba46e1e5de15ap-53,
+	0x1.73f9a48a58173p+0, 0x1.7ab9349cd1562p-53,
+	0x1.75feb564267c8p+0, 0x1.7edd354674916p-53,
+	0x1.780694fde5d3fp+0, 0x1.866b80a02162dp-54,
+	0x1.7a11473eb0186p+0, 0x1.afaa2047ed9b4p-53,
+	0x1.7c1ed0130c132p+0, 0x1.f124cd1164dd6p-54,
+	0x1.7e2f336cf4e62p+0, 0x1.05d02ba15797ep-56,
+	0x1.80427543e1a11p+0, 0x1.6c1bccec9346bp-53,
+	0x1.82589994cce12p+0, 0x1.159f115f56694p-53,
+	0x1.8471a4623c7acp+0, 0x1.9ca5ed72f8c81p-53,
+	0x1.868d99b4492ecp+0, 0x1.01c83b21584a3p-53,
+	0x1.88ac7d98a6699p+0, 0x1.994c2f37cb53ap-54,
+	0x1.8ace5422aa0dbp+0, 0x1.6e9f156864b27p-54,
+	0x1.8cf3216b5448bp+0, 0x1.de55439a2c38bp-53,
+	0x1.8f1ae99157736p+0, 0x1.5cc13a2e3976cp-55,
+	0x1.9145b0b91ffc5p+0, 0x1.114c368d3ed6ep-53,
+	0x1.93737b0cdc5e4p+0, 0x1.e8a0387e4a814p-53,
+	0x1.95a44cbc8520ep+0, 0x1.d36906d2b41f9p-53,
+	0x1.97d829fde4e4fp+0, 0x1.173d241f23d18p-53,
+	0x1.9a0f170ca07b9p+0, 0x1.7462137188ce7p-53,
+	0x1.9c49182a3f090p+0, 0x1.c7c46b071f2bep-56,
+	0x1.9e86319e32323p+0, 0x1.824ca78e64c6ep-56,
+	0x1.a0c667b5de564p+0, 0x1.6535b51719567p-53,
+	0x1.a309bec4a2d33p+0, 0x1.6305c7ddc36abp-54,
+	0x1.a5503b23e255cp+0, 0x1.1684892395f0fp-53,
+	0x1.a799e1330b358p+0, 0x1.bcb7ecac563c7p-54,
+	0x1.a9e6b5579fdbfp+0, 0x1.0fac90ef7fd31p-54,
+	0x1.ac36bbfd3f379p+0, 0x1.81b72cd4624ccp-53,
+	0x1.ae89f995ad3adp+0, 0x1.7a1cd345dcc81p-54,
+	0x1.b0e07298db665p+0, 0x1.2108559bf8deep-53,
+	0x1.b33a2b84f15fap+0, 0x1.ed7fa1cf7b290p-53,
+	0x1.b59728de55939p+0, 0x1.1c7102222c90ep-53,
+	0x1.b7f76f2fb5e46p+0, 0x1.d54f610356a79p-53,
+	0x1.ba5b030a10649p+0, 0x1.0819678d5eb69p-53,
+	0x1.bcc1e904bc1d2p+0, 0x1.23dd07a2d9e84p-55,
+	0x1.bf2c25bd71e08p+0, 0x1.0811ae04a31c7p-53,
+	0x1.c199bdd85529cp+0, 0x1.11065895048ddp-55,
+	0x1.c40ab5fffd07ap+0, 0x1.b4537e083c60ap-54,
+	0x1.c67f12e57d14bp+0, 0x1.2884dff483cadp-54,
+	0x1.c8f6d9406e7b5p+0, 0x1.1acbc48805c44p-56,
+	0x1.cb720dcef9069p+0, 0x1.503cbd1e949dbp-56,
+	0x1.cdf0b555dc3f9p+0, 0x1.889f12b1f58a3p-53,
+	0x1.d072d4a07897bp+0, 0x1.1a1e45e4342b2p-53,
+	0x1.d2f87080d89f1p+0, 0x1.15bc247313d44p-53,
+	0x1.d5818dcfba487p+0, 0x1.2ed02d75b3707p-55,
+	0x1.d80e316c98397p+0, 0x1.7709f3a09100cp-53,
+	0x1.da9e603db3285p+0, 0x1.c2300696db532p-54,
+	0x1.dd321f301b460p+0, 0x1.2da5778f018c3p-54,
+	0x1.dfc97337b9b5ep+0, 0x1.72d195873da52p-53,
+	0x1.e264614f5a128p+0, 0x1.424ec3f42f5b5p-53,
+	0x1.e502ee78b3ff6p+0, 0x1.39e8980a9cc8fp-55,
+	0x1.e7a51fbc74c83p+0, 0x1.2d522ca0c8de2p-54,
+	0x1.ea4afa2a490d9p+0, 0x1.0b1ee7431ebb6p-53,
+	0x1.ecf482d8e67f0p+0, 0x1.1b60625f7293ap-53,
+	0x1.efa1bee615a27p+0, 0x1.dc7f486a4b6b0p-54,
+	0x1.f252b376bba97p+0, 0x1.3a1a5bf0d8e43p-54,
+	0x1.f50765b6e4540p+0, 0x1.9d3e12dd8a18bp-54,
+	0x1.f7bfdad9cbe13p+0, 0x1.1227697fce57bp-53,
+	0x1.fa7c1819e90d8p+0, 0x1.74853f3a5931ep-55,
+	0x1.fd3c22b8f71f1p+0, 0x1.2eb74966579e7p-57
+};
+
+/*
+ * Kernel for expl(x).  x must be finite and not tiny or huge.
+ * "tiny" is anything that would make us underflow (|A6*x^6| < ~LDBL_MIN).
+ * "huge" is anything that would make fn*L1 inexact (|x| > ~2**17*ln2).
+ */
+static inline void
+__k_expl(long double x, long double *hip, long double *lop, int *kp)
+{
+	long double fn, q, r, r1, r2, t, z;
+	int n, n2;
+
+	/* 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. */
+#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;
+	/* Depend on the sign bit being propagated: */
+	*kp = n >> LOG2_INTERVALS;
+	r1 = x - fn * L1;
+	r2 = fn * -L2;
+
+	/* Evaluate expl(endpoint[n2] + r1 + r2) = tbl[n2] * expl(r1 + r2). */
+	z = r * r;
+#if 0
+	q = r2 + z * (A2 + r * A3) + z * z * (A4 + r * A5) + z * z * z * A6;
+#else
+	q = r2 + z * A2 + z * r * (A3 + r * A4 + z * (A5 + r * A6));
+#endif
+	t = (long double)tbl[n2].lo + tbl[n2].hi;
+	*hip = tbl[n2].hi;
+	*lop = tbl[n2].lo + t * (q + r1);
+}
+
+static inline void
+k_hexpl(long double x, long double *hip, long double *lop)
+{
+	float twopkm1;
+	int k;
+
+	__k_expl(x, hip, lop, &k);
+	SET_FLOAT_WORD(twopkm1, 0x3f800000 + ((k - 1) << 23));
+	*hip *= twopkm1;
+	*lop *= twopkm1;
+}
+
+static inline long double
+hexpl(long double x)
+{
+	long double hi, lo, twopkm2;
+	int k;
+
+	twopkm2 = 1;
+	__k_expl(x, &hi, &lo, &k);
+	SET_LDBL_EXPSIGN(twopkm2, BIAS + k - 2);
+	return (lo + hi) * 2 * twopkm2;
+}
+
+#ifdef _COMPLEX_H
+/*
+ * See ../src/k_exp.c for details.
+ */
+static inline long double complex
+__ldexp_cexpl(long double complex z, int expt)
+{
+	long double exp_x, hi, lo;
+	long double x, y, scale1, scale2;
+	int half_expt, k;
+
+	x = creall(z);
+	y = cimagl(z);
+	__k_expl(x, &hi, &lo, &k);
+
+	exp_x = (lo + hi) * 0x1p16382;
+	expt += k - 16382;
+
+	scale1 = 1;
+	half_expt = expt / 2;
+	SET_LDBL_EXPSIGN(scale1, BIAS + half_expt);
+	scale2 = 1;
+	SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
+
+	return (cpackl(cos(y) * exp_x * scale1 * scale2,
+	    sinl(y) * exp_x * scale1 * scale2));
+}
+#endif /* _COMPLEX_H */
diff --git a/lib/msun/ld80/s_erfl.c b/lib/msun/ld80/s_erfl.c
new file mode 100644
index 0000000..1ae2f90
--- /dev/null
+++ b/lib/msun/ld80/s_erfl.c
@@ -0,0 +1,337 @@
+/* @(#)s_erf.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See s_erf.c for complete comments.
+ *
+ * Converted to long double by Steven G. Kargl.
+ */
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile long double tiny = 0x1p-10000L;
+
+static const double
+half= 0.5,
+one = 1,
+two = 2;
+/*
+ * In the domain [0, 2**-34], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-102.
+ */
+static const union IEEEl2bits
+efxu  = LD80C(0x8375d410a6db446c, -3,  1.28379167095512573902e-1L),
+efx8u = LD80C(0x8375d410a6db446c,  0,  1.02703333676410059122e+0L),
+/*
+ * Domain [0, 0.84375], range ~[-1.423e-22, 1.423e-22]:
+ * |(erf(x) - x)/x - pp(x)/qq(x)| < 2**-72.573
+ */
+pp0u  = LD80C(0x8375d410a6db446c, -3,   1.28379167095512573902e-1L),
+pp1u  = LD80C(0xa46c7d09ec3d0cec, -2,  -3.21140201054840180596e-1L),
+pp2u  = LD80C(0x9b31e66325576f86, -5,  -3.78893851760347812082e-2L),
+pp3u  = LD80C(0x804ac72c9a0b97dd, -7,  -7.83032847030604679616e-3L),
+pp4u  = LD80C(0x9f42bcbc3d5a601d, -12, -3.03765663857082048459e-4L),
+pp5u  = LD80C(0x9ec4ad6193470693, -16, -1.89266527398167917502e-5L),
+qq1u  = LD80C(0xdb4b8eb713188d6b, -2,   4.28310832832310510579e-1L),
+qq2u  = LD80C(0xa5750835b2459bd1, -4,   8.07896272074540216658e-2L),
+qq3u  = LD80C(0x8b85d6bd6a90b51c, -7,   8.51579638189385354266e-3L),
+qq4u  = LD80C(0x87332f82cff4ff96, -11,  5.15746855583604912827e-4L),
+qq5u  = LD80C(0x83466cb6bf9dca00, -16,  1.56492109706256700009e-5L),
+qq6u  = LD80C(0xf5bf98c2f996bf63, -24,  1.14435527803073879724e-7L);
+#define	efx	(efxu.e)
+#define	efx8	(efx8u.e)
+#define	pp0	(pp0u.e)
+#define	pp1	(pp1u.e)
+#define	pp2	(pp2u.e)
+#define	pp3	(pp3u.e)
+#define	pp4	(pp4u.e)
+#define	pp5	(pp5u.e)
+#define	qq1	(qq1u.e)
+#define	qq2	(qq2u.e)
+#define	qq3	(qq3u.e)
+#define	qq4	(qq4u.e)
+#define	qq5	(qq5u.e)
+#define	qq6	(qq6u.e)
+static const union IEEEl2bits
+erxu  = LD80C(0xd7bb3d0000000000, -1,  8.42700779438018798828e-1L),
+/*
+ * Domain [0.84375, 1.25], range ~[-8.132e-22, 8.113e-22]:
+ * |(erf(x) - erx) - pa(x)/qa(x)| < 2**-71.762
+ */
+pa0u  = LD80C(0xe8211158da02c692, -27,  1.35116960705131296711e-8L),
+pa1u  = LD80C(0xd488f89f36988618, -2,   4.15107507167065612570e-1L),
+pa2u  = LD80C(0xece74f8c63fa3942, -4,  -1.15675565215949226989e-1L),
+pa3u  = LD80C(0xc8d31e020727c006, -4,   9.80589241379624665791e-2L),
+pa4u  = LD80C(0x985d5d5fafb0551f, -5,   3.71984145558422368847e-2L),
+pa5u  = LD80C(0xa5b6c4854d2f5452, -8,  -5.05718799340957673661e-3L),
+pa6u  = LD80C(0x85c8d58fe3993a47, -8,   4.08277919612202243721e-3L),
+pa7u  = LD80C(0xddbfbc23677b35cf, -13,  2.11476292145347530794e-4L),
+qa1u  = LD80C(0xb8a977896f5eff3f, -1,   7.21335860303380361298e-1L),
+qa2u  = LD80C(0x9fcd662c3d4eac86, -1,   6.24227891731886593333e-1L),
+qa3u  = LD80C(0x9d0b618eac67ba07, -2,   3.06727455774491855801e-1L),
+qa4u  = LD80C(0x881a4293f6d6c92d, -3,   1.32912674218195890535e-1L),
+qa5u  = LD80C(0xbab144f07dea45bf, -5,   4.55792134233613027584e-2L),
+qa6u  = LD80C(0xa6c34ba438bdc900, -7,   1.01783980070527682680e-2L),
+qa7u  = LD80C(0x8fa866dc20717a91, -9,   2.19204436518951438183e-3L);
+#define erx	(erxu.e)
+#define pa0	(pa0u.e)
+#define pa1	(pa1u.e)
+#define pa2	(pa2u.e)
+#define pa3	(pa3u.e)
+#define pa4	(pa4u.e)
+#define pa5	(pa5u.e)
+#define pa6	(pa6u.e)
+#define pa7	(pa7u.e)
+#define qa1	(qa1u.e)
+#define qa2	(qa2u.e)
+#define qa3	(qa3u.e)
+#define qa4	(qa4u.e)
+#define qa5	(qa5u.e)
+#define qa6	(qa6u.e)
+#define qa7	(qa7u.e)
+static const union IEEEl2bits
+/*
+ * Domain [1.25,2.85715], range ~[-2.334e-22,2.334e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - ra(x)/sa(x)| < 2**-71.860
+ */
+ra0u  = LD80C(0xa1a091e0fb4f335a, -7, -9.86494298915814308249e-3L),
+ra1u  = LD80C(0xc2b0d045ae37df6b, -1, -7.60510460864878271275e-1L),
+ra2u  = LD80C(0xf2cec3ee7da636c5, 3,  -1.51754798236892278250e+1L),
+ra3u  = LD80C(0x813cc205395adc7d, 7,  -1.29237335516455333420e+2L),
+ra4u  = LD80C(0x8737c8b7b4062c2f, 9,  -5.40871625829510494776e+2L),
+ra5u  = LD80C(0x8ffe5383c08d4943, 10, -1.15194769466026108551e+3L),
+ra6u  = LD80C(0x983573e64d5015a9, 10, -1.21767039790249025544e+3L),
+ra7u  = LD80C(0x92a794e763a6d4db, 9,  -5.86618463370624636688e+2L),
+ra8u  = LD80C(0xd5ad1fae77c3d9a3, 6,  -1.06838132335777049840e+2L),
+ra9u  = LD80C(0x934c1a247807bb9c, 2,  -4.60303980944467334806e+0L),
+sa1u  = LD80C(0xd342f90012bb1189, 4,   2.64077014928547064865e+1L),
+sa2u  = LD80C(0x839be13d9d5da883, 8,   2.63217811300123973067e+2L),
+sa3u  = LD80C(0x9f8cba6d1ae1b24b, 10,  1.27639775710344617587e+3L),
+sa4u  = LD80C(0xcaa83f403713e33e, 11,  3.24251544209971162003e+3L),
+sa5u  = LD80C(0x8796aff2f3c47968, 12,  4.33883591261332837874e+3L),
+sa6u  = LD80C(0xb6ef97f9c753157b, 11,  2.92697460344182158454e+3L),
+sa7u  = LD80C(0xe02aee5f83773d1c, 9,   8.96670799139389559818e+2L),
+sa8u  = LD80C(0xc82b83855b88e07e, 6,   1.00084987800048510018e+2L),
+sa9u  = LD80C(0x92f030aefadf28ad, 1,   2.29591004455459083843e+0L);
+#define ra0	(ra0u.e)
+#define ra1	(ra1u.e)
+#define ra2	(ra2u.e)
+#define ra3	(ra3u.e)
+#define ra4	(ra4u.e)
+#define ra5	(ra5u.e)
+#define ra6	(ra6u.e)
+#define ra7	(ra7u.e)
+#define ra8	(ra8u.e)
+#define ra9	(ra9u.e)
+#define sa1	(sa1u.e)
+#define sa2	(sa2u.e)
+#define sa3	(sa3u.e)
+#define sa4	(sa4u.e)
+#define sa5	(sa5u.e)
+#define sa6	(sa6u.e)
+#define sa7	(sa7u.e)
+#define sa8	(sa8u.e)
+#define sa9	(sa9u.e)
+/*
+ * Domain [2.85715,7], range ~[-8.323e-22,8.390e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rb(x)/sb(x)| < 2**-70.326
+ */
+static const union IEEEl2bits
+rb0u = LD80C(0xa1a091cf43abcd26, -7, -9.86494292470284646962e-3L),
+rb1u = LD80C(0xd19d2df1cbb8da0a, -1, -8.18804618389296662837e-1L),
+rb2u = LD80C(0x9a4dd1383e5daf5b, 4,  -1.92879967111618594779e+1L),
+rb3u = LD80C(0xbff0ae9fc0751de6, 7,  -1.91940164551245394969e+2L),
+rb4u = LD80C(0xdde08465310b472b, 9,  -8.87508080766577324539e+2L),
+rb5u = LD80C(0xe796e1d38c8c70a9, 10, -1.85271506669474503781e+3L),
+rb6u = LD80C(0xbaf655a76e0ab3b5, 10, -1.49569795581333675349e+3L),
+rb7u = LD80C(0x95d21e3e75503c21, 8,  -2.99641547972948019157e+2L),
+sb1u = LD80C(0x814487ed823c8cbd, 5,   3.23169247732868256569e+1L),
+sb2u = LD80C(0xbe4bfbb1301304be, 8,   3.80593618534539961773e+2L),
+sb3u = LD80C(0x809c4ade46b927c7, 11,  2.05776827838541292848e+3L),
+sb4u = LD80C(0xa55284359f3395a8, 12,  5.29031455540062116327e+3L),
+sb5u = LD80C(0xbcfa72da9b820874, 12,  6.04730608102312640462e+3L),
+sb6u = LD80C(0x9d09a35988934631, 11,  2.51260238030767176221e+3L),
+sb7u = LD80C(0xd675bbe542c159fa, 7,   2.14459898308561015684e+2L);
+#define rb0	(rb0u.e)
+#define rb1	(rb1u.e)
+#define rb2	(rb2u.e)
+#define rb3	(rb3u.e)
+#define rb4	(rb4u.e)
+#define rb5	(rb5u.e)
+#define rb6	(rb6u.e)
+#define rb7	(rb7u.e)
+#define sb1	(sb1u.e)
+#define sb2	(sb2u.e)
+#define sb3	(sb3u.e)
+#define sb4	(sb4u.e)
+#define sb5	(sb5u.e)
+#define sb6	(sb6u.e)
+#define sb7	(sb7u.e)
+/*
+ * Domain [7,108], range ~[-4.422e-22,4.422e-22]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rc(x)/sc(x)| < 2**-70.938
+ */
+static const union IEEEl2bits
+/* err = -4.422092275318925082e-22 -70.937689 */
+rc0u = LD80C(0xa1a091cf437a17ad, -7, -9.86494292470008707260e-3L),
+rc1u = LD80C(0xbe79c5a978122b00, -1, -7.44045595049165939261e-1L),
+rc2u = LD80C(0xdb26f9bbe31a2794, 3,  -1.36970155085888424425e+1L),
+rc3u = LD80C(0xb5f69a38f5747ac8, 6,  -9.09816453742625888546e+1L),
+rc4u = LD80C(0xd79676d970d0a21a, 7,  -2.15587750997584074147e+2L),
+rc5u = LD80C(0xfe528153c45ec97c, 6,  -1.27161142938347796666e+2L),
+sc1u = LD80C(0xc5e8cd46d5604a96, 4,   2.47386727842204312937e+1L),
+sc2u = LD80C(0xc5f0f5a5484520eb, 7,   1.97941248254913378865e+2L),
+sc3u = LD80C(0x964e3c7b34db9170, 9,   6.01222441484087787522e+2L),
+sc4u = LD80C(0x99be1b89faa0596a, 9,   6.14970430845978077827e+2L),
+sc5u = LD80C(0xf80dfcbf37ffc5ea, 6,   1.24027318931184605891e+2L);
+#define rc0	(rc0u.e)
+#define rc1	(rc1u.e)
+#define rc2	(rc2u.e)
+#define rc3	(rc3u.e)
+#define rc4	(rc4u.e)
+#define rc5	(rc5u.e)
+#define sc1	(sc1u.e)
+#define sc2	(sc2u.e)
+#define sc3	(sc3u.e)
+#define sc4	(sc4u.e)
+#define sc5	(sc5u.e)
+
+long double
+erfl(long double x)
+{
+	long double ax,R,S,P,Q,s,y,z,r;
+	uint64_t lx;
+	int32_t i;
+	uint16_t hx;
+
+	EXTRACT_LDBL80_WORDS(hx, lx, x);
+
+	if((hx & 0x7fff) == 0x7fff) {	/* erfl(nan)=nan */
+		i = (hx>>15)<<1;
+		return (1-i)+one/x;	/* erfl(+-inf)=+-1 */
+	}
+
+	ENTERI();
+
+	ax = fabsl(x);
+	if(ax < 0.84375) {
+	    if(ax < 0x1p-34L) {
+	        if(ax < 0x1p-16373L)	
+		    RETURNI((8*x+efx8*x)/8);	/* avoid spurious underflow */
+		RETURNI(x + efx*x);
+	    }
+	    z = x*x;
+	    r = pp0+z*(pp1+z*(pp2+z*(pp3+z*(pp4+z*pp5))));
+	    s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*(qq5+z*qq6)))));
+	    y = r/s;
+	    RETURNI(x + x*y);
+	}
+	if(ax < 1.25) {
+	    s = ax-one;
+	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*(pa6+s*pa7))))));
+	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*(qa6+s*qa7))))));
+	    if(x>=0) RETURNI(erx + P/Q); else RETURNI(-erx - P/Q);
+	}
+	if(ax >= 7) {			/* inf>|x|>= 7 */
+	    if(x>=0) RETURNI(one-tiny); else RETURNI(tiny-one);
+	}
+	s = one/(ax*ax);
+	if(ax < 2.85715) {	/* |x| < 2.85715 */
+	    R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*(ra7+
+		s*(ra8+s*ra9))))))));
+	    S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		s*(sa8+s*sa9))))))));
+	} else {	/* |x| >= 2.85715 */
+	    R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*(rb6+s*rb7))))));
+	    S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
+	}
+	z=(float)ax;
+	r=expl(-z*z-0.5625)*expl((z-ax)*(z+ax)+R/S);
+	if(x>=0) RETURNI(one-r/ax); else RETURNI(r/ax-one);
+}
+
+long double
+erfcl(long double x)
+{
+	long double ax,R,S,P,Q,s,y,z,r;
+	uint64_t lx;
+	uint16_t hx;
+
+	EXTRACT_LDBL80_WORDS(hx, lx, x);
+
+	if((hx & 0x7fff) == 0x7fff) {	/* erfcl(nan)=nan */
+					/* erfcl(+-inf)=0,2 */
+	    return ((hx>>15)<<1)+one/x;
+	}
+
+	ENTERI();
+
+	ax = fabsl(x);
+	if(ax < 0.84375L) {
+	    if(ax < 0x1p-34L)
+		RETURNI(one-x);
+	    z = x*x;
+	    r = pp0+z*(pp1+z*(pp2+z*(pp3+z*(pp4+z*pp5))));
+	    s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*(qq5+z*qq6)))));
+	    y = r/s;
+	    if(ax < 0.25L) {  	/* x<1/4 */
+		RETURNI(one-(x+x*y));
+	    } else {
+		r = x*y;
+		r += (x-half);
+	       RETURNI(half - r);
+	    }
+	}
+	if(ax < 1.25L) {
+	    s = ax-one;
+	    P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*(pa6+s*pa7))))));
+	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*(qa6+s*qa7))))));
+	    if(x>=0) {
+	        z  = one-erx; RETURNI(z - P/Q);
+	    } else {
+		z = (erx+P/Q); RETURNI(one+z);
+	    }
+	}
+
+	if(ax < 108) {			/* |x| < 108 */
+ 	    s = one/(ax*ax);
+	    if(ax < 2.85715) {		/* |x| < 2.85715 */
+		R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*(ra7+
+		    s*(ra8+s*ra9))))))));
+		S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		    s*(sa8+s*sa9))))))));
+	    } else if(ax < 7) {		/* | |x| < 7 */
+		R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*(rb6+s*rb7))))));
+		S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
+	    } else {
+		if(x < -7) RETURNI(two-tiny);/* x < -7 */
+		R=rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*rc5))));
+		S=one+s*(sc1+s*(sc2+s*(sc3+s*(sc4+s*sc5))));
+	    }
+	    z = (float)ax;
+	    r = expl(-z*z-0.5625)*expl((z-ax)*(z+ax)+R/S);
+	    if(x>0) RETURNI(r/ax); else RETURNI(two-r/ax);
+	} else {
+	    if(x>0) RETURNI(tiny*tiny); else RETURNI(two-tiny);
+	}
+}
diff --git a/lib/msun/ld80/s_expl.c b/lib/msun/ld80/s_expl.c
index ec748d3..3147d35 100644
--- a/lib/msun/ld80/s_expl.c
+++ b/lib/msun/ld80/s_expl.c
@@ -48,16 +48,15 @@ __FBSDID("$FreeBSD$");
 #include "fpmath.h"
 #include "math.h"
 #include "math_private.h"
+#include "k_expl.h"
 
-#define	INTERVALS	128
-#define	LOG2_INTERVALS	7
-#define	BIAS	(LDBL_MAX_EXP - 1)
+/* XXX Prevent compilers from erroneously constant folding these: */
+static const volatile long double
+huge = 0x1p10000L,
+tiny = 0x1p-10000L;
 
 static const long double
-huge = 0x1p10000L,
 twom10000 = 0x1p-10000L;
-/* XXX Prevent gcc from erroneously constant folding this: */
-static volatile const long double tiny = 0x1p-10000L;
 
 static const union IEEEl2bits
 /* log(2**16384 - 0.5) rounded towards zero: */
@@ -68,178 +67,16 @@ o_thresholdu = LD80C(0xb17217f7d1cf79ab, 13,  11356.5234062941439488L),
 u_thresholdu = LD80C(0xb21dfe7f09e2baa9, 13, -11399.4985314888605581L);
 #define u_threshold	 (u_thresholdu.e)
 
-static const double
-/*
- * ln2/INTERVALS = L1+L2 (hi+lo decomposition for multiplication).  L1 must
- * have at least 22 (= log2(|LDBL_MIN_EXP-extras|) + log2(INTERVALS)) lowest
- * bits zero so that multiplication of it by n is exact.
- */
-INV_L = 1.8466496523378731e+2,		/*  0x171547652b82fe.0p-45 */
-L1 =  5.4152123484527692e-3,		/*  0x162e42ff000000.0p-60 */
-L2 = -3.2819649005320973e-13,		/* -0x1718432a1b0e26.0p-94 */
-/*
- * Domain [-0.002708, 0.002708], range ~[-5.7136e-24, 5.7110e-24]:
- * |exp(x) - p(x)| < 2**-77.2
- * (0.002708 is ln2/(2*INTERVALS) rounded up a little).
- */
-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
- * the first 53 bits of the significand are stored in hi and the next 53
- * bits are in lo.  Tang's paper states that the trailing 6 bits of hi must
- * be zero for his algorithm in both single and double precision, because
- * the table is re-used in the implementation of expm1() where a floating
- * point addition involving hi must be exact.  Here hi is double, so
- * converting it to long double gives 11 trailing zero bits.
- */
-static const struct {
-	double	hi;
-	double	lo;
-} tbl[INTERVALS] = {
-	0x1p+0, 0x0p+0,
-	0x1.0163da9fb3335p+0, 0x1.b61299ab8cdb7p-54,
-	0x1.02c9a3e778060p+0, 0x1.dcdef95949ef4p-53,
-	0x1.04315e86e7f84p+0, 0x1.7ae71f3441b49p-53,
-	0x1.059b0d3158574p+0, 0x1.d73e2a475b465p-55,
-	0x1.0706b29ddf6ddp+0, 0x1.8db880753b0f6p-53,
-	0x1.0874518759bc8p+0, 0x1.186be4bb284ffp-57,
-	0x1.09e3ecac6f383p+0, 0x1.1487818316136p-54,
-	0x1.0b5586cf9890fp+0, 0x1.8a62e4adc610bp-54,
-	0x1.0cc922b7247f7p+0, 0x1.01edc16e24f71p-54,
-	0x1.0e3ec32d3d1a2p+0, 0x1.03a1727c57b53p-59,
-	0x1.0fb66affed31ap+0, 0x1.e464123bb1428p-53,
-	0x1.11301d0125b50p+0, 0x1.49d77e35db263p-53,
-	0x1.12abdc06c31cbp+0, 0x1.f72575a649ad2p-53,
-	0x1.1429aaea92ddfp+0, 0x1.66820328764b1p-53,
-	0x1.15a98c8a58e51p+0, 0x1.2406ab9eeab0ap-55,
-	0x1.172b83c7d517ap+0, 0x1.b9bef918a1d63p-53,
-	0x1.18af9388c8de9p+0, 0x1.777ee1734784ap-53,
-	0x1.1a35beb6fcb75p+0, 0x1.e5b4c7b4968e4p-55,
-	0x1.1bbe084045cd3p+0, 0x1.3563ce56884fcp-53,
-	0x1.1d4873168b9aap+0, 0x1.e016e00a2643cp-54,
-	0x1.1ed5022fcd91cp+0, 0x1.71033fec2243ap-53,
-	0x1.2063b88628cd6p+0, 0x1.dc775814a8495p-55,
-	0x1.21f49917ddc96p+0, 0x1.2a97e9494a5eep-55,
-	0x1.2387a6e756238p+0, 0x1.9b07eb6c70573p-54,
-	0x1.251ce4fb2a63fp+0, 0x1.ac155bef4f4a4p-55,
-	0x1.26b4565e27cddp+0, 0x1.2bd339940e9d9p-55,
-	0x1.284dfe1f56380p+0, 0x1.2d9e2b9e07941p-53,
-	0x1.29e9df51fdee1p+0, 0x1.612e8afad1255p-55,
-	0x1.2b87fd0dad98fp+0, 0x1.fbbd48ca71f95p-53,
-	0x1.2d285a6e4030bp+0, 0x1.0024754db41d5p-54,
-	0x1.2ecafa93e2f56p+0, 0x1.1ca0f45d52383p-56,
-	0x1.306fe0a31b715p+0, 0x1.6f46ad23182e4p-55,
-	0x1.32170fc4cd831p+0, 0x1.a9ce78e18047cp-55,
-	0x1.33c08b26416ffp+0, 0x1.32721843659a6p-54,
-	0x1.356c55f929ff0p+0, 0x1.928c468ec6e76p-53,
-	0x1.371a7373aa9cap+0, 0x1.4e28aa05e8a8fp-53,
-	0x1.38cae6d05d865p+0, 0x1.0b53961b37da2p-53,
-	0x1.3a7db34e59ff6p+0, 0x1.d43792533c144p-53,
-	0x1.3c32dc313a8e4p+0, 0x1.08003e4516b1ep-53,
-	0x1.3dea64c123422p+0, 0x1.ada0911f09ebcp-55,
-	0x1.3fa4504ac801bp+0, 0x1.417ee03548306p-53,
-	0x1.4160a21f72e29p+0, 0x1.f0864b71e7b6cp-53,
-	0x1.431f5d950a896p+0, 0x1.b8e088728219ap-53,
-	0x1.44e086061892dp+0, 0x1.89b7a04ef80d0p-59,
-	0x1.46a41ed1d0057p+0, 0x1.c944bd1648a76p-54,
-	0x1.486a2b5c13cd0p+0, 0x1.3c1a3b69062f0p-56,
-	0x1.4a32af0d7d3dep+0, 0x1.9cb62f3d1be56p-54,
-	0x1.4bfdad5362a27p+0, 0x1.d4397afec42e2p-56,
-	0x1.4dcb299fddd0dp+0, 0x1.8ecdbbc6a7833p-54,
-	0x1.4f9b2769d2ca6p+0, 0x1.5a67b16d3540ep-53,
-	0x1.516daa2cf6641p+0, 0x1.8225ea5909b04p-53,
-	0x1.5342b569d4f81p+0, 0x1.be1507893b0d5p-53,
-	0x1.551a4ca5d920ep+0, 0x1.8a5d8c4048699p-53,
-	0x1.56f4736b527dap+0, 0x1.9bb2c011d93adp-54,
-	0x1.58d12d497c7fdp+0, 0x1.295e15b9a1de8p-55,
-	0x1.5ab07dd485429p+0, 0x1.6324c054647adp-54,
-	0x1.5c9268a5946b7p+0, 0x1.c4b1b816986a2p-60,
-	0x1.5e76f15ad2148p+0, 0x1.ba6f93080e65ep-54,
-	0x1.605e1b976dc08p+0, 0x1.60edeb25490dcp-53,
-	0x1.6247eb03a5584p+0, 0x1.63e1f40dfa5b5p-53,
-	0x1.6434634ccc31fp+0, 0x1.8edf0e2989db3p-53,
-	0x1.6623882552224p+0, 0x1.224fb3c5371e6p-53,
-	0x1.68155d44ca973p+0, 0x1.038ae44f73e65p-57,
-	0x1.6a09e667f3bccp+0, 0x1.21165f626cdd5p-53,
-	0x1.6c012750bdabep+0, 0x1.daed533001e9ep-53,
-	0x1.6dfb23c651a2ep+0, 0x1.e441c597c3775p-53,
-	0x1.6ff7df9519483p+0, 0x1.9f0fc369e7c42p-53,
-	0x1.71f75e8ec5f73p+0, 0x1.ba46e1e5de15ap-53,
-	0x1.73f9a48a58173p+0, 0x1.7ab9349cd1562p-53,
-	0x1.75feb564267c8p+0, 0x1.7edd354674916p-53,
-	0x1.780694fde5d3fp+0, 0x1.866b80a02162dp-54,
-	0x1.7a11473eb0186p+0, 0x1.afaa2047ed9b4p-53,
-	0x1.7c1ed0130c132p+0, 0x1.f124cd1164dd6p-54,
-	0x1.7e2f336cf4e62p+0, 0x1.05d02ba15797ep-56,
-	0x1.80427543e1a11p+0, 0x1.6c1bccec9346bp-53,
-	0x1.82589994cce12p+0, 0x1.159f115f56694p-53,
-	0x1.8471a4623c7acp+0, 0x1.9ca5ed72f8c81p-53,
-	0x1.868d99b4492ecp+0, 0x1.01c83b21584a3p-53,
-	0x1.88ac7d98a6699p+0, 0x1.994c2f37cb53ap-54,
-	0x1.8ace5422aa0dbp+0, 0x1.6e9f156864b27p-54,
-	0x1.8cf3216b5448bp+0, 0x1.de55439a2c38bp-53,
-	0x1.8f1ae99157736p+0, 0x1.5cc13a2e3976cp-55,
-	0x1.9145b0b91ffc5p+0, 0x1.114c368d3ed6ep-53,
-	0x1.93737b0cdc5e4p+0, 0x1.e8a0387e4a814p-53,
-	0x1.95a44cbc8520ep+0, 0x1.d36906d2b41f9p-53,
-	0x1.97d829fde4e4fp+0, 0x1.173d241f23d18p-53,
-	0x1.9a0f170ca07b9p+0, 0x1.7462137188ce7p-53,
-	0x1.9c49182a3f090p+0, 0x1.c7c46b071f2bep-56,
-	0x1.9e86319e32323p+0, 0x1.824ca78e64c6ep-56,
-	0x1.a0c667b5de564p+0, 0x1.6535b51719567p-53,
-	0x1.a309bec4a2d33p+0, 0x1.6305c7ddc36abp-54,
-	0x1.a5503b23e255cp+0, 0x1.1684892395f0fp-53,
-	0x1.a799e1330b358p+0, 0x1.bcb7ecac563c7p-54,
-	0x1.a9e6b5579fdbfp+0, 0x1.0fac90ef7fd31p-54,
-	0x1.ac36bbfd3f379p+0, 0x1.81b72cd4624ccp-53,
-	0x1.ae89f995ad3adp+0, 0x1.7a1cd345dcc81p-54,
-	0x1.b0e07298db665p+0, 0x1.2108559bf8deep-53,
-	0x1.b33a2b84f15fap+0, 0x1.ed7fa1cf7b290p-53,
-	0x1.b59728de55939p+0, 0x1.1c7102222c90ep-53,
-	0x1.b7f76f2fb5e46p+0, 0x1.d54f610356a79p-53,
-	0x1.ba5b030a10649p+0, 0x1.0819678d5eb69p-53,
-	0x1.bcc1e904bc1d2p+0, 0x1.23dd07a2d9e84p-55,
-	0x1.bf2c25bd71e08p+0, 0x1.0811ae04a31c7p-53,
-	0x1.c199bdd85529cp+0, 0x1.11065895048ddp-55,
-	0x1.c40ab5fffd07ap+0, 0x1.b4537e083c60ap-54,
-	0x1.c67f12e57d14bp+0, 0x1.2884dff483cadp-54,
-	0x1.c8f6d9406e7b5p+0, 0x1.1acbc48805c44p-56,
-	0x1.cb720dcef9069p+0, 0x1.503cbd1e949dbp-56,
-	0x1.cdf0b555dc3f9p+0, 0x1.889f12b1f58a3p-53,
-	0x1.d072d4a07897bp+0, 0x1.1a1e45e4342b2p-53,
-	0x1.d2f87080d89f1p+0, 0x1.15bc247313d44p-53,
-	0x1.d5818dcfba487p+0, 0x1.2ed02d75b3707p-55,
-	0x1.d80e316c98397p+0, 0x1.7709f3a09100cp-53,
-	0x1.da9e603db3285p+0, 0x1.c2300696db532p-54,
-	0x1.dd321f301b460p+0, 0x1.2da5778f018c3p-54,
-	0x1.dfc97337b9b5ep+0, 0x1.72d195873da52p-53,
-	0x1.e264614f5a128p+0, 0x1.424ec3f42f5b5p-53,
-	0x1.e502ee78b3ff6p+0, 0x1.39e8980a9cc8fp-55,
-	0x1.e7a51fbc74c83p+0, 0x1.2d522ca0c8de2p-54,
-	0x1.ea4afa2a490d9p+0, 0x1.0b1ee7431ebb6p-53,
-	0x1.ecf482d8e67f0p+0, 0x1.1b60625f7293ap-53,
-	0x1.efa1bee615a27p+0, 0x1.dc7f486a4b6b0p-54,
-	0x1.f252b376bba97p+0, 0x1.3a1a5bf0d8e43p-54,
-	0x1.f50765b6e4540p+0, 0x1.9d3e12dd8a18bp-54,
-	0x1.f7bfdad9cbe13p+0, 0x1.1227697fce57bp-53,
-	0x1.fa7c1819e90d8p+0, 0x1.74853f3a5931ep-55,
-	0x1.fd3c22b8f71f1p+0, 0x1.2eb74966579e7p-57
-};
-
 long double
 expl(long double x)
 {
-	union IEEEl2bits u, v;
-	long double fn, q, r, r1, r2, t, twopk, twopkp10000;
-	long double z;
-	int k, n, n2;
+	union IEEEl2bits u;
+	long double hi, lo, t, twopk;
+	int k;
 	uint16_t hx, ix;
 
+	DOPRINT_START(&x);
+
 	/* Filter out exceptional cases. */
 	u.e = x;
 	hx = u.xbits.expsign;
@@ -247,59 +84,32 @@ expl(long double x)
 	if (ix >= BIAS + 13) {		/* |x| >= 8192 or x is NaN */
 		if (ix == BIAS + LDBL_MAX_EXP) {
 			if (hx & 0x8000)  /* x is -Inf, -NaN or unsupported */
-				return (-1 / x);
- 			return (x + x);	/* x is +Inf, +NaN or unsupported */
+				RETURNP(-1 / x);
+			RETURNP(x + x);	/* x is +Inf, +NaN or unsupported */
 		}
 		if (x > o_threshold)
-			return (huge * huge);
+			RETURNP(huge * huge);
 		if (x < u_threshold)
-			return (tiny * tiny);
-	} else if (ix < BIAS - 65) {	/* |x| < 0x1p-65 (includes pseudos) */
-		return (1 + x);		/* 1 with inexact iff x != 0 */
+			RETURNP(tiny * tiny);
+	} else if (ix < BIAS - 75) {	/* |x| < 0x1p-75 (includes pseudos) */
+		RETURN2P(1, x);		/* 1 with inexact iff x != 0 */
 	}
 
 	ENTERI();
 
-	/* 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. */
-#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;
-	/* Depend on the sign bit being propagated: */
-	k = n >> LOG2_INTERVALS;
-	r1 = x - fn * L1;
-	r2 = fn * -L2;
-
-	/* Prepare scale factors. */
-	v.e = 1;
-	if (k >= LDBL_MIN_EXP) {
-		v.xbits.expsign = BIAS + k;
-		twopk = v.e;
-	} else {
-		v.xbits.expsign = BIAS + k + 10000;
-		twopkp10000 = v.e;
-	}
-
-	/* Evaluate 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;
-	t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+	twopk = 1;
+	__k_expl(x, &hi, &lo, &k);
+	t = SUM2P(hi, lo);
 
 	/* Scale by 2**k. */
 	if (k >= LDBL_MIN_EXP) {
 		if (k == LDBL_MAX_EXP)
 			RETURNI(t * 2 * 0x1p16383L);
+		SET_LDBL_EXPSIGN(twopk, BIAS + k);
 		RETURNI(t * twopk);
 	} else {
-		RETURNI(t * twopkp10000 * twom10000);
+		SET_LDBL_EXPSIGN(twopk, BIAS + k + 10000);
+		RETURNI(t * twopk * twom10000);
 	}
 }
 
@@ -326,8 +136,11 @@ 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
+ * Domain [-0.1659, 0.1659], range ~[-2.6155e-22, 2.5507e-23]:
+ * |(exp(x)-1-x-x**2/2)/x - p(x)| < 2**-71.6
+ *
+ * XXX the coeffs aren't very carefully rounded, and I get 2.8 more bits,
+ * but unlike for ld128 we can't drop any terms.
  */
 static const union IEEEl2bits
 B3 = LD80C(0xaaaaaaaaaaaaaaab, -3,  1.66666666666666666671e-1L),
@@ -353,6 +166,8 @@ expm1l(long double x)
 	int k, n, n2;
 	uint16_t hx, ix;
 
+	DOPRINT_START(&x);
+
 	/* Filter out exceptional cases. */
 	u.e = x;
 	hx = u.xbits.expsign;
@@ -360,11 +175,11 @@ expm1l(long double x)
 	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 */
+				RETURNP(-1 / x - 1);
+			RETURNP(x + x);	/* x is +Inf, +NaN or unsupported */
 		}
 		if (x > o_threshold)
-			return (huge * huge);
+			RETURNP(huge * huge);
 		/*
 		 * expm1l() never underflows, but it must avoid
 		 * unrepresentable large negative exponents.  We used a
@@ -373,15 +188,15 @@ expm1l(long double x)
 		 * in the same way as large ones here.
 		 */
 		if (hx & 0x8000)	/* x <= -64 */
-			return (tiny - 1);	/* good for x < -65ln2 - eps */
+			RETURN2P(tiny, -1);	/* good for x < -65ln2 - eps */
 	}
 
 	ENTERI();
 
 	if (T1 < x && x < T2) {
-		if (ix < BIAS - 64) {	/* |x| < 0x1p-64 (includes pseudos) */
+		if (ix < BIAS - 74) {	/* |x| < 0x1p-74 (includes pseudos) */
 			/* x (rounded) with inexact if x != 0: */
-			RETURNI(x == 0 ? x :
+			RETURNPI(x == 0 ? x :
 			    (0x1p100 * x + fabsl(x)) * 0x1p-100);
 		}
 
@@ -402,9 +217,9 @@ expm1l(long double x)
 		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));
+			RETURN2PI(hx2_hi + x_hi, hx2_lo + x_lo + q);
 		else
-			RETURNI(hx2_lo + q + hx2_hi + x);
+			RETURN2PI(x, hx2_lo + q + hx2_hi);
 	}
 
 	/* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */
@@ -438,21 +253,21 @@ expm1l(long double x)
 	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);
+		t = SUM2P(tbl[n2].hi - 1, tbl[n2].lo * (r1 + 1) + t * q +
+		    tbl[n2].hi * r1);
 		RETURNI(t);
 	}
 	if (k == -1) {
-		t = tbl[n2].lo * (r1 + 1) + t * q + tbl[n2].hi * r1 + 
-		    (tbl[n2].hi - 2);
+		t = SUM2P(tbl[n2].hi - 2, tbl[n2].lo * (r1 + 1) + t * q +
+		    tbl[n2].hi * r1);
 		RETURNI(t / 2);
 	}
 	if (k < -7) {
-		t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+		t = SUM2P(tbl[n2].hi, tbl[n2].lo + t * (q + r1));
 		RETURNI(t * twopk - 1);
 	}
 	if (k > 2 * LDBL_MANT_DIG - 1) {
-		t = tbl[n2].lo + t * (q + r1) + tbl[n2].hi;
+		t = SUM2P(tbl[n2].hi, tbl[n2].lo + t * (q + r1));
 		if (k == LDBL_MAX_EXP)
 			RETURNI(t * 2 * 0x1p16383L - 1);
 		RETURNI(t * twopk - 1);
@@ -462,8 +277,8 @@ expm1l(long double x)
 	twomk = v.e;
 
 	if (k > LDBL_MANT_DIG - 1)
-		t = tbl[n2].lo - twomk + t * (q + r1) + tbl[n2].hi;
+		t = SUM2P(tbl[n2].hi, tbl[n2].lo - twomk + t * (q + r1));
 	else
-		t = tbl[n2].lo + t * (q + r1) + (tbl[n2].hi - twomk);
+		t = SUM2P(tbl[n2].hi - twomk, tbl[n2].lo + t * (q + r1));
 	RETURNI(t * twopk);
 }
diff --git a/lib/msun/man/cosh.3 b/lib/msun/man/cosh.3
index 96abc5a..334b564 100644
--- a/lib/msun/man/cosh.3
+++ b/lib/msun/man/cosh.3
@@ -28,12 +28,13 @@
 .\"     from: @(#)cosh.3	5.1 (Berkeley) 5/2/91
 .\" $FreeBSD$
 .\"
-.Dd January 14, 2005
+.Dd August 17, 2013
 .Dt COSH 3
 .Os
 .Sh NAME
 .Nm cosh ,
-.Nm coshf
+.Nm coshf ,
+.Nm coshl
 .Nd hyperbolic cosine functions
 .Sh LIBRARY
 .Lb libm
@@ -43,11 +44,14 @@
 .Fn cosh "double x"
 .Ft float
 .Fn coshf "float x"
+.Ft long double
+.Fn coshl "long double x"
 .Sh DESCRIPTION
 The
-.Fn cosh
-and the
-.Fn coshf
+.Fn cosh ,
+.Fn coshf ,
+and
+.Fn coshl
 functions compute the hyperbolic cosine of
 .Fa x .
 .Sh SEE ALSO
diff --git a/lib/msun/man/erf.3 b/lib/msun/man/erf.3
index 46886b0..a9a3e0a 100644
--- a/lib/msun/man/erf.3
+++ b/lib/msun/man/erf.3
@@ -28,14 +28,16 @@
 .\"     from: @(#)erf.3	6.4 (Berkeley) 4/20/91
 .\" $FreeBSD$
 .\"
-.Dd April 20, 1991
+.Dd July 13, 2014
 .Dt ERF 3
 .Os
 .Sh NAME
 .Nm erf ,
 .Nm erff ,
+.Nm erfl ,
 .Nm erfc ,
-.Nm erfcf
+.Nm erfcf ,
+.Nm erfcl
 .Nd error function operators
 .Sh LIBRARY
 .Lb libm
@@ -45,18 +47,23 @@
 .Fn erf "double x"
 .Ft float
 .Fn erff "float x"
+.Ft "long double"
+.Fn erfl "long double x"
 .Ft double
 .Fn erfc "double x"
 .Ft float
 .Fn erfcf "float x"
+.Ft "long double"
+.Fn erfcl "long double x"
 .Sh DESCRIPTION
 These functions calculate the error function of
 .Fa x .
 .Pp
 The
-.Fn erf
-and the
-.Fn erff
+.Fn erf ,
+.Fn erff ,
+and
+.Fn erfl
 functions calculate the error function of x; where
 .Bd -ragged -offset indent
 .if n \{\
@@ -69,9 +76,10 @@ erf\|(x) :=
 .Ed
 .Pp
 The
-.Fn erfc
-and the
-.Fn erfcf
+.Fn erfc ,
+.Fn erfcf ,
+and
+.Fn erfcl
 functions calculate the complementary error function of
 .Fa x ;
 that is
@@ -79,9 +87,6 @@ that is
 subtracts the result of the error function
 .Fn erf x
 from 1.0.
-This is useful, since for large
-.Fa x
-places disappear.
 .Sh SEE ALSO
 .Xr math 3
 .Sh HISTORY
diff --git a/lib/msun/man/sinh.3 b/lib/msun/man/sinh.3
index 02944cc..b34cc38 100644
--- a/lib/msun/man/sinh.3
+++ b/lib/msun/man/sinh.3
@@ -27,12 +27,14 @@
 .\"
 .\"	from: @(#)sinh.3	6.6 (Berkeley) 4/19/91
 .\" $FreeBSD$
-.Dd January 14, 2005
+.\"
+.Dd August 17, 2013
 .Dt SINH 3
 .Os
 .Sh NAME
 .Nm sinh ,
-.Nm sinhf
+.Nm sinhf ,
+.Nm sinhl
 .Nd hyperbolic sine function
 .Sh LIBRARY
 .Lb libm
@@ -42,11 +44,14 @@
 .Fn sinh "double x"
 .Ft float
 .Fn sinhf "float x"
+.Ft long double
+.Fn sinhl "long double x"
 .Sh DESCRIPTION
 The
-.Fn sinh
-and the
-.Fn sinhf
+.Fn sinh ,
+.Fn sinhf ,
+and
+.Fn sinhl
 functions compute the hyperbolic sine of
 .Fa x .
 .Sh SEE ALSO
diff --git a/lib/msun/man/tanh.3 b/lib/msun/man/tanh.3
index 6fb185c..ea2468f 100644
--- a/lib/msun/man/tanh.3
+++ b/lib/msun/man/tanh.3
@@ -28,12 +28,13 @@
 .\"     from: @(#)tanh.3	5.1 (Berkeley) 5/2/91
 .\" $FreeBSD$
 .\"
-.Dd May 2, 1991
+.Dd August 17, 2013
 .Dt TANH 3
 .Os
 .Sh NAME
 .Nm tanh ,
-.Nm tanhf
+.Nm tanhf ,
+.Nm tanhl
 .Nd hyperbolic tangent functions
 .Sh LIBRARY
 .Lb libm
@@ -43,20 +44,24 @@
 .Fn tanh "double x"
 .Ft float
 .Fn tanhf "float x"
+.Ft long double
+.Fn tanhl "long double x"
 .Sh DESCRIPTION
 The
-.Fn tanh
-and the
-.Fn tanhf
+.Fn tanh ,
+.Fn tanhf ,
+and
+.Fn tanhl
 functions compute the hyperbolic tangent of
 .Fa x .
 For a discussion of error due to roundoff, see
 .Xr math 3 .
 .Sh RETURN VALUES
 The
-.Fn tanh
+.Fn tanh ,
+.Fn tanhf ,
 and the
-.Fn tanhf
+.Fn tanhl
 functions return the hyperbolic tangent value.
 .Sh SEE ALSO
 .Xr acos 3 ,
diff --git a/lib/msun/src/e_cosh.c b/lib/msun/src/e_cosh.c
index a363695..246b5fb 100644
--- a/lib/msun/src/e_cosh.c
+++ b/lib/msun/src/e_cosh.c
@@ -35,6 +35,8 @@ __FBSDID("$FreeBSD$");
  *	only cosh(0)=1 is exact for finite x.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -77,3 +79,7 @@ __ieee754_cosh(double x)
     /* |x| > overflowthresold, cosh(x) overflow */
 	return huge*huge;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(cosh, coshl);
+#endif
diff --git a/lib/msun/src/e_coshl.c b/lib/msun/src/e_coshl.c
new file mode 100644
index 0000000..0a21277
--- /dev/null
+++ b/lib/msun/src/e_coshl.c
@@ -0,0 +1,130 @@
+/* from: FreeBSD: head/lib/msun/src/e_coshl.c XXX */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_cosh.c for complete comments.
+ *
+ * Converted to long double by Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+#include "k_expl.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+static const volatile long double huge = 0x1p10000L, tiny = 0x1p-10000L;
+#if LDBL_MANT_DIG == 64
+/*
+ * Domain [-1, 1], range ~[-1.8211e-21, 1.8211e-21]:
+ * |cosh(x) - c(x)| < 2**-68.8
+ */
+static const union IEEEl2bits
+C4u = LD80C(0xaaaaaaaaaaaaac78, -5,  4.16666666666666682297e-2L);
+#define	C4	C4u.e
+static const double
+C2  =  0.5,
+C6  =  1.3888888888888616e-3,		/*  0x16c16c16c16b99.0p-62 */
+C8  =  2.4801587301767953e-5,		/*  0x1a01a01a027061.0p-68 */
+C10 =  2.7557319163300398e-7,		/*  0x127e4fb6c9b55f.0p-74 */
+C12 =  2.0876768371393075e-9,		/*  0x11eed99406a3f4.0p-81 */
+C14 =  1.1469537039374480e-11,		/*  0x1938c67cd18c48.0p-89 */
+C16 =  4.8473490896852041e-14;		/*  0x1b49c429701e45.0p-97 */
+#elif LDBL_MANT_DIG == 113
+/*
+ * Domain [-1, 1], range ~[-2.3194e-37, 2.3194e-37]:
+ * |cosh(x) - c(x)| < 2**-121.69
+ */
+static const long double
+C4  =  4.16666666666666666666666666666666225e-2L,	/*  0x1555555555555555555555555554e.0p-117L */
+C6  =  1.38888888888888888888888888889434831e-3L,	/*  0x16c16c16c16c16c16c16c16c1dd7a.0p-122L */
+C8  =  2.48015873015873015873015871870962089e-5L,	/*  0x1a01a01a01a01a01a01a017af2756.0p-128L */
+C10 =  2.75573192239858906525574318600800201e-7L,	/*  0x127e4fb7789f5c72ef01c8a040640.0p-134L */
+C12 =  2.08767569878680989791444691755468269e-9L,	/*  0x11eed8eff8d897b543d0679607399.0p-141L */
+C14=  1.14707455977297247387801189650495351e-11L,	/*  0x193974a8c07c9d24ae169a7fa9b54.0p-149L */
+C16 =  4.77947733238737883626416876486279985e-14L;	/*  0x1ae7f3e733b814d4e1b90f5727fe4.0p-157L */
+static const double
+C2  =  0.5,
+C18 =  1.5619206968597871e-16,		/*  0x16827863b9900b.0p-105 */
+C20 =  4.1103176218528049e-19,		/*  0x1e542ba3d3c269.0p-114 */
+C22 =  8.8967926401641701e-22,		/*  0x10ce399542a014.0p-122 */
+C24 =  1.6116681626523904e-24,		/*  0x1f2c981d1f0cb7.0p-132 */
+C26 =  2.5022374732804632e-27;		/*  0x18c7ecf8b2c4a0.0p-141 */
+#else
+#error "Unsupported long double format"
+#endif /* LDBL_MANT_DIG == 64 */
+
+/* log(2**16385 - 0.5) rounded up: */
+static const float
+o_threshold =  1.13572168e4;		/*  0xb174de.0p-10 */
+
+long double
+coshl(long double x)
+{
+	long double hi,lo,x2,x4;
+	double dx2;
+	uint16_t ix;
+
+	GET_LDBL_EXPSIGN(ix,x);
+	ix &= 0x7fff;
+
+    /* x is INF or NaN */
+	if(ix>=0x7fff) return x*x;
+
+	ENTERI();
+
+    /* |x| < 1, return 1 or c(x) */
+	if(ix<0x3fff) {
+	    if (ix<BIAS-(LDBL_MANT_DIG+1)/2) 	/* |x| < TINY */
+		RETURNI(1+tiny);	/* cosh(tiny) = 1(+) with inexact */
+	    x2 = x*x;
+#if LDBL_MANT_DIG == 64
+	    x4 = x2*x2;
+	    RETURNI(((C16*x2 + C14)*x4 + (C12*x2 + C10))*(x4*x4*x2) +
+		((C8*x2 + C6)*x2 + C4)*x4 + C2*x2 + 1);
+#elif LDBL_MANT_DIG == 113
+	    dx2 = x2;
+	    RETURNI((((((((((((C26*dx2 + C24)*dx2 + C22)*dx2 +
+		C20)*x2 + C18)*x2 +
+		C16)*x2 + C14)*x2 + C12)*x2 + C10)*x2 + C8)*x2 + C6)*x2 +
+		C4)*(x2*x2) + C2*x2 + 1);
+#endif
+	}
+
+    /* |x| in [1, 64), return accurate exp(|x|)/2+1/exp(|x|)/2 */
+	if (ix < 0x4005) {
+	    k_hexpl(fabsl(x), &hi, &lo);
+	    RETURNI(lo + 0.25/(hi + lo) + hi);
+	}
+
+    /* |x| in [64, o_threshold], return correctly-overflowing exp(|x|)/2 */
+	if (fabsl(x) <= o_threshold)
+	    RETURNI(hexpl(fabsl(x)));
+
+    /* |x| > o_threshold, cosh(x) overflow */
+	RETURNI(huge*huge);
+}
diff --git a/lib/msun/src/e_lgamma_r.c b/lib/msun/src/e_lgamma_r.c
index 1cff592..7a95ea4 100644
--- a/lib/msun/src/e_lgamma_r.c
+++ b/lib/msun/src/e_lgamma_r.c
@@ -86,8 +86,10 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+static const volatile double vzero = 0;
+
 static const double
-two52=  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
+zero=  0.00000000000000000000e+00,
 half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 pi  =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
@@ -154,39 +156,35 @@ w4  = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
 w5  =  8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
 w6  = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
 
-static const double zero=  0.00000000000000000000e+00;
-
-	static double sin_pi(double x)
+/*
+ * Compute sin(pi*x) without actually doing the pi*x multiplication.
+ * sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is
+ * the precision of x.
+ */
+static double
+sin_pi(double x)
 {
+	volatile double vz;
 	double y,z;
-	int n,ix;
+	int n;
+
+	y = -x;
 
-	GET_HIGH_WORD(ix,x);
-	ix &= 0x7fffffff;
+	vz = y+0x1p52;			/* depend on 0 <= y < 0x1p52 */
+	z = vz-0x1p52;			/* rint(y) for the above range */
+	if (z == y)
+	    return zero;
 
-	if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0);
-	y = -x;		/* x is assume negative */
+	vz = y+0x1p50;
+	GET_LOW_WORD(n,vz);		/* bits for rounded y (units 0.25) */
+	z = vz-0x1p50;			/* y rounded to a multiple of 0.25 */
+	if (z > y) {
+	    z -= 0.25;			/* adjust to round down */
+	    n--;
+	}
+	n &= 7;				/* octant of y mod 2 */
+	y = y - z + n * 0.25;		/* y mod 2 */
 
-    /*
-     * argument reduction, make sure inexact flag not raised if input
-     * is an integer
-     */
-	z = floor(y);
-	if(z!=y) {				/* inexact anyway */
-	    y  *= 0.5;
-	    y   = 2.0*(y - floor(y));		/* y = |x| mod 2.0 */
-	    n   = (int) (y*4.0);
-	} else {
-            if(ix>=0x43400000) {
-                y = zero; n = 0;                 /* y must be even */
-            } else {
-                if(ix<0x43300000) z = y+two52;	/* exact */
-		GET_LOW_WORD(n,z);
-		n &= 1;
-                y  = n;
-                n<<= 2;
-            }
-        }
 	switch (n) {
 	    case 0:   y =  __kernel_sin(pi*y,zero,0); break;
 	    case 1:   
@@ -206,7 +204,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
 {
 	double t,y,z,nadj,p,p1,p2,p3,q,r,w;
 	int32_t hx;
-	int i,lx,ix;
+	int i,ix,lx;
 
 	EXTRACT_WORDS(hx,lx,x);
 
@@ -214,7 +212,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
 	if(ix>=0x7ff00000) return x*x;
-	if((ix|lx)==0) return one/zero;
+	if((ix|lx)==0) return one/vzero;
 	if(ix<0x3b900000) {	/* |x|<2**-70, return -log(|x|) */
 	    if(hx<0) {
 	        *signgamp = -1;
@@ -223,9 +221,9 @@ __ieee754_lgamma_r(double x, int *signgamp)
 	}
 	if(hx<0) {
 	    if(ix>=0x43300000) 	/* |x|>=2**52, must be -integer */
-		return one/zero;
+		return one/vzero;
 	    t = sin_pi(x);
-	    if(t==zero) return one/zero; /* -integer */
+	    if(t==zero) return one/vzero; /* -integer */
 	    nadj = __ieee754_log(pi/fabs(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
diff --git a/lib/msun/src/e_lgammaf_r.c b/lib/msun/src/e_lgammaf_r.c
index e2d90ef..9a7ab39 100644
--- a/lib/msun/src/e_lgammaf_r.c
+++ b/lib/msun/src/e_lgammaf_r.c
@@ -19,8 +19,10 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+static const volatile float vzero = 0;
+
 static const float
-two23=  8.3886080000e+06, /* 0x4b000000 */
+zero=  0.0000000000e+00,
 half=  5.0000000000e-01, /* 0x3f000000 */
 one =  1.0000000000e+00, /* 0x3f800000 */
 pi  =  3.1415927410e+00, /* 0x40490fdb */
@@ -87,39 +89,30 @@ w4  = -5.9518753551e-04, /* 0xba1c065c */
 w5  =  8.3633989561e-04, /* 0x3a5b3dd2 */
 w6  = -1.6309292987e-03; /* 0xbad5c4e8 */
 
-static const float zero=  0.0000000000e+00;
-
-	static float sin_pif(float x)
+static float
+sin_pif(float x)
 {
+	volatile float vz;
 	float y,z;
-	int n,ix;
+	int n;
 
-	GET_FLOAT_WORD(ix,x);
-	ix &= 0x7fffffff;
+	y = -x;
 
-	if(ix<0x3e800000) return __kernel_sindf(pi*x);
-	y = -x;		/* x is assume negative */
+	vz = y+0x1p23F;			/* depend on 0 <= y < 0x1p23 */
+	z = vz-0x1p23F;			/* rintf(y) for the above range */
+	if (z == y)
+	    return zero;
+
+	vz = y+0x1p21F;
+	GET_FLOAT_WORD(n,vz);		/* bits for rounded y (units 0.25) */
+	z = vz-0x1p21F;			/* y rounded to a multiple of 0.25 */
+	if (z > y) {
+	    z -= 0.25F;			/* adjust to round down */
+	    n--;
+	}
+	n &= 7;				/* octant of y mod 2 */
+	y = y - z + n * 0.25F;		/* y mod 2 */
 
-    /*
-     * argument reduction, make sure inexact flag not raised if input
-     * is an integer
-     */
-	z = floorf(y);
-	if(z!=y) {				/* inexact anyway */
-	    y  *= (float)0.5;
-	    y   = (float)2.0*(y - floorf(y));	/* y = |x| mod 2.0 */
-	    n   = (int) (y*(float)4.0);
-	} else {
-            if(ix>=0x4b800000) {
-                y = zero; n = 0;                 /* y must be even */
-            } else {
-                if(ix<0x4b000000) z = y+two23;	/* exact */
-		GET_FLOAT_WORD(n,z);
-		n &= 1;
-                y  = n;
-                n<<= 2;
-            }
-        }
 	switch (n) {
 	    case 0:   y =  __kernel_sindf(pi*y); break;
 	    case 1:
@@ -147,7 +140,7 @@ __ieee754_lgammaf_r(float x, int *signgamp)
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
 	if(ix>=0x7f800000) return x*x;
-	if(ix==0) return one/zero;
+	if(ix==0) return one/vzero;
 	if(ix<0x35000000) {	/* |x|<2**-21, return -log(|x|) */
 	    if(hx<0) {
 	        *signgamp = -1;
@@ -156,9 +149,9 @@ __ieee754_lgammaf_r(float x, int *signgamp)
 	}
 	if(hx<0) {
 	    if(ix>=0x4b000000) 	/* |x|>=2**23, must be -integer */
-		return one/zero;
+		return one/vzero;
 	    t = sin_pif(x);
-	    if(t==zero) return one/zero; /* -integer */
+	    if(t==zero) return one/vzero; /* -integer */
 	    nadj = __ieee754_logf(pi/fabsf(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
diff --git a/lib/msun/src/e_pow.c b/lib/msun/src/e_pow.c
index 7607a4a..d54af9d 100644
--- a/lib/msun/src/e_pow.c
+++ b/lib/msun/src/e_pow.c
@@ -19,20 +19,20 @@ __FBSDID("$FreeBSD$");
  *	1. Compute and return log2(x) in two pieces:
  *		log2(x) = w1 + w2,
  *	   where w1 has 53-24 = 29 bit trailing zeros.
- *	2. Perform y*log2(x) = n+y' by simulating muti-precision 
+ *	2. Perform y*log2(x) = n+y' by simulating multi-precision 
  *	   arithmetic, where |y'|<=0.5.
  *	3. Return x**y = 2**n*exp(y'*log2)
  *
  * Special cases:
  *	1.  (anything) ** 0  is 1
  *	2.  (anything) ** 1  is itself
- *	3.  (anything) ** NAN is NAN
+ *	3.  (anything) ** NAN is NAN except 1 ** NAN = 1
  *	4.  NAN ** (anything except 0) is NAN
  *	5.  +-(|x| > 1) **  +INF is +INF
  *	6.  +-(|x| > 1) **  -INF is +0
  *	7.  +-(|x| < 1) **  +INF is +0
  *	8.  +-(|x| < 1) **  -INF is +INF
- *	9.  +-1         ** +-INF is NAN
+ *	9.  +-1         ** +-INF is 1
  *	10. +0 ** (+anything except 0, NAN)               is +0
  *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
  *	12. +0 ** (-anything except 0, NAN)               is +INF
@@ -141,7 +141,7 @@ __ieee754_pow(double x, double y)
 	if(ly==0) { 	
 	    if (iy==0x7ff00000) {	/* y is +-inf */
 	        if(((ix-0x3ff00000)|lx)==0)
-		    return  one;	/* (-1)**+-inf is NaN */
+		    return  one;	/* (-1)**+-inf is 1 */
 	        else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
 		    return (hy>=0)? y: zero;
 	        else			/* (|x|<1)**-,+inf = inf,0 */
diff --git a/lib/msun/src/e_sinh.c b/lib/msun/src/e_sinh.c
index 17442d0..6c01f4a 100644
--- a/lib/msun/src/e_sinh.c
+++ b/lib/msun/src/e_sinh.c
@@ -32,6 +32,8 @@ __FBSDID("$FreeBSD$");
  *	only sinh(0)=0 is exact for finite x.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -71,3 +73,7 @@ __ieee754_sinh(double x)
     /* |x| > overflowthresold, sinh(x) overflow */
 	return x*shuge;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(sinh, sinhl);
+#endif
diff --git a/lib/msun/src/e_sinhl.c b/lib/msun/src/e_sinhl.c
new file mode 100644
index 0000000..ce7e333
--- /dev/null
+++ b/lib/msun/src/e_sinhl.c
@@ -0,0 +1,131 @@
+/* from: FreeBSD: head/lib/msun/src/e_sinhl.c XXX */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_sinh.c for complete comments.
+ *
+ * Converted to long double by Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+#include "k_expl.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+static const long double shuge = 0x1p16383L;
+#if LDBL_MANT_DIG == 64
+/*
+ * Domain [-1, 1], range ~[-6.6749e-22, 6.6749e-22]:
+ * |sinh(x)/x - s(x)| < 2**-70.3
+ */
+static const union IEEEl2bits
+S3u = LD80C(0xaaaaaaaaaaaaaaaa, -3,  1.66666666666666666658e-1L);
+#define	S3	S3u.e
+static const double
+S5  =  8.3333333333333332e-3,		/*  0x11111111111111.0p-59 */
+S7  =  1.9841269841270074e-4,		/*  0x1a01a01a01a070.0p-65 */
+S9  =  2.7557319223873889e-6,		/*  0x171de3a5565fe6.0p-71 */
+S11 =  2.5052108406704084e-8,		/*  0x1ae6456857530f.0p-78 */
+S13 =  1.6059042748655297e-10,		/*  0x161245fa910697.0p-85 */
+S15 =  7.6470006914396920e-13,		/*  0x1ae7ce4eff2792.0p-93 */
+S17 =  2.8346142308424267e-15;		/*  0x19882ce789ffc6.0p-101 */
+#elif LDBL_MANT_DIG == 113
+/*
+ * Domain [-1, 1], range ~[-2.9673e-36, 2.9673e-36]:
+ * |sinh(x)/x - s(x)| < 2**-118.0
+ */
+static const long double
+S3  =  1.66666666666666666666666666666666033e-1L,	/*  0x1555555555555555555555555553b.0p-115L */
+S5  =  8.33333333333333333333333333337643193e-3L,	/*  0x111111111111111111111111180f5.0p-119L */
+S7  =  1.98412698412698412698412697391263199e-4L,	/*  0x1a01a01a01a01a01a01a0176aad11.0p-125L */
+S9  =  2.75573192239858906525574406205464218e-6L,	/*  0x171de3a556c7338faac243aaa9592.0p-131L */
+S11 =  2.50521083854417187749675637460977997e-8L,	/*  0x1ae64567f544e38fe59b3380d7413.0p-138L */
+S13 =  1.60590438368216146368737762431552702e-10L,	/*  0x16124613a86d098059c7620850fc2.0p-145L */
+S15 =  7.64716373181980539786802470969096440e-13L,	/*  0x1ae7f3e733b814193af09ce723043.0p-153L */
+S17 =  2.81145725434775409870584280722701574e-15L;	/*  0x1952c77030c36898c3fd0b6dfc562.0p-161L */
+static const double
+S19=  8.2206352435411005e-18,		/*  0x12f49b4662b86d.0p-109 */
+S21=  1.9572943931418891e-20,		/*  0x171b8f2fab9628.0p-118 */
+S23 =  3.8679983530666939e-23,		/*  0x17617002b73afc.0p-127 */
+S25 =  6.5067867911512749e-26;		/*  0x1423352626048a.0p-136 */
+#else
+#error "Unsupported long double format"
+#endif /* LDBL_MANT_DIG == 64 */
+
+/* log(2**16385 - 0.5) rounded up: */
+static const float
+o_threshold =  1.13572168e4;		/*  0xb174de.0p-10 */
+
+long double
+sinhl(long double x)
+{
+	long double hi,lo,x2,x4;
+	double dx2,s;
+	int16_t ix,jx;
+
+	GET_LDBL_EXPSIGN(jx,x);
+	ix = jx&0x7fff;
+
+    /* x is INF or NaN */
+	if(ix>=0x7fff) return x+x;
+
+	ENTERI();
+
+	s = 1;
+	if (jx<0) s = -1;
+
+    /* |x| < 64, return x, s(x), or accurate s*(exp(|x|)/2-1/exp(|x|)/2) */
+	if (ix<0x4005) {		/* |x|<64 */
+	    if (ix<BIAS-(LDBL_MANT_DIG+1)/2) 	/* |x|<TINY */
+		if(shuge+x>1) RETURNI(x);  /* sinh(tiny) = tiny with inexact */
+	    if (ix<0x3fff) {		/* |x|<1 */
+		x2 = x*x;
+#if LDBL_MANT_DIG == 64
+		x4 = x2*x2;
+		RETURNI(((S17*x2 + S15)*x4 + (S13*x2 + S11))*(x2*x*x4*x4) +
+		    ((S9*x2 + S7)*x2 + S5)*(x2*x*x2) + S3*(x2*x) + x);
+#elif LDBL_MANT_DIG == 113
+		dx2 = x2;
+		RETURNI(((((((((((S25*dx2 + S23)*dx2 +
+		    S21)*x2 + S19)*x2 +
+		    S17)*x2 + S15)*x2 + S13)*x2 + S11)*x2 + S9)*x2 + S7)*x2 +
+		    S5)* (x2*x*x2) +
+		    S3*(x2*x) + x);
+#endif
+	    }
+	    k_hexpl(fabsl(x), &hi, &lo);
+	    RETURNI(s*(lo - 0.25/(hi + lo) + hi));
+	}
+
+    /* |x| in [64, o_threshold], return correctly-overflowing s*exp(|x|)/2 */
+	if (fabsl(x) <= o_threshold)
+	    RETURNI(s*hexpl(fabsl(x)));
+
+    /* |x| > o_threshold, sinh(x) overflow */
+	return x*shuge;
+}
diff --git a/lib/msun/src/imprecise.c b/lib/msun/src/imprecise.c
index a7503bf..92fb2d0 100644
--- a/lib/msun/src/imprecise.c
+++ b/lib/msun/src/imprecise.c
@@ -60,10 +60,5 @@ DECLARE_WEAK(powl);
 	long double imprecise_ ## f ## l(long double v) { return f(v); }\
 	DECLARE_WEAK(f ## l)
 
-DECLARE_IMPRECISE(cosh);
-DECLARE_IMPRECISE(erfc);
-DECLARE_IMPRECISE(erf);
 DECLARE_IMPRECISE(lgamma);
-DECLARE_IMPRECISE(sinh);
-DECLARE_IMPRECISE(tanh);
 DECLARE_IMPRECISE(tgamma);
diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h
index 1bd931c..3ab76f8 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -451,7 +451,10 @@ long double	atanl(long double);
 long double	cbrtl(long double);
 long double	ceill(long double);
 long double	copysignl(long double, long double) __pure2;
+long double	coshl(long double);
 long double	cosl(long double);
+long double	erfcl(long double);
+long double	erfl(long double);
 long double	exp2l(long double);
 long double	expl(long double);
 long double	expm1l(long double);
@@ -466,6 +469,7 @@ long double	frexpl(long double value, int *); /* fundamentally !__pure2 */
 long double	hypotl(long double, long double);
 int		ilogbl(long double) __pure2;
 long double	ldexpl(long double, int);
+long double	lgammal(long double);
 long long	llrintl(long double);
 long long	llroundl(long double);
 long double	log10l(long double);
@@ -482,45 +486,22 @@ long double	nextafterl(long double, long double);
 double		nexttoward(double, long double);
 float		nexttowardf(float, long double);
 long double	nexttowardl(long double, long double);
+long double	powl(long double, long double);
 long double	remainderl(long double, long double);
 long double	remquol(long double, long double, int *);
 long double	rintl(long double);
 long double	roundl(long double);
 long double	scalblnl(long double, long);
 long double	scalbnl(long double, int);
+long double	sinhl(long double);
 long double	sinl(long double);
 long double	sqrtl(long double);
+long double	tanhl(long double);
 long double	tanl(long double);
+long double	tgammal(long double);
 long double	truncl(long double);
 
 #endif /* __ISO_C_VISIBLE >= 1999 */
 __END_DECLS
 
 #endif /* !_MATH_H_ */
-
-/* separate header for cmath */
-#ifndef _MATH_EXTRA_H_
-#if __ISO_C_VISIBLE >= 1999
-#if _DECLARE_C99_LDBL_MATH
-
-#define _MATH_EXTRA_H_
-
-/*
- * extra long double versions of math functions for C99 and cmath
- */
-__BEGIN_DECLS
-
-long double	coshl(long double);
-long double	erfcl(long double);
-long double	erfl(long double);
-long double	lgammal(long double);
-long double	powl(long double, long double);
-long double	sinhl(long double);
-long double	tanhl(long double);
-long double	tgammal(long double);
-
-__END_DECLS
-
-#endif /* !_DECLARE_C99_LDBL_MATH */
-#endif /* __ISO_C_VISIBLE >= 1999 */
-#endif /* !_MATH_EXTRA_H_ */
diff --git a/lib/msun/src/s_erf.c b/lib/msun/src/s_erf.c
index 854767b..e1d63bc 100644
--- a/lib/msun/src/s_erf.c
+++ b/lib/msun/src/s_erf.c
@@ -111,18 +111,25 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile double tiny= 1e-300;
+
 static const double
-tiny	    = 1e-300,
-half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
-one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
-two =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
-	/* c = (float)0.84506291151 */
+half= 0.5,
+one = 1,
+two = 2,
+/* c = (float)0.84506291151 */
 erx =  8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
 /*
- * Coefficients for approximation to  erf on [0,0.84375]
+ * In the domain [0, 2**-28], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-84.
  */
 efx =  1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
 efx8=  1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
+/*
+ * Coefficients for approximation to erf on [0,0.84375]
+ */
 pp0  =  1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
 pp1  = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
 pp2  = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
@@ -134,7 +141,7 @@ qq3  =  5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
 qq4  =  1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
 qq5  = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
 /*
- * Coefficients for approximation to  erf  in [0.84375,1.25]
+ * Coefficients for approximation to erf in [0.84375,1.25]
  */
 pa0  = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
 pa1  =  4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
@@ -150,7 +157,7 @@ qa4  =  1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
 qa5  =  1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
 qa6  =  1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
 /*
- * Coefficients for approximation to  erfc in [1.25,1/0.35]
+ * Coefficients for approximation to erfc in [1.25,1/0.35]
  */
 ra0  = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
 ra1  = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
@@ -169,7 +176,7 @@ sa6  =  1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
 sa7  =  6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
 sa8  = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
 /*
- * Coefficients for approximation to  erfc in [1/.35,28]
+ * Coefficients for approximation to erfc in [1/.35,28]
  */
 rb0  = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
 rb1  = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
@@ -222,15 +229,12 @@ erf(double x)
 	x = fabs(x);
  	s = one/(x*x);
 	if(ix< 0x4006DB6E) {	/* |x| < 1/0.35 */
-	    R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
-				ra5+s*(ra6+s*ra7))))));
-	    S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
-				sa5+s*(sa6+s*(sa7+s*sa8)))))));
+	    R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
+	    S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		s*sa8)))))));
 	} else {	/* |x| >= 1/0.35 */
-	    R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
-				rb5+s*rb6)))));
-	    S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
-				sb5+s*(sb6+s*sb7))))));
+	    R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
+	    S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
 	}
 	z  = x;
 	SET_LOW_WORD(z,0);
@@ -238,6 +242,10 @@ erf(double x)
 	if(hx>=0) return one-r/x; else return  r/x-one;
 }
 
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(erf, erfl);
+#endif
+
 double
 erfc(double x)
 {
@@ -279,23 +287,23 @@ erfc(double x)
 	    x = fabs(x);
  	    s = one/(x*x);
 	    if(ix< 0x4006DB6D) {	/* |x| < 1/.35 ~ 2.857143*/
-	        R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
-				ra5+s*(ra6+s*ra7))))));
-	        S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
-				sa5+s*(sa6+s*(sa7+s*sa8)))))));
+		R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))));
+		S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+
+		    s*sa8)))))));
 	    } else {			/* |x| >= 1/.35 ~ 2.857143 */
 		if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
-	        R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
-				rb5+s*rb6)))));
-	        S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
-				sb5+s*(sb6+s*sb7))))));
+		R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))));
+		S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))));
 	    }
 	    z  = x;
 	    SET_LOW_WORD(z,0);
-	    r  =  __ieee754_exp(-z*z-0.5625)*
-			__ieee754_exp((z-x)*(z+x)+R/S);
+	    r  =  __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
 	    if(hx>0) return r/x; else return two-r/x;
 	} else {
 	    if(hx>0) return tiny*tiny; else return two-tiny;
 	}
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(erfc, erfcl);
+#endif
diff --git a/lib/msun/src/s_erff.c b/lib/msun/src/s_erff.c
index b97ca1d..d6cfbd2 100644
--- a/lib/msun/src/s_erff.c
+++ b/lib/msun/src/s_erff.c
@@ -19,64 +19,63 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile float tiny = 1e-30;
+
 static const float
-tiny	    = 1e-30,
-half=  5.0000000000e-01, /* 0x3F000000 */
-one =  1.0000000000e+00, /* 0x3F800000 */
-two =  2.0000000000e+00, /* 0x40000000 */
+half= 0.5,
+one = 1,
+two = 2,
+erx = 8.42697144e-01,			/* 0x3f57bb00 */
 /*
- * Coefficients for approximation to erf on [0,0.84375]
+ * In the domain [0, 2**-14], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-42.
  */
-efx =  1.2837916613e-01, /* 0x3e0375d4 */
-efx8=  1.0270333290e+00, /* 0x3f8375d4 */
+efx = 1.28379166e-01, /* 0x3e0375d4 */
+efx8= 1.02703333e+00, /* 0x3f8375d4 */
 /*
- *  Domain [0, 0.84375], range ~[-5.4446e-10,5.5197e-10]:
- *  |(erf(x) - x)/x - p(x)/q(x)| < 2**-31.
+ * Domain [0, 0.84375], range ~[-5.4419e-10, 5.5179e-10]:
+ * |(erf(x) - x)/x - pp(x)/qq(x)| < 2**-31
  */
-pp0  =  1.28379166e-01F, /*  0x1.06eba8p-3 */
-pp1  = -3.36030394e-01F, /* -0x1.58185ap-2 */
-pp2  = -1.86260219e-03F, /* -0x1.e8451ep-10 */
-qq1  =  3.12324286e-01F, /*  0x1.3fd1f0p-2 */
-qq2  =  2.16070302e-02F, /*  0x1.620274p-6 */
-qq3  = -1.98859419e-03F, /* -0x1.04a626p-9 */
+pp0  =  1.28379166e-01, /* 0x3e0375d4 */
+pp1  = -3.36030394e-01, /* 0xbeac0c2d */
+pp2  = -1.86261395e-03, /* 0xbaf422f4 */
+qq1  =  3.12324315e-01, /* 0x3e9fe8f9 */
+qq2  =  2.16070414e-02, /* 0x3cb10140 */
+qq3  = -1.98859372e-03, /* 0xbb025311 */
 /*
- * Domain [0.84375, 1.25], range ~[-1.953e-11,1.940e-11]:
- * |(erf(x) - erx) - p(x)/q(x)| < 2**-36.
+ * Domain [0.84375, 1.25], range ~[-1.023e-9, 1.023e-9]:
+ * |(erf(x) - erx) - pa(x)/qa(x)| < 2**-31
  */
-erx  =  8.42697144e-01F, /*  0x1.af7600p-1.  erf(1) rounded to 16 bits. */
-pa0  =  3.64939137e-06F, /*  0x1.e9d022p-19 */
-pa1  =  4.15109694e-01F, /*  0x1.a91284p-2 */
-pa2  = -1.65179938e-01F, /* -0x1.5249dcp-3 */
-pa3  =  1.10914491e-01F, /*  0x1.c64e46p-4 */
-qa1  =  6.02074385e-01F, /*  0x1.344318p-1 */
-qa2  =  5.35934687e-01F, /*  0x1.126608p-1 */
-qa3  =  1.68576106e-01F, /*  0x1.593e6ep-3 */
-qa4  =  5.62181212e-02F, /*  0x1.cc89f2p-5 */
+pa0  =  3.65041046e-06, /* 0x3674f993 */
+pa1  =  4.15109307e-01, /* 0x3ed48935 */
+pa2  = -2.09395722e-01, /* 0xbe566bd5 */
+pa3  =  8.67677554e-02, /* 0x3db1b34b */
+qa1  =  4.95560974e-01, /* 0x3efdba2b */
+qa2  =  3.71248513e-01, /* 0x3ebe1449 */
+qa3  =  3.92478965e-02, /* 0x3d20c267 */
 /*
- * Domain [1.25,1/0.35], range ~[-7.043e-10,7.457e-10]:
- * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-30
+ * Domain [1.25,1/0.35], range ~[-4.821e-9, 4.927e-9]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - ra(x)/sa(x)| < 2**-28
  */
-ra0  = -9.87132732e-03F, /* -0x1.4376b2p-7 */
-ra1  = -5.53605914e-01F, /* -0x1.1b723cp-1 */
-ra2  = -2.17589188e+00F, /* -0x1.1683a0p+1 */
-ra3  = -1.43268085e+00F, /* -0x1.6ec42cp+0 */
-sa1  =  5.45995426e+00F, /*  0x1.5d6fe4p+2 */
-sa2  =  6.69798088e+00F, /*  0x1.acabb8p+2 */
-sa3  =  1.43113089e+00F, /*  0x1.6e5e98p+0 */
-sa4  = -5.77397496e-02F, /* -0x1.d90108p-5 */
+ra0  = -9.88156721e-03, /* 0xbc21e64c */
+ra1  = -5.43658376e-01, /* 0xbf0b2d32 */
+ra2  = -1.66828310e+00, /* 0xbfd58a4d */
+ra3  = -6.91554189e-01, /* 0xbf3109b2 */
+sa1  =  4.48581553e+00, /* 0x408f8bcd */
+sa2  =  4.10799170e+00, /* 0x408374ab */
+sa3  =  5.53855181e-01, /* 0x3f0dc974 */
 /*
- * Domain [1/0.35, 11], range ~[-2.264e-13,2.336e-13]:
- * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-42
+ * Domain [2.85715, 11], range ~[-1.484e-9, 1.505e-9]:
+ * |log(x*erfc(x)) + x**2 + 0.5625 - rb(x)/sb(x)| < 2**-30
  */
-rb0  = -9.86494310e-03F, /* -0x1.434124p-7 */
-rb1  = -6.25171244e-01F, /* -0x1.401672p-1 */
-rb2  = -6.16498327e+00F, /* -0x1.8a8f16p+2 */
-rb3  = -1.66696873e+01F, /* -0x1.0ab70ap+4 */
-rb4  = -9.53764343e+00F, /* -0x1.313460p+3 */
-sb1  =  1.26884899e+01F, /*  0x1.96081cp+3 */
-sb2  =  4.51839523e+01F, /*  0x1.6978bcp+5 */
-sb3  =  4.72810211e+01F, /*  0x1.7a3f88p+5 */
-sb4  =  8.93033314e+00F; /*  0x1.1dc54ap+3 */
+rb0  = -9.86496918e-03, /* 0xbc21a0ae */
+rb1  = -5.48049808e-01, /* 0xbf0c4cfe */
+rb2  = -1.84115684e+00, /* 0xbfebab07 */
+sb1  =  4.87132740e+00, /* 0x409be1ea */
+sb2  =  3.04982710e+00, /* 0x4043305e */
+sb3  = -7.61900663e-01; /* 0xbf430bec */
 
 float
 erff(float x)
@@ -85,9 +84,9 @@ erff(float x)
 	float R,S,P,Q,s,y,z,r;
 	GET_FLOAT_WORD(hx,x);
 	ix = hx&0x7fffffff;
-	if(ix>=0x7f800000) {		/* erf(nan)=nan */
+	if(ix>=0x7f800000) {		/* erff(nan)=nan */
 	    i = ((u_int32_t)hx>>31)<<1;
-	    return (float)(1-i)+one/x;	/* erf(+-inf)=+-1 */
+	    return (float)(1-i)+one/x;	/* erff(+-inf)=+-1 */
 	}
 
 	if(ix < 0x3f580000) {		/* |x|<0.84375 */
@@ -105,7 +104,7 @@ erff(float x)
 	if(ix < 0x3fa00000) {		/* 0.84375 <= |x| < 1.25 */
 	    s = fabsf(x)-one;
 	    P = pa0+s*(pa1+s*(pa2+s*pa3));
-	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*qa4)));
+	    Q = one+s*(qa1+s*(qa2+s*qa3));
 	    if(hx>=0) return erx + P/Q; else return -erx - P/Q;
 	}
 	if (ix >= 0x40800000) {		/* inf>|x|>=4 */
@@ -113,12 +112,12 @@ erff(float x)
 	}
 	x = fabsf(x);
  	s = one/(x*x);
-	if(ix< 0x4036DB6E) {	/* |x| < 1/0.35 */
+	if(ix< 0x4036db8c) {	/* |x| < 2.85715 ~ 1/0.35 */
 	    R=ra0+s*(ra1+s*(ra2+s*ra3));
-	    S=one+s*(sa1+s*(sa2+s*(sa3+s*sa4)));
-	} else {	/* |x| >= 1/0.35 */
-	    R=rb0+s*(rb1+s*(rb2+s*(rb3+s*rb4)));
-	    S=one+s*(sb1+s*(sb2+s*(sb3+s*sb4)));
+	    S=one+s*(sa1+s*(sa2+s*sa3));
+	} else {	/* |x| >= 2.85715 ~ 1/0.35 */
+	    R=rb0+s*(rb1+s*rb2);
+	    S=one+s*(sb1+s*(sb2+s*sb3));
 	}
 	SET_FLOAT_WORD(z,hx&0xffffe000);
 	r  = expf(-z*z-0.5625F)*expf((z-x)*(z+x)+R/S);
@@ -132,8 +131,8 @@ erfcf(float x)
 	float R,S,P,Q,s,y,z,r;
 	GET_FLOAT_WORD(hx,x);
 	ix = hx&0x7fffffff;
-	if(ix>=0x7f800000) {			/* erfc(nan)=nan */
-						/* erfc(+-inf)=0,2 */
+	if(ix>=0x7f800000) {			/* erfcf(nan)=nan */
+						/* erfcf(+-inf)=0,2 */
 	    return (float)(((u_int32_t)hx>>31)<<1)+one/x;
 	}
 
@@ -155,7 +154,7 @@ erfcf(float x)
 	if(ix < 0x3fa00000) {		/* 0.84375 <= |x| < 1.25 */
 	    s = fabsf(x)-one;
 	    P = pa0+s*(pa1+s*(pa2+s*pa3));
-	    Q = one+s*(qa1+s*(qa2+s*(qa3+s*qa4)));
+	    Q = one+s*(qa1+s*(qa2+s*qa3));
 	    if(hx>=0) {
 	        z  = one-erx; return z - P/Q;
 	    } else {
@@ -165,13 +164,13 @@ erfcf(float x)
 	if (ix < 0x41300000) {		/* |x|<11 */
 	    x = fabsf(x);
  	    s = one/(x*x);
-	    if(ix< 0x4036DB6D) {	/* |x| < 1/.35 ~ 2.857143*/
-	        R=ra0+s*(ra1+s*(ra2+s*ra3));
-	        S=one+s*(sa1+s*(sa2+s*(sa3+s*sa4)));
-	    } else {			/* |x| >= 1/.35 ~ 2.857143 */
+	    if(ix< 0x4036db8c) {	/* |x| < 2.85715 ~ 1/.35 */
+		R=ra0+s*(ra1+s*(ra2+s*ra3));
+		S=one+s*(sa1+s*(sa2+s*sa3));
+	    } else {			/* |x| >= 2.85715 ~ 1/.35 */
 		if(hx<0&&ix>=0x40a00000) return two-tiny;/* x < -5 */
-	        R=rb0+s*(rb1+s*(rb2+s*(rb3+s*rb4)));
-		S=one+s*(sb1+s*(sb2+s*(sb3+s*sb4)));
+		R=rb0+s*(rb1+s*rb2);
+		S=one+s*(sb1+s*(sb2+s*sb3));
 	    }
 	    SET_FLOAT_WORD(z,hx&0xffffe000);
 	    r  = expf(-z*z-0.5625F)*expf((z-x)*(z+x)+R/S);
diff --git a/lib/msun/src/s_round.c b/lib/msun/src/s_round.c
index 65de31b..fab3019 100644
--- a/lib/msun/src/s_round.c
+++ b/lib/msun/src/s_round.c
@@ -27,25 +27,34 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-#include <math.h>
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
 
 double
 round(double x)
 {
 	double t;
+	uint32_t hx;
 
-	if (!isfinite(x))
-		return (x);
+	GET_HIGH_WORD(hx, x);
+	if ((hx & 0x7fffffff) == 0x7ff00000)
+		return (x + x);
 
-	if (x >= 0.0) {
+	if (!(hx & 0x80000000)) {
 		t = floor(x);
 		if (t - x <= -0.5)
-			t += 1.0;
+			t += 1;
 		return (t);
 	} else {
 		t = floor(-x);
 		if (t + x <= -0.5)
-			t += 1.0;
+			t += 1;
 		return (-t);
 	}
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(round, roundl);
+#endif
diff --git a/lib/msun/src/s_roundf.c b/lib/msun/src/s_roundf.c
index 952e8e7..e7e2eb9 100644
--- a/lib/msun/src/s_roundf.c
+++ b/lib/msun/src/s_roundf.c
@@ -27,25 +27,28 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-#include <math.h>
+#include "math.h"
+#include "math_private.h"
 
 float
 roundf(float x)
 {
 	float t;
+	uint32_t hx;
 
-	if (!isfinite(x))
-		return (x);
+	GET_FLOAT_WORD(hx, x);
+	if ((hx & 0x7fffffff) == 0x7f800000)
+		return (x + x);
 
-	if (x >= 0.0) {
+	if (!(hx & 0x80000000)) {
 		t = floorf(x);
-		if (t - x <= -0.5)
-			t += 1.0;
+		if (t - x <= -0.5F)
+			t += 1;
 		return (t);
 	} else {
 		t = floorf(-x);
-		if (t + x <= -0.5)
-			t += 1.0;
+		if (t + x <= -0.5F)
+			t += 1;
 		return (-t);
 	}
 }
diff --git a/lib/msun/src/s_roundl.c b/lib/msun/src/s_roundl.c
index a70b617..2d15e13 100644
--- a/lib/msun/src/s_roundl.c
+++ b/lib/msun/src/s_roundl.c
@@ -27,25 +27,36 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-#include <math.h>
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
 
 long double
 roundl(long double x)
 {
 	long double t;
+	uint16_t hx;
+
+	GET_LDBL_EXPSIGN(hx, x);
+	if ((hx & 0x7fff) == 0x7fff)
+		return (x + x);
 
-	if (!isfinite(x))
-		return (x);
+	ENTERI();
 
-	if (x >= 0.0) {
+	if (!(hx & 0x8000)) {
 		t = floorl(x);
-		if (t - x <= -0.5)
-			t += 1.0;
-		return (t);
+		if (t - x <= -0.5L)
+			t += 1;
+		RETURNI(t);
 	} else {
 		t = floorl(-x);
-		if (t + x <= -0.5)
-			t += 1.0;
-		return (-t);
+		if (t + x <= -0.5L)
+			t += 1;
+		RETURNI(-t);
 	}
 }
diff --git a/lib/msun/src/s_tanh.c b/lib/msun/src/s_tanh.c
index 96e3565..6d26c69 100644
--- a/lib/msun/src/s_tanh.c
+++ b/lib/msun/src/s_tanh.c
@@ -37,10 +37,13 @@ __FBSDID("$FreeBSD$");
  *	only tanh(0)=0 is exact for finite argument.
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
-static const double one = 1.0, two = 2.0, tiny = 1.0e-300, huge = 1.0e300;
+static const volatile double tiny = 1.0e-300;
+static const double one = 1.0, two = 2.0, huge = 1.0e300;
 
 double
 tanh(double x)
@@ -75,3 +78,7 @@ tanh(double x)
 	}
 	return (jx>=0)? z: -z;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(tanh, tanhl);
+#endif
diff --git a/lib/msun/src/s_tanhf.c b/lib/msun/src/s_tanhf.c
index 04f09c6..f537be4 100644
--- a/lib/msun/src/s_tanhf.c
+++ b/lib/msun/src/s_tanhf.c
@@ -19,7 +19,9 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
-static const float one=1.0, two=2.0, tiny = 1.0e-30, huge = 1.0e30;
+static const volatile float tiny = 1.0e-30;
+static const float one=1.0, two=2.0, huge = 1.0e30;
+
 float
 tanhf(float x)
 {
diff --git a/lib/msun/src/s_tanhl.c b/lib/msun/src/s_tanhl.c
new file mode 100644
index 0000000..886158b
--- /dev/null
+++ b/lib/msun/src/s_tanhl.c
@@ -0,0 +1,172 @@
+/* from: FreeBSD: head/lib/msun/src/s_tanhl.c XXX */
+
+/* @(#)s_tanh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See s_tanh.c for complete comments.
+ *
+ * Converted to long double by Bruce D. Evans.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "math.h"
+#include "math_private.h"
+#include "fpmath.h"
+#include "k_expl.h"
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual expsign encoding. */
+#error "Unsupported long double format"
+#endif
+
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+static const volatile double tiny = 1.0e-300;
+static const double one = 1.0;
+#if LDBL_MANT_DIG == 64
+/*
+ * Domain [-0.25, 0.25], range ~[-1.6304e-22, 1.6304e-22]:
+ * |tanh(x)/x - t(x)| < 2**-72.3
+ */
+static const union IEEEl2bits
+T3u = LD80C(0xaaaaaaaaaaaaaa9f, -2, -3.33333333333333333017e-1L);
+#define	T3	T3u.e
+static const double
+T5  =  1.3333333333333314e-1,		/*  0x1111111111110a.0p-55 */
+T7  = -5.3968253968210485e-2,		/* -0x1ba1ba1ba1a1a1.0p-57 */
+T9  =  2.1869488531393817e-2,		/*  0x1664f488172022.0p-58 */
+T11 = -8.8632352345964591e-3,		/* -0x1226e34bc138d5.0p-59 */
+T13 =  3.5921169709993771e-3,		/*  0x1d6d371d3e400f.0p-61 */
+T15 = -1.4555786415756001e-3,		/* -0x17d923aa63814d.0p-62 */
+T17 =  5.8645267876296793e-4,		/*  0x13378589b85aa7.0p-63 */
+T19 = -2.1121033571392224e-4;		/* -0x1baf0af80c4090.0p-65 */
+#elif LDBL_MANT_DIG == 113
+/*
+ * Domain [-0.25, 0.25], range ~[-2.4211e-37, 2.4211e-37]:
+ * |tanh(x)/x - t(x)| < 2**121.6
+ */
+static const long double
+T3 = -3.33333333333333333333333333333332980e-1L,	/* -0x1555555555555555555555555554e.0p-114L */
+T5  =  1.33333333333333333333333333332707260e-1L,	/*  0x1111111111111111111111110ab7b.0p-115L */
+T7  = -5.39682539682539682539682535723482314e-2L,	/* -0x1ba1ba1ba1ba1ba1ba1ba17b5fc98.0p-117L */
+T9  =  2.18694885361552028218693591149061717e-2L,	/*  0x1664f4882c10f9f32d6b1a12a25e5.0p-118L */
+T11 = -8.86323552990219656883762347736381851e-3L,	/* -0x1226e355e6c23c8f5a5a0f386cb4d.0p-119L */
+T13 =  3.59212803657248101358314398220822722e-3L,	/*  0x1d6d3d0e157ddfb403ad3637442c6.0p-121L */
+T15 = -1.45583438705131796512568010348874662e-3L;	/* -0x17da36452b75e150c44cc34253b34.0p-122L */
+static const double
+T17 =  5.9002744094556621e-4,		/*  0x1355824803668e.0p-63 */
+T19 = -2.3912911424260516e-4,		/* -0x1f57d7734c8dde.0p-65 */
+T21 =  9.6915379535512898e-5,		/*  0x1967e18ad6a6ca.0p-66 */
+T23 = -3.9278322983156353e-5,		/* -0x1497d8e6b75729.0p-67 */
+T25 =  1.5918887220143869e-5,		/*  0x10b1319998cafa.0p-68 */
+T27 = -6.4514295231630956e-6,		/* -0x1b0f2b71b218eb.0p-70 */
+T29 =  2.6120754043964365e-6,		/*  0x15e963a3cf3a39.0p-71 */
+T31 = -1.0407567231003314e-6,		/* -0x1176041e656869.0p-72 */
+T33 =  3.4744117554063574e-7;		/*  0x1750fe732cab9c.0p-74 */
+#endif /* LDBL_MANT_DIG == 64 */
+
+static inline long double
+divl(long double a, long double b, long double c, long double d,
+    long double e, long double f)
+{
+	long double inv, r;
+	float fr, fw;
+
+	_2sumF(a, c);
+	b = b + c;
+	_2sumF(d, f);
+	e = e + f;
+
+	inv = 1 / (d + e);
+
+	r = (a + b) * inv;
+	fr = r;
+	r = fr;
+
+	fw = d + e;
+	e = d - fw + e;
+	d = fw;
+
+	r = r + (a - d * r + b - e * r) * inv;
+
+	return r;
+}
+
+long double
+tanhl(long double x)
+{
+	long double hi,lo,s,x2,x4,z;
+	double dx2;
+	int16_t jx,ix;
+
+	GET_LDBL_EXPSIGN(jx,x);
+	ix = jx&0x7fff;
+
+    /* x is INF or NaN */
+	if(ix>=0x7fff) {
+	    if (jx>=0) return one/x+one;    /* tanh(+-inf)=+-1 */
+	    else       return one/x-one;    /* tanh(NaN) = NaN */
+	}
+
+	ENTERI();
+
+    /* |x| < 40 */
+	if (ix < 0x4004 || fabsl(x) < 40) {	/* |x|<40 */
+	    if (__predict_false(ix<BIAS-(LDBL_MANT_DIG+1)/2)) {	/* |x|<TINY */
+		/* tanh(+-0) = +0; tanh(tiny) = tiny(-+) with inexact: */
+		return (x == 0 ? x : (0x1p200 * x - x) * 0x1p-200);
+	    }
+	    if (ix<0x3ffd) {		/* |x|<0.25 */
+		x2 = x*x;
+#if LDBL_MANT_DIG == 64
+		x4 = x2*x2;
+		RETURNI(((T19*x2 + T17)*x4 + (T15*x2 + T13))*(x2*x*x2*x4*x4) +
+		    ((T11*x2 + T9)*x4 + (T7*x2 + T5))*(x2*x*x2) +
+		    T3*(x2*x) + x);
+#elif LDBL_MANT_DIG == 113
+		dx2 = x2;
+#if 0
+		RETURNI(((((((((((((((T33*dx2 + T31)*dx2 + T29)*dx2 + T27)*dx2 +
+		    T25)*x2 + T23)*x2 + T21)*x2 + T19)*x2 + T17)*x2 +
+		    T15)*x2 + T13)*x2 + T11)*x2 + T9)*x2 + T7)*x2 + T5)*
+		    (x2*x*x2) +
+		    T3*(x2*x) + x);
+#else
+		long double q = ((((((((((((((T33*dx2 + T31)*dx2 + T29)*dx2 + T27)*dx2 +
+		    T25)*x2 + T23)*x2 + T21)*x2 + T19)*x2 + T17)*x2 +
+		    T15)*x2 + T13)*x2 + T11)*x2 + T9)*x2 + T7)*x2 + T5)*
+		    (x2*x*x2);
+		RETURNI(q + T3*(x2*x) + x);
+#endif
+#endif
+	    }
+	    k_hexpl(2*fabsl(x), &hi, &lo);
+	    if (ix<0x4001 && fabsl(x) < 1.5)	/* |x|<1.5 */
+		z = divl(hi, lo, -0.5, hi, lo, 0.5);
+	    else
+		z = one - one/(lo+0.5+hi);
+    /* |x| >= 40, return +-1 */
+	} else {
+	    z = one - tiny;		/* raise inexact flag */
+	}
+	s = 1;
+	if (jx<0) s = -1;
+	RETURNI(s*z);
+}