- 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@

83 files changed, 4004 insertions, 304 deletions

diff --git a/lib/libm/Makefile b/lib/libm/Makefile
index 90e669abc91..453d53147e5 100644
--- a/lib/libm/Makefile
+++ b/lib/libm/Makefile
@@ -1,5 +1,5 @@
+#  $OpenBSD: Makefile,v 1.59 2008/12/09 20:00:35 martynas Exp $
 #  $NetBSD: Makefile,v 1.28 1995/11/20 22:06:19 jtc Exp $
-#  $OpenBSD: Makefile,v 1.58 2008/10/07 22:25:53 martynas Exp $
 #
 #  @(#)Makefile 5.1beta 93/09/24 
 # 
@@ -26,8 +26,9 @@ ARCH_SRCS = s_copysign.S s_copysignf.S
 .PATH:	${.CURDIR}/arch/i387
 ARCH_SRCS = e_acos.S e_asin.S e_atan2.S e_exp.S e_fmod.S e_log.S e_log10.S \
 	    e_remainder.S e_remainderf.S e_scalb.S e_sqrt.S e_sqrtf.S \
+	    invtrig.c \
 	    s_atan.S s_atanf.S s_ceil.S s_ceilf.S s_copysign.S s_copysignf.S \
-	    s_cos.S s_cosf.S s_finite.S s_finitef.S s_floor.S s_floorf.S \
+	    s_cos.S s_cosf.S s_floor.S s_floorf.S \
 	    s_ilogb.S s_ilogbf.S s_log1p.S s_log1pf.S s_logb.S s_logbf.S \
 	    s_llrint.S s_llrintf.S s_lrint.S s_lrintf.S s_rint.S s_rintf.S\
 	    s_scalbn.S s_scalbnf.S s_significand.S s_significandf.S \
@@ -37,8 +38,9 @@ ARCH_SRCS = e_acos.S e_asin.S e_atan2.S e_exp.S e_fmod.S e_log.S e_log10.S \
 CPPFLAGS+=-I${.CURDIR}/arch/amd64
 ARCH_SRCS = e_acos.S e_asin.S e_atan2.S e_exp.S e_fmod.S e_log.S e_log10.S \
 	    e_remainder.S e_remainderf.S e_scalb.S e_sqrt.S e_sqrtf.S \
+	    invtrig.c \
 	    s_atan.S s_atanf.S s_ceil.S s_ceilf.S s_copysign.S s_copysignf.S \
-	    s_cos.S s_cosf.S s_finite.S s_finitef.S s_floor.S s_floorf.S \
+	    s_cos.S s_cosf.S s_floor.S s_floorf.S \
 	    s_ilogb.S s_ilogbf.S s_log1p.S s_log1pf.S s_logb.S s_logbf.S \
 	    s_rint.S s_rintf.S s_scalbn.S s_scalbnf.S s_significand.S \
 	    s_significandf.S s_sin.S s_sinf.S s_tan.S s_tanf.S
@@ -46,7 +48,7 @@ ARCH_SRCS = e_acos.S e_asin.S e_atan2.S e_exp.S e_fmod.S e_log.S e_log10.S \
 .PATH:	${.CURDIR}/arch/mc68881
 ARCH_SRCS = e_acos.S e_asin.S e_atanh.S e_cosh.S e_exp.S e_log.S e_log10.S \
 	    e_remainder.S e_scalb.S e_sinh.S e_sqrt.S s_atan.S s_ceil.S \
-	    s_copysign.S s_cos.S s_expm1.S s_finite.S s_floor.S s_log1p.S \
+	    s_copysign.S s_cos.S s_expm1.S s_floor.S s_log1p.S \
 	    s_logb.S s_rint.S s_scalbn.S s_sin.S s_tan.S s_tanh.S
 .elif (${MACHINE_ARCH} == "hppa")
 .PATH:	${.CURDIR}/arch/hppa
@@ -86,8 +88,7 @@ COMMON_SRCS = b_exp__D.c b_log__D.c b_tgamma.c \
 	s_cproj.c s_cprojf.c s_creal.c s_crealf.c s_csin.c s_csinf.c s_csinh.c \
 	s_csinhf.c s_csqrt.c s_csqrtf.c s_ctan.c s_ctanf.c s_ctanh.c \
 	s_ctanhf.c s_erf.c s_erff.c s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c \
-	s_fabsf.c s_fdim.c s_fmax.c s_fmaxf.c s_fmin.c s_fminf.c s_finite.c \
-	s_finitef.c \
+	s_fabsf.c s_fdim.c s_fmax.c s_fmaxf.c s_fmin.c s_fminf.c \
 	s_floor.c s_floorf.c s_frexpf.c s_ilogb.c s_ilogbf.c \
 	s_ldexpf.c s_log1p.c \
 	s_log1pf.c s_logb.c s_logbf.c s_llrint.c s_llrintf.c s_lrint.c \
@@ -99,6 +100,13 @@ COMMON_SRCS = b_exp__D.c b_log__D.c b_tgamma.c \
 	s_truncf.c w_drem.c w_dremf.c w_gamma.c w_gamma_r.c w_gammaf.c \
 	w_gammaf_r.c w_lgamma.c w_lgammaf.c
 
+LONG_SRCS = e_acosl.c e_asinl.c e_atan2l.c e_sqrtl.c \
+	invtrig.c \
+	k_cosl.c k_sinl.c k_tanl.c \
+	s_atanl.c s_copysignl.c s_cosl.c s_exp2l.c s_fabsl.c s_fmaxl.c \
+	s_fminl.c s_frexpl.c s_ilogbl.c s_logbl.c s_nanl.c s_rintl.c \
+	s_scalbnl.c s_sinl.c s_tanl.c
+
 # math routines for non-IEEE architectures.
 NOIEEE_SRCS = n_acosh.c n_argred.c n_asincos.c n_asinh.c n_atan.c \
 	n_atan2.c n_atanh.c n_cabs.c n_cacos.c n_cacosh.c n_carg.c \
@@ -120,6 +128,18 @@ SRCS=	${NOIEEE_SRCS} ${NOIEEE_ARCH}
 MAN+=	infnan.3
 .else
 SRCS=	${COMMON_SRCS}
+.if (${MACHINE_ARCH} == "amd64") || (${MACHINE_ARCH} == "i386") || \
+	(${MACHINE_ARCH} == "m68k") || (${MACHINE_ARCH} == "m88k")
+.PATH:		${.CURDIR}/src/ld80
+CPPFLAGS+=	-I${.CURDIR}/src -I${.CURDIR}/src/ld80
+SRCS+=		${LONG_SRCS}
+.endif
+.if (${MACHINE_ARCH} == "hppa64") || (${MACHINE_ARCH} == "mips64") || \
+	(${MACHINE_ARCH} == "sparc64")
+.PATH:		${.CURDIR}/src/ld128
+CPPFLAGS+=	-I${.CURDIR}/src -I${.CURDIR}/src/ld128
+SRCS+=		${LONG_SRCS}
+.endif
 .endif
 
 # Substitute common sources with any arch specific sources 
@@ -231,5 +251,4 @@ MLINKS+=trunc.3 truncf.3
 #	cpp \
 #	/usr/src/lib/libm/arch/mc68881/e_acos.S |  as -o e_acos.o
 
-
 .include <bsd.lib.mk>
diff --git a/lib/libm/arch/alpha/s_copysign.S b/lib/libm/arch/alpha/s_copysign.S
index 5d9797a09c3..aebab3e3eb6 100644
--- a/lib/libm/arch/alpha/s_copysign.S
+++ b/lib/libm/arch/alpha/s_copysign.S
@@ -1,4 +1,4 @@
-/*	$OpenBSD: s_copysign.S,v 1.2 2008/06/26 05:42:05 ray Exp $	*/
+/*	$OpenBSD: s_copysign.S,v 1.3 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: s_copysign.S,v 1.4 1999/07/02 15:37:34 simonb Exp $	*/
 
 /*-
@@ -32,6 +32,9 @@
 
 #include <machine/asm.h>
 
+.weak copysignl
+	copysignl = copysign
+
 LEAF(copysign, 2)
 	cpys	fa1, fa0, fv0
 	RET
diff --git a/lib/libm/arch/amd64/invtrig.c b/lib/libm/arch/amd64/invtrig.c
new file mode 100644
index 00000000000..dc32f953c2d
--- /dev/null
+++ b/lib/libm/arch/amd64/invtrig.c
@@ -0,0 +1,84 @@
+/*	$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 <stdint.h>
+
+#define STRUCT_DECLS
+#include "invtrig.h"
+
+/*
+ * asinl() and acosl()
+ */
+const LONGDOUBLE
+pS0 = { 0xaaaaaaaaaaaaaaa8ULL, 0x3ffcU }, /*  1.66666666666666666631e-01L */
+pS1 = { 0xd5271b6699b48bfaULL, 0xbffdU }, /* -4.16313987993683104320e-01L */
+pS2 = { 0xbcf67ca9e9f669cfULL, 0x3ffdU }, /*  3.69068046323246813704e-01L */
+pS3 = { 0x8b7baa3d15f9830dULL, 0xbffcU }, /* -1.36213932016738603108e-01L */
+pS4 = { 0x92154b093a3bff1cULL, 0x3ff9U }, /*  1.78324189708471965733e-02L */
+pS5 = { 0xe5dd76401964508cULL, 0xbff2U }, /* -2.19216428382605211588e-04L */
+pS6 = { 0xee69c5b0fdb76951ULL, 0xbfedU }, /* -7.10526623669075243183e-06L */
+qS1 = { 0xbcaa2159c01436a0ULL, 0xc000U }, /* -2.94788392796209867269e+00L */
+qS2 = { 0xd17a73d1e1564c29ULL, 0x4000U }, /*  3.27309890266528636716e+00L */
+qS3 = { 0xd767e411c9cf4c2cULL, 0xbfffU }, /* -1.68285799854822427013e+00L */
+qS4 = { 0xc809c0dfb9b0d0b7ULL, 0x3ffdU }, /*  3.90699412641738801874e-01L */
+qS5 = { 0x80c3a2197c8ced57ULL, 0xbffaU }; /* -3.14365703596053263322e-02L */
+
+/*
+ * atanl()
+ */
+const LONGDOUBLE atanhi[] = {
+	{ 0xed63382b0dda7b45ULL, 0x3ffdU }, /*  4.63647609000806116202e-01L */
+	{ 0xc90fdaa22168c235ULL, 0x3ffeU }, /*  7.85398163397448309628e-01L */
+	{ 0xfb985e940fb4d900ULL, 0x3ffeU }, /*  9.82793723247329067960e-01L */
+	{ 0xc90fdaa22168c235ULL, 0x3fffU }, /*  1.57079632679489661926e+00L */
+};
+
+const LONGDOUBLE atanlo[] = {
+	{ 0xdfc88bd978751a07ULL, 0x3fbcU }, /*  1.18469937025062860669e-20L */
+	{ 0xece675d1fc8f8cbbULL, 0xbfbcU }, /* -1.25413940316708300586e-20L */
+	{ 0xf10f5e197793c283ULL, 0x3fbdU }, /*  2.55232234165405176172e-20L */
+	{ 0xece675d1fc8f8cbbULL, 0xbfbdU }, /* -2.50827880633416601173e-20L */
+};
+
+const LONGDOUBLE aT[] = {
+	{ 0xaaaaaaaaaaaaaa9fULL, 0x3ffdU }, /*  3.33333333333333333017e-01L */
+	{ 0xcccccccccccc62bcULL, 0xbffcU }, /* -1.99999999999999632011e-01L */
+	{ 0x9249249248b81e3fULL, 0x3ffcU }, /*  1.42857142857046531280e-01L */
+	{ 0xe38e38e3316f3de5ULL, 0xbffbU }, /* -1.11111111100562372733e-01L */
+	{ 0xba2e8b8dc280726aULL, 0x3ffbU }, /*  9.09090902935647302252e-02L */
+	{ 0x9d89d5b4c6847ec4ULL, 0xbffbU }, /* -7.69230552476207730353e-02L */
+	{ 0x8888461d3099c677ULL, 0x3ffbU }, /*  6.66661718042406260546e-02L */
+	{ 0xf0e8ee0f5328dc29ULL, 0xbffaU }, /* -5.88158892835030888692e-02L */
+	{ 0xd73ea84d24bae54aULL, 0x3ffaU }, /*  5.25499891539726639379e-02L */
+	{ 0xc08fa381dcd9213aULL, 0xbffaU }, /* -4.70119845393155721494e-02L */
+	{ 0xa54a26f4095f2a3aULL, 0x3ffaU }, /*  4.03539201366454414072e-02L */
+	{ 0xeea2d8d059ef3ad6ULL, 0xbff9U }, /* -2.91303858419364158725e-02L */
+	{ 0xcc82292ab894b051ULL, 0x3ff8U }, /*  1.24822046299269234080e-02L */
+};
+
+const LONGDOUBLE
+pi_lo = { 0xece675d1fc8f8cbbULL, 0xbfbeU }; /* -5.01655761266833202345e-20L */
diff --git a/lib/libm/arch/amd64/s_finite.S b/lib/libm/arch/amd64/s_finite.S
deleted file mode 100644
index 662e5aa4ba0..00000000000
--- a/lib/libm/arch/amd64/s_finite.S
+++ /dev/null
@@ -1,25 +0,0 @@
-/*	$OpenBSD: s_finite.S,v 1.2 2005/08/02 11:17:31 espie Exp $ */
-/*
- * Written by J.T. Conklin <jtc@NetBSD.org>.
- * Public domain.
- */
-
-#include <machine/asm.h>
-
-ENTRY(finite)
-#ifdef __i386__
-	movl	8(%esp),%eax
-	andl	$0x7ff00000, %eax
-	cmpl	$0x7ff00000, %eax
-	setne	%al
-	andl	$0x000000ff, %eax
-#else
-	xorl	%eax,%eax
-	movq	$0x7ff0000000000000,%rsi
-	movq	%rsi,%rdi
-	movsd	%xmm0,-8(%rsp)
-	andq	-8(%rsp),%rsi
-	cmpq	%rdi,%rsi
-	setne	%al
-#endif
-	ret
diff --git a/lib/libm/arch/amd64/s_finitef.S b/lib/libm/arch/amd64/s_finitef.S
deleted file mode 100644
index 2a30c540e98..00000000000
--- a/lib/libm/arch/amd64/s_finitef.S
+++ /dev/null
@@ -1,24 +0,0 @@
-/*	$OpenBSD: s_finitef.S,v 1.2 2005/08/02 11:17:31 espie Exp $ */
-/*
- * Written by J.T. Conklin <jtc@NetBSD.org>.
- * Public domain.
- */
-
-#include <machine/asm.h>
-
-ENTRY(finitef)
-#ifdef __i386__
-	movl	4(%esp),%eax
-	andl	$0x7f800000, %eax
-	cmpl	$0x7f800000, %eax
-	setne	%al
-	andl	$0x000000ff, %eax
-#else
-	xorl	%eax,%eax
-	movl	$0x7ff00000,%esi
-	movss	%xmm0,-4(%rsp)
-	andl	-4(%rsp),%esi
-	cmpl	$0x7ff00000,%esi
-	setne	%al
-#endif
-	ret
diff --git a/lib/libm/arch/hppa/e_sqrt.c b/lib/libm/arch/hppa/e_sqrt.c
index 5c2663a6934..a8622877b0d 100644
--- a/lib/libm/arch/hppa/e_sqrt.c
+++ b/lib/libm/arch/hppa/e_sqrt.c
@@ -3,10 +3,12 @@
  */
 
 #if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$OpenBSD: e_sqrt.c,v 1.2 2008/09/07 20:36:08 martynas Exp $";
+static char rcsid[] = "$OpenBSD: e_sqrt.c,v 1.3 2008/12/09 20:00:35 martynas Exp $";
 #endif
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
 
 double
 sqrt(double x)
@@ -14,3 +16,9 @@ sqrt(double x)
 	__asm__ __volatile__ ("fsqrt,dbl %0, %0" : "+f" (x));
 	return (x);
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(sqrtl, sqrt);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/arch/hppa/s_rint.c b/lib/libm/arch/hppa/s_rint.c
index e25f7291e6b..1f06be47f36 100644
--- a/lib/libm/arch/hppa/s_rint.c
+++ b/lib/libm/arch/hppa/s_rint.c
@@ -3,10 +3,12 @@
  */
 
 #if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$OpenBSD: s_rint.c,v 1.2 2002/09/11 15:16:52 mickey Exp $";
+static char rcsid[] = "$OpenBSD: s_rint.c,v 1.3 2008/12/09 20:00:35 martynas Exp $";
 #endif
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
 
 double
 rint(double x)
@@ -15,3 +17,9 @@ rint(double x)
 
 	return (x);
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias   
+__weak_alias(rintl, rint);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/arch/i387/invtrig.c b/lib/libm/arch/i387/invtrig.c
new file mode 100644
index 00000000000..dc32f953c2d
--- /dev/null
+++ b/lib/libm/arch/i387/invtrig.c
@@ -0,0 +1,84 @@
+/*	$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 <stdint.h>
+
+#define STRUCT_DECLS
+#include "invtrig.h"
+
+/*
+ * asinl() and acosl()
+ */
+const LONGDOUBLE
+pS0 = { 0xaaaaaaaaaaaaaaa8ULL, 0x3ffcU }, /*  1.66666666666666666631e-01L */
+pS1 = { 0xd5271b6699b48bfaULL, 0xbffdU }, /* -4.16313987993683104320e-01L */
+pS2 = { 0xbcf67ca9e9f669cfULL, 0x3ffdU }, /*  3.69068046323246813704e-01L */
+pS3 = { 0x8b7baa3d15f9830dULL, 0xbffcU }, /* -1.36213932016738603108e-01L */
+pS4 = { 0x92154b093a3bff1cULL, 0x3ff9U }, /*  1.78324189708471965733e-02L */
+pS5 = { 0xe5dd76401964508cULL, 0xbff2U }, /* -2.19216428382605211588e-04L */
+pS6 = { 0xee69c5b0fdb76951ULL, 0xbfedU }, /* -7.10526623669075243183e-06L */
+qS1 = { 0xbcaa2159c01436a0ULL, 0xc000U }, /* -2.94788392796209867269e+00L */
+qS2 = { 0xd17a73d1e1564c29ULL, 0x4000U }, /*  3.27309890266528636716e+00L */
+qS3 = { 0xd767e411c9cf4c2cULL, 0xbfffU }, /* -1.68285799854822427013e+00L */
+qS4 = { 0xc809c0dfb9b0d0b7ULL, 0x3ffdU }, /*  3.90699412641738801874e-01L */
+qS5 = { 0x80c3a2197c8ced57ULL, 0xbffaU }; /* -3.14365703596053263322e-02L */
+
+/*
+ * atanl()
+ */
+const LONGDOUBLE atanhi[] = {
+	{ 0xed63382b0dda7b45ULL, 0x3ffdU }, /*  4.63647609000806116202e-01L */
+	{ 0xc90fdaa22168c235ULL, 0x3ffeU }, /*  7.85398163397448309628e-01L */
+	{ 0xfb985e940fb4d900ULL, 0x3ffeU }, /*  9.82793723247329067960e-01L */
+	{ 0xc90fdaa22168c235ULL, 0x3fffU }, /*  1.57079632679489661926e+00L */
+};
+
+const LONGDOUBLE atanlo[] = {
+	{ 0xdfc88bd978751a07ULL, 0x3fbcU }, /*  1.18469937025062860669e-20L */
+	{ 0xece675d1fc8f8cbbULL, 0xbfbcU }, /* -1.25413940316708300586e-20L */
+	{ 0xf10f5e197793c283ULL, 0x3fbdU }, /*  2.55232234165405176172e-20L */
+	{ 0xece675d1fc8f8cbbULL, 0xbfbdU }, /* -2.50827880633416601173e-20L */
+};
+
+const LONGDOUBLE aT[] = {
+	{ 0xaaaaaaaaaaaaaa9fULL, 0x3ffdU }, /*  3.33333333333333333017e-01L */
+	{ 0xcccccccccccc62bcULL, 0xbffcU }, /* -1.99999999999999632011e-01L */
+	{ 0x9249249248b81e3fULL, 0x3ffcU }, /*  1.42857142857046531280e-01L */
+	{ 0xe38e38e3316f3de5ULL, 0xbffbU }, /* -1.11111111100562372733e-01L */
+	{ 0xba2e8b8dc280726aULL, 0x3ffbU }, /*  9.09090902935647302252e-02L */
+	{ 0x9d89d5b4c6847ec4ULL, 0xbffbU }, /* -7.69230552476207730353e-02L */
+	{ 0x8888461d3099c677ULL, 0x3ffbU }, /*  6.66661718042406260546e-02L */
+	{ 0xf0e8ee0f5328dc29ULL, 0xbffaU }, /* -5.88158892835030888692e-02L */
+	{ 0xd73ea84d24bae54aULL, 0x3ffaU }, /*  5.25499891539726639379e-02L */
+	{ 0xc08fa381dcd9213aULL, 0xbffaU }, /* -4.70119845393155721494e-02L */
+	{ 0xa54a26f4095f2a3aULL, 0x3ffaU }, /*  4.03539201366454414072e-02L */
+	{ 0xeea2d8d059ef3ad6ULL, 0xbff9U }, /* -2.91303858419364158725e-02L */
+	{ 0xcc82292ab894b051ULL, 0x3ff8U }, /*  1.24822046299269234080e-02L */
+};
+
+const LONGDOUBLE
+pi_lo = { 0xece675d1fc8f8cbbULL, 0xbfbeU }; /* -5.01655761266833202345e-20L */
diff --git a/lib/libm/arch/i387/s_finite.S b/lib/libm/arch/i387/s_finite.S
deleted file mode 100644
index a6ee62c9295..00000000000
--- a/lib/libm/arch/i387/s_finite.S
+++ /dev/null
@@ -1,15 +0,0 @@
-/*	$OpenBSD: s_finite.S,v 1.3 2005/08/02 11:17:31 espie Exp $ */
-/*
- * Written by J.T. Conklin <jtc@netbsd.org>.
- * Public domain.
- */
-
-#include <machine/asm.h>
-
-ENTRY(finite)
-	movl	8(%esp),%eax
-	andl	$0x7ff00000, %eax
-	cmpl	$0x7ff00000, %eax
-	setne	%al
-	andl	$0x000000ff, %eax
-	ret
diff --git a/lib/libm/arch/i387/s_finitef.S b/lib/libm/arch/i387/s_finitef.S
deleted file mode 100644
index 111a148611c..00000000000
--- a/lib/libm/arch/i387/s_finitef.S
+++ /dev/null
@@ -1,15 +0,0 @@
-/*	$OpenBSD: s_finitef.S,v 1.3 2005/08/02 11:17:31 espie Exp $ */
-/*
- * Written by J.T. Conklin <jtc@netbsd.org>.
- * Public domain.
- */
-
-#include <machine/asm.h>
-
-ENTRY(finitef)
-	movl	4(%esp),%eax
-	andl	$0x7f800000, %eax
-	cmpl	$0x7f800000, %eax
-	setne	%al
-	andl	$0x000000ff, %eax
-	ret
diff --git a/lib/libm/arch/vax/n_atan2.S b/lib/libm/arch/vax/n_atan2.S
index 7892e352a3b..2cf39beeb90 100644
--- a/lib/libm/arch/vax/n_atan2.S
+++ b/lib/libm/arch/vax/n_atan2.S
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_atan2.S,v 1.3 2008/05/21 20:37:10 miod Exp $	*/
+/*	$OpenBSD: n_atan2.S,v 1.4 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_atan2.S,v 1.1 1995/10/10 23:40:25 ragge Exp $	*/
 /*
  * Copyright (c) 1985, 1993
@@ -73,6 +73,9 @@
  *	atan2(y,x) returns the exact ARG(x+iy) nearly rounded. 
  */
 
+.weak atan2l
+	atan2l = atan2
+
 ENTRY(atan2, R2|R3|R4|R5|R6|R7|R8|R9|R10|R11)
 	movq	4(ap),r2		# r2 = y
 	movq	12(ap),r4		# r4 = x
diff --git a/lib/libm/arch/vax/n_sincos.S b/lib/libm/arch/vax/n_sincos.S
index a9f0bcc35b8..e0487088765 100644
--- a/lib/libm/arch/vax/n_sincos.S
+++ b/lib/libm/arch/vax/n_sincos.S
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_sincos.S,v 1.4 2008/06/21 08:26:19 martynas Exp $	*/
+/*	$OpenBSD: n_sincos.S,v 1.5 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_sincos.S,v 1.1 1995/10/10 23:40:28 ragge Exp $	*/
 /*
  * Copyright (c) 1985, 1993
@@ -48,6 +48,9 @@
  * S. McDonald, April 4,  1985
  */
 
+.weak sinl
+	sinl = sin
+
 ENTRY(sin, R2|R3|R4|R5|R6|R7|R8|R9|R10|R11)
 	movq	4(ap),r0
 	bicw3	$0x807f,r0,r2
@@ -77,6 +80,9 @@ ENTRY(sin, R2|R3|R4|R5|R6|R7|R8|R9|R10|R11)
  * S. McDonald, April 4,  1985
  */
 
+.weak cosl
+	cosl = cos
+
 ENTRY(cos, R2|R3|R4|R5|R6|R7|R8|R9|R10|R11)
 	movq	4(ap),r0
 	bicw3	$0x7f,r0,r2
diff --git a/lib/libm/arch/vax/n_sqrt.S b/lib/libm/arch/vax/n_sqrt.S
index 5050c9fefc5..a1fd2ba61fb 100644
--- a/lib/libm/arch/vax/n_sqrt.S
+++ b/lib/libm/arch/vax/n_sqrt.S
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_sqrt.S,v 1.5 2008/09/16 22:13:12 martynas Exp $	*/
+/*	$OpenBSD: n_sqrt.S,v 1.6 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_sqrt.S,v 1.1 1995/10/10 23:40:29 ragge Exp $	*/
 /*
  * Copyright (c) 1985, 1993
@@ -44,6 +44,10 @@
  *
  * entry points: sqrt		double arg is on the stack
  */
+
+.weak sqrtl
+	sqrtl = sqrt
+
 ENTRY(sqrt, R2|R3|R4|R5)
 	movq    4(ap),r0
 dsqrt2:	bicw3	$0x807f,r0,r2	# check exponent of input
diff --git a/lib/libm/arch/vax/n_support.S b/lib/libm/arch/vax/n_support.S
index 2f06891b861..f2572586526 100644
--- a/lib/libm/arch/vax/n_support.S
+++ b/lib/libm/arch/vax/n_support.S
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_support.S,v 1.13 2008/11/20 23:21:37 martynas Exp $	*/
+/*	$OpenBSD: n_support.S,v 1.14 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_support.S,v 1.1 1995/10/10 23:40:30 ragge Exp $	*/
 /*
  * Copyright (c) 1985, 1993
@@ -40,7 +40,6 @@
  * logb(x),
  * logbf(x),
  * scalbn(x,N),
- * finite(x),
  * remainder(x,y),
  * Coded in vax assembly language by K.C. Ng,  3/14/85.
  * Revised by K.C. Ng on 4/9/85.
@@ -52,6 +51,9 @@
  * copysign(double x, double y)
  */
 
+.weak copysignl
+	copysignl = copysign
+
 ENTRY(copysign, R2)
 	movq	4(ap),r0		# load x into r0
 	bicw3	$0x807f,r0,r2		# mask off the exponent of x
@@ -80,6 +82,9 @@ Fz:	ret
  * logb(double x)
  */
 
+.weak logbl
+	logbl = logb
+
 ENTRY(logb, 0)
 	bicl3	$0xffff807f,4(ap),r0	# mask off the exponent of x
 	beql    Ln
@@ -110,24 +115,13 @@ Fn:	movl	4(ap),r0		# r0:1 = x (zero or reserved op)
 1:	ret
 
 /*
- * long
- * finite(double x)
- */
-
-ENTRY(finite, 0)
-	bicw3	$0x7f,4(ap),r0		# mask off the mantissa
-	cmpw	r0,$0x8000		# to see if x is the reserved op
-	beql	1f			# if so, return FALSE (0)
-	movl	$1,r0			# else return TRUE (1)
-	ret
-1:	clrl	r0
-	ret
-
-/*
  * double
  * scalbn(double x, int N)
  */
 
+.weak scalbnl
+	scalbnl = scalbn
+
 ENTRY(scalbn, R2|R3)
 	movq	4(ap),r0
 	bicl3	$0xffff807f,r0,r3
diff --git a/lib/libm/arch/vax/n_tan.S b/lib/libm/arch/vax/n_tan.S
index 0f499a943fc..13ad1554429 100644
--- a/lib/libm/arch/vax/n_tan.S
+++ b/lib/libm/arch/vax/n_tan.S
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_tan.S,v 1.4 2008/06/21 08:26:19 martynas Exp $	*/
+/*	$OpenBSD: n_tan.S,v 1.5 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_tan.S,v 1.1 1995/10/10 23:40:31 ragge Exp $	*/
 /*
  * Copyright (c) 1985, 1993
@@ -46,6 +46,9 @@
  * S. McDonald, April 4,  1985
  */
 
+.weak tanl
+	tanl = tan
+
 ENTRY(tan, R2|R3|R4|R5|R6|R7|R8|R9|R10|R11)
 	movq	4(ap),r0
 	bicw3	$0x807f,r0,r2
diff --git a/lib/libm/noieee_src/n_asincos.c b/lib/libm/noieee_src/n_asincos.c
index 40cc86bcabb..6b720021291 100644
--- a/lib/libm/noieee_src/n_asincos.c
+++ b/lib/libm/noieee_src/n_asincos.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_asincos.c,v 1.7 2008/06/21 08:26:19 martynas Exp $	*/
+/*	$OpenBSD: n_asincos.c,v 1.8 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_asincos.c,v 1.1 1995/10/10 23:36:34 ragge Exp $	*/
 /*
  * Copyright (c) 1985, 1993
@@ -85,7 +85,9 @@ static char sccsid[] = "@(#)asincos.c	8.1 (Berkeley) 6/4/93";
  *      1.99 ulps.
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <math.h>
+
 #include "mathimpl.h"
 
 double
@@ -104,6 +106,10 @@ asin(double x)
 
 }
 
+#ifdef __weak_alias
+__weak_alias(asinl, asin);
+#endif /* __weak_alias */
+
 /* ACOS(X)
  * RETURNS ARC COS OF X
  * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits)
@@ -170,3 +176,7 @@ acos(double x)
 	    t=atan2(one,0.0);	/* t = PI/2 */
 	return(t+t);
 }
+
+#ifdef __weak_alias
+__weak_alias(acosl, acos);
+#endif /* __weak_alias */
diff --git a/lib/libm/noieee_src/n_atan.c b/lib/libm/noieee_src/n_atan.c
index f2dd0a684bf..9a0df75b1d2 100644
--- a/lib/libm/noieee_src/n_atan.c
+++ b/lib/libm/noieee_src/n_atan.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_atan.c,v 1.5 2008/06/21 08:26:19 martynas Exp $	*/
+/*	$OpenBSD: n_atan.c,v 1.6 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_atan.c,v 1.1 1995/10/10 23:36:36 ragge Exp $	*/
 /*
  * Copyright (c) 1985, 1993
@@ -77,7 +77,8 @@ static char sccsid[] = "@(#)atan.c	8.1 (Berkeley) 6/4/93";
  *	0.85 ulps.
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <math.h>
 
 double
 atan(double x)
@@ -85,3 +86,7 @@ atan(double x)
 	double one=1.0;
 	return(atan2(x,one));
 }
+
+#ifdef __weak_alias
+__weak_alias(atanl, atan);
+#endif /* __weak_alias */
diff --git a/lib/libm/noieee_src/n_atan2.c b/lib/libm/noieee_src/n_atan2.c
index 873946aabf2..919e2c662da 100644
--- a/lib/libm/noieee_src/n_atan2.c
+++ b/lib/libm/noieee_src/n_atan2.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_atan2.c,v 1.9 2008/07/17 15:36:28 martynas Exp $	*/
+/*	$OpenBSD: n_atan2.c,v 1.10 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_atan2.c,v 1.1 1995/10/10 23:36:37 ragge Exp $	*/
 /*
  * Copyright (c) 1985, 1993
@@ -107,7 +107,9 @@ static char sccsid[] = "@(#)atan2.c	8.1 (Berkeley) 6/4/93";
  * shown.
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <math.h>
+
 #include "mathimpl.h"
 
 vc(athfhi, 4.6364760900080611433E-1  ,6338,3fed,da7b,2b0d,  -1, .ED63382B0DDA7B)
@@ -281,3 +283,7 @@ begin:
 
 	return(copysign((signx>zero)?z:PI-z,signy));
 }
+
+#ifdef __weak_alias
+__weak_alias(atan2l, atan2);
+#endif /* __weak_alias */
diff --git a/lib/libm/noieee_src/n_fdim.c b/lib/libm/noieee_src/n_fdim.c
index 09644101820..409816672d0 100644
--- a/lib/libm/noieee_src/n_fdim.c
+++ b/lib/libm/noieee_src/n_fdim.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_fdim.c,v 1.1 2008/09/11 19:19:34 martynas Exp $	*/
+/*	$OpenBSD: n_fdim.c,v 1.2 2008/12/09 20:00:35 martynas Exp $	*/
 /*-
  * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
  * All rights reserved.
@@ -25,6 +25,7 @@
  * SUCH DAMAGE.
  */
 
+#include <machine/cdefs.h>
 #include <math.h>
 
 #define	DECL(type, fn)			\
@@ -41,4 +42,4 @@ fn(type x, type y)			\
 
 DECL(double, fdim)
 DECL(float, fdimf)
-DECL(long double, fdiml)
+__weak_alias(fdiml, fdim);
diff --git a/lib/libm/noieee_src/n_floor.c b/lib/libm/noieee_src/n_floor.c
index 3e45d05e159..7a9f3d66399 100644
--- a/lib/libm/noieee_src/n_floor.c
+++ b/lib/libm/noieee_src/n_floor.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_floor.c,v 1.8 2008/06/25 17:49:31 martynas Exp $	*/
+/*	$OpenBSD: n_floor.c,v 1.9 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_floor.c,v 1.1 1995/10/10 23:36:48 ragge Exp $	*/
 /*
  * Copyright (c) 1985, 1993
@@ -33,7 +33,9 @@
 static char sccsid[] = "@(#)floor.c	8.1 (Berkeley) 6/4/93";
 #endif /* not lint */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <math.h>
+
 #include "mathimpl.h"
 
 vc(L, 4503599627370496.0E0 ,0000,5c00,0000,0000, 55, 1.0) /* 2**55 */
@@ -121,3 +123,7 @@ rint(double x)
 	t = x + s;				/* x+s rounded to integer */
 	return (t - s);
 }
+
+#ifdef __weak_alias    
+__weak_alias(rintl, rint);
+#endif /* __weak_alias */
diff --git a/lib/libm/noieee_src/n_fmax.c b/lib/libm/noieee_src/n_fmax.c
index 057e09faa18..f1f2c28d908 100644
--- a/lib/libm/noieee_src/n_fmax.c
+++ b/lib/libm/noieee_src/n_fmax.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_fmax.c,v 1.1 2008/09/11 19:19:34 martynas Exp $	*/
+/*	$OpenBSD: n_fmax.c,v 1.2 2008/12/09 20:00:35 martynas Exp $	*/
 /*-
  * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
  * All rights reserved.
@@ -25,6 +25,7 @@
  * SUCH DAMAGE.
  */
 
+#include <machine/cdefs.h>
 #include <math.h>
 
 double
@@ -45,3 +46,5 @@ fmax(double x, double y)
 
 	return (x > y ? x : y);
 }
+
+__weak_alias(fmaxl, fmax);
diff --git a/lib/libm/noieee_src/n_fmin.c b/lib/libm/noieee_src/n_fmin.c
index 3ecaf8d1da4..7929997da0b 100644
--- a/lib/libm/noieee_src/n_fmin.c
+++ b/lib/libm/noieee_src/n_fmin.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_fmin.c,v 1.1 2008/09/11 19:19:34 martynas Exp $	*/
+/*	$OpenBSD: n_fmin.c,v 1.2 2008/12/09 20:00:35 martynas Exp $	*/
 /*-
  * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
  * All rights reserved.
@@ -25,6 +25,7 @@
  * SUCH DAMAGE.
  */
 
+#include <machine/cdefs.h>
 #include <math.h>
 
 double
@@ -45,3 +46,7 @@ fmin(double x, double y)
 
 	return (x < y ? x : y);
 }
+
+#ifdef __weak_alias    
+__weak_alias(fminl, fmin);
+#endif /* __weak_alias */
diff --git a/lib/libm/noieee_src/n_sincos.c b/lib/libm/noieee_src/n_sincos.c
index 44627cd2bc5..2da035851ca 100644
--- a/lib/libm/noieee_src/n_sincos.c
+++ b/lib/libm/noieee_src/n_sincos.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_sincos.c,v 1.6 2008/06/25 17:49:31 martynas Exp $	*/
+/*	$OpenBSD: n_sincos.c,v 1.7 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_sincos.c,v 1.1 1995/10/10 23:37:04 ragge Exp $	*/
 /*
  * Copyright (c) 1987, 1993
@@ -33,7 +33,9 @@
 static char sccsid[] = "@(#)sincos.c	8.1 (Berkeley) 6/4/93";
 #endif /* not lint */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <math.h>
+
 #include "mathimpl.h"
 
 double
@@ -65,6 +67,10 @@ sin(double x)
 	return x+x*sin__S(x*x);
 }
 
+#ifdef __weak_alias
+__weak_alias(sinl, sin);
+#endif /* __weak_alias */
+
 double
 cos(double x)
 {
@@ -94,3 +100,7 @@ cos(double x)
 	a = (z >= thresh ? half-((z-half)-c) : one-(z-c));
 	return copysign(a,s);
 }
+
+#ifdef __weak_alias
+__weak_alias(cosl, cos);
+#endif /* __weak_alias */
diff --git a/lib/libm/noieee_src/n_support.c b/lib/libm/noieee_src/n_support.c
index 63c6706d588..36123fc055e 100644
--- a/lib/libm/noieee_src/n_support.c
+++ b/lib/libm/noieee_src/n_support.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_support.c,v 1.15 2008/07/24 09:40:16 martynas Exp $	*/
+/*	$OpenBSD: n_support.c,v 1.16 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_support.c,v 1.1 1995/10/10 23:37:06 ragge Exp $	*/
 /*
  * Copyright (c) 1985, 1993
@@ -62,16 +62,15 @@ static char sccsid[] = "@(#)support.c	8.1 (Berkeley) 6/4/93";
  *              returns the unbiased exponent of x, a signed integer in
  *              double precision, except that logb(0) is -INF, logb(INF)
  *              is +INF, and logb(NAN) is that NAN.
- * (d) finite(x)
- *              returns the value TRUE if -INF < x < +INF and returns
- *              FALSE otherwise.
  *
  *
  * CODED IN C BY K.C. NG, 11/25/84;
  * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <math.h>
+
 #include "mathimpl.h"
 
 #if defined(__vax__)      /* VAX D format */
@@ -122,6 +121,9 @@ scalbn(double x, int N)
         return(x);
 }
 
+#ifdef __weak_alias    
+__weak_alias(scalbnl, scalbn);
+#endif /* __weak_alias */
 
 double
 copysign(double x, double y)
@@ -137,6 +139,10 @@ copysign(double x, double y)
         return(x);
 }
 
+#ifdef __weak_alias    
+__weak_alias(copysignl, copysign);
+#endif /* __weak_alias */
+
 double
 logb(double x)
 {
@@ -160,15 +166,9 @@ logb(double x)
 #endif	/* defined(__vax__) */
 }
 
-int
-finite(double x)
-{
-#if defined(__vax__)
-        return(1);
-#else	/* defined(__vax__) */
-        return( (*((short *) &x ) & mexp ) != mexp );
-#endif	/* defined(__vax__) */
-}
+#ifdef __weak_alias    
+__weak_alias(logbl, logb);
+#endif /* __weak_alias */
 
 double
 remainder(double x, double p)
@@ -319,6 +319,10 @@ sqrt(double x)
 end:        return(scalbn(q,n));
 }
 
+#ifdef __weak_alias    
+__weak_alias(sqrtl, sqrt);
+#endif /* __weak_alias */
+
 #if 0
 /* REMAINDER(X,Y)
  * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
diff --git a/lib/libm/noieee_src/n_tan.c b/lib/libm/noieee_src/n_tan.c
index bd3d97b96ef..e344a84b013 100644
--- a/lib/libm/noieee_src/n_tan.c
+++ b/lib/libm/noieee_src/n_tan.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: n_tan.c,v 1.6 2008/06/25 17:49:31 martynas Exp $	*/
+/*	$OpenBSD: n_tan.c,v 1.7 2008/12/09 20:00:35 martynas Exp $	*/
 /*	$NetBSD: n_tan.c,v 1.1 1995/10/10 23:37:07 ragge Exp $	*/
 /*
  * Copyright (c) 1987, 1993
@@ -33,7 +33,9 @@
 static char sccsid[] = "@(#)tan.c	8.1 (Berkeley) 6/4/93";
 #endif /* not lint */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <math.h>
+
 #include "mathimpl.h"
 
 double
@@ -67,3 +69,7 @@ tan(double x)
 	else
 		return c/(x+x*ss);		/* ... cos/sin */
 }
+
+#ifdef __weak_alias
+__weak_alias(tanl, tan);
+#endif /* __weak_alias */
diff --git a/lib/libm/shlib_version b/lib/libm/shlib_version
index d9961ea9fef..3066b9771e7 100644
--- a/lib/libm/shlib_version
+++ b/lib/libm/shlib_version
@@ -1,2 +1,2 @@
-major=4
+major=5
 minor=0
diff --git a/lib/libm/src/e_acos.c b/lib/libm/src/e_acos.c
index 2ae24586cc1..6f459cbf536 100644
--- a/lib/libm/src/e_acos.c
+++ b/lib/libm/src/e_acos.c
@@ -38,7 +38,10 @@ static char rcsid[] = "$NetBSD: e_acos.c,v 1.9 1995/05/12 04:57:13 jtc Exp $";
  * Function needed: sqrt
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 static const double 
@@ -101,3 +104,9 @@ acos(double x)
 	    return 2.0*(df+w);
 	}
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(acosl, acos);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/e_acosl.c b/lib/libm/src/e_acosl.c
new file mode 100644
index 00000000000..7e00b504dc5
--- /dev/null
+++ b/lib/libm/src/e_acosl.c
@@ -0,0 +1,104 @@
+/*	$OpenBSD: e_acosl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/* @(#)e_acos.c 1.3 95/01/18 */
+/* FreeBSD: head/lib/msun/src/e_acos.c 176451 2008-02-22 02:30:36Z das */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * See comments in e_acos.c.
+ * Converted to long double by David Schultz <das@FreeBSD.ORG>.
+ * Adapted for OpenBSD by Martynas Venckus <martynas@openbsd.org>.
+ */
+
+#include <float.h>
+#include <math.h>
+
+#include "invtrig.h"
+#include "math_private.h"
+
+#ifdef EXT_IMPLICIT_NBIT
+#define	LDBL_NBIT	0
+#else /* EXT_IMPLICIT_NBIT */
+#define	LDBL_NBIT	0x80000000
+#endif /* EXT_IMPLICIT_NBIT */
+
+static const long double
+one=  1.00000000000000000000e+00;
+
+#if defined(__amd64__) || defined(__i386__)
+/* XXX Work around the fact that gcc truncates long double constants on i386 */
+static volatile double
+pi1 =  3.14159265358979311600e+00,	/*  0x1.921fb54442d18p+1  */
+pi2 =  1.22514845490862001043e-16;	/*  0x1.1a80000000000p-53 */
+#define	pi	((long double)pi1 + pi2)
+#else
+static const long double
+pi =  3.14159265358979323846264338327950280e+00L;
+#endif
+
+long double
+acosl(long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} u;
+	long double z,p,q,r,w,s,c,df;
+	int16_t expsign, expt;
+	u.e = x;
+	expsign = (u.bits.ext_sign << 15) | u.bits.ext_exp;
+	expt = expsign & 0x7fff;
+	if(expt >= BIAS) {	/* |x| >= 1 */
+	    if(expt==BIAS && ((u.bits.ext_frach&~LDBL_NBIT)
+#ifdef EXT_FRACHMBITS
+		| u.bits.ext_frachm
+#endif /* EXT_FRACHMBITS */
+#ifdef EXT_FRACLMBITS
+		| u.bits.ext_fraclm
+#endif /* EXT_FRACLMBITS */
+		| u.bits.ext_fracl)==0) {
+		if (expsign>0) return 0.0;	/* acos(1) = 0  */
+		else return pi+2.0*pio2_lo;	/* acos(-1)= pi */
+	    }
+	    return (x-x)/(x-x);		/* acos(|x|>1) is NaN */
+	}
+	if(expt<BIAS-1) {	/* |x| < 0.5 */
+	    if(expt<ACOS_CONST) return pio2_hi+pio2_lo;/*x tiny: acosl=pi/2*/
+	    z = x*x;
+	    p = P(z);
+	    q = Q(z);
+	    r = p/q;
+	    return pio2_hi - (x - (pio2_lo-x*r));
+	} else  if (expsign<0) {	/* x < -0.5 */
+	    z = (one+x)*0.5;
+	    p = P(z);
+	    q = Q(z);
+	    s = sqrtl(z);
+	    r = p/q;
+	    w = r*s-pio2_lo;
+	    return pi - 2.0*(s+w);
+	} else {			/* x > 0.5 */
+	    z = (one-x)*0.5;
+	    s = sqrtl(z);
+	    u.e = s;
+	    u.bits.ext_fracl = 0;
+#ifdef EXT_FRACLMBITS
+	    u.bits.ext_fraclm = 0;
+#endif /* EXT_FRACLMBITS */
+	    df = u.e;
+	    c  = (z-df*df)/(s+df);
+	    p = P(z);
+	    q = Q(z);
+	    r = p/q;
+	    w = r*s+c;
+	    return 2.0*(df+w);
+	}
+}
diff --git a/lib/libm/src/e_asin.c b/lib/libm/src/e_asin.c
index 4c732971e40..5348ea0e510 100644
--- a/lib/libm/src/e_asin.c
+++ b/lib/libm/src/e_asin.c
@@ -44,8 +44,10 @@ static char rcsid[] = "$NetBSD: e_asin.c,v 1.9 1995/05/12 04:57:22 jtc Exp $";
  *
  */
 
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
 
-#include "math.h"
 #include "math_private.h"
 
 static const double 
@@ -110,3 +112,9 @@ asin(double x)
 	}    
 	if(hx>0) return t; else return -t;    
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(asinl, asin);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/e_asinl.c b/lib/libm/src/e_asinl.c
new file mode 100644
index 00000000000..316eed55ff7
--- /dev/null
+++ b/lib/libm/src/e_asinl.c
@@ -0,0 +1,100 @@
+/*	$OpenBSD: e_asinl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/* @(#)e_asin.c 1.3 95/01/18 */
+/* FreeBSD: head/lib/msun/src/e_asin.c 176451 2008-02-22 02:30:36Z das */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * See comments in e_asin.c.
+ * Converted to long double by David Schultz <das@FreeBSD.ORG>.
+ * Adapted for OpenBSD by Martynas Venckus <martynas@openbsd.org>.
+ */
+
+#include <float.h>
+#include <math.h>
+
+#include "invtrig.h"
+#include "math_private.h"
+
+#ifdef EXT_IMPLICIT_NBIT
+#define	LDBL_NBIT	0
+#else /* EXT_IMPLICIT_NBIT */
+#define	LDBL_NBIT	0x80000000
+#endif /* EXT_IMPLICIT_NBIT */
+
+static const long double
+one =  1.00000000000000000000e+00,
+huge = 1.000e+300;
+
+long double
+asinl(long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} u;
+	long double t=0.0,w,p,q,c,r,s;
+	int16_t expsign, expt;
+	u.e = x;
+	expsign = (u.bits.ext_sign << 15) | u.bits.ext_exp;
+	expt = expsign & 0x7fff;
+	if(expt >= BIAS) {		/* |x|>= 1 */
+		if(expt==BIAS && ((u.bits.ext_frach&~LDBL_NBIT)
+#ifdef EXT_FRACHMBITS
+			| u.bits.ext_frachm
+#endif /* EXT_FRACHMBITS */
+#ifdef EXT_FRACLMBITS
+			| u.bits.ext_fraclm
+#endif /* EXT_FRACLMBITS */
+			| u.bits.ext_fracl)==0)
+		    /* asin(1)=+-pi/2 with inexact */
+		    return x*pio2_hi+x*pio2_lo;	
+	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */   
+	} else if (expt<BIAS-1) {	/* |x|<0.5 */
+	    if(expt<ASIN_LINEAR) {	/* if |x| is small, asinl(x)=x */
+		if(huge+x>one) return x;/* return x with inexact if x!=0*/
+	    }
+	    t = x*x;
+	    p = P(t);
+	    q = Q(t);
+	    w = p/q;
+	    return x+x*w;
+	}
+	/* 1> |x|>= 0.5 */
+	w = one-fabsl(x);
+	t = w*0.5;
+	p = P(t);
+	q = Q(t);
+	s = sqrtl(t);
+#ifdef EXT_FRACHMBITS
+	if((((uint64_t)u.bits.ext_frach << EXT_FRACHMBITS)
+		| u.bits.ext_frachm) >= THRESH) {
+						/* if |x| is close to 1 */
+#else /* EXT_FRACHMBITS */
+	if(u.bits.ext_frach>=THRESH) {		/* if |x| is close to 1 */
+#endif /* EXT_FRACHMBITS */
+	    w = p/q;
+	    t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
+	} else {
+	    u.e = s;
+	    u.bits.ext_fracl = 0;
+#ifdef EXT_FRACLMBITS
+	    u.bits.ext_fraclm = 0;
+#endif /* EXT_FRACLMBITS */
+	    w = u.e;
+	    c  = (t-w*w)/(s+w);
+	    r  = p/q;
+	    p  = 2.0*s*r-(pio2_lo-2.0*c);
+	    q  = pio4_hi-2.0*w;
+	    t  = pio4_hi-(p-q);
+	}    
+	if(expsign>0) return t; else return -t;    
+}
diff --git a/lib/libm/src/e_atan2.c b/lib/libm/src/e_atan2.c
index f8a0330dcb8..7dc66c9650e 100644
--- a/lib/libm/src/e_atan2.c
+++ b/lib/libm/src/e_atan2.c
@@ -41,7 +41,10 @@ static char rcsid[] = "$NetBSD: e_atan2.c,v 1.8 1995/05/10 20:44:51 jtc Exp $";
  * to produce the hexadecimal values shown.
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 static const double 
@@ -120,3 +123,9 @@ atan2(double y, double x)
 	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
 	}
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(atan2l, atan2);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/e_atan2l.c b/lib/libm/src/e_atan2l.c
new file mode 100644
index 00000000000..94bdfa8e2f3
--- /dev/null
+++ b/lib/libm/src/e_atan2l.c
@@ -0,0 +1,163 @@
+/*	$OpenBSD: e_atan2l.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/* @(#)e_atan2.c 1.3 95/01/18 */
+/* FreeBSD: head/lib/msun/src/e_atan2.c 176451 2008-02-22 02:30:36Z das */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ */
+
+/*
+ * See comments in e_atan2.c.
+ * Converted to long double by David Schultz <das@FreeBSD.ORG>.
+ * Adapted for OpenBSD by Martynas Venckus <martynas@openbsd.org>.
+ */
+
+#include <float.h>
+#include <math.h>
+
+#include "invtrig.h"
+#include "math.h"
+#include "math_private.h"
+
+#ifdef EXT_IMPLICIT_NBIT
+#define	LDBL_NBIT	0
+#else /* EXT_IMPLICIT_NBIT */
+#define	LDBL_NBIT	0x80000000
+#endif /* EXT_IMPLICIT_NBIT */
+
+static volatile long double
+tiny  = 1.0e-300;
+static const long double
+zero  = 0.0;
+
+#if defined(__amd64__) || defined(__i386__)
+/* XXX Work around the fact that gcc truncates long double constants on i386 */
+static volatile double
+pi1 =  3.14159265358979311600e+00,	/*  0x1.921fb54442d18p+1  */
+pi2 =  1.22514845490862001043e-16;	/*  0x1.1a80000000000p-53 */
+#define	pi	((long double)pi1 + pi2)
+#else
+static const long double
+pi =  3.14159265358979323846264338327950280e+00L;
+#endif
+
+long double
+atan2l(long double y, long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} ux, uy;
+	long double z;
+	int32_t k,m;
+	int16_t exptx, expsignx, expty, expsigny;
+
+	uy.e = y;
+	expsigny = (uy.bits.ext_sign << 15) | uy.bits.ext_exp;
+	expty = expsigny & 0x7fff;
+	ux.e = x;
+	expsignx = (ux.bits.ext_sign << 15) | ux.bits.ext_exp;
+	exptx = expsignx & 0x7fff;
+
+	if ((exptx==BIAS+LDBL_MAX_EXP &&
+	     ((ux.bits.ext_frach&~LDBL_NBIT)
+#ifdef EXT_FRACHMBITS
+		| ux.bits.ext_frachm
+#endif /* EXT_FRACHMBITS */
+#ifdef EXT_FRACLMBITS
+		| ux.bits.ext_fraclm
+#endif /* EXT_FRACLMBITS */
+		| ux.bits.ext_fracl)!=0) ||	/* x is NaN */
+	    (expty==BIAS+LDBL_MAX_EXP &&
+	     ((uy.bits.ext_frach&~LDBL_NBIT)
+#ifdef EXT_FRACHMBITS
+		| uy.bits.ext_frachm
+#endif /* EXT_FRACHMBITS */
+#ifdef EXT_FRACLMBITS
+		| uy.bits.ext_fraclm
+#endif /* EXT_FRACLMBITS */
+		| uy.bits.ext_fracl)!=0))	/* y is NaN */
+	    return x+y;
+	if (expsignx==BIAS && ((ux.bits.ext_frach&~LDBL_NBIT)
+#ifdef EXT_FRACHMBITS
+		| ux.bits.ext_frachm
+#endif /* EXT_FRACHMBITS */
+#ifdef EXT_FRACLMBITS
+		| ux.bits.ext_fraclm
+#endif /* EXT_FRACLMBITS */
+		| ux.bits.ext_fracl)==0)
+	    return atanl(y);					/* x=1.0 */
+	m = ((expsigny>>15)&1)|((expsignx>>14)&2);	/* 2*sign(x)+sign(y) */
+
+    /* when y = 0 */
+	if(expty==0 && ((uy.bits.ext_frach&~LDBL_NBIT)
+#ifdef EXT_FRACHMBITS
+		| uy.bits.ext_frachm
+#endif /* EXT_FRACHMBITS */
+#ifdef EXT_FRACLMBITS
+		| uy.bits.ext_fraclm
+#endif /* EXT_FRACLMBITS */
+		| uy.bits.ext_fracl)==0) {
+	    switch(m) {
+		case 0: 
+		case 1: return y; 	/* atan(+-0,+anything)=+-0 */
+		case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
+		case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
+	    }
+	}
+    /* when x = 0 */
+	if(exptx==0 && ((ux.bits.ext_frach&~LDBL_NBIT)
+#ifdef EXT_FRACHMBITS
+		| ux.bits.ext_frachm
+#endif /* EXT_FRACHMBITS */
+#ifdef EXT_FRACLMBITS
+		| ux.bits.ext_fraclm
+#endif /* EXT_FRACLMBITS */
+		| ux.bits.ext_fracl)==0)
+	    return (expsigny<0)?  -pio2_hi-tiny: pio2_hi+tiny;
+
+    /* when x is INF */
+	if(exptx==BIAS+LDBL_MAX_EXP) {
+	    if(expty==BIAS+LDBL_MAX_EXP) {
+		switch(m) {
+		    case 0: return  pio2_hi*0.5+tiny;/* atan(+INF,+INF) */
+		    case 1: return -pio2_hi*0.5-tiny;/* atan(-INF,+INF) */
+		    case 2: return  1.5*pio2_hi+tiny;/*atan(+INF,-INF)*/
+		    case 3: return -1.5*pio2_hi-tiny;/*atan(-INF,-INF)*/
+		}
+	    } else {
+		switch(m) {
+		    case 0: return  zero  ;	/* atan(+...,+INF) */
+		    case 1: return -zero  ;	/* atan(-...,+INF) */
+		    case 2: return  pi+tiny  ;	/* atan(+...,-INF) */
+		    case 3: return -pi-tiny  ;	/* atan(-...,-INF) */
+		}
+	    }
+	}
+    /* when y is INF */
+	if(expty==BIAS+LDBL_MAX_EXP)
+	    return (expsigny<0)? -pio2_hi-tiny: pio2_hi+tiny;
+
+    /* compute y/x */
+	k = expty-exptx;
+	if(k > LDBL_MANT_DIG+2) {			/* |y/x| huge */
+	    z=pio2_hi+pio2_lo;
+	    m&=1;
+	}
+	else if(expsignx<0&&k<-LDBL_MANT_DIG-2) z=0.0; 	/* |y/x| tiny, x<0 */
+	else z=atanl(fabsl(y/x));		/* safe to do y/x */
+	switch (m) {
+	    case 0: return       z  ;	/* atan(+,+) */
+	    case 1: return      -z  ;	/* atan(-,+) */
+	    case 2: return  pi-(z-pi_lo);/* atan(+,-) */
+	    default: /* case 3 */
+	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
+	}
+}
diff --git a/lib/libm/src/e_rem_pio2.c b/lib/libm/src/e_rem_pio2.c
index 4de859150c8..c1eefcd96b5 100644
--- a/lib/libm/src/e_rem_pio2.c
+++ b/lib/libm/src/e_rem_pio2.c
@@ -23,23 +23,6 @@ static char rcsid[] = "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $
 #include "math.h"
 #include "math_private.h"
 
-/*
- * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 
- */
-static const int32_t two_over_pi[] = {
-0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 
-0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 
-0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 
-0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 
-0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 
-0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 
-0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 
-0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 
-0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 
-0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 
-0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 
-};
-
 static const int32_t npio2_hw[] = {
 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
@@ -161,7 +144,7 @@ __ieee754_rem_pio2(double x, double *y)
 	tx[2] = z;
 	nx = 3;
 	while(tx[nx-1]==zero) nx--;	/* skip zero term */
-	n  =  __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
+	n  =  __kernel_rem_pio2(tx,y,e0,nx,2);
 	if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
 	return n;
 }
diff --git a/lib/libm/src/e_sqrt.c b/lib/libm/src/e_sqrt.c
index 05355a33291..04c94bdeec1 100644
--- a/lib/libm/src/e_sqrt.c
+++ b/lib/libm/src/e_sqrt.c
@@ -84,7 +84,10 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $";
  *---------------
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 static	const double	one	= 1.0, tiny=1.0e-300;
@@ -442,4 +445,9 @@ B.  sqrt(x) by Reciproot Iteration
     (4)	Special cases (see (4) of Section A).	
  
  */
- 
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(sqrtl, sqrt);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/e_sqrtl.c b/lib/libm/src/e_sqrtl.c
new file mode 100644
index 00000000000..7ecc83381d7
--- /dev/null
+++ b/lib/libm/src/e_sqrtl.c
@@ -0,0 +1,230 @@
+/*	$OpenBSD: e_sqrtl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*-
+ * Copyright (c) 2007 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.
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+#include <ieeefp.h>
+#include <math.h>
+
+#ifdef EXT_IMPLICIT_NBIT
+#define	LDBL_NBIT	0
+#else /* EXT_IMPLICIT_NBIT */
+#define	LDBL_NBIT	0x80000000
+#endif /* EXT_IMPLICIT_NBIT */
+
+/* Return (x + ulp) for normal positive x. Assumes no overflow. */
+static inline long double
+inc(long double x)
+{
+	struct ieee_ext *p = (struct ieee_ext *)&x;
+
+#ifdef EXT_FRACHMBITS
+	uint64_t frach;
+
+	frach = ((uint64_t)p->ext_frach << EXT_FRACHMBITS) | p->ext_frachm;
+	frach++;
+	p->ext_frach = frach >> EXT_FRACHMBITS;
+	p->ext_frachm = frach & 0x00000000ffffffff;
+#else /* EXT_FRACHMBITS */
+	uint32_t frach;
+
+	p->ext_frach++;
+	frach = p->ext_frach;
+#endif /* EXT_FRACHMBITS */
+
+	if (frach == 0) {
+
+#ifdef EXT_FRACLMBITS
+		uint64_t fracl;
+
+		fracl = ((uint64_t)p->ext_fraclm << EXT_FRACLBITS) |
+			p->ext_fracl;
+		fracl++;
+		p->ext_fraclm = fracl >> EXT_FRACLBITS;
+		p->ext_fracl = fracl & 0x00000000ffffffff;
+#else /* EXT_FRACLMBITS */
+		uint32_t fracl;
+
+		p->ext_fracl++;
+		fracl = p->ext_fracl;
+#endif /* EXT_FRACLMBITS */
+
+		if (fracl == 0) {
+			p->ext_exp++;
+			p->ext_frach |= LDBL_NBIT;
+		}
+	}
+
+	return x;
+}
+
+/* Return (x - ulp) for normal positive x. Assumes no underflow. */
+static inline long double
+dec(long double x)
+{
+	struct ieee_ext *p = (struct ieee_ext *)&x;
+
+#ifdef EXT_FRACLMBITS
+	uint64_t fracl;
+
+	fracl = ((uint64_t)p->ext_fraclm << EXT_FRACLBITS) | p->ext_fracl;
+	fracl--;
+	p->ext_fraclm = fracl >> EXT_FRACLBITS;
+	p->ext_fracl = fracl & 0x00000000ffffffff;
+#else /* EXT_FRACLMBITS */
+	uint32_t fracl;
+
+	p->ext_fracl--;
+	fracl = p->ext_fracl;
+#endif /* EXT_FRACLMBITS */
+
+	if (fracl == 0) {
+
+#ifdef EXT_FRACHMBITS
+		uint64_t frach;
+
+		frach = ((uint64_t)p->ext_frach << EXT_FRACHMBITS) |
+			p->ext_frachm;
+		frach--;
+		p->ext_frach = frach >> EXT_FRACHMBITS;
+		p->ext_frachm = frach & 0x00000000ffffffff;
+#else /* EXT_FRACHMBITS */
+		uint32_t frach;
+
+		p->ext_frach--;
+		frach = p->ext_frach;
+#endif /* EXT_FRACHMBITS */
+
+		if (frach == LDBL_NBIT) {
+			p->ext_exp--;
+			p->ext_frach |= LDBL_NBIT;
+		}
+	}
+
+	return x;
+}
+
+/*
+ * This is slow, but simple and portable. You should use hardware sqrt
+ * if possible.
+ */
+
+long double
+sqrtl(long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} u;
+	int k, r;
+	long double lo, xn;
+
+	u.e = x;
+
+	/* If x = NaN, then sqrt(x) = NaN. */
+	/* If x = Inf, then sqrt(x) = Inf. */
+	/* If x = -Inf, then sqrt(x) = NaN. */
+	if (u.bits.ext_exp == LDBL_MAX_EXP * 2 - 1)
+		return (x * x + x);
+
+	/* If x = +-0, then sqrt(x) = +-0. */
+	if ((u.bits.ext_frach
+#ifdef EXT_FRACHMBITS
+		| u.bits.ext_frachm
+#endif /* EXT_FRACHMBITS */
+#ifdef EXT_FRACLMBITS
+		| u.bits.ext_fraclm
+#endif /* EXT_FRACLMBITS */
+		| u.bits.ext_fracl | u.bits.ext_exp) == 0)
+		return (x);
+
+	/* If x < 0, then raise invalid and return NaN */
+	if (u.bits.ext_sign)
+		return ((x - x) / (x - x));
+
+	if (u.bits.ext_exp == 0) {
+		/* Adjust subnormal numbers. */
+		u.e *= 0x1.0p514;
+		k = -514;
+	} else {
+		k = 0;
+	}
+	/*
+	 * u.e is a normal number, so break it into u.e = e*2^n where
+	 * u.e = (2*e)*2^2k for odd n and u.e = (4*e)*2^2k for even n.
+	 */
+	if ((u.bits.ext_exp - 0x3ffe) & 1) {	/* n is odd.     */
+		k += u.bits.ext_exp - 0x3fff;	/* 2k = n - 1.   */
+		u.bits.ext_exp = 0x3fff;	/* u.e in [1,2). */
+	} else {
+		k += u.bits.ext_exp - 0x4000;	/* 2k = n - 2.   */
+		u.bits.ext_exp = 0x4000;	/* u.e in [2,4). */
+	}
+
+	/*
+	 * Newton's iteration.
+	 * Split u.e into a high and low part to achieve additional precision.
+	 */
+	xn = sqrt(u.e);			/* 53-bit estimate of sqrtl(x). */
+#if LDBL_MANT_DIG > 100
+	xn = (xn + (u.e / xn)) * 0.5;	/* 106-bit estimate. */
+#endif
+	lo = u.e;
+	u.bits.ext_fracl = 0;		/* Zero out lower bits. */
+#ifdef EXT_FRACLMBITS
+	u.bits.ext_fraclm = 0;
+#endif /* EXT_FRACLMBITS */
+	lo = (lo - u.e) / xn;		/* Low bits divided by xn. */
+	xn = xn + (u.e / xn);		/* High portion of estimate. */
+	u.e = xn + lo;			/* Combine everything. */
+	u.bits.ext_exp += (k >> 1) - 1;
+
+	fpsetsticky(fpgetsticky() & ~FP_X_IMP);
+	r = fpsetround(FP_RZ);		/* Set to round-toward-zero. */
+	xn = x / u.e;			/* Chopped quotient (inexact?). */
+
+	if (!(fpgetsticky() & FP_X_IMP)) { /* Quotient is exact. */
+		if (xn == u.e) {
+			fpsetround(r);
+			return (u.e);
+		}
+		/* Round correctly for inputs like x = y**2 - ulp. */
+		xn = dec(xn);		/* xn = xn - ulp. */
+	}
+
+	if (r == FP_RN) {
+		xn = inc(xn);		/* xn = xn + ulp. */
+	} else if (r == FP_RP) {
+		u.e = inc(u.e);		/* u.e = u.e + ulp. */
+		xn = inc(xn);		/* xn  = xn + ulp. */
+	}
+	u.e = u.e + xn;			/* Chopped sum. */
+	fpsetround(r);			/* Restore env and raise inexact */
+	u.bits.ext_exp--;
+	return (u.e);
+}
diff --git a/lib/libm/src/k_rem_pio2.c b/lib/libm/src/k_rem_pio2.c
index ad651919258..6949d6ce783 100644
--- a/lib/libm/src/k_rem_pio2.c
+++ b/lib/libm/src/k_rem_pio2.c
@@ -15,8 +15,8 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $
 #endif
 
 /*
- * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
- * double x[],y[]; int e0,nx,prec; int ipio2[];
+ * __kernel_rem_pio2(x,y,e0,nx,prec)
+ * double x[],y[]; int e0,nx,prec;
  * 
  * __kernel_rem_pio2 return the last three digits of N with 
  *		y = x - N*pi/2
@@ -60,7 +60,8 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $
  *		r_head = t+w;
  *		r_tail = w - (r_head - t);
  *
- *	e0	The exponent of x[0]
+ *	e0	The exponent of x[0]. Must be <= 16360 or you need to
+ *              expand the ipio2 table.
  *
  *	nx	dimension of x[]
  *
@@ -70,13 +71,6 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $
  *			2	64  bits (extended)
  *			3	113 bits (quad)
  *
- *	ipio2[]
- *		integer array, contains the (24*i)-th to (24*i+23)-th 
- *		bit of 2/pi after binary point. The corresponding 
- *		floating value is
- *
- *			ipio2[i] * 2^(-24(i+1)).
- *
  * External function:
  *	double scalbn(), floor();
  *
@@ -130,11 +124,150 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $
  * to produce the hexadecimal values shown.
  */
 
-#include "math.h"
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
 
+/*
+ * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
+ *
+ *		integer array, contains the (24*i)-th to (24*i+23)-th 
+ *		bit of 2/pi after binary point. The corresponding 
+ *		floating value is
+ *
+ *			ipio2[i] * 2^(-24(i+1)).
+ *
+ * NB: This table must have at least (e0-3)/24 + jk terms.
+ *     For quad precision (e0 <= 16360, jk = 6), this is 686.
+ */
+static const int32_t ipio2[] = {
+0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 
+0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 
+0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 
+0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 
+0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 
+0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 
+0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 
+0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 
+0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 
+0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 
+0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 
+
+#if LDBL_MAX_EXP > 1024
+#if LDBL_MAX_EXP > 16384
+#error "ipio2 table needs to be expanded"
+#endif
+0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
+0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2,
+0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35,
+0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
+0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C,
+0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4,
+0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
+0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7,
+0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19,
+0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
+0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16,
+0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6,
+0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
+0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48,
+0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3,
+0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
+0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55,
+0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612,
+0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
+0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC,
+0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B,
+0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
+0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4,
+0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB,
+0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
+0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C,
+0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F,
+0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
+0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437,
+0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B,
+0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
+0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD,
+0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3,
+0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
+0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717,
+0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F,
+0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
+0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB,
+0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51,
+0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
+0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C,
+0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6,
+0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
+0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED,
+0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328,
+0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
+0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0,
+0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B,
+0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
+0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3,
+0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F,
+0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
+0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B,
+0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4,
+0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
+0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31,
+0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30,
+0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
+0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E,
+0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1,
+0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
+0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4,
+0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08,
+0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
+0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9,
+0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4,
+0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
+0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C,
+0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0,
+0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
+0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0,
+0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC,
+0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
+0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893,
+0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7,
+0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
+0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F,
+0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4,
+0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
+0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B,
+0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2,
+0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
+0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E,
+0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569,
+0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
+0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9,
+0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D,
+0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
+0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855,
+0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569,
+0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
+0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE,
+0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41,
+0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
+0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F,
+0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110,
+0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
+0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365,
+0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A,
+0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
+0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5,
+0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616,
+0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
+0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0,
+#endif
+
+};
+
 static const double PIo2[] = {
   1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
   7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
@@ -153,8 +286,7 @@ two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
 twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
 
 int
-__kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec,
-    const int32_t *ipio2) 
+__kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec)
 {
 	int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
 	double z,fw,f[20],fq[20],q[20];
@@ -278,6 +410,7 @@ recompute:
 	    case 2:
 		fw = 0.0;
 		for (i=jz;i>=0;i--) fw += fq[i]; 
+		STRICT_ASSIGN(double,fw,fw);
 		y[0] = (ih==0)? fw: -fw; 
 		fw = fq[0]-fw;
 		for (i=1;i<=jz;i++) fw += fq[i];
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/arch/mc68881/s_finite.S b/lib/libm/src/ld128/s_nanl.c
index dcc0fb457c5..1edcac513e7 100644
--- a/lib/libm/arch/mc68881/s_finite.S
+++ b/lib/libm/src/ld128/s_nanl.c
@@ -1,12 +1,8 @@
-/*	$OpenBSD: s_finite.S,v 1.3 2005/08/02 11:17:32 espie Exp $ */
+/*	$OpenBSD: s_nanl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
 /*-
- * Copyright (c) 1990 The Regents of the University of California.
+ * Copyright (c) 2007 David Schultz
  * All rights reserved.
  *
- * This code is derived from software contributed to Berkeley by
- * the Systems Programming Group of the University of Utah Computer
- * Science Department.
- *
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions
  * are met:
@@ -15,14 +11,11 @@
  * 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.
- * 3. Neither the name of the University nor the names of its contributors
- *    may be used to endorse or promote products derived from this software
- *    without specific prior written permission.
  *
- * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * 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 THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * 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)
@@ -30,20 +23,28 @@
  * 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 <machine/asm.h>
+#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 */
 
-| finite(x)
-| returns the value TRUE if -INF < x < +INF and returns FALSE otherwise.
-ENTRY(finite)
-	movw	#0x7FF0,d0
-	movw	sp@(4),d1
-	andw	d0,d1
-	cmpw	d0,d1
-	beq	Lnotfin
-	moveq	#1,d0
-	rts
-Lnotfin:
-	clrl	d0
-	rts
+	return (u.e);
+}
diff --git a/lib/libm/src/ld80/invtrig.c b/lib/libm/src/ld80/invtrig.c
new file mode 100644
index 00000000000..31d261a4f56
--- /dev/null
+++ b/lib/libm/src/ld80/invtrig.c
@@ -0,0 +1,80 @@
+/*	$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.66666666666666666631e-01L,
+pS1 = -4.16313987993683104320e-01L,
+pS2 =  3.69068046323246813704e-01L,
+pS3 = -1.36213932016738603108e-01L,
+pS4 =  1.78324189708471965733e-02L,
+pS5 = -2.19216428382605211588e-04L,
+pS6 = -7.10526623669075243183e-06L,
+qS1 = -2.94788392796209867269e+00L,
+qS2 =  3.27309890266528636716e+00L,
+qS3 = -1.68285799854822427013e+00L,
+qS4 =  3.90699412641738801874e-01L,
+qS5 = -3.14365703596053263322e-02L;
+
+/*
+ * atanl()
+ */
+const long double atanhi[] = {
+	 4.63647609000806116202e-01L,
+	 7.85398163397448309628e-01L,
+	 9.82793723247329067960e-01L,
+	 1.57079632679489661926e+00L,
+};
+
+const long double atanlo[] = {
+	 1.18469937025062860669e-20L,
+	-1.25413940316708300586e-20L,
+	 2.55232234165405176172e-20L,
+	-2.50827880633416601173e-20L,
+};
+
+const long double aT[] = {
+	 3.33333333333333333017e-01L,
+	-1.99999999999999632011e-01L,
+	 1.42857142857046531280e-01L,
+	-1.11111111100562372733e-01L,
+	 9.09090902935647302252e-02L,
+	-7.69230552476207730353e-02L,
+	 6.66661718042406260546e-02L,
+	-5.88158892835030888692e-02L,
+	 5.25499891539726639379e-02L,
+	-4.70119845393155721494e-02L,
+	 4.03539201366454414072e-02L,
+	-2.91303858419364158725e-02L,
+	 1.24822046299269234080e-02L,
+};
+
+const long double pi_lo = -5.01655761266833202345e-20L;
diff --git a/lib/libm/src/ld80/invtrig.h b/lib/libm/src/ld80/invtrig.h
new file mode 100644
index 00000000000..0dc77e01e34
--- /dev/null
+++ b/lib/libm/src/ld80/invtrig.h
@@ -0,0 +1,119 @@
+/*	$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/ld80/invtrig.h,v 1.2 2008/08/02 03:56:22 das Exp $
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+
+#define	BIAS		(LDBL_MAX_EXP - 1)
+#define	MANH_SIZE	EXT_FRACHBITS
+
+/* Approximation thresholds. */
+#define	ASIN_LINEAR	(BIAS - 32)	/* 2**-32 */
+#define	ACOS_CONST	(BIAS - 65)	/* 2**-65 */
+#define	ATAN_CONST	(BIAS + 65)	/* 2**65 */
+#define	ATAN_LINEAR	(BIAS - 32)	/* 2**-32 */
+
+/* 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	qS1	_ItL_qS1
+#define	qS2	_ItL_qS2
+#define	qS3	_ItL_qS3
+#define	qS4	_ItL_qS4
+#define	qS5	_ItL_qS5
+#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]
+
+#ifdef STRUCT_DECLS
+typedef struct longdouble {
+	uint64_t mant;
+	uint16_t expsign;
+} LONGDOUBLE;
+#else
+typedef long double LONGDOUBLE;
+#endif
+
+extern const LONGDOUBLE pS0, pS1, pS2, pS3, pS4, pS5, pS6;
+extern const LONGDOUBLE qS1, qS2, qS3, qS4, qS5;
+extern const LONGDOUBLE atanhi[], atanlo[], aT[];
+extern const LONGDOUBLE pi_lo;
+
+#ifndef STRUCT_DECLS
+
+static inline long double
+P(long double x)
+{
+
+	return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \
+		(pS4 + x * (pS5 + x * pS6)))))));
+}
+
+static inline long double
+Q(long double x)
+{
+
+	return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5)))));
+}
+
+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]))))));
+}
+
+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])))));
+}
+
+#endif
diff --git a/lib/libm/src/ld80/k_cosl.c b/lib/libm/src/ld80/k_cosl.c
new file mode 100644
index 00000000000..c9cec2fe372
--- /dev/null
+++ b/lib/libm/src/ld80/k_cosl.c
@@ -0,0 +1,76 @@
+/*	$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.
+ * ====================================================
+ */
+
+/*
+ * ld80 version of k_cos.c.  See ../k_cos.c for most comments.
+ */
+
+#include "math_private.h"
+
+/*
+ * Domain [-0.7854, 0.7854], range ~[-2.43e-23, 2.425e-23]:
+ * |cos(x) - c(x)| < 2**-75.1
+ *
+ * The coefficients of c(x) were generated by a pari-gp script using
+ * a Remez algorithm that searches for the best higher coefficients
+ * after rounding leading coefficients to a specified precision.
+ *
+ * Simpler methods like Chebyshev or basic Remez barely suffice for
+ * cos() in 64-bit precision, because we want the coefficient of x^2
+ * to be precisely -0.5 so that multiplying by it is exact, and plain
+ * rounding of the coefficients of a good polynomial approximation only
+ * gives this up to about 64-bit precision.  Plain rounding also gives
+ * a mediocre approximation for the coefficient of x^4, but a rounding
+ * error of 0.5 ulps for this coefficient would only contribute ~0.01
+ * ulps to the final error, so this is unimportant.  Rounding errors in
+ * higher coefficients are even less important.
+ *
+ * In fact, coefficients above the x^4 one only need to have 53-bit
+ * precision, and this is more efficient.  We get this optimization
+ * almost for free from the complications needed to search for the best
+ * higher coefficients.
+ */
+static const double
+one = 1.0;
+
+#if defined(__amd64__) || defined(__i386__)
+/* Long double constants are slow on these arches, and broken on i386. */
+static const volatile double
+C1hi = 0.041666666666666664,		/*  0x15555555555555.0p-57 */
+C1lo = 2.2598839032744733e-18;		/*  0x14d80000000000.0p-111 */
+#define	C1	((long double)C1hi + C1lo)
+#else
+static const long double
+C1 =  0.0416666666666666666136L;	/*  0xaaaaaaaaaaaaaa9b.0p-68 */
+#endif
+
+static const double
+C2 = -0.0013888888888888874,		/* -0x16c16c16c16c10.0p-62 */
+C3 =  0.000024801587301571716,		/*  0x1a01a01a018e22.0p-68 */
+C4 = -0.00000027557319215507120,	/* -0x127e4fb7602f22.0p-74 */
+C5 =  0.0000000020876754400407278,	/*  0x11eed8caaeccf1.0p-81 */
+C6 = -1.1470297442401303e-11,		/* -0x19393412bd1529.0p-89 */
+C7 =  4.7383039476436467e-14;		/*  0x1aac9d9af5c43e.0p-97 */
+
+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))))));
+	hz = 0.5*z;
+	w  = one-hz;
+	return w + (((one-w)-hz) + (z*r-x*y));
+}
diff --git a/lib/libm/src/ld80/k_sinl.c b/lib/libm/src/ld80/k_sinl.c
new file mode 100644
index 00000000000..1bcff39607e
--- /dev/null
+++ b/lib/libm/src/ld80/k_sinl.c
@@ -0,0 +1,60 @@
+/*	$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.
+ * ====================================================
+ */
+
+/*
+ * ld80 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.89e-22, 1.915e-22]
+ * |sin(x)/x - s(x)| < 2**-72.1
+ *
+ * See ../ld80/k_cosl.c for more details about the polynomial.
+ */
+#if defined(__amd64__) || defined(__i386__)
+/* Long double constants are slow on these arches, and broken on i386. */
+static const volatile double
+S1hi = -0.16666666666666666,		/* -0x15555555555555.0p-55 */
+S1lo = -9.2563760475949941e-18;		/* -0x15580000000000.0p-109 */
+#define	S1	((long double)S1hi + S1lo)
+#else
+static const long double
+S1 = -0.166666666666666666671L;		/* -0xaaaaaaaaaaaaaaab.0p-66 */
+#endif
+
+static const double
+S2 =  0.0083333333333333332,		/*  0x11111111111111.0p-59 */
+S3 = -0.00019841269841269427,		/* -0x1a01a01a019f81.0p-65 */
+S4 =  0.0000027557319223597490,		/*  0x171de3a55560f7.0p-71 */
+S5 = -0.000000025052108218074604,	/* -0x1ae64564f16cad.0p-78 */
+S6 =  1.6059006598854211e-10,		/*  0x161242b90243b5.0p-85 */
+S7 = -7.6429779983024564e-13,		/* -0x1ae42ebd1b2e00.0p-93 */
+S8 =  2.6174587166648325e-15;		/*  0x179372ea0b3f64.0p-101 */
+
+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)))));
+	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/ld80/k_tanl.c b/lib/libm/src/ld80/k_tanl.c
new file mode 100644
index 00000000000..0054cd2b22b
--- /dev/null
+++ b/lib/libm/src/ld80/k_tanl.c
@@ -0,0 +1,122 @@
+/*	$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.
+ * ====================================================
+ */
+
+/*
+ * ld80 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 ~[-2.25e-22, 1.921e-22]
+ * |tan(x)/x - t(x)| < 2**-71.9
+ *
+ * See k_cosl.c for more details about the polynomial.
+ */
+#if defined(__amd64__) || defined(__i386__)
+/* Long double constants are slow on these arches, and broken on i386. */
+static const volatile double
+T3hi =  0.33333333333333331,		/*  0x15555555555555.0p-54 */
+T3lo =  1.8350121769317163e-17,		/*  0x15280000000000.0p-108 */
+T5hi =  0.13333333333333336,		/*  0x11111111111112.0p-55 */
+T5lo =  1.3051083651294260e-17,		/*  0x1e180000000000.0p-109 */
+T7hi =  0.053968253968250494,		/*  0x1ba1ba1ba1b827.0p-57 */
+T7lo =  3.1509625637859973e-18,		/*  0x1d100000000000.0p-111 */
+pio4_hi =  0.78539816339744828,		/*  0x1921fb54442d18.0p-53 */
+pio4_lo =  3.0628711372715500e-17,	/*  0x11a80000000000.0p-107 */
+pio4lo_hi = -1.2541394031670831e-20,	/* -0x1d9cceba3f91f2.0p-119 */
+pio4lo_lo =  6.1493048227390915e-37;	/*  0x1a280000000000.0p-173 */
+#define	T3	((long double)T3hi + T3lo)
+#define	T5	((long double)T5hi + T5lo)
+#define	T7	((long double)T7hi + T7lo)
+#define	pio4	((long double)pio4_hi + pio4_lo)
+#define	pio4lo	((long double)pio4lo_hi + pio4lo_lo)
+#else
+static const long double
+T3 =   0.333333333333333333180L,	/*  0xaaaaaaaaaaaaaaa5.0p-65 */
+T5 =   0.133333333333333372290L,	/*  0x88888888888893c3.0p-66 */
+T7 =   0.0539682539682504975744L,	/*  0xdd0dd0dd0dc13ba2.0p-68 */
+pio4 = 0.785398163397448309628L,	/*  0xc90fdaa22168c235.0p-64 */
+pio4lo = -1.25413940316708300586e-20L;	/* -0xece675d1fc8f8cbb.0p-130 */
+#endif
+
+static const double
+T9  =  0.021869488536312216,		/*  0x1664f4882cc1c2.0p-58 */
+T11 =  0.0088632355256619590,		/*  0x1226e355c17612.0p-59 */
+T13 =  0.0035921281113786528,		/*  0x1d6d3d185d7ff8.0p-61 */
+T15 =  0.0014558334756312418,		/*  0x17da354aa3f96b.0p-62 */
+T17 =  0.00059003538700862256,		/*  0x13559358685b83.0p-63 */
+T19 =  0.00023907843576635544,		/*  0x1f56242026b5be.0p-65 */
+T21 =  0.000097154625656538905,		/*  0x1977efc26806f4.0p-66 */
+T23 =  0.000038440165747303162,		/*  0x14275a09b3ceac.0p-67 */
+T25 =  0.000018082171885432524,		/*  0x12f5e563e5487e.0p-68 */
+T27 =  0.0000024196006108814377,	/*  0x144c0d80cc6896.0p-71 */
+T29 =  0.0000078293456938132840,	/*  0x106b59141a6cb3.0p-69 */
+T31 = -0.0000032609076735050182,	/* -0x1b5abef3ba4b59.0p-71 */
+T33 =  0.0000023261313142559411;	/*  0x13835436c0c87f.0p-71 */
+
+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))))));
+	v = z * (T7 + w * (T11 + w * (T15 + w * (T19 + w * (T23 +
+	    w * (T27 + w * T31))))));
+	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/ld80/s_exp2l.c b/lib/libm/src/ld80/s_exp2l.c
new file mode 100644
index 00000000000..8e6cea6d183
--- /dev/null
+++ b/lib/libm/src/ld80/s_exp2l.c
@@ -0,0 +1,291 @@
+/*	$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)
+
+static const long double huge = 0x1p10000L;
+#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 double
+    redux     = 0x1.8p63 / TBLSIZE,
+    P1        = 0x1.62e42fefa39efp-1,
+    P2        = 0x1.ebfbdff82c58fp-3,
+    P3        = 0x1.c6b08d7049fap-5,
+    P4        = 0x1.3b2ab6fba4da5p-7,
+    P5        = 0x1.5d8804780a736p-10,
+    P6        = 0x1.430918835e33dp-13;
+
+static const double tbl[TBLSIZE * 2] = {
+	0x1.6a09e667f3bcdp-1,	-0x1.bdd3413b2648p-55,
+	0x1.6c012750bdabfp-1,	-0x1.2895667ff0cp-57,
+	0x1.6dfb23c651a2fp-1,	-0x1.bbe3a683c88p-58,
+	0x1.6ff7df9519484p-1,	-0x1.83c0f25860fp-56,
+	0x1.71f75e8ec5f74p-1,	-0x1.16e4786887bp-56,
+	0x1.73f9a48a58174p-1,	-0x1.0a8d96c65d5p-55,
+	0x1.75feb564267c9p-1,	-0x1.0245957316ep-55,
+	0x1.780694fde5d3fp-1,	 0x1.866b80a0216p-55,
+	0x1.7a11473eb0187p-1,	-0x1.41577ee0499p-56,
+	0x1.7c1ed0130c132p-1,	 0x1.f124cd1164ep-55,
+	0x1.7e2f336cf4e62p-1,	 0x1.05d02ba157ap-57,
+	0x1.80427543e1a12p-1,	-0x1.27c86626d97p-55,
+	0x1.82589994cce13p-1,	-0x1.d4c1dd41533p-55,
+	0x1.8471a4623c7adp-1,	-0x1.8d684a341cep-56,
+	0x1.868d99b4492edp-1,	-0x1.fc6f89bd4f68p-55,
+	0x1.88ac7d98a6699p-1,	 0x1.994c2f37cb5p-55,
+	0x1.8ace5422aa0dbp-1,	 0x1.6e9f156864bp-55,
+	0x1.8cf3216b5448cp-1,	-0x1.0d55e32e9e4p-57,
+	0x1.8f1ae99157736p-1,	 0x1.5cc13a2e397p-56,
+	0x1.9145b0b91ffc6p-1,	-0x1.dd6792e5825p-55,
+	0x1.93737b0cdc5e5p-1,	-0x1.75fc781b58p-58,
+	0x1.95a44cbc8520fp-1,	-0x1.64b7c96a5fp-57,
+	0x1.97d829fde4e5p-1,	-0x1.d185b7c1b86p-55,
+	0x1.9a0f170ca07bap-1,	-0x1.173bd91cee6p-55,
+	0x1.9c49182a3f09p-1,	 0x1.c7c46b071f2p-57,
+	0x1.9e86319e32323p-1,	 0x1.824ca78e64cp-57,
+	0x1.a0c667b5de565p-1,	-0x1.359495d1cd5p-55,
+	0x1.a309bec4a2d33p-1,	 0x1.6305c7ddc368p-55,
+	0x1.a5503b23e255dp-1,	-0x1.d2f6edb8d42p-55,
+	0x1.a799e1330b358p-1,	 0x1.bcb7ecac564p-55,
+	0x1.a9e6b5579fdbfp-1,	 0x1.0fac90ef7fdp-55,
+	0x1.ac36bbfd3f37ap-1,	-0x1.f9234cae76dp-56,
+	0x1.ae89f995ad3adp-1,	 0x1.7a1cd345dcc8p-55,
+	0x1.b0e07298db666p-1,	-0x1.bdef54c80e4p-55,
+	0x1.b33a2b84f15fbp-1,	-0x1.2805e3084d8p-58,
+	0x1.b59728de5593ap-1,	-0x1.c71dfbbba6ep-55,
+	0x1.b7f76f2fb5e47p-1,	-0x1.5584f7e54acp-57,
+	0x1.ba5b030a1064ap-1,	-0x1.efcd30e5429p-55,
+	0x1.bcc1e904bc1d2p-1,	 0x1.23dd07a2d9fp-56,
+	0x1.bf2c25bd71e09p-1,	-0x1.efdca3f6b9c8p-55,
+	0x1.c199bdd85529cp-1,	 0x1.11065895049p-56,
+	0x1.c40ab5fffd07ap-1,	 0x1.b4537e083c6p-55,
+	0x1.c67f12e57d14bp-1,	 0x1.2884dff483c8p-55,
+	0x1.c8f6d9406e7b5p-1,	 0x1.1acbc48805cp-57,
+	0x1.cb720dcef9069p-1,	 0x1.503cbd1e94ap-57,
+	0x1.cdf0b555dc3fap-1,	-0x1.dd83b53829dp-56,
+	0x1.d072d4a07897cp-1,	-0x1.cbc3743797a8p-55,
+	0x1.d2f87080d89f2p-1,	-0x1.d487b719d858p-55,
+	0x1.d5818dcfba487p-1,	 0x1.2ed02d75b37p-56,
+	0x1.d80e316c98398p-1,	-0x1.11ec18bedep-55,
+	0x1.da9e603db3285p-1,	 0x1.c2300696db5p-55,
+	0x1.dd321f301b46p-1,	 0x1.2da5778f019p-55,
+	0x1.dfc97337b9b5fp-1,	-0x1.1a5cd4f184b8p-55,
+	0x1.e264614f5a129p-1,	-0x1.7b627817a148p-55,
+	0x1.e502ee78b3ff6p-1,	 0x1.39e8980a9cdp-56,
+	0x1.e7a51fbc74c83p-1,	 0x1.2d522ca0c8ep-55,
+	0x1.ea4afa2a490dap-1,	-0x1.e9c23179c288p-55,
+	0x1.ecf482d8e67f1p-1,	-0x1.c93f3b411ad8p-55,
+	0x1.efa1bee615a27p-1,	 0x1.dc7f486a4b68p-55,
+	0x1.f252b376bba97p-1,	 0x1.3a1a5bf0d8e8p-55,
+	0x1.f50765b6e454p-1,	 0x1.9d3e12dd8a18p-55,
+	0x1.f7bfdad9cbe14p-1,	-0x1.dbb12d00635p-55,
+	0x1.fa7c1819e90d8p-1,	 0x1.74853f3a593p-56,
+	0x1.fd3c22b8f71f1p-1,	 0x1.2eb74966578p-58,
+	0x1p+0,	 0x0p+0,
+	0x1.0163da9fb3335p+0,	 0x1.b61299ab8cd8p-54,
+	0x1.02c9a3e778061p+0,	-0x1.19083535b08p-56,
+	0x1.04315e86e7f85p+0,	-0x1.0a31c1977c98p-54,
+	0x1.059b0d3158574p+0,	 0x1.d73e2a475b4p-55,
+	0x1.0706b29ddf6dep+0,	-0x1.c91dfe2b13cp-55,
+	0x1.0874518759bc8p+0,	 0x1.186be4bb284p-57,
+	0x1.09e3ecac6f383p+0,	 0x1.14878183161p-54,
+	0x1.0b5586cf9890fp+0,	 0x1.8a62e4adc61p-54,
+	0x1.0cc922b7247f7p+0,	 0x1.01edc16e24f8p-54,
+	0x1.0e3ec32d3d1a2p+0,	 0x1.03a1727c58p-59,
+	0x1.0fb66affed31bp+0,	-0x1.b9bedc44ebcp-57,
+	0x1.11301d0125b51p+0,	-0x1.6c51039449bp-54,
+	0x1.12abdc06c31ccp+0,	-0x1.1b514b36ca8p-58,
+	0x1.1429aaea92dep+0,	-0x1.32fbf9af1368p-54,
+	0x1.15a98c8a58e51p+0,	 0x1.2406ab9eeabp-55,
+	0x1.172b83c7d517bp+0,	-0x1.19041b9d78ap-55,
+	0x1.18af9388c8deap+0,	-0x1.11023d1970f8p-54,
+	0x1.1a35beb6fcb75p+0,	 0x1.e5b4c7b4969p-55,
+	0x1.1bbe084045cd4p+0,	-0x1.95386352ef6p-54,
+	0x1.1d4873168b9aap+0,	 0x1.e016e00a264p-54,
+	0x1.1ed5022fcd91dp+0,	-0x1.1df98027bb78p-54,
+	0x1.2063b88628cd6p+0,	 0x1.dc775814a85p-55,
+	0x1.21f49917ddc96p+0,	 0x1.2a97e9494a6p-55,
+	0x1.2387a6e756238p+0,	 0x1.9b07eb6c7058p-54,
+	0x1.251ce4fb2a63fp+0,	 0x1.ac155bef4f5p-55,
+	0x1.26b4565e27cddp+0,	 0x1.2bd339940eap-55,
+	0x1.284dfe1f56381p+0,	-0x1.a4c3a8c3f0d8p-54,
+	0x1.29e9df51fdee1p+0,	 0x1.612e8afad12p-55,
+	0x1.2b87fd0dad99p+0,	-0x1.10adcd6382p-59,
+	0x1.2d285a6e4030bp+0,	 0x1.0024754db42p-54,
+	0x1.2ecafa93e2f56p+0,	 0x1.1ca0f45d524p-56,
+	0x1.306fe0a31b715p+0,	 0x1.6f46ad23183p-55,
+	0x1.32170fc4cd831p+0,	 0x1.a9ce78e1804p-55,
+	0x1.33c08b26416ffp+0,	 0x1.327218436598p-54,
+	0x1.356c55f929ff1p+0,	-0x1.b5cee5c4e46p-55,
+	0x1.371a7373aa9cbp+0,	-0x1.63aeabf42ebp-54,
+	0x1.38cae6d05d866p+0,	-0x1.e958d3c99048p-54,
+	0x1.3a7db34e59ff7p+0,	-0x1.5e436d661f6p-56,
+	0x1.3c32dc313a8e5p+0,	-0x1.efff8375d2ap-54,
+	0x1.3dea64c123422p+0,	 0x1.ada0911f09fp-55,
+	0x1.3fa4504ac801cp+0,	-0x1.7d023f956fap-54,
+	0x1.4160a21f72e2ap+0,	-0x1.ef3691c309p-58,
+	0x1.431f5d950a897p+0,	-0x1.1c7dde35f7ap-55,
+	0x1.44e086061892dp+0,	 0x1.89b7a04ef8p-59,
+	0x1.46a41ed1d0057p+0,	 0x1.c944bd1648a8p-54,
+	0x1.486a2b5c13cdp+0,	 0x1.3c1a3b69062p-56,
+	0x1.4a32af0d7d3dep+0,	 0x1.9cb62f3d1be8p-54,
+	0x1.4bfdad5362a27p+0,	 0x1.d4397afec42p-56,
+	0x1.4dcb299fddd0dp+0,	 0x1.8ecdbbc6a78p-54,
+	0x1.4f9b2769d2ca7p+0,	-0x1.4b309d25958p-54,
+	0x1.516daa2cf6642p+0,	-0x1.f768569bd94p-55,
+	0x1.5342b569d4f82p+0,	-0x1.07abe1db13dp-55,
+	0x1.551a4ca5d920fp+0,	-0x1.d689cefede6p-55,
+	0x1.56f4736b527dap+0,	 0x1.9bb2c011d938p-54,
+	0x1.58d12d497c7fdp+0,	 0x1.295e15b9a1ep-55,
+	0x1.5ab07dd485429p+0,	 0x1.6324c0546478p-54,
+	0x1.5c9268a5946b7p+0,	 0x1.c4b1b81698p-60,
+	0x1.5e76f15ad2148p+0,	 0x1.ba6f93080e68p-54,
+	0x1.605e1b976dc09p+0,	-0x1.3e2429b56de8p-54,
+	0x1.6247eb03a5585p+0,	-0x1.383c17e40b48p-54,
+	0x1.6434634ccc32p+0,	-0x1.c483c759d89p-55,
+	0x1.6623882552225p+0,	-0x1.bb60987591cp-54,
+	0x1.68155d44ca973p+0,	 0x1.038ae44f74p-57,
+};
+
+/*
+ * exp2l(x): compute the base 2 exponential of x
+ *
+ * Accuracy: Peak error < 0.511 ulp.
+ *
+ * Method: (equally-spaced tables)
+ *
+ *   Reduce x:
+ *     x = 2**k + y, for integer k and |y| <= 1/2.
+ *     Thus we have exp2l(x) = 2**k * exp2(y).
+ *
+ *   Reduce y:
+ *     y = i/TBLSIZE + z for integer i near y * TBLSIZE.
+ *     Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z),
+ *     with |z| <= 2**-(TBLBITS+1).
+ *
+ *   We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a
+ *   degree-6 minimax polynomial with maximum error under 2**-69.
+ *   The table entries each have 104 bits of accuracy, encoded as
+ *   a pair of double precision values.
+ */
+long double
+exp2l(long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} u, v;
+	long double r, twopk, twopkp10000, z;
+	uint32_t es, hx, ix, i0;
+	int k;
+
+	/* Filter out exceptional cases. */
+	u.e = x;
+	hx = (u.bits.ext_sign << 15) | u.bits.ext_exp;
+	ix = hx & EXPMASK;
+	if (ix >= BIAS + 14) {		/* |x| >= 16384 or x is NaN */
+		if (ix == BIAS + LDBL_MAX_EXP) {
+			if ((u.bits.ext_frach != 1U << 31 &&
+			    u.bits.ext_fracl != 0) || (hx & 0x8000) == 0)
+				return (x + x);	/* x is +Inf or NaN */
+			else
+				return (0.0);	/* x is -Inf */
+		}
+		if (x >= 16384)
+			return (huge * huge);	/* overflow */
+		if (x <= -16446)
+			return (twom10000 * twom10000);	/* underflow */
+	} else if (ix <= BIAS - 66) {		/* |x| < 0x1p-66 */
+		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 = u.bits.ext_fracl + TBLSIZE / 2;
+	k = (int)i0 >> TBLBITS;
+	i0 = (i0 & (TBLSIZE - 1)) << 1;
+	u.e -= redux;
+	z = x - u.e;
+
+	v.bits.ext_frach = 1U << 31;
+	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 = exp2l(y) = exp2lt[i0] * p(z). */
+	long double t_hi = tbl[i0];
+	long double t_lo = tbl[i0 + 1];
+	/* XXX This gives > 1 ulp errors outside of FE_TONEAREST mode */
+	r = t_lo + (t_hi + t_lo) * z * (P1 + z * (P2 + z * (P3 + z * (P4
+	    + z * (P5 + z * P6))))) + t_hi;
+
+	/* 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/ld80/s_nanl.c b/lib/libm/src/ld80/s_nanl.c
new file mode 100644
index 00000000000..157e34e8be0
--- /dev/null
+++ b/lib/libm/src/ld80/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/ld80/s_nanl.c,v 1.2 2007/12/18 23:46:31 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[3];
+	} u;
+
+	_scan_nan(u.bits, 3, s);
+	u.ieee.ext_exp = 0x7fff;
+	u.ieee.ext_frach |= 0xc0000000;	/* make it a quiet NaN */
+
+	return (u.e);
+}
diff --git a/lib/libm/src/math_private.h b/lib/libm/src/math_private.h
index 0fc4f6a78d3..324e825e760 100644
--- a/lib/libm/src/math_private.h
+++ b/lib/libm/src/math_private.h
@@ -1,4 +1,4 @@
-/*	$OpenBSD: math_private.h,v 1.10 2008/09/07 20:36:09 martynas Exp $	*/
+/*	$OpenBSD: math_private.h,v 1.11 2008/12/09 20:00:35 martynas Exp $	*/
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -199,7 +199,7 @@ do {								\
  */
 #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
 #define	STRICT_ASSIGN(type, lval, rval)	((lval) = (rval))
-#else
+#else /* FLT_EVAL_METHOD == 0 || __GNUC__ == 0 */
 #define	STRICT_ASSIGN(type, lval, rval) do {	\
 	volatile type __lval;			\
 						\
@@ -210,15 +210,15 @@ do {								\
 		(lval) = __lval;		\
 	}					\
 } while (0)
-#endif
-#endif
+#endif /* FLT_EVAL_METHOD == 0 || __GNUC__ == 0 */
+#endif /* FLT_EVAL_METHOD */
 
 /* fdlibm kernel function */
 extern int    __ieee754_rem_pio2(double,double*);
 extern double __kernel_sin(double,double,int);
 extern double __kernel_cos(double,double);
 extern double __kernel_tan(double,double,int);
-extern int    __kernel_rem_pio2(double*,double*,int,int,int,const int*);
+extern int    __kernel_rem_pio2(double*,double*,int,int,int);
 
 /* float versions of fdlibm kernel functions */
 extern int   __ieee754_rem_pio2f(float,float*);
@@ -227,6 +227,16 @@ extern float __kernel_cosf(float,float);
 extern float __kernel_tanf(float,float,int);
 extern int   __kernel_rem_pio2f(float*,float*,int,int,int,const int*);
 
+/* long double precision kernel functions */
+long double __kernel_sinl(long double, long double, int);
+long double __kernel_cosl(long double, long double);
+long double __kernel_tanl(long double, long double, int);
+
+/*
+ * Common routine to process the arguments to nan(), nanf(), and nanl().
+ */
+void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
+
 /*
  * TRUNC() is a macro that sets the trailing 27 bits in the mantissa
  * of an IEEE double variable to zero.  It must be expression-like
diff --git a/lib/libm/src/s_atan.c b/lib/libm/src/s_atan.c
index 1a88d3225d3..ccc678ae55f 100644
--- a/lib/libm/src/s_atan.c
+++ b/lib/libm/src/s_atan.c
@@ -34,7 +34,10 @@ static char rcsid[] = "$NetBSD: s_atan.c,v 1.8 1995/05/10 20:46:45 jtc Exp $";
  * to produce the hexadecimal values shown.
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 static const double atanhi[] = {
@@ -117,3 +120,9 @@ atan(double x)
 	    return (hx<0)? -z:z;
 	}
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(atanl, atan);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_atanl.c b/lib/libm/src/s_atanl.c
new file mode 100644
index 00000000000..6ce61533a18
--- /dev/null
+++ b/lib/libm/src/s_atanl.c
@@ -0,0 +1,101 @@
+/*	$OpenBSD: s_atanl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/* @(#)s_atan.c 5.1 93/09/24 */
+/* FreeBSD: head/lib/msun/src/s_atan.c 176451 2008-02-22 02:30:36Z das */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/*
+ * See comments in s_atan.c.
+ * Converted to long double by David Schultz <das@FreeBSD.ORG>.
+ * Adapted for OpenBSD by Martynas Venckus <martynas@openbsd.org>.
+ */
+
+#include <float.h>
+#include <math.h>
+
+#include "invtrig.h"
+#include "math_private.h"
+
+#ifdef EXT_IMPLICIT_NBIT
+#define	LDBL_NBIT	0
+#else /* EXT_IMPLICIT_NBIT */
+#define	LDBL_NBIT	0x80000000
+#endif /* EXT_IMPLICIT_NBIT */
+
+static const long double
+one   = 1.0,
+huge   = 1.0e300;
+
+long double
+atanl(long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} u;
+	long double w,s1,s2,z;
+	int id;
+	int16_t expsign, expt;
+	int32_t expman;
+
+	u.e = x;
+	expsign = (u.bits.ext_sign << 15) | u.bits.ext_exp;
+	expt = expsign & 0x7fff;
+	if(expt >= ATAN_CONST) {	/* if |x| is large, atan(x)~=pi/2 */
+	    if(expt == BIAS + LDBL_MAX_EXP &&
+	       ((u.bits.ext_frach&~LDBL_NBIT)
+#ifdef EXT_FRACHMBITS
+		| u.bits.ext_frachm
+#endif /* EXT_FRACHMBITS */
+#ifdef EXT_FRACLMBITS
+		| u.bits.ext_fraclm
+#endif /* EXT_FRACLMBITS */
+		| u.bits.ext_fracl)!=0)
+		return x+x;		/* NaN */
+	    if(expsign>0) return  atanhi[3]+atanlo[3];
+	    else     return -atanhi[3]-atanlo[3];
+	}
+	/* Extract the exponent and the first few bits of the mantissa. */
+	/* XXX There should be a more convenient way to do this. */
+	expman = (expt << 8) |
+		((u.bits.ext_frach >> (EXT_FRACHBITS - 9)) & 0xff);
+	if (expman < ((BIAS - 2) << 8) + 0xc0) {	/* |x| < 0.4375 */
+	    if (expt < ATAN_LINEAR) {	/* if |x| is small, atanl(x)~=x */
+		if(huge+x>one) return x;	/* raise inexact */
+	    }
+	    id = -1;
+	} else {
+	x = fabsl(x);
+	if (expman < (BIAS << 8) + 0x30) {		/* |x| < 1.1875 */
+	    if (expman < ((BIAS - 1) << 8) + 0x60) {	/* 7/16 <=|x|<11/16 */
+		id = 0; x = (2.0*x-one)/(2.0+x);
+	    } else {			/* 11/16<=|x|< 19/16 */
+		id = 1; x  = (x-one)/(x+one);
+	    }
+	} else {
+	    if (expman < ((BIAS + 1) << 8) + 0x38) {	/* |x| < 2.4375 */
+		id = 2; x  = (x-1.5)/(one+1.5*x);
+	    } else {			/* 2.4375 <= |x| < 2^ATAN_CONST */
+		id = 3; x  = -1.0/x;
+	    }
+	}}
+    /* end of argument reduction */
+	z = x*x;
+	w = z*z;
+    /* break sum aT[i]z**(i+1) into odd and even poly */
+	s1 = z*T_even(w);
+	s2 = w*T_odd(w);
+	if (id<0) return x - x*(s1+s2);
+	else {
+	    z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
+	    return (expsign<0)? -z:z;
+	}
+}
diff --git a/lib/libm/src/s_copysign.c b/lib/libm/src/s_copysign.c
index ac060dc4977..20be8edb540 100644
--- a/lib/libm/src/s_copysign.c
+++ b/lib/libm/src/s_copysign.c
@@ -20,7 +20,10 @@ static char rcsid[] = "$NetBSD: s_copysign.c,v 1.8 1995/05/10 20:46:57 jtc Exp $
  * with the sign bit of y.
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 double
@@ -32,3 +35,9 @@ copysign(double x, double y)
 	SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
         return x;
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(copysignl, copysign);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_copysignl.c b/lib/libm/src/s_copysignl.c
new file mode 100644
index 00000000000..3d50a4118c0
--- /dev/null
+++ b/lib/libm/src/s_copysignl.c
@@ -0,0 +1,31 @@
+/*	$OpenBSD: s_copysignl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <math.h>
+
+long double
+copysignl(long double x, long double y)
+{
+	struct ieee_ext *px = (struct ieee_ext *)&x;
+	struct ieee_ext *py = (struct ieee_ext *)&y;
+
+	px->ext_sign = py->ext_sign;
+
+	return x;
+}
diff --git a/lib/libm/src/s_cos.c b/lib/libm/src/s_cos.c
index 54bf0e07d39..fa46074cfe5 100644
--- a/lib/libm/src/s_cos.c
+++ b/lib/libm/src/s_cos.c
@@ -45,7 +45,10 @@ static char rcsid[] = "$NetBSD: s_cos.c,v 1.7 1995/05/10 20:47:02 jtc Exp $";
  *	TRIG(x) returns trig(x) nearly rounded 
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 double
@@ -76,3 +79,9 @@ cos(double x)
 	    }
 	}
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(cosl, cos);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_cosl.c b/lib/libm/src/s_cosl.c
new file mode 100644
index 00000000000..6104751b6d5
--- /dev/null
+++ b/lib/libm/src/s_cosl.c
@@ -0,0 +1,116 @@
+/*	$OpenBSD: s_cosl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*-
+ * Copyright (c) 2007 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.
+ */
+
+/*
+ * Compute cos(x) for x where x is reduced to y = x - k * pi / 2.
+ * Limited testing on pseudorandom numbers drawn within [-2e8:4e8] shows
+ * an accuracy of <= 0.7412 ULP.
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+#include <math.h>
+
+#include "math_private.h"
+
+#if LDBL_MANT_DIG == 64
+#define	NX	3
+#define	PREC	2
+#elif LDBL_MANT_DIG == 113
+#define	NX	5
+#define	PREC	3
+#else
+#error "Unsupported long double format"
+#endif
+
+static const long double two24 = 1.67772160000000000000e+07L;
+
+long double
+cosl(long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} z;
+	int i, e0;
+	double xd[NX], yd[PREC];
+	long double hi, lo;
+
+	z.e = x;
+	z.bits.ext_sign = 0;
+
+	/* If x = +-0 or x is a subnormal number, then cos(x) = 1 */
+	if (z.bits.ext_exp == 0)
+		return (1.0);
+
+	/* If x = NaN or Inf, then cos(x) = NaN. */
+	if (z.bits.ext_exp == 32767)
+		return ((x - x) / (x - x));
+
+	/* Optimize the case where x is already within range. */
+	if (z.e < M_PI_4)
+		return (__kernel_cosl(z.e, 0));
+
+	/* Split z.e into a 24-bit representation. */
+	e0 = ilogbl(z.e) - 23;
+	z.e = scalbnl(z.e, -e0);
+	for (i = 0; i < NX; i++) {
+		xd[i] = (double)((int32_t)z.e);
+		z.e = (z.e - xd[i]) * two24;
+	}
+	
+	/* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */
+	e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC);
+
+#if PREC == 2
+	hi = (long double)yd[0] + yd[1];
+	lo = yd[1] - (hi - yd[0]);
+#else /* PREC == 3 */
+	long double t;
+	t = (long double)yd[2] + yd[1];
+	hi = t + yd[0];
+	lo = yd[0] - (hi - t);
+#endif
+
+	switch (e0 & 3) {
+	case 0:
+	    hi = __kernel_cosl(hi, lo);
+	    break;
+	case 1:
+	    hi = - __kernel_sinl(hi, lo, 1);
+	    break;
+	case 2:
+	    hi = - __kernel_cosl(hi, lo);
+	    break;
+	case 3:
+	    hi = __kernel_sinl(hi, lo, 1);
+	    break;
+	}
+	
+	return (hi);
+}
diff --git a/lib/libm/src/s_exp2.c b/lib/libm/src/s_exp2.c
index 0a864a49940..7359265f86d 100644
--- a/lib/libm/src/s_exp2.c
+++ b/lib/libm/src/s_exp2.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: s_exp2.c,v 1.1 2008/07/24 09:40:16 martynas Exp $	*/
+/*	$OpenBSD: s_exp2.c,v 1.2 2008/12/09 20:00:35 martynas Exp $	*/
 /*-
  * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
  * All rights reserved.
@@ -27,8 +27,8 @@
 
 #include <sys/cdefs.h>
 #include <float.h>
+#include <math.h>
 
-#include "math.h"
 #include "math_private.h"
 
 #define	TBLBITS	8
@@ -389,3 +389,9 @@ exp2(double x)
 		return (r * twopkp1000 * twom1000);
 	}
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(exp2l, exp2);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_fabs.c b/lib/libm/src/s_fabs.c
index c5332e60d2c..8cbd5f2d4f2 100644
--- a/lib/libm/src/s_fabs.c
+++ b/lib/libm/src/s_fabs.c
@@ -18,7 +18,10 @@ static char rcsid[] = "$NetBSD: s_fabs.c,v 1.7 1995/05/10 20:47:13 jtc Exp $";
  * fabs(x) returns the absolute value of x.
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 double
@@ -29,3 +32,9 @@ fabs(double x)
 	SET_HIGH_WORD(x,high&0x7fffffff);
         return x;
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(fabsl, fabs);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_fabsl.c b/lib/libm/src/s_fabsl.c
new file mode 100644
index 00000000000..f4cdeecf51c
--- /dev/null
+++ b/lib/libm/src/s_fabsl.c
@@ -0,0 +1,30 @@
+/*	$OpenBSD: s_fabsl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*
+ * Copyright (c) 2008 Martynas Venckus <martynas@openbsd.org>
+ *
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose with or without fee is hereby granted, provided that the above
+ * copyright notice and this permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+ * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+ * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+ * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <math.h>
+
+long double
+fabsl(long double x)
+{
+	struct ieee_ext *p = (struct ieee_ext *)&x;
+
+	p->ext_sign = 0;
+
+	return x;
+}
diff --git a/lib/libm/src/s_fdim.c b/lib/libm/src/s_fdim.c
index 6af4ba038b9..1e071a76dd9 100644
--- a/lib/libm/src/s_fdim.c
+++ b/lib/libm/src/s_fdim.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: s_fdim.c,v 1.1 2008/09/07 20:36:09 martynas Exp $	*/
+/*	$OpenBSD: s_fdim.c,v 1.2 2008/12/09 20:00:35 martynas Exp $	*/
 /*-
  * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
  * All rights reserved.
@@ -25,6 +25,8 @@
  * SUCH DAMAGE.
  */
 
+#include <machine/cdefs.h>
+#include <float.h>
 #include <math.h>
 
 #define	DECL(type, fn)			\
@@ -41,4 +43,9 @@ fn(type x, type y)			\
 
 DECL(double, fdim)
 DECL(float, fdimf)
-DECL(long double, fdiml)
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(fdiml, fdim);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_finite.c b/lib/libm/src/s_finite.c
deleted file mode 100644
index a1eea501193..00000000000
--- a/lib/libm/src/s_finite.c
+++ /dev/null
@@ -1,31 +0,0 @@
-/* @(#)s_finite.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.
- * ====================================================
- */
-
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_finite.c,v 1.8 1995/05/10 20:47:17 jtc Exp $";
-#endif
-
-/*
- * finite(x) returns 1 if x is finite, else 0;
- * no branching!
- */
-
-#include "math.h"
-#include "math_private.h"
-
-int
-finite(double x)
-{
-	int32_t hx;
-	GET_HIGH_WORD(hx,x);
-	return (int)((u_int32_t)((hx&0x7fffffff)-0x7ff00000)>>31);
-}
diff --git a/lib/libm/src/s_finitef.c b/lib/libm/src/s_finitef.c
deleted file mode 100644
index 09295830d60..00000000000
--- a/lib/libm/src/s_finitef.c
+++ /dev/null
@@ -1,34 +0,0 @@
-/* s_finitef.c -- float version of s_finite.c.
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-
-/*
- * ====================================================
- * 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.
- * ====================================================
- */
-
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_finitef.c,v 1.4 1995/05/10 20:47:18 jtc Exp $";
-#endif
-
-/*
- * finitef(x) returns 1 if x is finite, else 0;
- * no branching!
- */
-
-#include "math.h"
-#include "math_private.h"
-
-int
-finitef(float x)
-{
-	int32_t ix;
-	GET_FLOAT_WORD(ix,x);
-	return (int)((u_int32_t)((ix&0x7fffffff)-0x7f800000)>>31);
-}
diff --git a/lib/libm/src/s_fmax.c b/lib/libm/src/s_fmax.c
index 7ea5b8dc4db..436a898b143 100644
--- a/lib/libm/src/s_fmax.c
+++ b/lib/libm/src/s_fmax.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: s_fmax.c,v 1.2 2008/09/11 19:18:12 martynas Exp $	*/
+/*	$OpenBSD: s_fmax.c,v 1.3 2008/12/09 20:00:35 martynas Exp $	*/
 /*-
  * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
  * All rights reserved.
@@ -25,6 +25,8 @@
  * SUCH DAMAGE.
  */
 
+#include <machine/cdefs.h>
+#include <float.h>
 #include <math.h>
 
 double
@@ -45,3 +47,9 @@ fmax(double x, double y)
 
 	return (x > y ? x : y);
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(fmaxl, fmax);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_fmaxl.c b/lib/libm/src/s_fmaxl.c
new file mode 100644
index 00000000000..0968f084d2b
--- /dev/null
+++ b/lib/libm/src/s_fmaxl.c
@@ -0,0 +1,47 @@
+/*	$OpenBSD: s_fmaxl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*-
+ * Copyright (c) 2004 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 <math.h>
+
+long double
+fmaxl(long double x, long double y)
+{
+	/* Check for NaNs to avoid raising spurious exceptions. */
+	if (isnan(x))
+		return (y);
+	if (isnan(y))
+		return (x);
+
+	/* Handle comparisons of signed zeroes. */
+	if (signbit(x) != signbit(y))
+		if (signbit(x))
+			return (y);
+		else
+			return (x);
+
+	return (x > y ? x : y);
+}
diff --git a/lib/libm/src/s_fmin.c b/lib/libm/src/s_fmin.c
index 7309a6ffd98..051f4e62d20 100644
--- a/lib/libm/src/s_fmin.c
+++ b/lib/libm/src/s_fmin.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: s_fmin.c,v 1.2 2008/09/11 19:18:12 martynas Exp $	*/
+/*	$OpenBSD: s_fmin.c,v 1.3 2008/12/09 20:00:35 martynas Exp $	*/
 /*-
  * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
  * All rights reserved.
@@ -25,6 +25,8 @@
  * SUCH DAMAGE.
  */
 
+#include <machine/cdefs.h>
+#include <float.h>
 #include <math.h>
 
 double
@@ -45,3 +47,9 @@ fmin(double x, double y)
 
 	return (x < y ? x : y);
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(fminl, fmin);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_fminl.c b/lib/libm/src/s_fminl.c
new file mode 100644
index 00000000000..4bdd92d94c6
--- /dev/null
+++ b/lib/libm/src/s_fminl.c
@@ -0,0 +1,47 @@
+/*	$OpenBSD: s_fminl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*-
+ * Copyright (c) 2004 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 <math.h>
+
+long double
+fminl(long double x, long double y)
+{
+	/* Check for NaNs to avoid raising spurious exceptions. */
+	if (isnan(x))
+		return (y);
+	if (isnan(y))
+		return (x);
+
+	/* Handle comparisons of signed zeroes. */
+	if (signbit(x) != signbit(y))
+		if (signbit(y))
+			return (y);
+		else
+			return (x);
+
+	return (x < y ? x : y);
+}
diff --git a/lib/libm/src/s_frexp.c b/lib/libm/src/s_frexp.c
index 3034855d4b3..f18515a5ba2 100644
--- a/lib/libm/src/s_frexp.c
+++ b/lib/libm/src/s_frexp.c
@@ -24,7 +24,10 @@ static char rcsid[] = "$NetBSD: s_frexp.c,v 1.9 1995/05/10 20:47:24 jtc Exp $";
  * with *exp=0. 
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 static const double
@@ -49,3 +52,9 @@ frexp(double x, int *eptr)
 	SET_HIGH_WORD(x,hx);
 	return x;
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(frexpl, frexp);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_frexpl.c b/lib/libm/src/s_frexpl.c
new file mode 100644
index 00000000000..9093af12bd8
--- /dev/null
+++ b/lib/libm/src/s_frexpl.c
@@ -0,0 +1,70 @@
+/*	$OpenBSD: s_frexpl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*-
+ * Copyright (c) 2004-2005 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/src/s_frexpl.c,v 1.1 2005/03/07 04:54:51 das Exp $
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+#include <math.h>
+
+#if LDBL_MAX_EXP != 0x4000
+#error "Unsupported long double format"
+#endif
+
+long double
+frexpl(long double x, int *ex)
+{
+	struct ieee_ext *p = (struct ieee_ext *)&x;
+
+	switch (p->ext_exp) {
+	case 0:		/* 0 or subnormal */
+		if ((p->ext_fracl
+#ifdef EXT_FRACLMBITS
+			| p->ext_fraclm
+#endif /* EXT_FRACLMBITS */
+#ifdef EXT_FRACHMBITS
+			| p->ext_frachm
+#endif /* EXT_FRACHMBITS */
+			| p->ext_frach) == 0) {
+			*ex = 0;
+		} else {
+			x *= 0x1.0p514;
+			*ex = p->ext_exp - 0x4200;
+			p->ext_exp = 0x3ffe;
+		}
+		break;
+	case 0x7fff:	/* infinity or NaN; value of *ex is unspecified */
+		break;
+	default:	/* normal */
+		*ex = p->ext_exp - 0x3ffe;
+		p->ext_exp = 0x3ffe;
+		break;
+	}
+
+	return x;
+}
diff --git a/lib/libm/src/s_ilogb.c b/lib/libm/src/s_ilogb.c
index 9c4ce207b0b..a843eff9480 100644
--- a/lib/libm/src/s_ilogb.c
+++ b/lib/libm/src/s_ilogb.c
@@ -20,7 +20,10 @@ static char rcsid[] = "$NetBSD: s_ilogb.c,v 1.9 1995/05/10 20:47:28 jtc Exp $";
  * ilogb(inf/NaN) = 0x7fffffff (no signal is raised)
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 int
@@ -45,3 +48,9 @@ ilogb(double x)
 	else if (hx<0x7ff00000) return (hx>>20)-1023;
 	else return 0x7fffffff;
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(ilogbl, ilogb);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_ilogbl.c b/lib/libm/src/s_ilogbl.c
new file mode 100644
index 00000000000..da06cc9e6d8
--- /dev/null
+++ b/lib/libm/src/s_ilogbl.c
@@ -0,0 +1,78 @@
+/*	$OpenBSD: s_ilogbl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*
+ * From: @(#)s_ilogb.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/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+
+int
+ilogbl(long double x)
+{
+	struct ieee_ext *p = (struct ieee_ext *)&x;
+	unsigned long m;
+	int b;
+
+	if (p->ext_exp == 0) {
+		if ((p->ext_fracl
+#ifdef EXT_FRACLMBITS
+			| p->ext_fraclm
+#endif /* EXT_FRACLMBITS */
+#ifdef EXT_FRACHMBITS
+			| p->ext_frachm
+#endif /* EXT_FRACHMBITS */
+			| p->ext_frach) == 0)
+			return (FP_ILOGB0);
+		/* denormalized */
+		if (p->ext_frach == 0
+#ifdef EXT_FRACHMBITS
+			&& p->ext_frachm == 0
+#endif
+			) {
+			m = 1lu << (EXT_FRACLBITS - 1);
+			for (b = EXT_FRACHBITS; !(p->ext_fracl & m); m >>= 1)
+				b++;
+#if defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS)
+			m = 1lu << (EXT_FRACLMBITS - 1);
+			for (b += EXT_FRACHMBITS; !(p->ext_fraclm & m); m >>= 1)
+				b++;
+#endif /* defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS) */
+		} else {
+			m = 1lu << (EXT_FRACHBITS - 1);
+			for (b = 0; !(p->ext_frach & m); m >>= 1)
+				b++;
+#ifdef EXT_FRACHMBITS
+			m = 1lu << (EXT_FRACHMBITS - 1);
+			for (; !(p->ext_frachm & m); m >>= 1)
+				b++;
+#endif /* EXT_FRACHMBITS */
+		}
+#ifdef EXT_IMPLICIT_NBIT
+		b++;
+#endif
+		return (LDBL_MIN_EXP - b - 1);
+	} else if (p->ext_exp < (LDBL_MAX_EXP << 1) - 1)
+		return (p->ext_exp - LDBL_MAX_EXP + 1);
+	else if (p->ext_fracl != 0
+#ifdef EXT_FRACLMBITS
+		|| p->ext_fraclm != 0
+#endif /* EXT_FRACLMBITS */
+#ifdef EXT_FRACHMBITS
+		|| p->ext_frachm != 0
+#endif /* EXT_FRACHMBITS */
+		|| p->ext_frach != 0)
+		return (FP_ILOGBNAN);
+	else
+		return (INT_MAX);
+}
diff --git a/lib/libm/src/s_infinity.c b/lib/libm/src/s_infinity.c
deleted file mode 100644
index a24c1b38c5f..00000000000
--- a/lib/libm/src/s_infinity.c
+++ /dev/null
@@ -1,12 +0,0 @@
-/*
- * Written by J.T. Conklin <jtc@netbsd.org>.
- * Public domain.
- */
-
-#include <machine/endian.h>
-
-#if BYTE_ORDER == LITTLE_ENDIAN
-char __infinity[] = { 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xf0, 0x7f };
-#else
-char __infinity[] = { 0x7f, 0xf0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00 };
-#endif
diff --git a/lib/libm/src/s_ldexp.c b/lib/libm/src/s_ldexp.c
index aeef0ba9b62..cfbe791f1a3 100644
--- a/lib/libm/src/s_ldexp.c
+++ b/lib/libm/src/s_ldexp.c
@@ -14,9 +14,12 @@
 static char rcsid[] = "$NetBSD: s_ldexp.c,v 1.6 1995/05/10 20:47:40 jtc Exp $";
 #endif
 
-#include "math.h"
-#include "math_private.h"
+#include <machine/cdefs.h>
 #include <errno.h>
+#include <float.h>
+#include <math.h>
+
+#include "math_private.h"
 
 double
 ldexp(double value, int exp)
@@ -26,3 +29,9 @@ ldexp(double value, int exp)
 	if(!finite(value)||value==0.0) errno = ERANGE;
 	return value;
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(ldexpl, ldexp);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_logbl.c b/lib/libm/src/s_logbl.c
new file mode 100644
index 00000000000..7de6ede4fd0
--- /dev/null
+++ b/lib/libm/src/s_logbl.c
@@ -0,0 +1,77 @@
+/*	$OpenBSD: s_logbl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*
+ * From: @(#)s_ilogb.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/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+
+long double
+logbl(long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} u;
+	unsigned long m;
+	int b;
+
+	u.e = x;
+	if (u.bits.ext_exp == 0) {
+		if ((u.bits.ext_fracl
+#ifdef EXT_FRACLMBITS
+			| u.bits.ext_fraclm
+#endif /* EXT_FRACLMBITS */
+#ifdef EXT_FRACHMBITS
+			| u.bits.ext_frachm
+#endif /* EXT_FRACHMBITS */
+			| u.bits.ext_frach) == 0) {	/* x == 0 */
+			u.bits.ext_sign = 1;
+			return (1.0L / u.e);
+		}
+		/* denormalized */
+		if (u.bits.ext_frach == 0
+#ifdef EXT_FRACHMBITS
+			&& u.bits.ext_frachm == 0
+#endif
+			) {
+			m = 1lu << (EXT_FRACLBITS - 1);
+			for (b = EXT_FRACHBITS; !(u.bits.ext_fracl & m); m >>= 1)
+				b++;
+#if defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS)
+			m = 1lu << (EXT_FRACLMBITS - 1);
+			for (b += EXT_FRACHMBITS; !(u.bits.ext_fraclm & m);
+				m >>= 1)
+				b++;
+#endif /* defined(EXT_FRACHMBITS) && defined(EXT_FRACLMBITS) */
+		} else {
+			m = 1lu << (EXT_FRACHBITS - 1);
+			for (b = 0; !(u.bits.ext_frach & m); m >>= 1)
+				b++;
+#ifdef EXT_FRACHMBITS
+			m = 1lu << (EXT_FRACHMBITS - 1);
+			for (; !(u.bits.ext_frachm & m); m >>= 1)
+				b++;
+#endif /* EXT_FRACHMBITS */
+		}
+#ifdef EXT_IMPLICIT_NBIT
+		b++;
+#endif
+		return ((long double)(LDBL_MIN_EXP - b - 1));
+	}
+	if (u.bits.ext_exp < (LDBL_MAX_EXP << 1) - 1)	/* normal */
+		return ((long double)(u.bits.ext_exp - LDBL_MAX_EXP + 1));
+	else						/* +/- inf or nan */
+		return (x * x);
+}
diff --git a/lib/libm/src/s_nan.c b/lib/libm/src/s_nan.c
index 3e99cfda129..99b55590313 100644
--- a/lib/libm/src/s_nan.c
+++ b/lib/libm/src/s_nan.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: s_nan.c,v 1.4 2008/09/07 20:36:09 martynas Exp $	*/
+/*	$OpenBSD: s_nan.c,v 1.5 2008/12/09 20:00:35 martynas Exp $	*/
 /*-
  * Copyright (c) 2007 David Schultz
  * All rights reserved.
@@ -27,6 +27,7 @@
  * $FreeBSD: src/lib/msun/src/s_nan.c,v 1.2 2007/12/18 23:46:32 das Exp $
  */
 
+#include <machine/cdefs.h>
 #include <sys/types.h>
 #include <sys/endian.h>
 #include <ctype.h>
@@ -121,3 +122,9 @@ nanf(const char *s)
 	u.bits[0] |= 0x7fc00000;
 	return (u.f);
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(nanl, nan);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_rint.c b/lib/libm/src/s_rint.c
index 192a765d7d0..262a8fab59e 100644
--- a/lib/libm/src/s_rint.c
+++ b/lib/libm/src/s_rint.c
@@ -24,7 +24,10 @@ static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $";
  *	Inexact flag raised if x not equal to rint(x).
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 static const double
@@ -76,3 +79,9 @@ rint(double x)
 	w = TWO52[sx]+x;
 	return w-TWO52[sx];
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(rintl, rint);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_rintl.c b/lib/libm/src/s_rintl.c
new file mode 100644
index 00000000000..55ebca5ba0e
--- /dev/null
+++ b/lib/libm/src/s_rintl.c
@@ -0,0 +1,91 @@
+/*	$OpenBSD: s_rintl.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 <sys/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+#include <math.h>
+
+#if LDBL_MAX_EXP != 0x4000
+/* We also require the usual bias, min exp and expsign packing. */
+#error "Unsupported long double format"
+#endif
+
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+static const float
+shift[2] = {
+#if LDBL_MANT_DIG == 64
+	0x1.0p63, -0x1.0p63
+#elif LDBL_MANT_DIG == 113
+	0x1.0p112, -0x1.0p112
+#else
+#error "Unsupported long double format"
+#endif
+};
+static const float zero[2] = { 0.0, -0.0 };
+
+long double
+rintl(long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} u;
+	uint32_t expsign;
+	int ex, sign;
+
+	u.e = x;
+	expsign = (u.bits.ext_sign << 15) | u.bits.ext_exp;
+	ex = expsign & 0x7fff;
+
+	if (ex >= BIAS + LDBL_MANT_DIG - 1) {
+		if (ex == BIAS + LDBL_MAX_EXP)
+			return (x + x);	/* Inf, NaN, or unsupported format */
+		return (x);		/* finite and already an integer */
+	}
+	sign = expsign >> 15;
+
+	/*
+	 * The following code assumes that intermediate results are
+	 * evaluated in long double precision. If they are evaluated in
+	 * greater precision, double rounding may occur, and if they are
+	 * evaluated in less precision (as on i386), results will be
+	 * wildly incorrect.
+	 */
+	x += shift[sign];
+	x -= shift[sign];
+
+	/*
+	 * If the result is +-0, then it must have the same sign as x, but
+	 * the above calculation doesn't always give this.  Fix up the sign.
+	 */
+	if (ex < BIAS && x == 0.0L)
+		return (zero[sign]);
+
+	return (x);
+}
diff --git a/lib/libm/src/s_scalbn.c b/lib/libm/src/s_scalbn.c
index 9444f51f0c0..4fbbe8f9c55 100644
--- a/lib/libm/src/s_scalbn.c
+++ b/lib/libm/src/s_scalbn.c
@@ -21,7 +21,10 @@ static char rcsid[] = "$NetBSD: s_scalbn.c,v 1.8 1995/05/10 20:48:08 jtc Exp $";
  * exponentiation or a multiplication.
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 static const double
@@ -56,3 +59,9 @@ scalbn (double x, int n)
 	SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
         return x*twom54;
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(scalbnl, scalbn);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_scalbnl.c b/lib/libm/src/s_scalbnl.c
new file mode 100644
index 00000000000..eda7a5e2333
--- /dev/null
+++ b/lib/libm/src/s_scalbnl.c
@@ -0,0 +1,82 @@
+/*	$OpenBSD: s_scalbnl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/* @(#)s_scalbn.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.
+ * ====================================================
+ */
+
+/*
+ * scalbnl (long double x, int n)
+ * scalbnl(x,n) returns x* 2**n  computed by  exponent
+ * manipulation rather than by actually performing an
+ * exponentiation or a multiplication.
+ */
+
+/*
+ * We assume that a long double has a 15-bit exponent.  On systems
+ * where long double is the same as double, scalbnl() is an alias
+ * for scalbn(), so we don't use this routine.
+ */
+
+#include <sys/cdefs.h>
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+#include <math.h>
+
+#if LDBL_MAX_EXP != 0x4000
+#error "Unsupported long double format"
+#endif
+
+static const long double
+huge = 0x1p16000L,
+tiny = 0x1p-16000L;
+
+long double
+scalbnl (long double x, int n)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} u;
+	int k;
+	u.e = x;
+        k = u.bits.ext_exp;			/* extract exponent */
+        if (k==0) {				/* 0 or subnormal x */
+            if ((u.bits.ext_frach
+#ifdef EXT_FRACHMBITS
+		| u.bits.ext_frachm
+#endif /* EXT_FRACHMBITS */
+#ifdef EXT_FRACLMBITS
+		| u.bits.ext_fraclm
+#endif /* EXT_FRACLMBITS */
+		| u.bits.ext_fracl)==0) return x;	/* +-0 */
+	    u.e *= 0x1p+128;
+	    k = u.bits.ext_exp - 128;
+            if (n< -50000) return tiny*x; 	/*underflow*/
+	    }
+        if (k==0x7fff) return x+x;		/* NaN or Inf */
+        k = k+n;
+        if (k >= 0x7fff) return huge*copysignl(huge,x); /* overflow  */
+        if (k > 0) 				/* normal result */
+	    {u.bits.ext_exp = k; return u.e;}
+        if (k <= -128)
+            if (n > 50000) 	/* in case integer overflow in n+k */
+		return huge*copysign(huge,x);	/*overflow*/
+	    else return tiny*copysign(tiny,x); 	/*underflow*/
+        k += 128;				/* subnormal result */
+	u.bits.ext_exp = k;
+        return u.e*0x1p-128;
+}
+
+long double
+ldexpl(long double x, int n)
+{
+	return scalbnl(x, n);
+}
diff --git a/lib/libm/src/s_sin.c b/lib/libm/src/s_sin.c
index 1366556e45a..9e4b729f13a 100644
--- a/lib/libm/src/s_sin.c
+++ b/lib/libm/src/s_sin.c
@@ -45,7 +45,10 @@ static char rcsid[] = "$NetBSD: s_sin.c,v 1.7 1995/05/10 20:48:15 jtc Exp $";
  *	TRIG(x) returns trig(x) nearly rounded 
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 double
@@ -76,3 +79,9 @@ sin(double x)
 	    }
 	}
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(sinl, sin);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_sinl.c b/lib/libm/src/s_sinl.c
new file mode 100644
index 00000000000..dc940dd7837
--- /dev/null
+++ b/lib/libm/src/s_sinl.c
@@ -0,0 +1,117 @@
+/*	$OpenBSD: s_sinl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*-
+ * Copyright (c) 2007 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.
+ */
+
+/*
+ * Compute sin(x) for x where x is reduced to y = x - k * pi / 2.
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+#include <math.h>
+
+#include "math_private.h"
+
+#if LDBL_MANT_DIG == 64
+#define	NX	3
+#define	PREC	2
+#elif LDBL_MANT_DIG == 113
+#define	NX	5
+#define	PREC	3
+#else
+#error "Unsupported long double format"
+#endif
+
+static const long double two24 = 1.67772160000000000000e+07L;
+
+long double
+sinl(long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} z;
+	int i, e0, s;
+	double xd[NX], yd[PREC];
+	long double hi, lo;
+
+	z.e = x;
+	s = z.bits.ext_sign;
+	z.bits.ext_sign = 0;
+
+	/* If x = +-0 or x is a subnormal number, then sin(x) = x */
+	if (z.bits.ext_exp == 0)
+		return (x);
+
+	/* If x = NaN or Inf, then sin(x) = NaN. */
+	if (z.bits.ext_exp == 32767)
+		return ((x - x) / (x - x));
+
+	/* Optimize the case where x is already within range. */
+	if (z.e < M_PI_4) {
+		hi = __kernel_sinl(z.e, 0, 0);
+		return  (s ? -hi : hi);
+	}
+
+	/* Split z.e into a 24-bit representation. */
+	e0 = ilogbl(z.e) - 23;
+	z.e = scalbnl(z.e, -e0);
+	for (i = 0; i < NX; i++) {
+		xd[i] = (double)((int32_t)z.e);
+		z.e = (z.e - xd[i]) * two24;
+	}
+
+	/* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */
+	e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC);
+
+#if PREC == 2
+	hi = (long double)yd[0] + yd[1];
+	lo = yd[1] - (hi - yd[0]);
+#else /* PREC == 3 */
+	long double t;
+	t = (long double)yd[2] + yd[1];
+	hi = t + yd[0];
+	lo = yd[0] - (hi - t);
+#endif
+
+	switch (e0 & 3) {
+	case 0:
+	    hi = __kernel_sinl(hi, lo, 1);
+	    break;
+	case 1:
+	    hi = __kernel_cosl(hi, lo);
+	    break;
+	case 2:
+	    hi = - __kernel_sinl(hi, lo, 1);
+	    break;
+	case 3:
+	    hi = - __kernel_cosl(hi, lo);
+	    break;
+	}
+	
+	return (s ? -hi : hi);
+}
diff --git a/lib/libm/src/s_tan.c b/lib/libm/src/s_tan.c
index 2ad099dd535..8855e79d71d 100644
--- a/lib/libm/src/s_tan.c
+++ b/lib/libm/src/s_tan.c
@@ -44,7 +44,10 @@ static char rcsid[] = "$NetBSD: s_tan.c,v 1.7 1995/05/10 20:48:18 jtc Exp $";
  *	TRIG(x) returns trig(x) nearly rounded 
  */
 
-#include "math.h"
+#include <machine/cdefs.h>
+#include <float.h>
+#include <math.h>
+
 #include "math_private.h"
 
 double
@@ -70,3 +73,9 @@ tan(double x)
 							-1 -- n odd */
 	}
 }
+
+#if LDBL_MANT_DIG == 53
+#ifdef __weak_alias
+__weak_alias(tanl, tan);
+#endif /* __weak_alias */
+#endif /* LDBL_MANT_DIG == 53 */
diff --git a/lib/libm/src/s_tanl.c b/lib/libm/src/s_tanl.c
new file mode 100644
index 00000000000..c71e8d5c7a6
--- /dev/null
+++ b/lib/libm/src/s_tanl.c
@@ -0,0 +1,116 @@
+/*	$OpenBSD: s_tanl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $	*/
+/*-
+ * Copyright (c) 2007 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.
+ */
+
+/*
+ * Compute tan(x) for x where x is reduced to y = x - k * pi / 2.
+ * Limited testing on pseudorandom numbers drawn within [0:4e8] shows
+ * an accuracy of <= 1.5 ULP where 247024 values of x out of 40 million
+ * possibles resulted in tan(x) that exceeded 0.5 ULP (ie., 0.6%).
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+#include <float.h>
+#include <math.h>
+
+#include "math_private.h"
+
+#if LDBL_MANT_DIG == 64
+#define	NX	3
+#define	PREC	2
+#elif LDBL_MANT_DIG == 113
+#define	NX	5
+#define	PREC	3
+#else
+#error "Unsupported long double format"
+#endif
+
+static const long double two24 = 1.67772160000000000000e+07L;
+
+long double
+tanl(long double x)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} z;
+	int i, e0, s;
+	double xd[NX], yd[PREC];
+	long double hi, lo;
+
+	z.e = x;
+	s = z.bits.ext_sign;
+	z.bits.ext_sign = 0;
+
+	/* If x = +-0 or x is subnormal, then tan(x) = x. */
+	if (z.bits.ext_exp == 0)
+		return (x);
+
+	/* If x = NaN or Inf, then tan(x) = NaN. */
+	if (z.bits.ext_exp == 32767)
+		return ((x - x) / (x - x));
+
+	/* Optimize the case where x is already within range. */
+	if (z.e < M_PI_4) {
+		hi = __kernel_tanl(z.e, 0, 0);
+		return (s ? -hi : hi);
+	}
+
+	/* Split z.e into a 24-bit representation. */
+	e0 = ilogbl(z.e) - 23;
+	z.e = scalbnl(z.e, -e0);
+	for (i = 0; i < NX; i++) {
+		xd[i] = (double)((int32_t)z.e);
+		z.e = (z.e - xd[i]) * two24;
+	}
+	
+	/* yd contains the pieces of xd rem pi/2 such that |yd| < pi/4. */
+	e0 = __kernel_rem_pio2(xd, yd, e0, NX, PREC);
+
+#if PREC == 2
+	hi = (long double)yd[0] + yd[1];
+	lo = yd[1] - (hi - yd[0]);
+#else /* PREC == 3 */
+	long double t;
+	t = (long double)yd[2] + yd[1];
+	hi = t + yd[0];
+	lo = yd[0] - (hi - t);
+#endif
+
+	switch (e0 & 3) {
+	case 0:
+	case 2:
+	    hi = __kernel_tanl(hi, lo, 0);
+	    break;
+	case 1:
+	case 3:
+	    hi = __kernel_tanl(hi, lo, 1);
+	    break;
+	}
+
+	return (s ? -hi : hi);
+}