- 80-bit and quad precision trigonometric and other most

important functions: acosl, asinl, atanl, atan2l, cosl,
sinl, tanl, exp2l, frexpl, ilogbl, ldexpl, logbl, scalbnl,
fabsl, hypotl, powl, sqrtl, rintl, copysignl, nanl, fdiml,
fmaxl, fminl.  mostly taken from freebsd, needed alot of
changes to adapt.  note, these are all c versions;  and are
quite slow when architectures have, e.g. sqrt.  assembly
versions will be added afterwards
- make them .weak/__weak_alias to the double precision
versions on other archs
- no need to have two finites.  finite() and finitef() are
non-standard 3BSD obsolete versions of isfinite.  remove
from libm.  make them weak_alias in libc to __isfinite and
__isfinitef instead.  similarly make 3BSD obsolete versions
of isinf, isinff, isnan, isnanf weak_aliases to C99's
__isinf, __isinff, __isnan, __isnanf
- remove unused infinity.c.  the c library has infinities
for each supported platform
- use STRICT_ASSIGN cast hack for _kernel_rem_pio2, so that
the double version has a chance of working on i386 with
extra precision
- avoid storing multiple copies of the pi/2 array, since
it won't vary
- bump major due to removed finite/finitef.  although they
will be in libc, which anything is linked to, minor bump
might be enough
ok millert@.   tested by sthen@, jsg@, ajacoutot@, kili@, naddy@

7 files changed, 940 insertions, 0 deletions

diff --git a/lib/libm/src/ld128/invtrig.c b/lib/libm/src/ld128/invtrig.c
new file mode 100644
index 00000000000..15db44a1804
--- /dev/null
+++ b/lib/libm/src/ld128/invtrig.c
@@ -0,0 +1,98 @@
+/*	$OpenBSD: invtrig.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*-
+ * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include "invtrig.h"
+
+/*
+ * asinl() and acosl()
+ */
+const long double
+pS0 =  1.66666666666666666666666666666700314e-01L,
+pS1 = -7.32816946414566252574527475428622708e-01L,
+pS2 =  1.34215708714992334609030036562143589e+00L,
+pS3 = -1.32483151677116409805070261790752040e+00L,
+pS4 =  7.61206183613632558824485341162121989e-01L,
+pS5 = -2.56165783329023486777386833928147375e-01L,
+pS6 =  4.80718586374448793411019434585413855e-02L,
+pS7 = -4.42523267167024279410230886239774718e-03L,
+pS8 =  1.44551535183911458253205638280410064e-04L,
+pS9 = -2.10558957916600254061591040482706179e-07L,
+qS1 = -4.84690167848739751544716485245697428e+00L,
+qS2 =  9.96619113536172610135016921140206980e+00L,
+qS3 = -1.13177895428973036660836798461641458e+01L,
+qS4 =  7.74004374389488266169304117714658761e+00L,
+qS5 = -3.25871986053534084709023539900339905e+00L,
+qS6 =  8.27830318881232209752469022352928864e-01L,
+qS7 = -1.18768052702942805423330715206348004e-01L,
+qS8 =  8.32600764660522313269101537926539470e-03L,
+qS9 = -1.99407384882605586705979504567947007e-04L;
+
+/*
+ * atanl()
+ */
+const long double atanhi[] = {
+	 4.63647609000806116214256231461214397e-01L,
+	 7.85398163397448309615660845819875699e-01L,       
+	 9.82793723247329067985710611014666038e-01L,       
+	 1.57079632679489661923132169163975140e+00L,
+};
+
+const long double atanlo[] = {
+	 4.89509642257333492668618435220297706e-36L,
+	 2.16795253253094525619926100651083806e-35L,
+	-2.31288434538183565909319952098066272e-35L,
+	 4.33590506506189051239852201302167613e-35L,
+};
+
+const long double aT[] = {
+	 3.33333333333333333333333333333333125e-01L,
+	-1.99999999999999999999999999999180430e-01L,
+	 1.42857142857142857142857142125269827e-01L,
+	-1.11111111111111111111110834490810169e-01L,
+	 9.09090909090909090908522355708623681e-02L,
+	-7.69230769230769230696553844935357021e-02L,
+	 6.66666666666666660390096773046256096e-02L,
+	-5.88235294117646671706582985209643694e-02L,
+	 5.26315789473666478515847092020327506e-02L,
+	-4.76190476189855517021024424991436144e-02L,
+	 4.34782608678695085948531993458097026e-02L,
+	-3.99999999632663469330634215991142368e-02L,
+	 3.70370363987423702891250829918659723e-02L,
+	-3.44827496515048090726669907612335954e-02L,
+	 3.22579620681420149871973710852268528e-02L,
+	-3.03020767654269261041647570626778067e-02L,
+	 2.85641979882534783223403715930946138e-02L,
+	-2.69824879726738568189929461383741323e-02L,
+	 2.54194698498808542954187110873675769e-02L,
+	-2.35083879708189059926183138130183215e-02L,
+	 2.04832358998165364349957325067131428e-02L,
+	-1.54489555488544397858507248612362957e-02L,
+	 8.64492360989278761493037861575248038e-03L,
+	-2.58521121597609872727919154569765469e-03L,
+};
+
+const long double pi_lo = 8.67181013012378102479704402604335225e-35L;
diff --git a/lib/libm/src/ld128/invtrig.h b/lib/libm/src/ld128/invtrig.h
new file mode 100644
index 00000000000..89b49766b86
--- /dev/null
+++ b/lib/libm/src/ld128/invtrig.h
@@ -0,0 +1,118 @@
+/*	$OpenBSD: invtrig.h,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*-
+ * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD: src/lib/msun/ld128/invtrig.h,v 1.1 2008/07/31 22:41:26 das Exp $
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+
+#define	BIAS		(LDBL_MAX_EXP - 1)
+#define	MANH_SIZE	(EXT_FRACHBITS + EXT_FRACHMBITS + 1)
+
+/* Approximation thresholds. */
+#define	ASIN_LINEAR	(BIAS - 56)	/* 2**-56 */
+#define	ACOS_CONST	(BIAS - 113)	/* 2**-113 */
+#define	ATAN_CONST	(BIAS + 113)	/* 2**113 */
+#define	ATAN_LINEAR	(BIAS - 56)	/* 2**-56 */
+
+/* 0.95 */
+#ifdef EXT_IMPLICIT_NBIT
+#define	THRESH	(0xe666666666666666ULL>>(64-(MANH_SIZE-1)))
+#else /* EXT_IMPLICIT_NBIT */
+#define	THRESH	((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|0x80000000)
+#endif /* EXT_IMPLICIT_NBIT */
+
+/* Constants shared by the long double inverse trig functions. */
+#define	pS0	_ItL_pS0
+#define	pS1	_ItL_pS1
+#define	pS2	_ItL_pS2
+#define	pS3	_ItL_pS3
+#define	pS4	_ItL_pS4
+#define	pS5	_ItL_pS5
+#define	pS6	_ItL_pS6
+#define	pS7	_ItL_pS7
+#define	pS8	_ItL_pS8
+#define	pS9	_ItL_pS9
+#define	qS1	_ItL_qS1
+#define	qS2	_ItL_qS2
+#define	qS3	_ItL_qS3
+#define	qS4	_ItL_qS4
+#define	qS5	_ItL_qS5
+#define	qS6	_ItL_qS6
+#define	qS7	_ItL_qS7
+#define	qS8	_ItL_qS8
+#define	qS9	_ItL_qS9
+#define	atanhi	_ItL_atanhi
+#define	atanlo	_ItL_atanlo
+#define	aT	_ItL_aT
+#define	pi_lo	_ItL_pi_lo
+
+#define	pio2_hi	atanhi[3]
+#define	pio2_lo	atanlo[3]
+#define	pio4_hi	atanhi[1]
+
+/* Constants shared by the long double inverse trig functions. */
+extern const long double pS0, pS1, pS2, pS3, pS4, pS5, pS6, pS7, pS8, pS9;
+extern const long double qS1, qS2, qS3, qS4, qS5, qS6, qS7, qS8, qS9;
+extern const long double atanhi[], atanlo[], aT[];
+extern const long double pi_lo;
+
+static inline long double
+P(long double x)
+{
+
+	return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \
+		(pS4 + x * (pS5 + x * (pS6 + x * (pS7 + x * (pS8 + x * \
+		pS9))))))))));
+}
+
+static inline long double
+Q(long double x)
+{
+
+	return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * \
+		(qS5 + x * (qS6 + x * (qS7 + x * (qS8 + x * qS9)))))))));
+}
+
+static inline long double
+T_even(long double x)
+{
+
+	return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * \
+		(aT[8] + x * (aT[10] + x * (aT[12] + x * (aT[14] + x * \
+		(aT[16] + x * (aT[18] + x * (aT[20] + x * aT[22])))))))))));
+}
+
+static inline long double
+T_odd(long double x)
+{
+
+	return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * \
+		(aT[9] + x * (aT[11] + x * (aT[13] + x * (aT[15] + x * \
+		(aT[17] + x * (aT[19] + x * (aT[21] + x * aT[23])))))))))));
+}
diff --git a/lib/libm/src/ld128/k_cosl.c b/lib/libm/src/ld128/k_cosl.c
new file mode 100644
index 00000000000..70966786aed
--- /dev/null
+++ b/lib/libm/src/ld128/k_cosl.c
@@ -0,0 +1,59 @@
+/*	$OpenBSD: k_cosl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/* From: @(#)k_cos.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * ld128 version of k_cos.c.  See ../k_cos.c for most comments.
+ */
+
+#include "math_private.h"
+
+/*
+ * Domain [-0.7854, 0.7854], range ~[-1.80e-37, 1.79e-37]:
+ * |cos(x) - c(x))| < 2**-122.0
+ *
+ * 113-bit precision requires more care than 64-bit precision, since
+ * simple methods give a minimax polynomial with coefficient for x^2
+ * that is 1 ulp below 0.5, but we want it to be precisely 0.5.  See
+ * ../ld80/k_cosl.c for more details.
+ */
+static const double
+one = 1.0;
+
+static const long double
+C1 =  0.04166666666666666666666666666666658424671L,
+C2 = -0.001388888888888888888888888888863490893732L,
+C3 =  0.00002480158730158730158730158600795304914210L,
+C4 = -0.2755731922398589065255474947078934284324e-6L,
+C5 =  0.2087675698786809897659225313136400793948e-8L,
+C6 = -0.1147074559772972315817149986812031204775e-10L,
+C7 =  0.4779477332386808976875457937252120293400e-13L;
+
+static const double
+C8 = -0.1561920696721507929516718307820958119868e-15,
+C9 =  0.4110317413744594971475941557607804508039e-18,
+C10 = -0.8896592467191938803288521958313920156409e-21,
+C11 =  0.1601061435794535138244346256065192782581e-23;
+
+long double
+__kernel_cosl(long double x, long double y)
+{
+	long double hz,z,r,w;
+
+	z  = x*x;
+	r  = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*(C7+
+	    z*(C8+z*(C9+z*(C10+z*C11))))))))));
+	hz = 0.5*z;
+	w  = one-hz;
+	return w + (((one-w)-hz) + (z*r-x*y));
+}
diff --git a/lib/libm/src/ld128/k_sinl.c b/lib/libm/src/ld128/k_sinl.c
new file mode 100644
index 00000000000..15a6cedc22f
--- /dev/null
+++ b/lib/libm/src/ld128/k_sinl.c
@@ -0,0 +1,57 @@
+/*	$OpenBSD: k_sinl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/* From: @(#)k_sin.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * ld128 version of k_sin.c.  See ../k_sin.c for most comments.
+ */
+
+#include "math_private.h"
+
+static const double
+half =  0.5;
+
+/*
+ * Domain [-0.7854, 0.7854], range ~[-1.53e-37, 1.659e-37]
+ * |sin(x)/x - s(x)| < 2**-122.1
+ *
+ * See ../ld80/k_cosl.c for more details about the polynomial.
+ */
+static const long double
+S1 = -0.16666666666666666666666666666666666606732416116558L,
+S2 =  0.0083333333333333333333333333333331135404851288270047L,
+S3 = -0.00019841269841269841269841269839935785325638310428717L,
+S4 =  0.27557319223985890652557316053039946268333231205686e-5L,
+S5 = -0.25052108385441718775048214826384312253862930064745e-7L,
+S6 =  0.16059043836821614596571832194524392581082444805729e-9L,
+S7 = -0.76471637318198151807063387954939213287488216303768e-12L,
+S8 =  0.28114572543451292625024967174638477283187397621303e-14L;
+
+static const double
+S9  = -0.82206352458348947812512122163446202498005154296863e-17,
+S10 =  0.19572940011906109418080609928334380560135358385256e-19,
+S11 = -0.38680813379701966970673724299207480965452616911420e-22,
+S12 =  0.64038150078671872796678569586315881020659912139412e-25;
+
+long double
+__kernel_sinl(long double x, long double y, int iy)
+{
+	long double z,r,v;
+
+	z	=  x*x;
+	v	=  z*x;
+	r	=  S2+z*(S3+z*(S4+z*(S5+z*(S6+z*(S7+z*(S8+
+	    z*(S9+z*(S10+z*(S11+z*S12)))))))));
+	if(iy==0) return x+v*(S1+z*r);
+	else      return x-((z*(half*y-v*r)-y)-v*S1);
+}
diff --git a/lib/libm/src/ld128/k_tanl.c b/lib/libm/src/ld128/k_tanl.c
new file mode 100644
index 00000000000..43b1ad87ce8
--- /dev/null
+++ b/lib/libm/src/ld128/k_tanl.c
@@ -0,0 +1,117 @@
+/*	$OpenBSD: k_tanl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/* From: @(#)k_tan.c 1.5 04/04/22 SMI */
+/*
+ * ====================================================
+ * Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
+ * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * ld128 version of k_tan.c.  See ../k_tan.c for most comments.
+ */
+
+#include <math.h>
+
+#include "math_private.h"
+
+/*
+ * Domain [-0.67434, 0.67434], range ~[-3.37e-36, 1.982e-37]
+ * |tan(x)/x - t(x)| < 2**-117.8 (XXX should be ~1e-37)
+ *
+ * See ../ld80/k_cosl.c for more details about the polynomial.
+ */
+static const long double
+T3 = 0x1.5555555555555555555555555553p-2L,
+T5 = 0x1.1111111111111111111111111eb5p-3L,
+T7 = 0x1.ba1ba1ba1ba1ba1ba1ba1b694cd6p-5L,
+T9 = 0x1.664f4882c10f9f32d6bbe09d8bcdp-6L,
+T11 = 0x1.226e355e6c23c8f5b4f5762322eep-7L,
+T13 = 0x1.d6d3d0e157ddfb5fed8e84e27b37p-9L,
+T15 = 0x1.7da36452b75e2b5fce9ee7c2c92ep-10L,
+T17 = 0x1.355824803674477dfcf726649efep-11L,
+T19 = 0x1.f57d7734d1656e0aceb716f614c2p-13L,
+T21 = 0x1.967e18afcb180ed942dfdc518d6cp-14L,
+T23 = 0x1.497d8eea21e95bc7e2aa79b9f2cdp-15L,
+T25 = 0x1.0b132d39f055c81be49eff7afd50p-16L,
+T27 = 0x1.b0f72d33eff7bfa2fbc1059d90b6p-18L,
+T29 = 0x1.5ef2daf21d1113df38d0fbc00267p-19L,
+T31 = 0x1.1c77d6eac0234988cdaa04c96626p-20L,
+T33 = 0x1.cd2a5a292b180e0bdd701057dfe3p-22L,
+T35 = 0x1.75c7357d0298c01a31d0a6f7d518p-23L,
+T37 = 0x1.2f3190f4718a9a520f98f50081fcp-24L,
+pio4 = 0x1.921fb54442d18469898cc51701b8p-1L,
+pio4lo = 0x1.cd129024e088a67cc74020bbea60p-116L;
+
+static const double
+T39 =  0.000000028443389121318352,	/*  0x1e8a7592977938.0p-78 */
+T41 =  0.000000011981013102001973,	/*  0x19baa1b1223219.0p-79 */
+T43 =  0.0000000038303578044958070,	/*  0x107385dfb24529.0p-80 */
+T45 =  0.0000000034664378216909893,	/*  0x1dc6c702a05262.0p-81 */
+T47 = -0.0000000015090641701997785,	/* -0x19ecef3569ebb6.0p-82 */
+T49 =  0.0000000029449552300483952,	/*  0x194c0668da786a.0p-81 */
+T51 = -0.0000000022006995706097711,	/* -0x12e763b8845268.0p-81 */
+T53 =  0.0000000015468200913196612,	/*  0x1a92fc98c29554.0p-82 */
+T55 = -0.00000000061311613386849674,	/* -0x151106cbc779a9.0p-83 */
+T57 =  1.4912469681508012e-10;		/*  0x147edbdba6f43a.0p-85 */
+
+long double
+__kernel_tanl(long double x, long double y, int iy) {
+	long double z, r, v, w, s;
+	long double osign;
+	int i;
+
+	iy = (iy == 1 ? -1 : 1);	/* XXX recover original interface */
+	osign = (x >= 0 ? 1.0 : -1.0);	/* XXX slow, probably wrong for -0 */
+	if (fabsl(x) >= 0.67434) {
+		if (x < 0) {
+			x = -x;
+			y = -y;
+		}
+		z = pio4 - x;
+		w = pio4lo - y;
+		x = z + w;
+		y = 0.0;
+		i = 1;
+	} else
+		i = 0;
+	z = x * x;
+	w = z * z;
+	r = T5 + w * (T9 + w * (T13 + w * (T17 + w * (T21 +
+	    w * (T25 + w * (T29 + w * (T33 +
+	    w * (T37 + w * (T41 + w * (T45 + w * (T49 + w * (T53 +
+	    w * T57))))))))))));
+	v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
+	    w * (T27 + w * (T31 + w * (T35 +
+	    w * (T39 + w * (T43 + w * (T47 + w * (T51 + w * T55))))))))))));
+	s = z * x;
+	r = y + z * (s * (r + v) + y);
+	r += T3 * s;
+	w = x + r;
+	if (i == 1) {
+		v = (long double) iy;
+		return osign *
+			(v - 2.0 * (x - (w * w / (w + v) - r)));
+	}
+	if (iy == 1)
+		return w;
+	else {
+		/*
+		 * if allow error up to 2 ulp, simply return
+		 * -1.0 / (x+r) here
+		 */
+		/* compute -1.0 / (x+r) accurately */
+		long double a, t;
+		z = w;
+		z = z + 0x1p32 - 0x1p32;
+		v = r - (z - x);	/* z+v = r+x */
+		t = a = -1.0 / w;	/* a = -1.0/w */
+		t = t + 0x1p32 - 0x1p32;
+		s = 1.0 + t * z;
+		return t + a * (s + t * v);
+	}
+}
diff --git a/lib/libm/src/ld128/s_exp2l.c b/lib/libm/src/ld128/s_exp2l.c
new file mode 100644
index 00000000000..1bba73bb688
--- /dev/null
+++ b/lib/libm/src/ld128/s_exp2l.c
@@ -0,0 +1,441 @@
+/*	$OpenBSD: s_exp2l.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*-
+ * Copyright (c) 2005-2008 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include <sys/types.h>
+#include <sys/endian.h>
+#include <machine/ieee.h>
+#include <float.h>
+#include <math.h>
+#include <stdint.h>
+
+#define	TBLBITS	7
+#define	TBLSIZE	(1 << TBLBITS)
+
+#define	BIAS	(LDBL_MAX_EXP - 1)
+#define	EXPMASK	(BIAS + LDBL_MAX_EXP)
+
+#if 0 /* XXX Prevent gcc from erroneously constant folding this. */
+static const long double twom10000 = 0x1p-10000L;
+#else
+static volatile long double twom10000 = 0x1p-10000L;
+#endif
+
+static const long double
+    huge      = 0x1p10000L,
+    P1        = 0x1.62e42fefa39ef35793c7673007e6p-1L,
+    P2	      = 0x1.ebfbdff82c58ea86f16b06ec9736p-3L,
+    P3        = 0x1.c6b08d704a0bf8b33a762bad3459p-5L,
+    P4        = 0x1.3b2ab6fba4e7729ccbbe0b4f3fc2p-7L,
+    P5        = 0x1.5d87fe78a67311071dee13fd11d9p-10L,
+    P6        = 0x1.430912f86c7876f4b663b23c5fe5p-13L;
+
+static const double
+    P7        = 0x1.ffcbfc588b041p-17,
+    P8        = 0x1.62c0223a5c7c7p-20,
+    P9        = 0x1.b52541ff59713p-24,
+    P10       = 0x1.e4cf56a391e22p-28,
+    redux     = 0x1.8p112 / TBLSIZE;
+
+static const long double tbl[TBLSIZE] = {
+	0x1.6a09e667f3bcc908b2fb1366dfeap-1L,
+	0x1.6c012750bdabeed76a99800f4edep-1L,
+	0x1.6dfb23c651a2ef220e2cbe1bc0d4p-1L,
+	0x1.6ff7df9519483cf87e1b4f3e1e98p-1L,
+	0x1.71f75e8ec5f73dd2370f2ef0b148p-1L,
+	0x1.73f9a48a58173bd5c9a4e68ab074p-1L,
+	0x1.75feb564267c8bf6e9aa33a489a8p-1L,
+	0x1.780694fde5d3f619ae02808592a4p-1L,
+	0x1.7a11473eb0186d7d51023f6ccb1ap-1L,
+	0x1.7c1ed0130c1327c49334459378dep-1L,
+	0x1.7e2f336cf4e62105d02ba1579756p-1L,
+	0x1.80427543e1a11b60de67649a3842p-1L,
+	0x1.82589994cce128acf88afab34928p-1L,
+	0x1.8471a4623c7acce52f6b97c6444cp-1L,
+	0x1.868d99b4492ec80e41d90ac2556ap-1L,
+	0x1.88ac7d98a669966530bcdf2d4cc0p-1L,
+	0x1.8ace5422aa0db5ba7c55a192c648p-1L,
+	0x1.8cf3216b5448bef2aa1cd161c57ap-1L,
+	0x1.8f1ae991577362b982745c72eddap-1L,
+	0x1.9145b0b91ffc588a61b469f6b6a0p-1L,
+	0x1.93737b0cdc5e4f4501c3f2540ae8p-1L,
+	0x1.95a44cbc8520ee9b483695a0e7fep-1L,
+	0x1.97d829fde4e4f8b9e920f91e8eb6p-1L,
+	0x1.9a0f170ca07b9ba3109b8c467844p-1L,
+	0x1.9c49182a3f0901c7c46b071f28dep-1L,
+	0x1.9e86319e323231824ca78e64c462p-1L,
+	0x1.a0c667b5de564b29ada8b8cabbacp-1L,
+	0x1.a309bec4a2d3358c171f770db1f4p-1L,
+	0x1.a5503b23e255c8b424491caf88ccp-1L,
+	0x1.a799e1330b3586f2dfb2b158f31ep-1L,
+	0x1.a9e6b5579fdbf43eb243bdff53a2p-1L,
+	0x1.ac36bbfd3f379c0db966a3126988p-1L,
+	0x1.ae89f995ad3ad5e8734d17731c80p-1L,
+	0x1.b0e07298db66590842acdfc6fb4ep-1L,
+	0x1.b33a2b84f15faf6bfd0e7bd941b0p-1L,
+	0x1.b59728de559398e3881111648738p-1L,
+	0x1.b7f76f2fb5e46eaa7b081ab53ff6p-1L,
+	0x1.ba5b030a10649840cb3c6af5b74cp-1L,
+	0x1.bcc1e904bc1d2247ba0f45b3d06cp-1L,
+	0x1.bf2c25bd71e088408d7025190cd0p-1L,
+	0x1.c199bdd85529c2220cb12a0916bap-1L,
+	0x1.c40ab5fffd07a6d14df820f17deap-1L,
+	0x1.c67f12e57d14b4a2137fd20f2a26p-1L,
+	0x1.c8f6d9406e7b511acbc48805c3f6p-1L,
+	0x1.cb720dcef90691503cbd1e949d0ap-1L,
+	0x1.cdf0b555dc3f9c44f8958fac4f12p-1L,
+	0x1.d072d4a07897b8d0f22f21a13792p-1L,
+	0x1.d2f87080d89f18ade123989ea50ep-1L,
+	0x1.d5818dcfba48725da05aeb66dff8p-1L,
+	0x1.d80e316c98397bb84f9d048807a0p-1L,
+	0x1.da9e603db3285708c01a5b6d480cp-1L,
+	0x1.dd321f301b4604b695de3c0630c0p-1L,
+	0x1.dfc97337b9b5eb968cac39ed284cp-1L,
+	0x1.e264614f5a128a12761fa17adc74p-1L,
+	0x1.e502ee78b3ff6273d130153992d0p-1L,
+	0x1.e7a51fbc74c834b548b2832378a4p-1L,
+	0x1.ea4afa2a490d9858f73a18f5dab4p-1L,
+	0x1.ecf482d8e67f08db0312fb949d50p-1L,
+	0x1.efa1bee615a27771fd21a92dabb6p-1L,
+	0x1.f252b376bba974e8696fc3638f24p-1L,
+	0x1.f50765b6e4540674f84b762861a6p-1L,
+	0x1.f7bfdad9cbe138913b4bfe72bd78p-1L,
+	0x1.fa7c1819e90d82e90a7e74b26360p-1L,
+	0x1.fd3c22b8f71f10975ba4b32bd006p-1L,
+	0x1.0000000000000000000000000000p+0L,
+	0x1.0163da9fb33356d84a66ae336e98p+0L,
+	0x1.02c9a3e778060ee6f7caca4f7a18p+0L,
+	0x1.04315e86e7f84bd738f9a20da442p+0L,
+	0x1.059b0d31585743ae7c548eb68c6ap+0L,
+	0x1.0706b29ddf6ddc6dc403a9d87b1ep+0L,
+	0x1.0874518759bc808c35f25d942856p+0L,
+	0x1.09e3ecac6f3834521e060c584d5cp+0L,
+	0x1.0b5586cf9890f6298b92b7184200p+0L,
+	0x1.0cc922b7247f7407b705b893dbdep+0L,
+	0x1.0e3ec32d3d1a2020742e4f8af794p+0L,
+	0x1.0fb66affed31af232091dd8a169ep+0L,
+	0x1.11301d0125b50a4ebbf1aed9321cp+0L,
+	0x1.12abdc06c31cbfb92bad324d6f84p+0L,
+	0x1.1429aaea92ddfb34101943b2588ep+0L,
+	0x1.15a98c8a58e512480d573dd562aep+0L,
+	0x1.172b83c7d517adcdf7c8c50eb162p+0L,
+	0x1.18af9388c8de9bbbf70b9a3c269cp+0L,
+	0x1.1a35beb6fcb753cb698f692d2038p+0L,
+	0x1.1bbe084045cd39ab1e72b442810ep+0L,
+	0x1.1d4873168b9aa7805b8028990be8p+0L,
+	0x1.1ed5022fcd91cb8819ff61121fbep+0L,
+	0x1.2063b88628cd63b8eeb0295093f6p+0L,
+	0x1.21f49917ddc962552fd29294bc20p+0L,
+	0x1.2387a6e75623866c1fadb1c159c0p+0L,
+	0x1.251ce4fb2a63f3582ab7de9e9562p+0L,
+	0x1.26b4565e27cdd257a673281d3068p+0L,
+	0x1.284dfe1f5638096cf15cf03c9fa0p+0L,
+	0x1.29e9df51fdee12c25d15f5a25022p+0L,
+	0x1.2b87fd0dad98ffddea46538fca24p+0L,
+	0x1.2d285a6e4030b40091d536d0733ep+0L,
+	0x1.2ecafa93e2f5611ca0f45d5239a4p+0L,
+	0x1.306fe0a31b7152de8d5a463063bep+0L,
+	0x1.32170fc4cd8313539cf1c3009330p+0L,
+	0x1.33c08b26416ff4c9c8610d96680ep+0L,
+	0x1.356c55f929ff0c94623476373be4p+0L,
+	0x1.371a7373aa9caa7145502f45452ap+0L,
+	0x1.38cae6d05d86585a9cb0d9bed530p+0L,
+	0x1.3a7db34e59ff6ea1bc9299e0a1fep+0L,
+	0x1.3c32dc313a8e484001f228b58cf0p+0L,
+	0x1.3dea64c12342235b41223e13d7eep+0L,
+	0x1.3fa4504ac801ba0bf701aa417b9cp+0L,
+	0x1.4160a21f72e29f84325b8f3dbacap+0L,
+	0x1.431f5d950a896dc704439410b628p+0L,
+	0x1.44e086061892d03136f409df0724p+0L,
+	0x1.46a41ed1d005772512f459229f0ap+0L,
+	0x1.486a2b5c13cd013c1a3b69062f26p+0L,
+	0x1.4a32af0d7d3de672d8bcf46f99b4p+0L,
+	0x1.4bfdad5362a271d4397afec42e36p+0L,
+	0x1.4dcb299fddd0d63b36ef1a9e19dep+0L,
+	0x1.4f9b2769d2ca6ad33d8b69aa0b8cp+0L,
+	0x1.516daa2cf6641c112f52c84d6066p+0L,
+	0x1.5342b569d4f81df0a83c49d86bf4p+0L,
+	0x1.551a4ca5d920ec52ec620243540cp+0L,
+	0x1.56f4736b527da66ecb004764e61ep+0L,
+	0x1.58d12d497c7fd252bc2b7343d554p+0L,
+	0x1.5ab07dd48542958c93015191e9a8p+0L,
+	0x1.5c9268a5946b701c4b1b81697ed4p+0L,
+	0x1.5e76f15ad21486e9be4c20399d12p+0L,
+	0x1.605e1b976dc08b076f592a487066p+0L,
+	0x1.6247eb03a5584b1f0fa06fd2d9eap+0L,
+	0x1.6434634ccc31fc76f8714c4ee122p+0L,
+	0x1.66238825522249127d9e29b92ea2p+0L,
+	0x1.68155d44ca973081c57227b9f69ep+0L,
+};
+
+static const float eps[TBLSIZE] = {
+	-0x1.5c50p-101,
+	-0x1.5d00p-106,
+	 0x1.8e90p-102,
+	-0x1.5340p-103,
+	 0x1.1bd0p-102,
+	-0x1.4600p-105,
+	-0x1.7a40p-104,
+	 0x1.d590p-102,
+	-0x1.d590p-101,
+	 0x1.b100p-103,
+	-0x1.0d80p-105,
+	 0x1.6b00p-103,
+	-0x1.9f00p-105,
+	 0x1.c400p-103,
+	 0x1.e120p-103,
+	-0x1.c100p-104,
+	-0x1.9d20p-103,
+	 0x1.a800p-108,
+	 0x1.4c00p-106,
+	-0x1.9500p-106,
+	 0x1.6900p-105,
+	-0x1.29d0p-100,
+	 0x1.4c60p-103,
+	 0x1.13a0p-102,
+	-0x1.5b60p-103,
+	-0x1.1c40p-103,
+	 0x1.db80p-102,
+	 0x1.91a0p-102,
+	 0x1.dc00p-105,
+	 0x1.44c0p-104,
+	 0x1.9710p-102,
+	 0x1.8760p-103,
+	-0x1.a720p-103,
+	 0x1.ed20p-103,
+	-0x1.49c0p-102,
+	-0x1.e000p-111,
+	 0x1.86a0p-103,
+	 0x1.2b40p-103,
+	-0x1.b400p-108,
+	 0x1.1280p-99,
+	-0x1.02d8p-102,
+	-0x1.e3d0p-103,
+	-0x1.b080p-105,
+	-0x1.f100p-107,
+	-0x1.16c0p-105,
+	-0x1.1190p-103,
+	-0x1.a7d2p-100,
+	 0x1.3450p-103,
+	-0x1.67c0p-105,
+	 0x1.4b80p-104,
+	-0x1.c4e0p-103,
+	 0x1.6000p-108,
+	-0x1.3f60p-105,
+	 0x1.93f0p-104,
+	 0x1.5fe0p-105,
+	 0x1.6f80p-107,
+	-0x1.7600p-106,
+	 0x1.21e0p-106,
+	-0x1.3a40p-106,
+	-0x1.40c0p-104,
+	-0x1.9860p-105,
+	-0x1.5d40p-108,
+	-0x1.1d70p-106,
+	 0x1.2760p-105,
+	 0x0.0000p+0,
+	 0x1.21e2p-104,
+	-0x1.9520p-108,
+	-0x1.5720p-106,
+	-0x1.4810p-106,
+	-0x1.be00p-109,
+	 0x1.0080p-105,
+	-0x1.5780p-108,
+	-0x1.d460p-105,
+	-0x1.6140p-105,
+	 0x1.4630p-104,
+	 0x1.ad50p-103,
+	 0x1.82e0p-105,
+	 0x1.1d3cp-101,
+	 0x1.6100p-107,
+	 0x1.ec30p-104,
+	 0x1.f200p-108,
+	 0x1.0b40p-103,
+	 0x1.3660p-102,
+	 0x1.d9d0p-103,
+	-0x1.02d0p-102,
+	 0x1.b070p-103,
+	 0x1.b9c0p-104,
+	-0x1.01c0p-103,
+	-0x1.dfe0p-103,
+	 0x1.1b60p-104,
+	-0x1.ae94p-101,
+	-0x1.3340p-104,
+	 0x1.b3d8p-102,
+	-0x1.6e40p-105,
+	-0x1.3670p-103,
+	 0x1.c140p-104,
+	 0x1.1840p-101,
+	 0x1.1ab0p-102,
+	-0x1.a400p-104,
+	 0x1.1f00p-104,
+	-0x1.7180p-103,
+	 0x1.4ce0p-102,
+	 0x1.9200p-107,
+	-0x1.54c0p-103,
+	 0x1.1b80p-105,
+	-0x1.1828p-101,
+	 0x1.5720p-102,
+	-0x1.a060p-100,
+	 0x1.9160p-102,
+	 0x1.a280p-104,
+	 0x1.3400p-107,
+	 0x1.2b20p-102,
+	 0x1.7800p-108,
+	 0x1.cfd0p-101,
+	 0x1.2ef0p-102,
+	-0x1.2760p-99,
+	 0x1.b380p-104,
+	 0x1.0048p-101,
+	-0x1.60b0p-102,
+	 0x1.a1ccp-100,
+	-0x1.a640p-104,
+	-0x1.08a0p-101,
+	 0x1.7e60p-102,
+	 0x1.22c0p-103,
+	-0x1.7200p-106,
+	 0x1.f0f0p-102,
+	 0x1.eb4ep-99,
+	 0x1.c6e0p-103,
+};
+
+/*
+ * exp2l(x): compute the base 2 exponential of x
+ *
+ * Accuracy: Peak error < 0.502 ulp.
+ *
+ * Method: (accurate tables)
+ *
+ *   Reduce x:
+ *     x = 2**k + y, for integer k and |y| <= 1/2.
+ *     Thus we have exp2(x) = 2**k * exp2(y).
+ *
+ *   Reduce y:
+ *     y = i/TBLSIZE + z - eps[i] for integer i near y * TBLSIZE.
+ *     Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z - eps[i]),
+ *     with |z - eps[i]| <= 2**-8 + 2**-98 for the table used.
+ *
+ *   We compute exp2(i/TBLSIZE) via table lookup and exp2(z - eps[i]) via
+ *   a degree-10 minimax polynomial with maximum error under 2**-120.
+ *   The values in exp2t[] and eps[] are chosen such that
+ *   exp2t[i] = exp2(i/TBLSIZE + eps[i]), and eps[i] is a small offset such
+ *   that exp2t[i] is accurate to 2**-122.
+ *
+ *   Note that the range of i is +-TBLSIZE/2, so we actually index the tables
+ *   by i0 = i + TBLSIZE/2.
+ *
+ *   This method is due to Gal, with many details due to Gal and Bachelis:
+ *
+ *	Gal, S. and Bachelis, B.  An Accurate Elementary Mathematical Library
+ *	for the IEEE Floating Point Standard.  TOMS 17(1), 26-46 (1991).
+ */
+long double
+exp2l(long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} u, v;
+	long double r, t, twopk, twopkp10000, z;
+	uint32_t es, hx, ix, i0;
+	int k;
+
+	u.e = x;
+
+	/* Filter out exceptional cases. */
+	hx = (u.bits.ext_sign << 15) | u.bits.ext_exp;
+	ix = hx & EXPMASK;
+	if (ix >= BIAS + 14) {		/* |x| >= 16384 */
+		if (ix == BIAS + LDBL_MAX_EXP) {
+			if (u.bits.ext_frach != 0
+			    || u.bits.ext_frachm != 0
+			    || u.bits.ext_fraclm != 0
+			    || u.bits.ext_fracl != 0
+			    || (hx & 0x8000) == 0)
+				return (x + x);	/* x is NaN or +Inf */
+			else
+				return (0.0);	/* x is -Inf */
+		}
+		if (x >= 16384)
+			return (huge * huge); /* overflow */
+		if (x <= -16495)
+			return (twom10000 * twom10000); /* underflow */
+	} else if (ix <= BIAS - 115) {		/* |x| < 0x1p-115 */
+		return (1.0 + x);
+	}
+
+	/*
+	 * Reduce x, computing z, i0, and k. The low bits of x + redux
+	 * contain the 16-bit integer part of the exponent (k) followed by
+	 * TBLBITS fractional bits (i0). We use bit tricks to extract these
+	 * as integers, then set z to the remainder.
+	 *
+	 * Example: Suppose x is 0xabc.123456p0 and TBLBITS is 8.
+	 * Then the low-order word of x + redux is 0x000abc12,
+	 * We split this into k = 0xabc and i0 = 0x12 (adjusted to
+	 * index into the table), then we compute z = 0x0.003456p0.
+	 *
+	 * XXX If the exponent is negative, the computation of k depends on
+	 *     '>>' doing sign extension.
+	 */
+	u.e = x + redux;
+	i0 = ((((uint64_t)u.bits.ext_fraclm << EXT_FRACLBITS)
+		| u.bits.ext_fracl) & 0xffffffff) + TBLSIZE / 2;
+	k = (int)i0 >> TBLBITS;
+	i0 = i0 & (TBLSIZE - 1);
+	u.e -= redux;
+	z = x - u.e;
+	v.bits.ext_frach = 0;
+	v.bits.ext_frachm = 0;
+	v.bits.ext_fraclm = 0;
+	v.bits.ext_fracl = 0;
+	if (k >= LDBL_MIN_EXP) {
+		es = LDBL_MAX_EXP - 1 + k;
+		v.bits.ext_exp = es & 0x7fffffff;
+		v.bits.ext_sign = u.bits.ext_sign >> 15;
+		twopk = v.e;
+	} else {
+		es = LDBL_MAX_EXP - 1 + k + 10000;
+		v.bits.ext_exp = es & 0x7fffffff;
+		v.bits.ext_sign = u.bits.ext_sign >> 15;
+		twopkp10000 = v.e;
+	}
+
+	/* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
+	t = tbl[i0];		/* exp2t[i0] */
+	z -= eps[i0];		/* eps[i0]   */
+	r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * (P5 + z * (P6
+	    + z * (P7 + z * (P8 + z * (P9 + z * P10)))))))));
+
+	/* Scale by 2**k. */
+	if(k >= LDBL_MIN_EXP) {
+		if (k == LDBL_MAX_EXP)
+			return (r * 2.0 * 0x1p16383L);
+		return (r * twopk);
+	} else {
+		return (r * twopkp10000 * twom10000);
+	}
+}
diff --git a/lib/libm/src/ld128/s_nanl.c b/lib/libm/src/ld128/s_nanl.c
new file mode 100644
index 00000000000..1edcac513e7
--- /dev/null
+++ b/lib/libm/src/ld128/s_nanl.c
@@ -0,0 +1,50 @@
+/*	$OpenBSD: s_nanl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*-
+ * Copyright (c) 2007 David Schultz
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ *
+ * $FreeBSD: src/lib/msun/ld128/s_nanl.c,v 1.3 2008/03/02 20:16:55 das Exp $
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <math.h>
+
+#include "math_private.h"
+
+long double
+nanl(const char *s)
+{
+	union {
+		long double e;
+		struct ieee_ext ieee;
+		uint32_t bits[4];
+	} u;
+
+	_scan_nan(u.bits, 4, s);
+	u.ieee.ext_exp = 0x7fff;
+	u.ieee.ext_frach |= 1U << 15;	/* make it a quiet NaN */
+
+	return (u.e);
+}