Meta: Add OpenLibm

2025-12-16 16:09:53 +02:00 · 2025-04-04 18:16:29 -04:00
parent 4a72f35e9b
commit 172b4ce7a3
724 changed files with 51312 additions and 0 deletions
--- a/src/openlibm/ld80/Make.files
+++ b/src/openlibm/ld80/Make.files
@@ -0,0 +1,13 @@
+$(CUR_SRCS) += 	invtrig.c \
+            e_acoshl.c     e_powl.c       k_tanl.c       s_exp2l.c \
+            e_atanhl.c     e_lgammal_r.c  e_sinhl.c      s_asinhl.c     s_expm1l.c \
+            e_coshl.c      e_log10l.c     e_tgammal.c \
+            e_expl.c       e_log2l.c      k_cosl.c       s_log1pl.c     s_tanhl.c \
+            e_logl.c       k_sinl.c       s_erfl.c
+
+#           s_remquol.c    e_fmodl.c      s_truncl.c
+#           e_hypotl.c     s_floorl.c     s_nextafterl.c s_ceill.c      s_modfl.c
+
+ifneq ($(OS), WINNT)
+$(CUR_SRCS) += s_nanl.c
+endif
--- a/src/openlibm/ld80/e_acoshl.c
+++ b/src/openlibm/ld80/e_acoshl.c
@@ -0,0 +1,57 @@
+/* @(#)e_acosh.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.
+ * ====================================================
+ */
+
+/* acoshl(x)
+ * Method :
+ *	Based on
+ *		acoshl(x) = logl [ x + sqrtl(x*x-1) ]
+ *	we have
+ *		acoshl(x) := logl(x)+ln2,	if x is large; else
+ *		acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else
+ *		acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1.
+ *
+ * Special cases:
+ *	acoshl(x) is NaN with signal if x<1.
+ *	acoshl(NaN) is NaN without signal.
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static const long double
+one	= 1.0,
+ln2	= 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
+
+long double
+acoshl(long double x)
+{
+	long double t;
+	u_int32_t se,i0,i1;
+	GET_LDOUBLE_WORDS(se,i0,i1,x);
+	if(se<0x3fff || se & 0x8000) {	/* x < 1 */
+	    return (x-x)/(x-x);
+	} else if(se >=0x401d) {	/* x > 2**30 */
+	    if(se >=0x7fff) {		/* x is inf of NaN */
+		return x+x;
+	    } else
+		return logl(x)+ln2;	/* acoshl(huge)=logl(2x) */
+	} else if(((se-0x3fff)|i0|i1)==0) {
+	    return 0.0;			/* acosh(1) = 0 */
+	} else if (se > 0x4000) {	/* 2**28 > x > 2 */
+	    t=x*x;
+	    return logl(2.0*x-one/(x+sqrtl(t-one)));
+	} else {			/* 1<x<2 */
+	    t = x-one;
+	    return log1pl(t+sqrtl(2.0*t+t*t));
+	}
+}
--- a/src/openlibm/ld80/e_acoshl.c.o
+++ b/src/openlibm/ld80/e_acoshl.c.o
--- a/src/openlibm/ld80/e_atanhl.c
+++ b/src/openlibm/ld80/e_atanhl.c
@@ -0,0 +1,60 @@
+/* @(#)e_atanh.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.
+ * ====================================================
+ */
+
+/* atanhl(x)
+ * Method :
+ *    1.Reduced x to positive by atanh(-x) = -atanh(x)
+ *    2.For x>=0.5
+ *                   1              2x                          x
+ *	atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
+ *                   2             1 - x                      1 - x
+ *
+ * 	For x<0.5
+ *	atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x))
+ *
+ * Special cases:
+ *	atanhl(x) is NaN if |x| > 1 with signal;
+ *	atanhl(NaN) is that NaN with no signal;
+ *	atanhl(+-1) is +-INF with signal.
+ *
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static const long double one = 1.0, huge = 1e4900L;
+
+static const long double zero = 0.0;
+
+long double
+atanhl(long double x)
+{
+	long double t;
+	int32_t ix;
+	u_int32_t se,i0,i1;
+	GET_LDOUBLE_WORDS(se,i0,i1,x);
+	ix = se&0x7fff;
+	if ((ix+((((i0&0x7fffffff)|i1)|(-((i0&0x7fffffff)|i1)))>>31))>0x3fff)
+	  /* |x|>1 */
+	    return (x-x)/(x-x);
+	if(ix==0x3fff)
+	    return x/zero;
+	if(ix<0x3fe3&&(huge+x)>zero) return x;	/* x<2**-28 */
+	SET_LDOUBLE_EXP(x,ix);
+	if(ix<0x3ffe) {		/* x < 0.5 */
+	    t = x+x;
+	    t = 0.5*log1pl(t+t*x/(one-x));
+	} else
+	    t = 0.5*log1pl((x+x)/(one-x));
+	if(se<=0x7fff) return t; else return -t;
+}
--- a/src/openlibm/ld80/e_atanhl.c.o
+++ b/src/openlibm/ld80/e_atanhl.c.o
--- a/src/openlibm/ld80/e_coshl.c
+++ b/src/openlibm/ld80/e_coshl.c
@@ -0,0 +1,83 @@
+/* @(#)e_cosh.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.
+ * ====================================================
+ */
+
+/* coshl(x)
+ * Method :
+ * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
+ *	1. Replace x by |x| (coshl(x) = coshl(-x)).
+ *	2.
+ *		                                        [ exp(x) - 1 ]^2
+ *	    0        <= x <= ln2/2  :  coshl(x) := 1 + -------------------
+ *			       			           2*exp(x)
+ *
+ *		                                   exp(x) +  1/exp(x)
+ *	    ln2/2    <= x <= 22     :  coshl(x) := -------------------
+ *			       			           2
+ *	    22       <= x <= lnovft :  coshl(x) := expl(x)/2
+ *	    lnovft   <= x <= ln2ovft:  coshl(x) := expl(x/2)/2 * expl(x/2)
+ *	    ln2ovft  <  x	    :  coshl(x) := huge*huge (overflow)
+ *
+ * Special cases:
+ *	coshl(x) is |x| if x is +INF, -INF, or NaN.
+ *	only coshl(0)=1 is exact for finite x.
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static const long double one = 1.0, half=0.5, huge = 1.0e4900L;
+
+long double
+coshl(long double x)
+{
+	long double t,w;
+	int32_t ex;
+	u_int32_t mx,lx;
+
+    /* High word of |x|. */
+	GET_LDOUBLE_WORDS(ex,mx,lx,x);
+	ex &= 0x7fff;
+
+    /* x is INF or NaN */
+	if(ex==0x7fff) return x*x;
+
+    /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
+	if(ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) {
+	    t = expm1l(fabsl(x));
+	    w = one+t;
+	    if (ex<0x3fbc) return w;	/* cosh(tiny) = 1 */
+	    return one+(t*t)/(w+w);
+	}
+
+    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
+	if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) {
+		t = expl(fabsl(x));
+		return half*t+half/t;
+	}
+
+    /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
+	if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u))
+		return half*expl(fabsl(x));
+
+    /* |x| in [log(maxdouble), log(2*maxdouble)) */
+	if (ex == 0x400c && (mx < 0xb174ddc0u
+			     || (mx == 0xb174ddc0u && lx < 0x31aec0ebu)))
+	{
+	    w = expl(half*fabsl(x));
+	    t = half*w;
+	    return t*w;
+	}
+
+    /* |x| >= log(2*maxdouble), cosh(x) overflow */
+	return huge*huge;
+}
--- a/src/openlibm/ld80/e_coshl.c.o
+++ b/src/openlibm/ld80/e_coshl.c.o
--- a/src/openlibm/ld80/e_expl.c
+++ b/src/openlibm/ld80/e_expl.c
@@ -0,0 +1,131 @@
+/*	$OpenBSD: e_expl.c,v 1.3 2013/11/12 20:35:19 martynas Exp $	*/
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * 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.
+ */
+
+/*							expl.c
+ *
+ *	Exponential function, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, expl();
+ *
+ * y = expl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns e (2.71828...) raised to the x power.
+ *
+ * Range reduction is accomplished by separating the argument
+ * into an integer k and fraction f such that
+ *
+ *     x    k  f
+ *    e  = 2  e.
+ *
+ * A Pade' form of degree 2/3 is used to approximate exp(f) - 1
+ * in the basic range [-0.5 ln 2, 0.5 ln 2].
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      +-10000     50000       1.12e-19    2.81e-20
+ *
+ *
+ * Error amplification in the exponential function can be
+ * a serious matter.  The error propagation involves
+ * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
+ * which shows that a 1 lsb error in representing X produces
+ * a relative error of X times 1 lsb in the function.
+ * While the routine gives an accurate result for arguments
+ * that are exactly represented by a long double precision
+ * computer number, the result contains amplified roundoff
+ * error for large arguments not exactly represented.
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * exp underflow    x < MINLOG         0.0
+ * exp overflow     x > MAXLOG         MAXNUM
+ *
+ */
+
+/*	Exponential function	*/
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static long double P[3] = {
+ 1.2617719307481059087798E-4L,
+ 3.0299440770744196129956E-2L,
+ 9.9999999999999999991025E-1L,
+};
+static long double Q[4] = {
+ 3.0019850513866445504159E-6L,
+ 2.5244834034968410419224E-3L,
+ 2.2726554820815502876593E-1L,
+ 2.0000000000000000000897E0L,
+};
+static const long double C1 = 6.9314575195312500000000E-1L;
+static const long double C2 = 1.4286068203094172321215E-6L;
+static const long double MAXLOGL = 1.1356523406294143949492E4L;
+static const long double MINLOGL = -1.13994985314888605586758E4L;
+static const long double LOG2EL = 1.4426950408889634073599E0L;
+
+long double
+expl(long double x)
+{
+long double px, xx;
+int n;
+
+if( isnan(x) )
+	return(x);
+if( x > MAXLOGL)
+	return( INFINITY );
+
+if( x < MINLOGL )
+	return(0.0L);
+
+/* Express e**x = e**g 2**n
+ *   = e**g e**( n loge(2) )
+ *   = e**( g + n loge(2) )
+ */
+px = floorl( LOG2EL * x + 0.5L ); /* floor() truncates toward -infinity. */
+n = px;
+x -= px * C1;
+x -= px * C2;
+
+
+/* rational approximation for exponential
+ * of the fractional part:
+ * e**x =  1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
+ */
+xx = x * x;
+px = x * __polevll( xx, P, 2 );
+x =  px/( __polevll( xx, Q, 3 ) - px );
+x = 1.0L + ldexpl( x, 1 );
+
+x = ldexpl( x, n );
+return(x);
+}
--- a/src/openlibm/ld80/e_expl.c.o
+++ b/src/openlibm/ld80/e_expl.c.o
--- a/src/openlibm/ld80/e_fmodl.c
+++ b/src/openlibm/ld80/e_fmodl.c
@@ -0,0 +1,142 @@
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/types.h>
+//#include <machine/ieee.h>
+
+#include <float.h>
+#include <openlibm_math.h>
+#include <stdint.h>
+
+#include "math_private.h"
+
+#define	BIAS (LDBL_MAX_EXP - 1)
+
+/*
+ * These macros add and remove an explicit integer bit in front of the
+ * fractional mantissa, if the architecture doesn't have such a bit by
+ * default already.
+ */
+#ifdef LDBL_IMPLICIT_NBIT
+#define	LDBL_NBIT	0
+#define	SET_NBIT(hx)	((hx) | (1ULL << LDBL_MANH_SIZE))
+#define	HFRAC_BITS	EXT_FRACHBITS
+#else
+#define	LDBL_NBIT	0x80000000
+#define	SET_NBIT(hx)	(hx)
+#define	HFRAC_BITS	(EXT_FRACHBITS - 1)
+#endif
+
+#define	MANL_SHIFT	(EXT_FRACLBITS - 1)
+
+static const long double one = 1.0, Zero[] = {0.0, -0.0,};
+
+/*
+ * fmodl(x,y)
+ * Return x mod y in exact arithmetic
+ * Method: shift and subtract
+ *
+ * Assumptions:
+ * - The low part of the mantissa fits in a manl_t exactly.
+ * - The high part of the mantissa fits in an int64_t with enough room
+ *   for an explicit integer bit in front of the fractional bits.
+ */
+long double
+fmodl(long double x, long double y)
+{
+	union {
+		long double e;
+		struct ieee_ext bits;
+	} ux, uy;
+	int64_t hx,hz;	/* We need a carry bit even if LDBL_MANH_SIZE is 32. */
+	uint32_t hy;
+	uint32_t lx,ly,lz;
+	int ix,iy,n,sx;
+
+	ux.e = x;
+	uy.e = y;
+	sx = ux.bits.ext_sign;
+
+    /* purge off exception values */
+	if((uy.bits.ext_exp|uy.bits.ext_frach|uy.bits.ext_fracl)==0 || /* y=0 */
+	   (ux.bits.ext_exp == BIAS + LDBL_MAX_EXP) ||	 /* or x not finite */
+	   (uy.bits.ext_exp == BIAS + LDBL_MAX_EXP &&
+	    ((uy.bits.ext_frach&~LDBL_NBIT)|uy.bits.ext_fracl)!=0)) /* or y is NaN */
+	    return (x*y)/(x*y);
+	if(ux.bits.ext_exp<=uy.bits.ext_exp) {
+	    if((ux.bits.ext_exp<uy.bits.ext_exp) ||
+	       (ux.bits.ext_frach<=uy.bits.ext_frach &&
+		(ux.bits.ext_frach<uy.bits.ext_frach ||
+		 ux.bits.ext_fracl<uy.bits.ext_fracl))) {
+		return x;		/* |x|<|y| return x or x-y */
+	    }
+	    if(ux.bits.ext_frach==uy.bits.ext_frach &&
+		ux.bits.ext_fracl==uy.bits.ext_fracl) {
+		return Zero[sx];	/* |x|=|y| return x*0*/
+	    }
+	}
+
+    /* determine ix = ilogb(x) */
+	if(ux.bits.ext_exp == 0) {	/* subnormal x */
+	    ux.e *= 0x1.0p512;
+	    ix = ux.bits.ext_exp - (BIAS + 512);
+	} else {
+	    ix = ux.bits.ext_exp - BIAS;
+	}
+
+    /* determine iy = ilogb(y) */
+	if(uy.bits.ext_exp == 0) {	/* subnormal y */
+	    uy.e *= 0x1.0p512;
+	    iy = uy.bits.ext_exp - (BIAS + 512);
+	} else {
+	    iy = uy.bits.ext_exp - BIAS;
+	}
+
+    /* set up {hx,lx}, {hy,ly} and align y to x */
+	hx = SET_NBIT(ux.bits.ext_frach);
+	hy = SET_NBIT(uy.bits.ext_frach);
+	lx = ux.bits.ext_fracl;
+	ly = uy.bits.ext_fracl;
+
+    /* fix point fmod */
+	n = ix - iy;
+
+	while(n--) {
+	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	    if(hz<0){hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;}
+	    else {
+		if ((hz|lz)==0)		/* return sign(x)*0 */
+		    return Zero[sx];
+		hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz;
+	    }
+	}
+	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	if(hz>=0) {hx=hz;lx=lz;}
+
+    /* convert back to floating value and restore the sign */
+	if((hx|lx)==0)			/* return sign(x)*0 */
+	    return Zero[sx];
+	while(hx<(1ULL<<HFRAC_BITS)) {	/* normalize x */
+	    hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;
+	    iy -= 1;
+	}
+	ux.bits.ext_frach = hx; /* The mantissa is truncated here if needed. */
+	ux.bits.ext_fracl = lx;
+	if (iy < LDBL_MIN_EXP) {
+	    ux.bits.ext_exp = iy + (BIAS + 512);
+	    ux.e *= 0x1p-512;
+	} else {
+	    ux.bits.ext_exp = iy + BIAS;
+	}
+	x = ux.e * one;		/* create necessary signal */
+	return x;		/* exact output */
+}
--- a/src/openlibm/ld80/e_hypotl.c
+++ b/src/openlibm/ld80/e_hypotl.c
@@ -0,0 +1,122 @@
+/* @(#)e_hypot.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.
+ * ====================================================
+ */
+
+/* hypotl(x,y)
+ *
+ * Method :
+ *	If (assume round-to-nearest) z=x*x+y*y
+ *	has error less than sqrt(2)/2 ulp, than
+ *	sqrt(z) has error less than 1 ulp (exercise).
+ *
+ *	So, compute sqrt(x*x+y*y) with some care as
+ *	follows to get the error below 1 ulp:
+ *
+ *	Assume x>y>0;
+ *	(if possible, set rounding to round-to-nearest)
+ *	1. if x > 2y  use
+ *		x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
+ *	where x1 = x with lower 32 bits cleared, x2 = x-x1; else
+ *	2. if x <= 2y use
+ *		t1*yy1+((x-y)*(x-y)+(t1*y2+t2*y))
+ *	where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
+ *	yy1= y with lower 32 bits chopped, y2 = y-yy1.
+ *
+ *	NOTE: scaling may be necessary if some argument is too
+ *	      large or too tiny
+ *
+ * Special cases:
+ *	hypot(x,y) is INF if x or y is +INF or -INF; else
+ *	hypot(x,y) is NAN if x or y is NAN.
+ *
+ * Accuracy:
+ * 	hypot(x,y) returns sqrt(x^2+y^2) with error less
+ * 	than 1 ulps (units in the last place)
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+long double
+hypotl(long double x, long double y)
+{
+	long double a,b,t1,t2,yy1,y2,w;
+	u_int32_t j,k,ea,eb;
+
+	GET_LDOUBLE_EXP(ea,x);
+	ea &= 0x7fff;
+	GET_LDOUBLE_EXP(eb,y);
+	eb &= 0x7fff;
+	if(eb > ea) {a=y;b=x;j=ea; ea=eb;eb=j;} else {a=x;b=y;}
+	SET_LDOUBLE_EXP(a,ea);	/* a <- |a| */
+	SET_LDOUBLE_EXP(b,eb);	/* b <- |b| */
+	if((ea-eb)>0x46) {return a+b;} /* x/y > 2**70 */
+	k=0;
+	if(ea > 0x5f3f) {	/* a>2**8000 */
+	   if(ea == 0x7fff) {	/* Inf or NaN */
+	       u_int32_t es,high,low;
+	       w = a+b;			/* for sNaN */
+	       GET_LDOUBLE_WORDS(es,high,low,a);
+	       if(((high&0x7fffffff)|low)==0) w = a;
+	       GET_LDOUBLE_WORDS(es,high,low,b);
+	       if(((eb^0x7fff)|(high&0x7fffffff)|low)==0) w = b;
+	       return w;
+	   }
+	   /* scale a and b by 2**-9600 */
+	   ea -= 0x2580; eb -= 0x2580;	k += 9600;
+	   SET_LDOUBLE_EXP(a,ea);
+	   SET_LDOUBLE_EXP(b,eb);
+	}
+	if(eb < 0x20bf) {	/* b < 2**-8000 */
+	    if(eb == 0) {	/* subnormal b or 0 */
+		u_int32_t es,high,low;
+		GET_LDOUBLE_WORDS(es,high,low,b);
+		if((high|low)==0) return a;
+		SET_LDOUBLE_WORDS(t1, 0x7ffd, 0, 0);	/* t1=2^16382 */
+		b *= t1;
+		a *= t1;
+		k -= 16382;
+	    } else {		/* scale a and b by 2^9600 */
+		ea += 0x2580;	/* a *= 2^9600 */
+		eb += 0x2580;	/* b *= 2^9600 */
+		k -= 9600;
+		SET_LDOUBLE_EXP(a,ea);
+		SET_LDOUBLE_EXP(b,eb);
+	    }
+	}
+    /* medium size a and b */
+	w = a-b;
+	if (w>b) {
+	    u_int32_t high;
+	    GET_LDOUBLE_MSW(high,a);
+	    SET_LDOUBLE_WORDS(t1,ea,high,0);
+	    t2 = a-t1;
+	    w  = sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));
+	} else {
+	    u_int32_t high;
+	    GET_LDOUBLE_MSW(high,b);
+	    a  = a+a;
+	    SET_LDOUBLE_WORDS(yy1,eb,high,0);
+	    y2 = b - yy1;
+	    GET_LDOUBLE_MSW(high,a);
+	    SET_LDOUBLE_WORDS(t1,ea+1,high,0);
+	    t2 = a - t1;
+	    w  = sqrtl(t1*yy1-(w*(-w)-(t1*y2+t2*b)));
+	}
+	if(k!=0) {
+	    u_int32_t es;
+	    t1 = 1.0;
+	    GET_LDOUBLE_EXP(es,t1);
+	    SET_LDOUBLE_EXP(t1,es+k);
+	    return t1*w;
+	} else return w;
+}
--- a/src/openlibm/ld80/e_lgammal_r.c
+++ b/src/openlibm/ld80/e_lgammal_r.c
@@ -0,0 +1,425 @@
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * 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.
+ */
+
+/* lgammal_r(x, signgamp)
+ * Reentrant version of the logarithm of the Gamma function
+ * with user provide pointer for the sign of Gamma(x).
+ *
+ * Method:
+ *   1. Argument Reduction for 0 < x <= 8
+ *	Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
+ *	reduce x to a number in [1.5,2.5] by
+ *		lgamma(1+s) = log(s) + lgamma(s)
+ *	for example,
+ *		lgamma(7.3) = log(6.3) + lgamma(6.3)
+ *			    = log(6.3*5.3) + lgamma(5.3)
+ *			    = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
+ *   2. Polynomial approximation of lgamma around its
+ *	minimun ymin=1.461632144968362245 to maintain monotonicity.
+ *	On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
+ *		Let z = x-ymin;
+ *		lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
+ *   2. Rational approximation in the primary interval [2,3]
+ *	We use the following approximation:
+ *		s = x-2.0;
+ *		lgamma(x) = 0.5*s + s*P(s)/Q(s)
+ *	Our algorithms are based on the following observation
+ *
+ *                             zeta(2)-1    2    zeta(3)-1    3
+ * lgamma(2+s) = s*(1-Euler) + --------- * s  -  --------- * s  + ...
+ *                                 2                 3
+ *
+ *	where Euler = 0.5771... is the Euler constant, which is very
+ *	close to 0.5.
+ *
+ *   3. For x>=8, we have
+ *	lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
+ *	(better formula:
+ *	   lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
+ *	Let z = 1/x, then we approximation
+ *		f(z) = lgamma(x) - (x-0.5)(log(x)-1)
+ *	by
+ *				    3       5             11
+ *		w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
+ *
+ *   4. For negative x, since (G is gamma function)
+ *		-x*G(-x)*G(x) = pi/sin(pi*x),
+ *	we have
+ *		G(x) = pi/(sin(pi*x)*(-x)*G(-x))
+ *	since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
+ *	Hence, for x<0, signgam = sign(sin(pi*x)) and
+ *		lgamma(x) = log(|Gamma(x)|)
+ *			  = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
+ *	Note: one should avoid compute pi*(-x) directly in the
+ *	      computation of sin(pi*(-x)).
+ *
+ *   5. Special Cases
+ *		lgamma(2+s) ~ s*(1-Euler) for tiny s
+ *		lgamma(1)=lgamma(2)=0
+ *		lgamma(x) ~ -log(x) for tiny x
+ *		lgamma(0) = lgamma(inf) = inf
+ *		lgamma(-integer) = +-inf
+ *
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static const long double
+  half = 0.5L,
+  one = 1.0L,
+  pi = 3.14159265358979323846264L,
+  two63 = 9.223372036854775808e18L,
+
+  /* lgam(1+x) = 0.5 x + x a(x)/b(x)
+     -0.268402099609375 <= x <= 0
+     peak relative error 6.6e-22 */
+  a0 = -6.343246574721079391729402781192128239938E2L,
+  a1 =  1.856560238672465796768677717168371401378E3L,
+  a2 =  2.404733102163746263689288466865843408429E3L,
+  a3 =  8.804188795790383497379532868917517596322E2L,
+  a4 =  1.135361354097447729740103745999661157426E2L,
+  a5 =  3.766956539107615557608581581190400021285E0L,
+
+  b0 =  8.214973713960928795704317259806842490498E3L,
+  b1 =  1.026343508841367384879065363925870888012E4L,
+  b2 =  4.553337477045763320522762343132210919277E3L,
+  b3 =  8.506975785032585797446253359230031874803E2L,
+  b4 =  6.042447899703295436820744186992189445813E1L,
+  /* b5 =  1.000000000000000000000000000000000000000E0 */
+
+
+  tc =  1.4616321449683623412626595423257213284682E0L,
+  tf = -1.2148629053584961146050602565082954242826E-1,/* double precision */
+/* tt = (tail of tf), i.e. tf + tt has extended precision. */
+  tt = 3.3649914684731379602768989080467587736363E-18L,
+  /* lgam ( 1.4616321449683623412626595423257213284682E0 ) =
+-1.2148629053584960809551455717769158215135617312999903886372437313313530E-1 */
+
+  /* lgam (x + tc) = tf + tt + x g(x)/h(x)
+     - 0.230003726999612341262659542325721328468 <= x
+     <= 0.2699962730003876587373404576742786715318
+     peak relative error 2.1e-21 */
+  g0 = 3.645529916721223331888305293534095553827E-18L,
+  g1 = 5.126654642791082497002594216163574795690E3L,
+  g2 = 8.828603575854624811911631336122070070327E3L,
+  g3 = 5.464186426932117031234820886525701595203E3L,
+  g4 = 1.455427403530884193180776558102868592293E3L,
+  g5 = 1.541735456969245924860307497029155838446E2L,
+  g6 = 4.335498275274822298341872707453445815118E0L,
+
+  h0 = 1.059584930106085509696730443974495979641E4L,
+  h1 =  2.147921653490043010629481226937850618860E4L,
+  h2 = 1.643014770044524804175197151958100656728E4L,
+  h3 =  5.869021995186925517228323497501767586078E3L,
+  h4 =  9.764244777714344488787381271643502742293E2L,
+  h5 =  6.442485441570592541741092969581997002349E1L,
+  /* h6 = 1.000000000000000000000000000000000000000E0 */
+
+
+  /* lgam (x+1) = -0.5 x + x u(x)/v(x)
+     -0.100006103515625 <= x <= 0.231639862060546875
+     peak relative error 1.3e-21 */
+  u0 = -8.886217500092090678492242071879342025627E1L,
+  u1 =  6.840109978129177639438792958320783599310E2L,
+  u2 =  2.042626104514127267855588786511809932433E3L,
+  u3 =  1.911723903442667422201651063009856064275E3L,
+  u4 =  7.447065275665887457628865263491667767695E2L,
+  u5 =  1.132256494121790736268471016493103952637E2L,
+  u6 =  4.484398885516614191003094714505960972894E0L,
+
+  v0 =  1.150830924194461522996462401210374632929E3L,
+  v1 =  3.399692260848747447377972081399737098610E3L,
+  v2 =  3.786631705644460255229513563657226008015E3L,
+  v3 =  1.966450123004478374557778781564114347876E3L,
+  v4 =  4.741359068914069299837355438370682773122E2L,
+  v5 =  4.508989649747184050907206782117647852364E1L,
+  /* v6 =  1.000000000000000000000000000000000000000E0 */
+
+
+  /* lgam (x+2) = .5 x + x s(x)/r(x)
+     0 <= x <= 1
+     peak relative error 7.2e-22 */
+  s0 =  1.454726263410661942989109455292824853344E6L,
+  s1 = -3.901428390086348447890408306153378922752E6L,
+  s2 = -6.573568698209374121847873064292963089438E6L,
+  s3 = -3.319055881485044417245964508099095984643E6L,
+  s4 = -7.094891568758439227560184618114707107977E5L,
+  s5 = -6.263426646464505837422314539808112478303E4L,
+  s6 = -1.684926520999477529949915657519454051529E3L,
+
+  r0 = -1.883978160734303518163008696712983134698E7L,
+  r1 = -2.815206082812062064902202753264922306830E7L,
+  r2 = -1.600245495251915899081846093343626358398E7L,
+  r3 = -4.310526301881305003489257052083370058799E6L,
+  r4 = -5.563807682263923279438235987186184968542E5L,
+  r5 = -3.027734654434169996032905158145259713083E4L,
+  r6 = -4.501995652861105629217250715790764371267E2L,
+  /* r6 =  1.000000000000000000000000000000000000000E0 */
+
+
+/* lgam(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x w(1/x^2)
+   x >= 8
+   Peak relative error 1.51e-21
+   w0 = LS2PI - 0.5 */
+  w0 =  4.189385332046727417803e-1L,
+  w1 =  8.333333333333331447505E-2L,
+  w2 = -2.777777777750349603440E-3L,
+  w3 =  7.936507795855070755671E-4L,
+  w4 = -5.952345851765688514613E-4L,
+  w5 =  8.412723297322498080632E-4L,
+  w6 = -1.880801938119376907179E-3L,
+  w7 =  4.885026142432270781165E-3L;
+
+static const long double zero = 0.0L;
+
+static long double
+sin_pi(long double x)
+{
+  long double y, z;
+  int n, ix;
+  u_int32_t se, i0, i1;
+
+  GET_LDOUBLE_WORDS (se, i0, i1, x);
+  ix = se & 0x7fff;
+  ix = (ix << 16) | (i0 >> 16);
+  if (ix < 0x3ffd8000) /* 0.25 */
+    return sinl (pi * x);
+  y = -x;			/* x is assume negative */
+
+  /*
+   * argument reduction, make sure inexact flag not raised if input
+   * is an integer
+   */
+  z = floorl (y);
+  if (z != y)
+    {				/* inexact anyway */
+      y  *= 0.5;
+      y = 2.0*(y - floorl(y));		/* y = |x| mod 2.0 */
+      n = (int) (y*4.0);
+    }
+  else
+    {
+      if (ix >= 0x403f8000)  /* 2^64 */
+	{
+	  y = zero; n = 0;		/* y must be even */
+	}
+      else
+	{
+	if (ix < 0x403e8000)  /* 2^63 */
+	  z = y + two63;	/* exact */
+	GET_LDOUBLE_WORDS (se, i0, i1, z);
+	n = i1 & 1;
+	y  = n;
+	n <<= 2;
+      }
+    }
+
+  switch (n)
+    {
+    case 0:
+      y = sinl (pi * y);
+      break;
+    case 1:
+    case 2:
+      y = cosl (pi * (half - y));
+      break;
+    case 3:
+    case 4:
+      y = sinl (pi * (one - y));
+      break;
+    case 5:
+    case 6:
+      y = -cosl (pi * (y - 1.5));
+      break;
+    default:
+      y = sinl (pi * (y - 2.0));
+      break;
+    }
+  return -y;
+}
+
+
+long double
+lgammal_r(long double x, int *signgamp)
+{
+  long double t, y, z, nadj, p, p1, p2, q, r, w;
+  int i, ix;
+  u_int32_t se, i0, i1;
+
+  *signgamp = 1;
+  GET_LDOUBLE_WORDS (se, i0, i1, x);
+  ix = se & 0x7fff;
+
+  if ((ix | i0 | i1) == 0)
+    {
+      if (se & 0x8000)
+	*signgamp = -1;
+      return one / fabsl (x);
+    }
+
+  ix = (ix << 16) | (i0 >> 16);
+
+  /* purge off +-inf, NaN, +-0, and negative arguments */
+  if (ix >= 0x7fff0000)
+    return x * x;
+
+  if (ix < 0x3fc08000) /* 2^-63 */
+    {				/* |x|<2**-63, return -log(|x|) */
+      if (se & 0x8000)
+	{
+	  *signgamp = -1;
+	  return -logl (-x);
+	}
+      else
+	return -logl (x);
+    }
+  if (se & 0x8000)
+    {
+      t = sin_pi (x);
+      if (t == zero)
+	return one / fabsl (t);	/* -integer */
+      nadj = logl (pi / fabsl (t * x));
+      if (t < zero)
+	*signgamp = -1;
+      x = -x;
+    }
+
+  /* purge off 1 and 2 */
+  if ((((ix - 0x3fff8000) | i0 | i1) == 0)
+      || (((ix - 0x40008000) | i0 | i1) == 0))
+    r = 0;
+  else if (ix < 0x40008000) /* 2.0 */
+    {
+      /* x < 2.0 */
+      if (ix <= 0x3ffee666) /* 8.99993896484375e-1 */
+	{
+	  /* lgamma(x) = lgamma(x+1) - log(x) */
+	  r = -logl (x);
+	  if (ix >= 0x3ffebb4a) /* 7.31597900390625e-1 */
+	    {
+	      y = x - one;
+	      i = 0;
+	    }
+	  else if (ix >= 0x3ffced33)/* 2.31639862060546875e-1 */
+	    {
+	      y = x - (tc - one);
+	      i = 1;
+	    }
+	  else
+	    {
+	      /* x < 0.23 */
+	      y = x;
+	      i = 2;
+	    }
+	}
+      else
+	{
+	  r = zero;
+	  if (ix >= 0x3fffdda6) /* 1.73162841796875 */
+	    {
+	      /* [1.7316,2] */
+	      y = x - 2.0;
+	      i = 0;
+	    }
+	  else if (ix >= 0x3fff9da6)/* 1.23162841796875 */
+	    {
+	      /* [1.23,1.73] */
+	      y = x - tc;
+	      i = 1;
+	    }
+	  else
+	    {
+	      /* [0.9, 1.23] */
+	      y = x - one;
+	      i = 2;
+	    }
+	}
+      switch (i)
+	{
+	case 0:
+	  p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))));
+	  p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y))));
+	  r += half * y + y * p1/p2;
+	  break;
+	case 1:
+    p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6)))));
+    p2 = h0 + y * (h1 + y * (h2 + y * (h3 + y * (h4 + y * (h5 + y)))));
+    p = tt + y * p1/p2;
+	  r += (tf + p);
+	  break;
+	case 2:
+ p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6))))));
+      p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y)))));
+	  r += (-half * y + p1 / p2);
+	}
+    }
+  else if (ix < 0x40028000) /* 8.0 */
+    {
+      /* x < 8.0 */
+      i = (int) x;
+      t = zero;
+      y = x - (double) i;
+  p = y *
+     (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
+  q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
+      r = half * y + p / q;
+      z = one;			/* lgamma(1+s) = log(s) + lgamma(s) */
+      switch (i)
+	{
+	case 7:
+	  z *= (y + 6.0);	/* FALLTHRU */
+	case 6:
+	  z *= (y + 5.0);	/* FALLTHRU */
+	case 5:
+	  z *= (y + 4.0);	/* FALLTHRU */
+	case 4:
+	  z *= (y + 3.0);	/* FALLTHRU */
+	case 3:
+	  z *= (y + 2.0);	/* FALLTHRU */
+	  r += logl (z);
+	  break;
+	}
+    }
+  else if (ix < 0x40418000) /* 2^66 */
+    {
+      /* 8.0 <= x < 2**66 */
+      t = logl (x);
+      z = one / x;
+      y = z * z;
+      w = w0 + z * (w1
+	  + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7))))));
+      r = (x - half) * (t - one) + w;
+    }
+  else
+    /* 2**66 <= x <= inf */
+    r = x * (logl (x) - one);
+  if (se & 0x8000)
+    r = nadj - r;
+  return r;
+}
--- a/src/openlibm/ld80/e_lgammal_r.c.o
+++ b/src/openlibm/ld80/e_lgammal_r.c.o
--- a/src/openlibm/ld80/e_log10l.c
+++ b/src/openlibm/ld80/e_log10l.c
@@ -0,0 +1,205 @@
+/*	$OpenBSD: e_log10l.c,v 1.2 2013/11/12 20:35:19 martynas Exp $	*/
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * 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.
+ */
+
+/*							log10l.c
+ *
+ *	Common logarithm, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, log10l();
+ *
+ * y = log10l( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base 10 logarithm of x.
+ *
+ * The argument is separated into its exponent and fractional
+ * parts.  If the exponent is between -1 and +1, the logarithm
+ * of the fraction is approximated by
+ *
+ *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
+ *
+ * Otherwise, setting  z = 2(x-1)/x+1),
+ *
+ *     log(x) = z + z**3 P(z)/Q(z).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0.5, 2.0     30000      9.0e-20     2.6e-20
+ *    IEEE     exp(+-10000)  30000      6.0e-20     2.3e-20
+ *
+ * In the tests over the interval exp(+-10000), the logarithms
+ * of the random arguments were uniformly distributed over
+ * [-10000, +10000].
+ *
+ * ERROR MESSAGES:
+ *
+ * log singularity:  x = 0; returns MINLOG
+ * log domain:       x < 0; returns MINLOG
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 6.2e-22
+ */
+static long double P[] = {
+ 4.9962495940332550844739E-1L,
+ 1.0767376367209449010438E1L,
+ 7.7671073698359539859595E1L,
+ 2.5620629828144409632571E2L,
+ 4.2401812743503691187826E2L,
+ 3.4258224542413922935104E2L,
+ 1.0747524399916215149070E2L,
+};
+static long double Q[] = {
+/* 1.0000000000000000000000E0,*/
+ 2.3479774160285863271658E1L,
+ 1.9444210022760132894510E2L,
+ 7.7952888181207260646090E2L,
+ 1.6911722418503949084863E3L,
+ 2.0307734695595183428202E3L,
+ 1.2695660352705325274404E3L,
+ 3.2242573199748645407652E2L,
+};
+
+/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
+ * where z = 2(x-1)/(x+1)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 6.16e-22
+ */
+
+static long double R[4] = {
+ 1.9757429581415468984296E-3L,
+-7.1990767473014147232598E-1L,
+ 1.0777257190312272158094E1L,
+-3.5717684488096787370998E1L,
+};
+static long double S[4] = {
+/* 1.00000000000000000000E0L,*/
+-2.6201045551331104417768E1L,
+ 1.9361891836232102174846E2L,
+-4.2861221385716144629696E2L,
+};
+/* log10(2) */
+#define L102A 0.3125L
+#define L102B -1.1470004336018804786261e-2L
+/* log10(e) */
+#define L10EA 0.5L
+#define L10EB -6.5705518096748172348871e-2L
+
+#define SQRTH 0.70710678118654752440L
+
+long double
+log10l(long double x)
+{
+long double y;
+volatile long double z;
+int e;
+
+if( isnan(x) )
+	return(x);
+/* Test for domain */
+if( x <= 0.0L )
+	{
+	if( x == 0.0L )
+		return (-1.0L / (x - x));
+	else
+		return (x - x) / (x - x);
+	}
+if( x == INFINITY )
+	return(INFINITY);
+/* separate mantissa from exponent */
+
+/* Note, frexp is used so that denormal numbers
+ * will be handled properly.
+ */
+x = frexpl( x, &e );
+
+
+/* logarithm using log(x) = z + z**3 P(z)/Q(z),
+ * where z = 2(x-1)/x+1)
+ */
+if( (e > 2) || (e < -2) )
+{
+if( x < SQRTH )
+	{ /* 2( 2x-1 )/( 2x+1 ) */
+	e -= 1;
+	z = x - 0.5L;
+	y = 0.5L * z + 0.5L;
+	}	
+else
+	{ /*  2 (x-1)/(x+1)   */
+	z = x - 0.5L;
+	z -= 0.5L;
+	y = 0.5L * x  + 0.5L;
+	}
+x = z / y;
+z = x*x;
+y = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) );
+goto done;
+}
+
+
+/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
+
+if( x < SQRTH )
+	{
+	e -= 1;
+	x = ldexpl( x, 1 ) - 1.0L; /*  2x - 1  */
+	}	
+else
+	{
+	x = x - 1.0L;
+	}
+z = x*x;
+y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 7 ) );
+y = y - ldexpl( z, -1 );   /* -0.5x^2 + ... */
+
+done:
+
+/* Multiply log of fraction by log10(e)
+ * and base 2 exponent by log10(2).
+ *
+ * ***CAUTION***
+ *
+ * This sequence of operations is critical and it may
+ * be horribly defeated by some compiler optimizers.
+ */
+z = y * (L10EB);
+z += x * (L10EB);
+z += e * (L102B);
+z += y * (L10EA);
+z += x * (L10EA);
+z += e * (L102A);
+
+return( z );
+}
--- a/src/openlibm/ld80/e_log10l.c.o
+++ b/src/openlibm/ld80/e_log10l.c.o
--- a/src/openlibm/ld80/e_log2l.c
+++ b/src/openlibm/ld80/e_log2l.c
@@ -0,0 +1,199 @@
+/*	$OpenBSD: e_log2l.c,v 1.2 2013/11/12 20:35:19 martynas Exp $	*/
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * 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.
+ */
+
+/*							log2l.c
+ *
+ *	Base 2 logarithm, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, log2l();
+ *
+ * y = log2l( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base 2 logarithm of x.
+ *
+ * The argument is separated into its exponent and fractional
+ * parts.  If the exponent is between -1 and +1, the (natural)
+ * logarithm of the fraction is approximated by
+ *
+ *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
+ *
+ * Otherwise, setting  z = 2(x-1)/x+1),
+ *
+ *     log(x) = z + z**3 P(z)/Q(z).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0.5, 2.0     30000      9.8e-20     2.7e-20
+ *    IEEE     exp(+-10000)  70000      5.4e-20     2.3e-20
+ *
+ * In the tests over the interval exp(+-10000), the logarithms
+ * of the random arguments were uniformly distributed over
+ * [-10000, +10000].
+ *
+ * ERROR MESSAGES:
+ *
+ * log singularity:  x = 0; returns -INFINITY
+ * log domain:       x < 0; returns NAN
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+/* Coefficients for ln(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 6.2e-22
+ */
+static long double P[] = {
+ 4.9962495940332550844739E-1L,
+ 1.0767376367209449010438E1L,
+ 7.7671073698359539859595E1L,
+ 2.5620629828144409632571E2L,
+ 4.2401812743503691187826E2L,
+ 3.4258224542413922935104E2L,
+ 1.0747524399916215149070E2L,
+};
+static long double Q[] = {
+/* 1.0000000000000000000000E0,*/
+ 2.3479774160285863271658E1L,
+ 1.9444210022760132894510E2L,
+ 7.7952888181207260646090E2L,
+ 1.6911722418503949084863E3L,
+ 2.0307734695595183428202E3L,
+ 1.2695660352705325274404E3L,
+ 3.2242573199748645407652E2L,
+};
+
+/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
+ * where z = 2(x-1)/(x+1)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 6.16e-22
+ */
+static long double R[4] = {
+ 1.9757429581415468984296E-3L,
+-7.1990767473014147232598E-1L,
+ 1.0777257190312272158094E1L,
+-3.5717684488096787370998E1L,
+};
+static long double S[4] = {
+/* 1.00000000000000000000E0L,*/
+-2.6201045551331104417768E1L,
+ 1.9361891836232102174846E2L,
+-4.2861221385716144629696E2L,
+};
+/* log2(e) - 1 */
+#define LOG2EA 4.4269504088896340735992e-1L
+
+#define SQRTH 0.70710678118654752440L
+
+long double
+log2l(long double x)
+{
+volatile long double z;
+long double y;
+int e;
+
+if( isnan(x) )
+	return(x);
+if( x == INFINITY )
+	return(x);
+/* Test for domain */
+if( x <= 0.0L )
+	{
+	if( x == 0.0L )
+		return( -INFINITY );
+	else
+		return( NAN );
+	}
+
+/* separate mantissa from exponent */
+
+/* Note, frexp is used so that denormal numbers
+ * will be handled properly.
+ */
+x = frexpl( x, &e );
+
+
+/* logarithm using log(x) = z + z**3 P(z)/Q(z),
+ * where z = 2(x-1)/x+1)
+ */
+if( (e > 2) || (e < -2) )
+{
+if( x < SQRTH )
+	{ /* 2( 2x-1 )/( 2x+1 ) */
+	e -= 1;
+	z = x - 0.5L;
+	y = 0.5L * z + 0.5L;
+	}	
+else
+	{ /*  2 (x-1)/(x+1)   */
+	z = x - 0.5L;
+	z -= 0.5L;
+	y = 0.5L * x  + 0.5L;
+	}
+x = z / y;
+z = x*x;
+y = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) );
+goto done;
+}
+
+
+/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
+
+if( x < SQRTH )
+	{
+	e -= 1;
+	x = ldexpl( x, 1 ) - 1.0L; /*  2x - 1  */
+	}	
+else
+	{
+	x = x - 1.0L;
+	}
+z = x*x;
+y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 7 ) );
+y = y - ldexpl( z, -1 );   /* -0.5x^2 + ... */
+
+done:
+
+/* Multiply log of fraction by log2(e)
+ * and base 2 exponent by 1
+ *
+ * ***CAUTION***
+ *
+ * This sequence of operations is critical and it may
+ * be horribly defeated by some compiler optimizers.
+ */
+z = y * LOG2EA;
+z += x * LOG2EA;
+z += y;
+z += x;
+z += e;
+return( z );
+}
--- a/src/openlibm/ld80/e_log2l.c.o
+++ b/src/openlibm/ld80/e_log2l.c.o
--- a/src/openlibm/ld80/e_logl.c
+++ b/src/openlibm/ld80/e_logl.c
@@ -0,0 +1,190 @@
+/*	$OpenBSD: e_logl.c,v 1.3 2013/11/12 20:35:19 martynas Exp $	*/
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * 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.
+ */
+
+/*							logl.c
+ *
+ *	Natural logarithm, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, logl();
+ *
+ * y = logl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base e (2.718...) logarithm of x.
+ *
+ * The argument is separated into its exponent and fractional
+ * parts.  If the exponent is between -1 and +1, the logarithm
+ * of the fraction is approximated by
+ *
+ *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
+ *
+ * Otherwise, setting  z = 2(x-1)/x+1),
+ *
+ *     log(x) = z + z**3 P(z)/Q(z).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0.5, 2.0    150000      8.71e-20    2.75e-20
+ *    IEEE     exp(+-10000) 100000      5.39e-20    2.34e-20
+ *
+ * In the tests over the interval exp(+-10000), the logarithms
+ * of the random arguments were uniformly distributed over
+ * [-10000, +10000].
+ *
+ * ERROR MESSAGES:
+ *
+ * log singularity:  x = 0; returns -INFINITY
+ * log domain:       x < 0; returns NAN
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 2.32e-20
+ */
+static long double P[] = {
+ 4.5270000862445199635215E-5L,
+ 4.9854102823193375972212E-1L,
+ 6.5787325942061044846969E0L,
+ 2.9911919328553073277375E1L,
+ 6.0949667980987787057556E1L,
+ 5.7112963590585538103336E1L,
+ 2.0039553499201281259648E1L,
+};
+static long double Q[] = {
+/* 1.0000000000000000000000E0,*/
+ 1.5062909083469192043167E1L,
+ 8.3047565967967209469434E1L,
+ 2.2176239823732856465394E2L,
+ 3.0909872225312059774938E2L,
+ 2.1642788614495947685003E2L,
+ 6.0118660497603843919306E1L,
+};
+
+/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
+ * where z = 2(x-1)/(x+1)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 6.16e-22
+ */
+
+static long double R[4] = {
+ 1.9757429581415468984296E-3L,
+-7.1990767473014147232598E-1L,
+ 1.0777257190312272158094E1L,
+-3.5717684488096787370998E1L,
+};
+static long double S[4] = {
+/* 1.00000000000000000000E0L,*/
+-2.6201045551331104417768E1L,
+ 1.9361891836232102174846E2L,
+-4.2861221385716144629696E2L,
+};
+static const long double C1 = 6.9314575195312500000000E-1L;
+static const long double C2 = 1.4286068203094172321215E-6L;
+
+#define SQRTH 0.70710678118654752440L
+
+long double
+logl(long double x)
+{
+long double y, z;
+int e;
+
+if( isnan(x) )
+	return(x);
+if( x == INFINITY )
+	return(x);
+/* Test for domain */
+if( x <= 0.0L )
+	{
+	if( x == 0.0L )
+		return( -INFINITY );
+	else
+		return( NAN );
+	}
+
+/* separate mantissa from exponent */
+
+/* Note, frexp is used so that denormal numbers
+ * will be handled properly.
+ */
+x = frexpl( x, &e );
+
+/* logarithm using log(x) = z + z**3 P(z)/Q(z),
+ * where z = 2(x-1)/x+1)
+ */
+if( (e > 2) || (e < -2) )
+{
+if( x < SQRTH )
+	{ /* 2( 2x-1 )/( 2x+1 ) */
+	e -= 1;
+	z = x - 0.5L;
+	y = 0.5L * z + 0.5L;
+	}	
+else
+	{ /*  2 (x-1)/(x+1)   */
+	z = x - 0.5L;
+	z -= 0.5L;
+	y = 0.5L * x  + 0.5L;
+	}
+x = z / y;
+z = x*x;
+z = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) );
+z = z + e * C2;
+z = z + x;
+z = z + e * C1;
+return( z );
+}
+
+
+/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
+
+if( x < SQRTH )
+	{
+	e -= 1;
+	x = ldexpl( x, 1 ) - 1.0L; /*  2x - 1  */
+	}	
+else
+	{
+	x = x - 1.0L;
+	}
+z = x*x;
+y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 6 ) );
+y = y + e * C2;
+z = y - ldexpl( z, -1 );   /*  y - 0.5 * z  */
+/* Note, the sum of above terms does not exceed x/4,
+ * so it contributes at most about 1/4 lsb to the error.
+ */
+z = z + x;
+z = z + e * C1; /* This sum has an error of 1/2 lsb. */
+return( z );
+}
--- a/src/openlibm/ld80/e_logl.c.o
+++ b/src/openlibm/ld80/e_logl.c.o
--- a/src/openlibm/ld80/e_powl.c
+++ b/src/openlibm/ld80/e_powl.c
@@ -0,0 +1,615 @@
+/*	$OpenBSD: e_powl.c,v 1.5 2013/11/12 20:35:19 martynas Exp $	*/
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * 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.
+ */
+
+/*							powl.c
+ *
+ *	Power function, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, z, powl();
+ *
+ * z = powl( x, y );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Computes x raised to the yth power.  Analytically,
+ *
+ *      x**y  =  exp( y log(x) ).
+ *
+ * Following Cody and Waite, this program uses a lookup table
+ * of 2**-i/32 and pseudo extended precision arithmetic to
+ * obtain several extra bits of accuracy in both the logarithm
+ * and the exponential.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * The relative error of pow(x,y) can be estimated
+ * by   y dl ln(2),   where dl is the absolute error of
+ * the internally computed base 2 logarithm.  At the ends
+ * of the approximation interval the logarithm equal 1/32
+ * and its relative error is about 1 lsb = 1.1e-19.  Hence
+ * the predicted relative error in the result is 2.3e-21 y .
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *
+ *    IEEE     +-1000       40000      2.8e-18      3.7e-19
+ * .001 < x < 1000, with log(x) uniformly distributed.
+ * -1000 < y < 1000, y uniformly distributed.
+ *
+ *    IEEE     0,8700       60000      6.5e-18      1.0e-18
+ * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * pow overflow     x**y > MAXNUM      INFINITY
+ * pow underflow   x**y < 1/MAXNUM       0.0
+ * pow domain      x<0 and y noninteger  0.0
+ *
+ */
+
+#include <float.h>
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+/* Table size */
+#define NXT 32
+/* log2(Table size) */
+#define LNXT 5
+
+/* log(1+x) =  x - .5x^2 + x^3 *  P(z)/Q(z)
+ * on the domain  2^(-1/32) - 1  <=  x  <=  2^(1/32) - 1
+ */
+static long double P[] = {
+ 8.3319510773868690346226E-4L,
+ 4.9000050881978028599627E-1L,
+ 1.7500123722550302671919E0L,
+ 1.4000100839971580279335E0L,
+};
+static long double Q[] = {
+/* 1.0000000000000000000000E0L,*/
+ 5.2500282295834889175431E0L,
+ 8.4000598057587009834666E0L,
+ 4.2000302519914740834728E0L,
+};
+/* A[i] = 2^(-i/32), rounded to IEEE long double precision.
+ * If i is even, A[i] + B[i/2] gives additional accuracy.
+ */
+static long double A[33] = {
+ 1.0000000000000000000000E0L,
+ 9.7857206208770013448287E-1L,
+ 9.5760328069857364691013E-1L,
+ 9.3708381705514995065011E-1L,
+ 9.1700404320467123175367E-1L,
+ 8.9735453750155359320742E-1L,
+ 8.7812608018664974155474E-1L,
+ 8.5930964906123895780165E-1L,
+ 8.4089641525371454301892E-1L,
+ 8.2287773907698242225554E-1L,
+ 8.0524516597462715409607E-1L,
+ 7.8799042255394324325455E-1L,
+ 7.7110541270397041179298E-1L,
+ 7.5458221379671136985669E-1L,
+ 7.3841307296974965571198E-1L,
+ 7.2259040348852331001267E-1L,
+ 7.0710678118654752438189E-1L,
+ 6.9195494098191597746178E-1L,
+ 6.7712777346844636413344E-1L,
+ 6.6261832157987064729696E-1L,
+ 6.4841977732550483296079E-1L,
+ 6.3452547859586661129850E-1L,
+ 6.2092890603674202431705E-1L,
+ 6.0762367999023443907803E-1L,
+ 5.9460355750136053334378E-1L,
+ 5.8186242938878875689693E-1L,
+ 5.6939431737834582684856E-1L,
+ 5.5719337129794626814472E-1L,
+ 5.4525386633262882960438E-1L,
+ 5.3357020033841180906486E-1L,
+ 5.2213689121370692017331E-1L,
+ 5.1094857432705833910408E-1L,
+ 5.0000000000000000000000E-1L,
+};
+static long double B[17] = {
+ 0.0000000000000000000000E0L,
+ 2.6176170809902549338711E-20L,
+-1.0126791927256478897086E-20L,
+ 1.3438228172316276937655E-21L,
+ 1.2207982955417546912101E-20L,
+-6.3084814358060867200133E-21L,
+ 1.3164426894366316434230E-20L,
+-1.8527916071632873716786E-20L,
+ 1.8950325588932570796551E-20L,
+ 1.5564775779538780478155E-20L,
+ 6.0859793637556860974380E-21L,
+-2.0208749253662532228949E-20L,
+ 1.4966292219224761844552E-20L,
+ 3.3540909728056476875639E-21L,
+-8.6987564101742849540743E-22L,
+-1.2327176863327626135542E-20L,
+ 0.0000000000000000000000E0L,
+};
+
+/* 2^x = 1 + x P(x),
+ * on the interval -1/32 <= x <= 0
+ */
+static long double R[] = {
+ 1.5089970579127659901157E-5L,
+ 1.5402715328927013076125E-4L,
+ 1.3333556028915671091390E-3L,
+ 9.6181291046036762031786E-3L,
+ 5.5504108664798463044015E-2L,
+ 2.4022650695910062854352E-1L,
+ 6.9314718055994530931447E-1L,
+};
+
+#define douba(k) A[k]
+#define doubb(k) B[k]
+#define MEXP (NXT*16384.0L)
+/* The following if denormal numbers are supported, else -MEXP: */
+#define MNEXP (-NXT*(16384.0L+64.0L))
+/* log2(e) - 1 */
+#define LOG2EA 0.44269504088896340735992L
+
+#define F W
+#define Fa Wa
+#define Fb Wb
+#define G W
+#define Ga Wa
+#define Gb u
+#define H W
+#define Ha Wb
+#define Hb Wb
+
+static const long double MAXLOGL = 1.1356523406294143949492E4L;
+static const long double MINLOGL = -1.13994985314888605586758E4L;
+static const long double LOGE2L = 6.9314718055994530941723E-1L;
+static volatile long double z;
+static long double w, W, Wa, Wb, ya, yb, u;
+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 long double reducl( long double );
+static long double powil ( long double, int );
+
+long double
+powl(long double x, long double y)
+{
+/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
+int i, nflg, iyflg, yoddint;
+long e;
+
+if( y == 0.0L )
+	return( 1.0L );
+
+if( x == 1.0L )
+	return( 1.0L );
+
+if( isnan(x) )
+	return( x );
+if( isnan(y) )
+	return( y );
+
+if( y == 1.0L )
+	return( x );
+
+if( !isfinite(y) && x == -1.0L )
+	return( 1.0L );
+
+if( y >= LDBL_MAX )
+	{
+	if( x > 1.0L )
+		return( INFINITY );
+	if( x > 0.0L && x < 1.0L )
+		return( 0.0L );
+	if( x < -1.0L )
+		return( INFINITY );
+	if( x > -1.0L && x < 0.0L )
+		return( 0.0L );
+	}
+if( y <= -LDBL_MAX )
+	{
+	if( x > 1.0L )
+		return( 0.0L );
+	if( x > 0.0L && x < 1.0L )
+		return( INFINITY );
+	if( x < -1.0L )
+		return( 0.0L );
+	if( x > -1.0L && x < 0.0L )
+		return( INFINITY );
+	}
+if( x >= LDBL_MAX )
+	{
+	if( y > 0.0L )
+		return( INFINITY );
+	return( 0.0L );
+	}
+
+w = floorl(y);
+/* Set iyflg to 1 if y is an integer.  */
+iyflg = 0;
+if( w == y )
+	iyflg = 1;
+
+/* Test for odd integer y.  */
+yoddint = 0;
+if( iyflg )
+	{
+	ya = fabsl(y);
+	ya = floorl(0.5L * ya);
+	yb = 0.5L * fabsl(w);
+	if( ya != yb )
+		yoddint = 1;
+	}
+
+if( x <= -LDBL_MAX )
+	{
+	if( y > 0.0L )
+		{
+		if( yoddint )
+			return( -INFINITY );
+		return( INFINITY );
+		}
+	if( y < 0.0L )
+		{
+		if( yoddint )
+			return( -0.0L );
+		return( 0.0 );
+		}
+	}
+
+
+nflg = 0;	/* flag = 1 if x<0 raised to integer power */
+if( x <= 0.0L )
+	{
+	if( x == 0.0L )
+		{
+		if( y < 0.0 )
+			{
+			if( signbit(x) && yoddint )
+				return( -INFINITY );
+			return( INFINITY );
+			}
+		if( y > 0.0 )
+			{
+			if( signbit(x) && yoddint )
+				return( -0.0L );
+			return( 0.0 );
+			}
+		if( y == 0.0L )
+			return( 1.0L );  /*   0**0   */
+		else
+			return( 0.0L );  /*   0**y   */
+		}
+	else
+		{
+		if( iyflg == 0 )
+			return (x - x) / (x - x); /* (x<0)**(non-int) is NaN */
+		nflg = 1;
+		}
+	}
+
+/* Integer power of an integer.  */
+
+if( iyflg )
+	{
+	i = w;
+	w = floorl(x);
+	if( (w == x) && (fabsl(y) < 32768.0) )
+		{
+		w = powil( x, (int) y );
+		return( w );
+		}
+	}
+
+
+if( nflg )
+	x = fabsl(x);
+
+/* separate significand from exponent */
+x = frexpl( x, &i );
+e = i;
+
+/* find significand in antilog table A[] */
+i = 1;
+if( x <= douba(17) )
+	i = 17;
+if( x <= douba(i+8) )
+	i += 8;
+if( x <= douba(i+4) )
+	i += 4;
+if( x <= douba(i+2) )
+	i += 2;
+if( x >= douba(1) )
+	i = -1;
+i += 1;
+
+
+/* Find (x - A[i])/A[i]
+ * in order to compute log(x/A[i]):
+ *
+ * log(x) = log( a x/a ) = log(a) + log(x/a)
+ *
+ * log(x/a) = log(1+v),  v = x/a - 1 = (x-a)/a
+ */
+x -= douba(i);
+x -= doubb(i/2);
+x /= douba(i);
+
+
+/* rational approximation for log(1+v):
+ *
+ * log(1+v)  =  v  -  v**2/2  +  v**3 P(v) / Q(v)
+ */
+z = x*x;
+w = x * ( z * __polevll( x, P, 3 ) / __p1evll( x, Q, 3 ) );
+w = w - ldexpl( z, -1 );   /*  w - 0.5 * z  */
+
+/* Convert to base 2 logarithm:
+ * multiply by log2(e) = 1 + LOG2EA
+ */
+z = LOG2EA * w;
+z += w;
+z += LOG2EA * x;
+z += x;
+
+/* Compute exponent term of the base 2 logarithm. */
+w = -i;
+w = ldexpl( w, -LNXT );	/* divide by NXT */
+w += e;
+/* Now base 2 log of x is w + z. */
+
+/* Multiply base 2 log by y, in extended precision. */
+
+/* separate y into large part ya
+ * and small part yb less than 1/NXT
+ */
+ya = reducl(y);
+yb = y - ya;
+
+/* (w+z)(ya+yb)
+ * = w*ya + w*yb + z*y
+ */
+F = z * y  +  w * yb;
+Fa = reducl(F);
+Fb = F - Fa;
+
+G = Fa + w * ya;
+Ga = reducl(G);
+Gb = G - Ga;
+
+H = Fb + Gb;
+Ha = reducl(H);
+w = ldexpl( Ga+Ha, LNXT );
+
+/* Test the power of 2 for overflow */
+if( w > MEXP )
+	return (huge * huge);		/* overflow */
+
+if( w < MNEXP )
+	return (twom10000 * twom10000);	/* underflow */
+
+e = w;
+Hb = H - Ha;
+
+if( Hb > 0.0L )
+	{
+	e += 1;
+	Hb -= (1.0L/NXT);  /*0.0625L;*/
+	}
+
+/* Now the product y * log2(x)  =  Hb + e/NXT.
+ *
+ * Compute base 2 exponential of Hb,
+ * where -0.0625 <= Hb <= 0.
+ */
+z = Hb * __polevll( Hb, R, 6 );  /*    z  =  2**Hb - 1    */
+
+/* Express e/NXT as an integer plus a negative number of (1/NXT)ths.
+ * Find lookup table entry for the fractional power of 2.
+ */
+if( e < 0 )
+	i = 0;
+else
+	i = 1;
+i = e/NXT + i;
+e = NXT*i - e;
+w = douba( e );
+z = w * z;      /*    2**-e * ( 1 + (2**Hb-1) )    */
+z = z + w;
+z = ldexpl( z, i );  /* multiply by integer power of 2 */
+
+if( nflg )
+	{
+/* For negative x,
+ * find out if the integer exponent
+ * is odd or even.
+ */
+	w = ldexpl( y, -1 );
+	w = floorl(w);
+	w = ldexpl( w, 1 );
+	if( w != y )
+		z = -z; /* odd exponent */
+	}
+
+return( z );
+}
+
+
+/* Find a multiple of 1/NXT that is within 1/NXT of x. */
+static long double
+reducl(long double x)
+{
+long double t;
+
+t = ldexpl( x, LNXT );
+t = floorl( t );
+t = ldexpl( t, -LNXT );
+return(t);
+}
+
+/*							powil.c
+ *
+ *	Real raised to integer power, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, powil();
+ * int n;
+ *
+ * y = powil( x, n );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns argument x raised to the nth power.
+ * The routine efficiently decomposes n as a sum of powers of
+ * two. The desired power is a product of two-to-the-kth
+ * powers of x.  Thus to compute the 32767 power of x requires
+ * 28 multiplications instead of 32767 multiplications.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *
+ *                      Relative error:
+ * arithmetic   x domain   n domain  # trials      peak         rms
+ *    IEEE     .001,1000  -1022,1023  50000       4.3e-17     7.8e-18
+ *    IEEE        1,2     -1022,1023  20000       3.9e-17     7.6e-18
+ *    IEEE     .99,1.01     0,8700    10000       3.6e-16     7.2e-17
+ *
+ * Returns MAXNUM on overflow, zero on underflow.
+ *
+ */
+
+static long double
+powil(long double x, int nn)
+{
+long double ww, y;
+long double s;
+int n, e, sign, asign, lx;
+
+if( x == 0.0L )
+	{
+	if( nn == 0 )
+		return( 1.0L );
+	else if( nn < 0 )
+		return( LDBL_MAX );
+	else
+		return( 0.0L );
+	}
+
+if( nn == 0 )
+	return( 1.0L );
+
+
+if( x < 0.0L )
+	{
+	asign = -1;
+	x = -x;
+	}
+else
+	asign = 0;
+
+
+if( nn < 0 )
+	{
+	sign = -1;
+	n = -nn;
+	}
+else
+	{
+	sign = 1;
+	n = nn;
+	}
+
+/* Overflow detection */
+
+/* Calculate approximate logarithm of answer */
+s = x;
+s = frexpl( s, &lx );
+e = (lx - 1)*n;
+if( (e == 0) || (e > 64) || (e < -64) )
+	{
+	s = (s - 7.0710678118654752e-1L) / (s +  7.0710678118654752e-1L);
+	s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;
+	}
+else
+	{
+	s = LOGE2L * e;
+	}
+
+if( s > MAXLOGL )
+	return (huge * huge);		/* overflow */
+
+if( s < MINLOGL )
+	return (twom10000 * twom10000);	/* underflow */
+/* Handle tiny denormal answer, but with less accuracy
+ * since roundoff error in 1.0/x will be amplified.
+ * The precise demarcation should be the gradual underflow threshold.
+ */
+if( s < (-MAXLOGL+2.0L) )
+	{
+	x = 1.0L/x;
+	sign = -sign;
+	}
+
+/* First bit of the power */
+if( n & 1 )
+	y = x;
+		
+else
+	{
+	y = 1.0L;
+	asign = 0;
+	}
+
+ww = x;
+n >>= 1;
+while( n )
+	{
+	ww = ww * ww;	/* arg to the 2-to-the-kth power */
+	if( n & 1 )	/* if that bit is set, then include in product */
+		y *= ww;
+	n >>= 1;
+	}
+
+if( asign )
+	y = -y; /* odd power of negative number */
+if( sign < 0 )
+	y = 1.0L/y;
+return(y);
+}
--- a/src/openlibm/ld80/e_powl.c.o
+++ b/src/openlibm/ld80/e_powl.c.o
--- a/src/openlibm/ld80/e_rem_pio2l.h
+++ b/src/openlibm/ld80/e_rem_pio2l.h
@@ -0,0 +1,152 @@
+/* From: @(#)e_rem_pio2.c 1.4 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.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+#include "cdefs-compat.h"
+//__FBSDID("$FreeBSD: src/lib/msun/ld80/e_rem_pio2l.h,v 1.3 2011/06/18 13:56:33 benl Exp $");
+
+/* ld80 version of __ieee754_rem_pio2l(x,y)
+ *
+ * return the remainder of x rem pi/2 in y[0]+y[1]
+ * use __kernel_rem_pio2()
+ */
+
+#include <float.h>
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+#define	BIAS	(LDBL_MAX_EXP - 1)
+
+/*
+ * invpio2:  64 bits of 2/pi
+ * pio2_1:   first  39 bits of pi/2
+ * pio2_1t:  pi/2 - pio2_1
+ * pio2_2:   second 39 bits of pi/2
+ * pio2_2t:  pi/2 - (pio2_1+pio2_2)
+ * pio2_3:   third  39 bits of pi/2
+ * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
+ */
+
+static const double
+zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+two24 =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
+pio2_1  =  1.57079632679597125389e+00,	/* 0x3FF921FB, 0x54444000 */
+pio2_2  = -1.07463465549783099519e-12,	/* -0x12e7b967674000.0p-92 */
+pio2_3  =  6.36831716351370313614e-25;	/*  0x18a2e037074000.0p-133 */
+
+#if defined(__amd64__) || defined(__i386__)
+/* Long double constants are slow on these arches, and broken on i386. */
+static const volatile double
+invpio2hi =  6.3661977236758138e-01,	/*  0x145f306dc9c883.0p-53 */
+invpio2lo = -3.9356538861223811e-17,	/* -0x16b00000000000.0p-107 */
+pio2_1thi = -1.0746346554971943e-12,	/* -0x12e7b9676733af.0p-92 */
+pio2_1tlo =  8.8451028997905949e-29,	/*  0x1c080000000000.0p-146 */
+pio2_2thi =  6.3683171635109499e-25,	/*  0x18a2e03707344a.0p-133 */
+pio2_2tlo =  2.3183081793789774e-41,	/*  0x10280000000000.0p-187 */
+pio2_3thi = -2.7529965190440717e-37,	/* -0x176b7ed8fbbacc.0p-174 */
+pio2_3tlo = -4.2006647512740502e-54;	/* -0x19c00000000000.0p-230 */
+#define	invpio2	((long double)invpio2hi + invpio2lo)
+#define	pio2_1t	((long double)pio2_1thi + pio2_1tlo)
+#define	pio2_2t	((long double)pio2_2thi + pio2_2tlo)
+#define	pio2_3t	((long double)pio2_3thi + pio2_3tlo)
+#else
+static const long double
+invpio2 =  6.36619772367581343076e-01L,	/*  0xa2f9836e4e44152a.0p-64 */
+pio2_1t = -1.07463465549719416346e-12L,	/* -0x973dcb3b399d747f.0p-103 */
+pio2_2t =  6.36831716351095013979e-25L,	/*  0xc51701b839a25205.0p-144 */
+pio2_3t = -2.75299651904407171810e-37L;	/* -0xbb5bf6c7ddd660ce.0p-185 */
+#endif
+
+//VBS
+//static inline __always_inline int
+//__ieee754_rem_pio2l(long double x, long double *y)
+
+static inline int
+__ieee754_rem_pio2l(long double x, long double *y)
+{
+	union IEEEl2bits u,u1;
+	long double z,w,t,r,fn;
+	double tx[3],ty[2];
+	int e0,ex,i,j,nx,n;
+	int16_t expsign;
+
+	u.e = x;
+	expsign = u.xbits.expsign;
+	ex = expsign & 0x7fff;
+	if (ex < BIAS + 25 || (ex == BIAS + 25 && u.bits.manh < 0xc90fdaa2)) {
+	    /* |x| ~< 2^25*(pi/2), medium size */
+	    /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
+	    fn = x*invpio2+0x1.8p63;
+	    fn = fn-0x1.8p63;
+#ifdef HAVE_EFFICIENT_IRINT
+	    n  = irint(fn);
+#else
+	    n  = fn;
+#endif
+	    r  = x-fn*pio2_1;
+	    w  = fn*pio2_1t;	/* 1st round good to 102 bit */
+	    {
+		union IEEEl2bits u2;
+	        int ex1;
+	        j  = ex;
+	        y[0] = r-w;
+		u2.e = y[0];
+		ex1 = u2.xbits.expsign & 0x7fff;
+	        i = j-ex1;
+	        if(i>22) {  /* 2nd iteration needed, good to 141 */
+		    t  = r;
+		    w  = fn*pio2_2;
+		    r  = t-w;
+		    w  = fn*pio2_2t-((t-r)-w);
+		    y[0] = r-w;
+		    u2.e = y[0];
+		    ex1 = u2.xbits.expsign & 0x7fff;
+		    i = j-ex1;
+		    if(i>61) {	/* 3rd iteration need, 180 bits acc */
+		    	t  = r;	/* will cover all possible cases */
+		    	w  = fn*pio2_3;
+		    	r  = t-w;
+		    	w  = fn*pio2_3t-((t-r)-w);
+		    	y[0] = r-w;
+		    }
+		}
+	    }
+	    y[1] = (r-y[0])-w;
+	    return n;
+	}
+    /*
+     * all other (large) arguments
+     */
+	if(ex==0x7fff) {		/* x is inf or NaN */
+	    y[0]=y[1]=x-x; return 0;
+	}
+    /* set z = scalbn(|x|,ilogb(x)-23) */
+	u1.e = x;
+	e0 = ex - BIAS - 23;		/* e0 = ilogb(|x|)-23; */
+	u1.xbits.expsign = ex - e0;
+	z = u1.e;
+	for(i=0;i<2;i++) {
+		tx[i] = (double)((int32_t)(z));
+		z     = (z-tx[i])*two24;
+	}
+	tx[2] = z;
+	nx = 3;
+	while(tx[nx-1]==zero) nx--;	/* skip zero term */
+	n  =  __kernel_rem_pio2(tx,ty,e0,nx,2);
+	r = (long double)ty[0] + ty[1];
+	w = ty[1] - (r - ty[0]);
+	if(expsign<0) {y[0] = -r; y[1] = -w; return -n;}
+	y[0] = r; y[1] = w; return n;
+}
--- a/src/openlibm/ld80/e_sinhl.c
+++ b/src/openlibm/ld80/e_sinhl.c
@@ -0,0 +1,76 @@
+/* @(#)e_sinh.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.
+ * ====================================================
+ */
+
+/* sinhl(x)
+ * Method :
+ * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
+ *	1. Replace x by |x| (sinhl(-x) = -sinhl(x)).
+ *	2.
+ *		                                     E + E/(E+1)
+ *	    0        <= x <= 25     :  sinhl(x) := --------------, E=expm1l(x)
+ *			       			         2
+ *
+ *	    25       <= x <= lnovft :  sinhl(x) := expl(x)/2
+ *	    lnovft   <= x <= ln2ovft:  sinhl(x) := expl(x/2)/2 * expl(x/2)
+ *	    ln2ovft  <  x	    :  sinhl(x) := x*shuge (overflow)
+ *
+ * Special cases:
+ *	sinhl(x) is |x| if x is +INF, -INF, or NaN.
+ *	only sinhl(0)=0 is exact for finite x.
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static const long double one = 1.0, shuge = 1.0e4931L;
+
+long double
+sinhl(long double x)
+{
+	long double t,w,h;
+	u_int32_t jx,ix,i0,i1;
+
+    /* Words of |x|. */
+	GET_LDOUBLE_WORDS(jx,i0,i1,x);
+	ix = jx&0x7fff;
+
+    /* x is INF or NaN */
+	if(ix==0x7fff) return x+x;
+
+	h = 0.5;
+	if (jx & 0x8000) h = -h;
+    /* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */
+	if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x|<25 */
+	    if (ix<0x3fdf)		 /* |x|<2**-32 */
+		if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
+	    t = expm1l(fabsl(x));
+	    if(ix<0x3fff) return h*(2.0*t-t*t/(t+one));
+	    return h*(t+t/(t+one));
+	}
+
+    /* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */
+	if (ix < 0x400c || (ix == 0x400c && i0 < 0xb17217f7))
+		return h*expl(fabsl(x));
+
+    /* |x| in [log(maxdouble), overflowthreshold] */
+	if (ix<0x400c || (ix == 0x400c && (i0 < 0xb174ddc0
+					   || (i0 == 0xb174ddc0
+					       && i1 <= 0x31aec0ea)))) {
+	    w = expl(0.5*fabsl(x));
+	    t = h*w;
+	    return t*w;
+	}
+
+    /* |x| > overflowthreshold, sinhl(x) overflow */
+	return x*shuge;
+}
--- a/src/openlibm/ld80/e_sinhl.c.o
+++ b/src/openlibm/ld80/e_sinhl.c.o
--- a/src/openlibm/ld80/e_tgammal.c
+++ b/src/openlibm/ld80/e_tgammal.c
@@ -0,0 +1,313 @@
+/*	$OpenBSD: e_tgammal.c,v 1.4 2013/11/12 20:35:19 martynas Exp $	*/
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * 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.
+ */
+
+/*							tgammal.c
+ *
+ *	Gamma function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, tgammal();
+ *
+ * y = tgammal( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns gamma function of the argument.  The result is correctly
+ * signed.  This variable is also filled in by the logarithmic gamma
+ * function lgamma().
+ *
+ * Arguments |x| <= 13 are reduced by recurrence and the function
+ * approximated by a rational function of degree 7/8 in the
+ * interval (2,3).  Large arguments are handled by Stirling's
+ * formula. Large negative arguments are made positive using
+ * a reflection formula.
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     -40,+40      10000       3.6e-19     7.9e-20
+ *    IEEE    -1755,+1755   10000       4.8e-18     6.5e-19
+ *
+ * Accuracy for large arguments is dominated by error in powl().
+ *
+ */
+
+#include <float.h>
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+/*
+tgamma(x+2)  = tgamma(x+2) P(x)/Q(x)
+0 <= x <= 1
+Relative error
+n=7, d=8
+Peak error =  1.83e-20
+Relative error spread =  8.4e-23
+*/
+
+static long double P[8] = {
+ 4.212760487471622013093E-5L,
+ 4.542931960608009155600E-4L,
+ 4.092666828394035500949E-3L,
+ 2.385363243461108252554E-2L,
+ 1.113062816019361559013E-1L,
+ 3.629515436640239168939E-1L,
+ 8.378004301573126728826E-1L,
+ 1.000000000000000000009E0L,
+};
+static long double Q[9] = {
+-1.397148517476170440917E-5L,
+ 2.346584059160635244282E-4L,
+-1.237799246653152231188E-3L,
+-7.955933682494738320586E-4L,
+ 2.773706565840072979165E-2L,
+-4.633887671244534213831E-2L,
+-2.243510905670329164562E-1L,
+ 4.150160950588455434583E-1L,
+ 9.999999999999999999908E-1L,
+};
+
+/*
+static long double P[] = {
+-3.01525602666895735709e0L,
+-3.25157411956062339893e1L,
+-2.92929976820724030353e2L,
+-1.70730828800510297666e3L,
+-7.96667499622741999770e3L,
+-2.59780216007146401957e4L,
+-5.99650230220855581642e4L,
+-7.15743521530849602425e4L
+};
+static long double Q[] = {
+ 1.00000000000000000000e0L,
+-1.67955233807178858919e1L,
+ 8.85946791747759881659e1L,
+ 5.69440799097468430177e1L,
+-1.98526250512761318471e3L,
+ 3.31667508019495079814e3L,
+ 1.60577839621734713377e4L,
+-2.97045081369399940529e4L,
+-7.15743521530849602412e4L
+};
+*/
+#define MAXGAML 1755.455L
+/*static const long double LOGPI = 1.14472988584940017414L;*/
+
+/* Stirling's formula for the gamma function
+tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
+z(x) = x
+13 <= x <= 1024
+Relative error
+n=8, d=0
+Peak error =  9.44e-21
+Relative error spread =  8.8e-4
+*/
+
+static long double STIR[9] = {
+ 7.147391378143610789273E-4L,
+-2.363848809501759061727E-5L,
+-5.950237554056330156018E-4L,
+ 6.989332260623193171870E-5L,
+ 7.840334842744753003862E-4L,
+-2.294719747873185405699E-4L,
+-2.681327161876304418288E-3L,
+ 3.472222222230075327854E-3L,
+ 8.333333333333331800504E-2L,
+};
+
+#define MAXSTIR 1024.0L
+static const long double SQTPI = 2.50662827463100050242E0L;
+
+/* 1/tgamma(x) = z P(z)
+ * z(x) = 1/x
+ * 0 < x < 0.03125
+ * Peak relative error 4.2e-23
+ */
+
+static long double S[9] = {
+-1.193945051381510095614E-3L,
+ 7.220599478036909672331E-3L,
+-9.622023360406271645744E-3L,
+-4.219773360705915470089E-2L,
+ 1.665386113720805206758E-1L,
+-4.200263503403344054473E-2L,
+-6.558780715202540684668E-1L,
+ 5.772156649015328608253E-1L,
+ 1.000000000000000000000E0L,
+};
+
+/* 1/tgamma(-x) = z P(z)
+ * z(x) = 1/x
+ * 0 < x < 0.03125
+ * Peak relative error 5.16e-23
+ * Relative error spread =  2.5e-24
+ */
+
+static long double SN[9] = {
+ 1.133374167243894382010E-3L,
+ 7.220837261893170325704E-3L,
+ 9.621911155035976733706E-3L,
+-4.219773343731191721664E-2L,
+-1.665386113944413519335E-1L,
+-4.200263503402112910504E-2L,
+ 6.558780715202536547116E-1L,
+ 5.772156649015328608727E-1L,
+-1.000000000000000000000E0L,
+};
+
+static const long double PIL = 3.1415926535897932384626L;
+
+static long double stirf ( long double );
+
+/* Gamma function computed by Stirling's formula.
+ */
+static long double stirf(long double x)
+{
+long double y, w, v;
+
+w = 1.0L/x;
+/* For large x, use rational coefficients from the analytical expansion.  */
+if( x > 1024.0L )
+	w = (((((6.97281375836585777429E-5L * w
+		+ 7.84039221720066627474E-4L) * w
+		- 2.29472093621399176955E-4L) * w
+		- 2.68132716049382716049E-3L) * w
+		+ 3.47222222222222222222E-3L) * w
+		+ 8.33333333333333333333E-2L) * w
+		+ 1.0L;
+else
+	w = 1.0L + w * __polevll( w, STIR, 8 );
+y = expl(x);
+if( x > MAXSTIR )
+	{ /* Avoid overflow in pow() */
+	v = powl( x, 0.5L * x - 0.25L );
+	y = v * (v / y);
+	}
+else
+	{
+	y = powl( x, x - 0.5L ) / y;
+	}
+y = SQTPI * y * w;
+return( y );
+}
+
+long double
+tgammal(long double x)
+{
+long double p, q, z;
+int i;
+
+if( isnan(x) )
+	return(NAN);
+if(x == INFINITY)
+	return(INFINITY);
+if(x == -INFINITY)
+	return(x - x);
+if( x == 0.0L )
+	return( 1.0L / x );
+q = fabsl(x);
+
+if( q > 13.0L )
+	{
+	int sign = 1;
+	if( q > MAXGAML )
+		goto goverf;
+	if( x < 0.0L )
+		{
+		p = floorl(q);
+		if( p == q )
+			return (x - x) / (x - x);
+		i = p;
+		if( (i & 1) == 0 )
+			sign = -1;
+		z = q - p;
+		if( z > 0.5L )
+			{
+			p += 1.0L;
+			z = q - p;
+			}
+		z = q * sinl( PIL * z );
+		z = fabsl(z) * stirf(q);
+		if( z <= PIL/LDBL_MAX )
+			{
+goverf:
+			return( sign * INFINITY);
+			}
+		z = PIL/z;
+		}
+	else
+		{
+		z = stirf(x);
+		}
+	return( sign * z );
+	}
+
+z = 1.0L;
+while( x >= 3.0L )
+	{
+	x -= 1.0L;
+	z *= x;
+	}
+
+while( x < -0.03125L )
+	{
+	z /= x;
+	x += 1.0L;
+	}
+
+if( x <= 0.03125L )
+	goto small;
+
+while( x < 2.0L )
+	{
+	z /= x;
+	x += 1.0L;
+	}
+
+if( x == 2.0L )
+	return(z);
+
+x -= 2.0L;
+p = __polevll( x, P, 7 );
+q = __polevll( x, Q, 8 );
+z = z * p / q;
+return z;
+
+small:
+if( x == 0.0L )
+	return (x - x) / (x - x);
+else
+	{
+	if( x < 0.0L )
+		{
+		x = -x;
+		q = z / (x * __polevll( x, SN, 8 ));
+		}
+	else
+		q = z / (x * __polevll( x, S, 8 ));
+	}
+return q;
+}
--- a/src/openlibm/ld80/e_tgammal.c.o
+++ b/src/openlibm/ld80/e_tgammal.c.o
--- a/src/openlibm/ld80/invtrig.c
+++ b/src/openlibm/ld80/invtrig.c
@@ -0,0 +1,82 @@
+/*-
+ * 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 "cdefs-compat.h"
+//__FBSDID("$FreeBSD: src/lib/msun/ld80/invtrig.c,v 1.1 2008/07/31 22:41:26 das Exp $");
+
+#include "ld80/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;
--- a/src/openlibm/ld80/invtrig.c.o
+++ b/src/openlibm/ld80/invtrig.c.o
--- a/src/openlibm/ld80/invtrig.h
+++ b/src/openlibm/ld80/invtrig.h
@@ -0,0 +1,114 @@
+/*-
+ * 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 <float.h>
+
+#include <openlibm_math.h>
+
+#define	BIAS		(LDBL_MAX_EXP - 1)
+#define	MANH_SIZE	LDBL_MANH_SIZE
+
+/* 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 */
+#define	THRESH	((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_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
--- a/src/openlibm/ld80/k_cosl.c
+++ b/src/openlibm/ld80/k_cosl.c
@@ -0,0 +1,78 @@
+/* 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.
+ * ====================================================
+ */
+
+#include "cdefs-compat.h"
+//__FBSDID("$FreeBSD: src/lib/msun/ld80/k_cosl.c,v 1.1 2008/02/17 07:32:14 das Exp $");
+
+/*
+ * ld80 version of k_cos.c.  See ../src/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));
+}
--- a/src/openlibm/ld80/k_cosl.c.o
+++ b/src/openlibm/ld80/k_cosl.c.o
--- a/src/openlibm/ld80/k_sinl.c
+++ b/src/openlibm/ld80/k_sinl.c
@@ -0,0 +1,62 @@
+/* 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.
+ * ====================================================
+ */
+
+#include "cdefs-compat.h"
+//__FBSDID("$FreeBSD: src/lib/msun/ld80/k_sinl.c,v 1.1 2008/02/17 07:32:14 das Exp $");
+
+/*
+ * ld80 version of k_sin.c.  See ../src/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);
+}
--- a/src/openlibm/ld80/k_sinl.c.o
+++ b/src/openlibm/ld80/k_sinl.c.o
--- a/src/openlibm/ld80/k_tanl.c
+++ b/src/openlibm/ld80/k_tanl.c
@@ -0,0 +1,125 @@
+/* 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.
+ * ====================================================
+ */
+
+#include "cdefs-compat.h"
+//__FBSDID("$FreeBSD: src/lib/msun/ld80/k_tanl.c,v 1.3 2008/02/18 15:39:52 bde Exp $");
+
+/*
+ * ld80 version of k_tan.c.  See ../src/k_tan.c for most comments.
+ */
+
+#include <openlibm_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);
+	}
+}
--- a/src/openlibm/ld80/k_tanl.c.o
+++ b/src/openlibm/ld80/k_tanl.c.o
--- a/src/openlibm/ld80/s_asinhl.c
+++ b/src/openlibm/ld80/s_asinhl.c
@@ -0,0 +1,54 @@
+/* @(#)s_asinh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* asinhl(x)
+ * Method :
+ *	Based on
+ *		asinhl(x) = signl(x) * logl [ |x| + sqrtl(x*x+1) ]
+ *	we have
+ *	asinhl(x) := x  if  1+x*x=1,
+ *		  := signl(x)*(logl(x)+ln2)) for large |x|, else
+ *		  := signl(x)*logl(2|x|+1/(|x|+sqrtl(x*x+1))) if|x|>2, else
+ *		  := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2)))
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static const long double
+one =  1.000000000000000000000e+00L, /* 0x3FFF, 0x00000000, 0x00000000 */
+ln2 =  6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */
+huge=  1.000000000000000000e+4900L;
+
+long double
+asinhl(long double x)
+{
+	long double t,w;
+	int32_t hx,ix;
+	GET_LDOUBLE_EXP(hx,x);
+	ix = hx&0x7fff;
+	if(ix==0x7fff) return x+x;	/* x is inf or NaN */
+	if(ix< 0x3fde) {	/* |x|<2**-34 */
+	    if(huge+x>one) return x;	/* return x inexact except 0 */
+	}
+	if(ix>0x4020) {		/* |x| > 2**34 */
+	    w = logl(fabsl(x))+ln2;
+	} else if (ix>0x4000) {	/* 2**34 > |x| > 2.0 */
+	    t = fabsl(x);
+	    w = logl(2.0*t+one/(sqrtl(x*x+one)+t));
+	} else {		/* 2.0 > |x| > 2**-28 */
+	    t = x*x;
+	    w =log1pl(fabsl(x)+t/(one+sqrtl(one+t)));
+	}
+	if(hx&0x8000) return -w; else return w;
+}
--- a/src/openlibm/ld80/s_asinhl.c.o
+++ b/src/openlibm/ld80/s_asinhl.c.o
--- a/src/openlibm/ld80/s_ceill.c
+++ b/src/openlibm/ld80/s_ceill.c
@@ -0,0 +1,78 @@
+/* @(#)s_ceil.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.
+ * ====================================================
+ */
+
+/*
+ * ceill(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to ceil(x).
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static const long double huge = 1.0e4930L;
+
+long double
+ceill(long double x)
+{
+	int32_t i1,jj0;
+	u_int32_t i,j,se,i0,sx;
+	GET_LDOUBLE_WORDS(se,i0,i1,x);
+	sx = (se>>15)&1;
+	jj0 = (se&0x7fff)-0x3fff;
+	if(jj0<31) {
+	    if(jj0<0) {	/* raise inexact if x != 0 */
+		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
+		    if(sx) {se=0x8000;i0=0;i1=0;}
+		    else if((i0|i1)!=0) { se=0x3fff;i0=0;i1=0;}
+		}
+	    } else {
+		i = (0x7fffffff)>>jj0;
+		if(((i0&i)|i1)==0) return x; /* x is integral */
+		if(huge+x>0.0) {	/* raise inexact flag */
+		    if(sx==0) {
+			if (jj0>0 && (i0+(0x80000000>>jj0))>i0)
+			  i0+=0x80000000>>jj0;
+			else
+			  {
+			    i = 0x7fffffff;
+			    ++se;
+			  }
+		    }
+		    i0 &= (~i); i1=0;
+		}
+	    }
+	} else if (jj0>62) {
+	    if(jj0==0x4000) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	} else {
+	    i = ((u_int32_t)(0xffffffff))>>(jj0-31);
+	    if((i1&i)==0) return x;	/* x is integral */
+	    if(huge+x>0.0) {		/* raise inexact flag */
+		if(sx==0) {
+		    if(jj0==31) i0+=1;
+		    else {
+			j = i1 + (1<<(63-jj0));
+			if(j<i1) i0+=1;	/* got a carry */
+			i1 = j;
+		    }
+		}
+		i1 &= (~i);
+	    }
+	}
+	SET_LDOUBLE_WORDS(x,se,i0,i1);
+	return x;
+}
--- a/src/openlibm/ld80/s_erfl.c
+++ b/src/openlibm/ld80/s_erfl.c
@@ -0,0 +1,430 @@
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * 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.
+ */
+
+/* double erf(double x)
+ * double erfc(double x)
+ *			     x
+ *		      2      |\
+ *     erf(x)  =  ---------  | exp(-t*t)dt
+ *		   sqrt(pi) \|
+ *			     0
+ *
+ *     erfc(x) =  1-erf(x)
+ *  Note that
+ *		erf(-x) = -erf(x)
+ *		erfc(-x) = 2 - erfc(x)
+ *
+ * Method:
+ *	1. For |x| in [0, 0.84375]
+ *	    erf(x)  = x + x*R(x^2)
+ *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
+ *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
+ *	   Remark. The formula is derived by noting
+ *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
+ *	   and that
+ *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
+ *	   is close to one. The interval is chosen because the fix
+ *	   point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
+ *	   near 0.6174), and by some experiment, 0.84375 is chosen to
+ *	   guarantee the error is less than one ulp for erf.
+ *
+ *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
+ *         c = 0.84506291151 rounded to single (24 bits)
+ *	erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
+ *	erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
+ *			  1+(c+P1(s)/Q1(s))    if x < 0
+ *	   Remark: here we use the taylor series expansion at x=1.
+ *		erf(1+s) = erf(1) + s*Poly(s)
+ *			 = 0.845.. + P1(s)/Q1(s)
+ *	   Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
+ *
+ *      3. For x in [1.25,1/0.35(~2.857143)],
+ *	erfc(x) = (1/x)*exp(-x*x-0.5625+R1(z)/S1(z))
+ *              z=1/x^2
+ *	erf(x)  = 1 - erfc(x)
+ *
+ *      4. For x in [1/0.35,107]
+ *	erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
+ *			= 2.0 - (1/x)*exp(-x*x-0.5625+R2(z)/S2(z))
+ *                             if -6.666<x<0
+ *			= 2.0 - tiny		(if x <= -6.666)
+ *              z=1/x^2
+ *	erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6.666, else
+ *	erf(x)  = sign(x)*(1.0 - tiny)
+ *      Note1:
+ *	   To compute exp(-x*x-0.5625+R/S), let s be a single
+ *	   precision number and s := x; then
+ *		-x*x = -s*s + (s-x)*(s+x)
+ *	        exp(-x*x-0.5626+R/S) =
+ *			exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
+ *      Note2:
+ *	   Here 4 and 5 make use of the asymptotic series
+ *			  exp(-x*x)
+ *		erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
+ *			  x*sqrt(pi)
+ *
+ *      5. For inf > x >= 107
+ *	erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
+ *	erfc(x) = tiny*tiny (raise underflow) if x > 0
+ *			= 2 - tiny if x<0
+ *
+ *      7. Special case:
+ *	erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
+ *	erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
+ *		erfc/erf(NaN) is NaN
+ */
+
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static const long double
+tiny = 1e-4931L,
+  half = 0.5L,
+  one = 1.0L,
+  two = 2.0L,
+	/* c = (float)0.84506291151 */
+  erx = 0.845062911510467529296875L,
+/*
+ * Coefficients for approximation to  erf on [0,0.84375]
+ */
+  /* 2/sqrt(pi) - 1 */
+  efx = 1.2837916709551257389615890312154517168810E-1L,
+  /* 8 * (2/sqrt(pi) - 1) */
+  efx8 = 1.0270333367641005911692712249723613735048E0L,
+
+  pp[6] = {
+    1.122751350964552113068262337278335028553E6L,
+    -2.808533301997696164408397079650699163276E6L,
+    -3.314325479115357458197119660818768924100E5L,
+    -6.848684465326256109712135497895525446398E4L,
+    -2.657817695110739185591505062971929859314E3L,
+    -1.655310302737837556654146291646499062882E2L,
+  },
+
+  qq[6] = {
+    8.745588372054466262548908189000448124232E6L,
+    3.746038264792471129367533128637019611485E6L,
+    7.066358783162407559861156173539693900031E5L,
+    7.448928604824620999413120955705448117056E4L,
+    4.511583986730994111992253980546131408924E3L,
+    1.368902937933296323345610240009071254014E2L,
+    /* 1.000000000000000000000000000000000000000E0 */
+  },
+
+/*
+ * Coefficients for approximation to  erf  in [0.84375,1.25]
+ */
+/* erf(x+1) = 0.845062911510467529296875 + pa(x)/qa(x)
+   -0.15625 <= x <= +.25
+   Peak relative error 8.5e-22  */
+
+  pa[8] = {
+    -1.076952146179812072156734957705102256059E0L,
+     1.884814957770385593365179835059971587220E2L,
+    -5.339153975012804282890066622962070115606E1L,
+     4.435910679869176625928504532109635632618E1L,
+     1.683219516032328828278557309642929135179E1L,
+    -2.360236618396952560064259585299045804293E0L,
+     1.852230047861891953244413872297940938041E0L,
+     9.394994446747752308256773044667843200719E-2L,
+  },
+
+  qa[7] =  {
+    4.559263722294508998149925774781887811255E2L,
+    3.289248982200800575749795055149780689738E2L,
+    2.846070965875643009598627918383314457912E2L,
+    1.398715859064535039433275722017479994465E2L,
+    6.060190733759793706299079050985358190726E1L,
+    2.078695677795422351040502569964299664233E1L,
+    4.641271134150895940966798357442234498546E0L,
+    /* 1.000000000000000000000000000000000000000E0 */
+  },
+
+/*
+ * Coefficients for approximation to  erfc in [1.25,1/0.35]
+ */
+/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + ra(x^2)/sa(x^2))
+   1/2.85711669921875 < 1/x < 1/1.25
+   Peak relative error 3.1e-21  */
+
+    ra[] = {
+      1.363566591833846324191000679620738857234E-1L,
+      1.018203167219873573808450274314658434507E1L,
+      1.862359362334248675526472871224778045594E2L,
+      1.411622588180721285284945138667933330348E3L,
+      5.088538459741511988784440103218342840478E3L,
+      8.928251553922176506858267311750789273656E3L,
+      7.264436000148052545243018622742770549982E3L,
+      2.387492459664548651671894725748959751119E3L,
+      2.220916652813908085449221282808458466556E2L,
+    },
+
+    sa[] = {
+      -1.382234625202480685182526402169222331847E1L,
+      -3.315638835627950255832519203687435946482E2L,
+      -2.949124863912936259747237164260785326692E3L,
+      -1.246622099070875940506391433635999693661E4L,
+      -2.673079795851665428695842853070996219632E4L,
+      -2.880269786660559337358397106518918220991E4L,
+      -1.450600228493968044773354186390390823713E4L,
+      -2.874539731125893533960680525192064277816E3L,
+      -1.402241261419067750237395034116942296027E2L,
+      /* 1.000000000000000000000000000000000000000E0 */
+    },
+/*
+ * Coefficients for approximation to  erfc in [1/.35,107]
+ */
+/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rb(x^2)/sb(x^2))
+   1/6.6666259765625 < 1/x < 1/2.85711669921875
+   Peak relative error 4.2e-22  */
+    rb[] = {
+      -4.869587348270494309550558460786501252369E-5L,
+      -4.030199390527997378549161722412466959403E-3L,
+      -9.434425866377037610206443566288917589122E-2L,
+      -9.319032754357658601200655161585539404155E-1L,
+      -4.273788174307459947350256581445442062291E0L,
+      -8.842289940696150508373541814064198259278E0L,
+      -7.069215249419887403187988144752613025255E0L,
+      -1.401228723639514787920274427443330704764E0L,
+    },
+
+    sb[] = {
+      4.936254964107175160157544545879293019085E-3L,
+      1.583457624037795744377163924895349412015E-1L,
+      1.850647991850328356622940552450636420484E0L,
+      9.927611557279019463768050710008450625415E0L,
+      2.531667257649436709617165336779212114570E1L,
+      2.869752886406743386458304052862814690045E1L,
+      1.182059497870819562441683560749192539345E1L,
+      /* 1.000000000000000000000000000000000000000E0 */
+    },
+/* erfc(1/x) = x exp (-1/x^2 - 0.5625 + rc(x^2)/sc(x^2))
+   1/107 <= 1/x <= 1/6.6666259765625
+   Peak relative error 1.1e-21  */
+    rc[] = {
+      -8.299617545269701963973537248996670806850E-5L,
+      -6.243845685115818513578933902532056244108E-3L,
+      -1.141667210620380223113693474478394397230E-1L,
+      -7.521343797212024245375240432734425789409E-1L,
+      -1.765321928311155824664963633786967602934E0L,
+      -1.029403473103215800456761180695263439188E0L,
+    },
+
+    sc[] = {
+      8.413244363014929493035952542677768808601E-3L,
+      2.065114333816877479753334599639158060979E-1L,
+      1.639064941530797583766364412782135680148E0L,
+      4.936788463787115555582319302981666347450E0L,
+      5.005177727208955487404729933261347679090E0L,
+      /* 1.000000000000000000000000000000000000000E0 */
+    };
+
+long double
+erfl(long double x)
+{
+  long double R, S, P, Q, s, y, z, r;
+  int32_t ix, i;
+  u_int32_t se, i0, i1;
+
+  GET_LDOUBLE_WORDS (se, i0, i1, x);
+  ix = se & 0x7fff;
+
+  if (ix >= 0x7fff)
+    {				/* erf(nan)=nan */
+      i = ((se & 0xffff) >> 15) << 1;
+      return (long double) (1 - i) + one / x;	/* erf(+-inf)=+-1 */
+    }
+
+  ix = (ix << 16) | (i0 >> 16);
+  if (ix < 0x3ffed800) /* |x|<0.84375 */
+    {
+      if (ix < 0x3fde8000) /* |x|<2**-33 */
+	{
+	  if (ix < 0x00080000)
+	    return 0.125 * (8.0 * x + efx8 * x);	/*avoid underflow */
+	  return x + efx * x;
+	}
+      z = x * x;
+      r = pp[0] + z * (pp[1]
+	+ z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
+      s = qq[0] + z * (qq[1]
+	+ z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
+      y = r / s;
+      return x + x * y;
+    }
+  if (ix < 0x3fffa000) /* 1.25 */
+    {				/* 0.84375 <= |x| < 1.25 */
+      s = fabsl (x) - one;
+      P = pa[0] + s * (pa[1] + s * (pa[2]
+	+ s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
+      Q = qa[0] + s * (qa[1] + s * (qa[2]
+	+ s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
+      if ((se & 0x8000) == 0)
+	return erx + P / Q;
+      else
+	return -erx - P / Q;
+    }
+  if (ix >= 0x4001d555) /* 6.6666259765625 */
+    {				/* inf>|x|>=6.666 */
+      if ((se & 0x8000) == 0)
+	return one - tiny;
+      else
+	return tiny - one;
+    }
+  x = fabsl (x);
+  s = one / (x * x);
+  if (ix < 0x4000b6db) /* 2.85711669921875 */
+    {
+      R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
+	s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
+      S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
+	s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
+    }
+  else
+    {				/* |x| >= 1/0.35 */
+      R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
+	s * (rb[5] + s * (rb[6] + s * rb[7]))))));
+      S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
+	s * (sb[5] + s * (sb[6] + s))))));
+    }
+  z = x;
+  GET_LDOUBLE_WORDS (i, i0, i1, z);
+  i1 = 0;
+  SET_LDOUBLE_WORDS (z, i, i0, i1);
+  r =
+    expl (-z * z - 0.5625) * expl ((z - x) * (z + x) + R / S);
+  if ((se & 0x8000) == 0)
+    return one - r / x;
+  else
+    return r / x - one;
+}
+
+long double
+erfcl(long double x)
+{
+  int32_t hx, ix;
+  long double R, S, P, Q, s, y, z, r;
+  u_int32_t se, i0, i1;
+
+  GET_LDOUBLE_WORDS (se, i0, i1, x);
+  ix = se & 0x7fff;
+  if (ix >= 0x7fff)
+    {				/* erfc(nan)=nan */
+      /* erfc(+-inf)=0,2 */
+      return (long double) (((se & 0xffff) >> 15) << 1) + one / x;
+    }
+
+  ix = (ix << 16) | (i0 >> 16);
+  if (ix < 0x3ffed800) /* |x|<0.84375 */
+    {
+      if (ix < 0x3fbe0000) /* |x|<2**-65 */
+	return one - x;
+      z = x * x;
+      r = pp[0] + z * (pp[1]
+	+ z * (pp[2] + z * (pp[3] + z * (pp[4] + z * pp[5]))));
+      s = qq[0] + z * (qq[1]
+	+ z * (qq[2] + z * (qq[3] + z * (qq[4] + z * (qq[5] + z)))));
+      y = r / s;
+      if (ix < 0x3ffd8000) /* x<1/4 */
+	{
+	  return one - (x + x * y);
+	}
+      else
+	{
+	  r = x * y;
+	  r += (x - half);
+	  return half - r;
+	}
+    }
+  if (ix < 0x3fffa000) /* 1.25 */
+    {				/* 0.84375 <= |x| < 1.25 */
+      s = fabsl (x) - one;
+      P = pa[0] + s * (pa[1] + s * (pa[2]
+	+ s * (pa[3] + s * (pa[4] + s * (pa[5] + s * (pa[6] + s * pa[7]))))));
+      Q = qa[0] + s * (qa[1] + s * (qa[2]
+	+ s * (qa[3] + s * (qa[4] + s * (qa[5] + s * (qa[6] + s))))));
+      if ((se & 0x8000) == 0)
+	{
+	  z = one - erx;
+	  return z - P / Q;
+	}
+      else
+	{
+	  z = erx + P / Q;
+	  return one + z;
+	}
+    }
+  if (ix < 0x4005d600) /* 107 */
+    {				/* |x|<107 */
+      x = fabsl (x);
+      s = one / (x * x);
+      if (ix < 0x4000b6db) /* 2.85711669921875 */
+	{			/* |x| < 1/.35 ~ 2.857143 */
+	  R = ra[0] + s * (ra[1] + s * (ra[2] + s * (ra[3] + s * (ra[4] +
+	    s * (ra[5] + s * (ra[6] + s * (ra[7] + s * ra[8])))))));
+	  S = sa[0] + s * (sa[1] + s * (sa[2] + s * (sa[3] + s * (sa[4] +
+	    s * (sa[5] + s * (sa[6] + s * (sa[7] + s * (sa[8] + s))))))));
+	}
+      else if (ix < 0x4001d555) /* 6.6666259765625 */
+	{			/* 6.666 > |x| >= 1/.35 ~ 2.857143 */
+	  R = rb[0] + s * (rb[1] + s * (rb[2] + s * (rb[3] + s * (rb[4] +
+	    s * (rb[5] + s * (rb[6] + s * rb[7]))))));
+	  S = sb[0] + s * (sb[1] + s * (sb[2] + s * (sb[3] + s * (sb[4] +
+	    s * (sb[5] + s * (sb[6] + s))))));
+	}
+      else
+	{			/* |x| >= 6.666 */
+	  if (se & 0x8000)
+	    return two - tiny;	/* x < -6.666 */
+
+	  R = rc[0] + s * (rc[1] + s * (rc[2] + s * (rc[3] +
+						    s * (rc[4] + s * rc[5]))));
+	  S = sc[0] + s * (sc[1] + s * (sc[2] + s * (sc[3] +
+						    s * (sc[4] + s))));
+	}
+      z = x;
+      GET_LDOUBLE_WORDS (hx, i0, i1, z);
+      i1 = 0;
+      i0 &= 0xffffff00;
+      SET_LDOUBLE_WORDS (z, hx, i0, i1);
+      r = expl (-z * z - 0.5625) *
+	expl ((z - x) * (z + x) + R / S);
+      if ((se & 0x8000) == 0)
+	return r / x;
+      else
+	return two - r / x;
+    }
+  else
+    {
+      if ((se & 0x8000) == 0)
+	return tiny * tiny;
+      else
+	return two - tiny;
+    }
+}
--- a/src/openlibm/ld80/s_erfl.c.o
+++ b/src/openlibm/ld80/s_erfl.c.o
--- a/src/openlibm/ld80/s_exp2l.c
+++ b/src/openlibm/ld80/s_exp2l.c
@@ -0,0 +1,295 @@
+/*-
+ * 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 "cdefs-compat.h"
+//__FBSDID("$FreeBSD: src/lib/msun/ld80/s_exp2l.c,v 1.3 2008/02/13 10:44:44 bde Exp $");
+
+#include <float.h>
+#include <stdint.h>
+
+#include "bsd_cdefs.h"
+#include "amd64/bsd_ieeefp.h"
+
+#include <openlibm_math.h>
+
+#include "math_private.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.
+ */
+OLM_DLLEXPORT long double
+exp2l(long double x)
+{
+	union IEEEl2bits u, v;
+	long double r, twopk, twopkp10000, z;
+	uint32_t hx, ix, i0;
+	int k;
+
+	/* Filter out exceptional cases. */
+	u.e = x;
+	hx = u.xbits.expsign;
+	ix = hx & EXPMASK;
+	if (ix >= BIAS + 14) {		/* |x| >= 16384 or x is NaN */
+		if (ix == BIAS + LDBL_MAX_EXP) {
+			if (u.xbits.man != 1ULL << 63 || (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);
+	}
+
+#ifdef __i386__
+	/*
+	 * The default precision on i386 is 53 bits, so long doubles are
+	 * broken. Call exp2() to get an accurate (double precision) result.
+	 */
+	if (__fpgetprec() != FP_PE)
+		return (exp2(x));
+#endif
+
+
+	/*
+	 * 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.manl + TBLSIZE / 2;
+	k = (int)i0 >> TBLBITS;
+	i0 = (i0 & (TBLSIZE - 1)) << 1;
+	u.e -= redux;
+	z = x - u.e;
+	v.xbits.man = 1ULL << 63;
+	if (k >= LDBL_MIN_EXP) {
+		v.xbits.expsign = LDBL_MAX_EXP - 1 + k;
+		twopk = v.e;
+	} else {
+		v.xbits.expsign = LDBL_MAX_EXP - 1 + k + 10000;
+		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);
+	}
+}
--- a/src/openlibm/ld80/s_exp2l.c.o
+++ b/src/openlibm/ld80/s_exp2l.c.o
--- a/src/openlibm/ld80/s_expm1l.c
+++ b/src/openlibm/ld80/s_expm1l.c
@@ -0,0 +1,138 @@
+/*	$OpenBSD: s_expm1l.c,v 1.2 2011/07/20 21:02:51 martynas Exp $	*/
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * 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.
+ */
+
+/*							expm1l.c
+ *
+ *	Exponential function, minus 1
+ *      Long double precision
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, expm1l();
+ *
+ * y = expm1l( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns e (2.71828...) raised to the x power, minus 1.
+ *
+ * Range reduction is accomplished by separating the argument
+ * into an integer k and fraction f such that
+ *
+ *     x    k  f
+ *    e  = 2  e.
+ *
+ * An expansion x + .5 x^2 + x^3 R(x) approximates exp(f) - 1
+ * in the basic range [-0.5 ln 2, 0.5 ln 2].
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE    -45,+MAXLOG   200,000     1.2e-19     2.5e-20
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * expm1l overflow   x > MAXLOG         MAXNUM
+ *
+ */
+
+#include <openlibm_math.h>
+
+static const long double MAXLOGL = 1.1356523406294143949492E4L;
+
+/* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
+   -.5 ln 2  <  x  <  .5 ln 2
+   Theoretical peak relative error = 3.4e-22  */
+
+static const long double
+  P0 = -1.586135578666346600772998894928250240826E4L,
+  P1 =  2.642771505685952966904660652518429479531E3L,
+  P2 = -3.423199068835684263987132888286791620673E2L,
+  P3 =  1.800826371455042224581246202420972737840E1L,
+  P4 = -5.238523121205561042771939008061958820811E-1L,
+
+  Q0 = -9.516813471998079611319047060563358064497E4L,
+  Q1 =  3.964866271411091674556850458227710004570E4L,
+  Q2 = -7.207678383830091850230366618190187434796E3L,
+  Q3 =  7.206038318724600171970199625081491823079E2L,
+  Q4 = -4.002027679107076077238836622982900945173E1L,
+  /* Q5 = 1.000000000000000000000000000000000000000E0 */
+
+/* C1 + C2 = ln 2 */
+C1 = 6.93145751953125E-1L,
+C2 = 1.428606820309417232121458176568075500134E-6L,
+/* ln 2^-65 */
+minarg = -4.5054566736396445112120088E1L;
+static const long double huge = 0x1p10000L;
+
+long double
+expm1l(long double x)
+{
+long double px, qx, xx;
+int k;
+
+/* Overflow.  */
+if (x > MAXLOGL)
+  return (huge*huge);	/* overflow */
+
+if (x == 0.0)
+  return x;
+
+/* Minimum value.  */
+if (x < minarg)
+  return -1.0L;
+
+xx = C1 + C2;
+
+/* Express x = ln 2 (k + remainder), remainder not exceeding 1/2. */
+px = floorl (0.5 + x / xx);
+k = px;
+/* remainder times ln 2 */
+x -= px * C1;
+x -= px * C2;
+
+/* Approximate exp(remainder ln 2).  */
+px = (((( P4 * x
+	 + P3) * x
+	+ P2) * x
+       + P1) * x
+      + P0) * x;
+
+qx = (((( x
+	 + Q4) * x
+	+ Q3) * x
+       + Q2) * x
+      + Q1) * x
+     + Q0;
+
+xx = x * x;
+qx = x + (0.5 * xx + xx * px / qx);
+
+/* exp(x) = exp(k ln 2) exp(remainder ln 2) = 2^k exp(remainder ln 2).
+   We have qx = exp(remainder ln 2) - 1, so
+   exp(x) - 1  =  2^k (qx + 1) - 1  =  2^k qx + 2^k - 1.  */
+px = ldexpl(1.0L, k);
+x = px * qx + (px - 1.0);
+return x;
+}
--- a/src/openlibm/ld80/s_expm1l.c.o
+++ b/src/openlibm/ld80/s_expm1l.c.o
--- a/src/openlibm/ld80/s_floorl.c
+++ b/src/openlibm/ld80/s_floorl.c
@@ -0,0 +1,80 @@
+/* @(#)s_floor.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.
+ * ====================================================
+ */
+
+/*
+ * floorl(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to floor(x).
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static const long double huge = 1.0e4930L;
+
+long double
+floorl(long double x)
+{
+	int32_t i1,jj0;
+	u_int32_t i,j,se,i0,sx;
+	GET_LDOUBLE_WORDS(se,i0,i1,x);
+	sx = (se>>15)&1;
+	jj0 = (se&0x7fff)-0x3fff;
+	if(jj0<31) {
+	    if(jj0<0) {	/* raise inexact if x != 0 */
+		if(huge+x>0.0) {
+		    if(sx==0)
+			return 0.0L;
+		    else if(((se&0x7fff)|i0|i1)!=0)
+			return -1.0L;
+		}
+	    } else {
+		i = (0x7fffffff)>>jj0;
+		if(((i0&i)|i1)==0) return x; /* x is integral */
+		if(huge+x>0.0) {	/* raise inexact flag */
+		    if(sx) {
+			if (jj0>0 && (i0+(0x80000000>>jj0))>i0)
+			  i0 += (0x80000000)>>jj0;
+			else
+			  {
+			    i = 0x7fffffff;
+			    ++se;
+			  }
+		    }
+		    i0 &= (~i); i1=0;
+		}
+	    }
+	} else if (jj0>62) {
+	    if(jj0==0x4000) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	} else {
+	    i = ((u_int32_t)(0xffffffff))>>(jj0-31);
+	    if((i1&i)==0) return x;	/* x is integral */
+	    if(huge+x>0.0) {		/* raise inexact flag */
+		if(sx) {
+		    if(jj0==31) i0+=1;
+		    else {
+			j = i1+(1<<(63-jj0));
+			if(j<i1) i0 +=1 ;	/* got a carry */
+			i1=j;
+		    }
+		}
+		i1 &= (~i);
+	    }
+	}
+	SET_LDOUBLE_WORDS(x,se,i0,i1);
+	return x;
+}
--- a/src/openlibm/ld80/s_log1pl.c
+++ b/src/openlibm/ld80/s_log1pl.c
@@ -0,0 +1,191 @@
+/*	$OpenBSD: s_log1pl.c,v 1.3 2013/11/12 20:35:19 martynas Exp $	*/
+
+/*
+ * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
+ *
+ * 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.
+ */
+
+/*							log1pl.c
+ *
+ *      Relative error logarithm
+ *	Natural logarithm of 1+x, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, log1pl();
+ *
+ * y = log1pl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base e (2.718...) logarithm of 1+x.
+ *
+ * The argument 1+x is separated into its exponent and fractional
+ * parts.  If the exponent is between -1 and +1, the logarithm
+ * of the fraction is approximated by
+ *
+ *     log(1+x) = x - 0.5 x^2 + x^3 P(x)/Q(x).
+ *
+ * Otherwise, setting  z = 2(x-1)/x+1),
+ *
+ *     log(x) = z + z^3 P(z)/Q(z).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     -1.0, 9.0    100000      8.2e-20    2.5e-20
+ *
+ * ERROR MESSAGES:
+ *
+ * log singularity:  x-1 = 0; returns -INFINITY
+ * log domain:       x-1 < 0; returns NAN
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+/* Coefficients for log(1+x) = x - x^2 / 2 + x^3 P(x)/Q(x)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 2.32e-20
+ */
+
+static long double P[] = {
+ 4.5270000862445199635215E-5L,
+ 4.9854102823193375972212E-1L,
+ 6.5787325942061044846969E0L,
+ 2.9911919328553073277375E1L,
+ 6.0949667980987787057556E1L,
+ 5.7112963590585538103336E1L,
+ 2.0039553499201281259648E1L,
+};
+static long double Q[] = {
+/* 1.0000000000000000000000E0,*/
+ 1.5062909083469192043167E1L,
+ 8.3047565967967209469434E1L,
+ 2.2176239823732856465394E2L,
+ 3.0909872225312059774938E2L,
+ 2.1642788614495947685003E2L,
+ 6.0118660497603843919306E1L,
+};
+
+/* Coefficients for log(x) = z + z^3 P(z^2)/Q(z^2),
+ * where z = 2(x-1)/(x+1)
+ * 1/sqrt(2) <= x < sqrt(2)
+ * Theoretical peak relative error = 6.16e-22
+ */
+
+static long double R[4] = {
+ 1.9757429581415468984296E-3L,
+-7.1990767473014147232598E-1L,
+ 1.0777257190312272158094E1L,
+-3.5717684488096787370998E1L,
+};
+static long double S[4] = {
+/* 1.00000000000000000000E0L,*/
+-2.6201045551331104417768E1L,
+ 1.9361891836232102174846E2L,
+-4.2861221385716144629696E2L,
+};
+static const long double C1 = 6.9314575195312500000000E-1L;
+static const long double C2 = 1.4286068203094172321215E-6L;
+
+#define SQRTH 0.70710678118654752440L
+
+long double
+log1pl(long double xm1)
+{
+long double x, y, z;
+int e;
+
+if( isnan(xm1) )
+	return(xm1);
+if( xm1 == INFINITY )
+	return(xm1);
+if(xm1 == 0.0)
+	return(xm1);
+
+x = xm1 + 1.0L;
+
+/* Test for domain errors.  */
+if( x <= 0.0L )
+	{
+	if( x == 0.0L )
+		return( -INFINITY );
+	else
+		return( NAN );
+	}
+
+/* Separate mantissa from exponent.
+   Use frexp so that denormal numbers will be handled properly.  */
+x = frexpl( x, &e );
+
+/* logarithm using log(x) = z + z^3 P(z)/Q(z),
+   where z = 2(x-1)/x+1)  */
+if( (e > 2) || (e < -2) )
+{
+if( x < SQRTH )
+	{ /* 2( 2x-1 )/( 2x+1 ) */
+	e -= 1;
+	z = x - 0.5L;
+	y = 0.5L * z + 0.5L;
+	}	
+else
+	{ /*  2 (x-1)/(x+1)   */
+	z = x - 0.5L;
+	z -= 0.5L;
+	y = 0.5L * x  + 0.5L;
+	}
+x = z / y;
+z = x*x;
+z = x * ( z * __polevll( z, R, 3 ) / __p1evll( z, S, 3 ) );
+z = z + e * C2;
+z = z + x;
+z = z + e * C1;
+return( z );
+}
+
+
+/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
+
+if( x < SQRTH )
+	{
+	e -= 1;
+	if (e != 0)
+	  x = 2.0 * x - 1.0L;
+	else
+	  x = xm1;
+	}	
+else
+	{
+	  if (e != 0)
+	    x = x - 1.0L;
+	  else
+	    x = xm1;
+	}
+z = x*x;
+y = x * ( z * __polevll( x, P, 6 ) / __p1evll( x, Q, 6 ) );
+y = y + e * C2;
+z = y - 0.5 * z;
+z = z + x;
+z = z + e * C1;
+return( z );
+}
--- a/src/openlibm/ld80/s_log1pl.c.o
+++ b/src/openlibm/ld80/s_log1pl.c.o
--- a/src/openlibm/ld80/s_modfl.c
+++ b/src/openlibm/ld80/s_modfl.c
@@ -0,0 +1,69 @@
+/* @(#)s_modf.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.
+ * ====================================================
+ */
+
+/*
+ * modfl(long double x, long double *iptr)
+ * return fraction part of x, and return x's integral part in *iptr.
+ * Method:
+ *	Bit twiddling.
+ *
+ * Exception:
+ *	No exception.
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static const long double one = 1.0;
+
+long double
+modfl(long double x, long double *iptr)
+{
+	int32_t i0,i1,jj0;
+	u_int32_t i,se;
+	GET_LDOUBLE_WORDS(se,i0,i1,x);
+	jj0 = (se&0x7fff)-0x3fff;	/* exponent of x */
+	if(jj0<32) {			/* integer part in high x */
+	    if(jj0<0) {			/* |x|<1 */
+		SET_LDOUBLE_WORDS(*iptr,se&0x8000,0,0);	/* *iptr = +-0 */
+		return x;
+	    } else {
+		i = (0x7fffffff)>>jj0;
+		if(((i0&i)|i1)==0) {		/* x is integral */
+		    *iptr = x;
+		    SET_LDOUBLE_WORDS(x,se&0x8000,0,0);	/* return +-0 */
+		    return x;
+		} else {
+		    SET_LDOUBLE_WORDS(*iptr,se,i0&(~i),0);
+		    return x - *iptr;
+		}
+	    }
+	} else if (jj0>63) {		/* no fraction part */
+	    *iptr = x*one;
+	    /* We must handle NaNs separately.  */
+	    if (jj0 == 0x4000 && ((i0 & 0x7fffffff) | i1))
+	      return x*one;
+	    SET_LDOUBLE_WORDS(x,se&0x8000,0,0);	/* return +-0 */
+	    return x;
+	} else {			/* fraction part in low x */
+	    i = ((u_int32_t)(0x7fffffff))>>(jj0-32);
+	    if((i1&i)==0) {		/* x is integral */
+		*iptr = x;
+		SET_LDOUBLE_WORDS(x,se&0x8000,0,0);	/* return +-0 */
+		return x;
+	    } else {
+		SET_LDOUBLE_WORDS(*iptr,se,i0,i1&(~i));
+		return x - *iptr;
+	    }
+	}
+}
--- a/src/openlibm/ld80/s_nanl.c
+++ b/src/openlibm/ld80/s_nanl.c
@@ -0,0 +1,45 @@
+/*-
+ * 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 <openlibm_math.h>
+
+#include "math_private.h"
+
+OLM_DLLEXPORT long double
+nanl(const char *s)
+{
+	union {
+		union IEEEl2bits ieee;
+		uint32_t bits[3];
+	} u;
+
+	__scan_nan(u.bits, 3, s);
+	u.ieee.bits.exp = 0x7fff;
+	u.ieee.bits.manh |= 0xc0000000;	/* make it a quiet NaN */
+	return (u.ieee.e);
+}
--- a/src/openlibm/ld80/s_nanl.c.o
+++ b/src/openlibm/ld80/s_nanl.c.o
--- a/src/openlibm/ld80/s_nextafterl.c
+++ b/src/openlibm/ld80/s_nextafterl.c
@@ -0,0 +1,90 @@
+/* @(#)s_nextafter.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.
+ * ====================================================
+ */
+
+/* IEEE functions
+ *	nextafterl(x,y)
+ *	return the next machine floating-point number of x in the
+ *	direction toward y.
+ *   Special cases:
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+long double
+nextafterl(long double x, long double y)
+{
+	int32_t hx,hy,ix,iy;
+	u_int32_t lx,ly,esx,esy;
+
+	GET_LDOUBLE_WORDS(esx,hx,lx,x);
+	GET_LDOUBLE_WORDS(esy,hy,ly,y);
+	ix = esx&0x7fff;		/* |x| */
+	iy = esy&0x7fff;		/* |y| */
+
+	if (((ix==0x7fff)&&((hx&0x7fffffff|lx)!=0)) ||   /* x is nan */
+	    ((iy==0x7fff)&&((hy&0x7fffffff|ly)!=0)))     /* y is nan */
+	   return x+y;
+	if(x==y) return y;		/* x=y, return y */
+	if((ix|hx|lx)==0) {			/* x == 0 */
+	    volatile long double u;
+	    SET_LDOUBLE_WORDS(x,esy&0x8000,0,1);/* return +-minsubnormal */
+	    u = x;
+	    u = u * u;				/* raise underflow flag */
+	    return x;
+	}
+	if(esx<0x8000) {			/* x > 0 */
+	    if(ix>iy||((ix==iy) && (hx>hy||((hx==hy)&&(lx>ly))))) {
+	      /* x > y, x -= ulp */
+		if(lx==0) {
+		    if ((hx&0x7fffffff)==0) esx -= 1;
+		    hx = (hx - 1) | (hx & 0x80000000);
+		}
+		lx -= 1;
+	    } else {				/* x < y, x += ulp */
+		lx += 1;
+		if(lx==0) {
+		    hx = (hx + 1) | (hx & 0x80000000);
+		    if ((hx&0x7fffffff)==0) esx += 1;
+		}
+	    }
+	} else {				/* x < 0 */
+	    if(esy>=0||(ix>iy||((ix==iy)&&(hx>hy||((hx==hy)&&(lx>ly)))))){
+	      /* x < y, x -= ulp */
+		if(lx==0) {
+		    if ((hx&0x7fffffff)==0) esx -= 1;
+		    hx = (hx - 1) | (hx & 0x80000000);
+		}
+		lx -= 1;
+	    } else {				/* x > y, x += ulp */
+		lx += 1;
+		if(lx==0) {
+		    hx = (hx + 1) | (hx & 0x80000000);
+		    if ((hx&0x7fffffff)==0) esx += 1;
+		}
+	    }
+	}
+	esy = esx&0x7fff;
+	if(esy==0x7fff) return x+x;		/* overflow  */
+	if(esy==0) {
+	    volatile long double u = x*x;	/* underflow */
+	    if(u==x) {
+		SET_LDOUBLE_WORDS(x,esx,hx,lx);
+		return x;
+	    }
+	}
+	SET_LDOUBLE_WORDS(x,esx,hx,lx);
+	return x;
+}
+
+//__strong_alias(nexttowardl, nextafterl);
--- a/src/openlibm/ld80/s_nexttoward.c
+++ b/src/openlibm/ld80/s_nexttoward.c
@@ -0,0 +1,86 @@
+/* @(#)s_nextafter.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.
+ * ====================================================
+ */
+
+/* IEEE functions
+ *	nexttoward(x,y)
+ *	return the next machine floating-point number of x in the
+ *	direction toward y.
+ *   Special cases:
+ */
+
+#include <float.h>
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+double
+nexttoward(double x, long double y)
+{
+	int32_t hx,ix,iy;
+	u_int32_t lx,hy,ly,esy;
+
+	EXTRACT_WORDS(hx,lx,x);
+	GET_LDOUBLE_WORDS(esy,hy,ly,y);
+	ix = hx&0x7fffffff;		/* |x| */
+	iy = esy&0x7fff;		/* |y| */
+
+	if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) ||   /* x is nan */
+	   ((iy>=0x7fff)&&(hy|ly)!=0))		/* y is nan */
+	   return x+y;
+	if((long double) x==y) return y;	/* x=y, return y */
+	if((ix|lx)==0) {			/* x == 0 */
+	    volatile double u;
+	    INSERT_WORDS(x,(esy&0x8000)<<16,1); /* return +-minsub */
+	    u = x;
+	    u = u * u;				/* raise underflow flag */
+	    return x;
+	}
+	if(hx>=0) {				/* x > 0 */
+	    if (esy>=0x8000||((ix>>20)&0x7ff)>iy-0x3c00
+		|| (((ix>>20)&0x7ff)==iy-0x3c00
+		    && (((hx<<11)|(lx>>21))>(hy&0x7fffffff)
+			|| (((hx<<11)|(lx>>21))==(hy&0x7fffffff)
+			    && (lx<<11)>ly)))) {	/* x > y, x -= ulp */
+		if(lx==0) hx -= 1;
+		lx -= 1;
+	    } else {				/* x < y, x += ulp */
+		lx += 1;
+		if(lx==0) hx += 1;
+	    }
+	} else {				/* x < 0 */
+	    if (esy<0x8000||((ix>>20)&0x7ff)>iy-0x3c00
+		|| (((ix>>20)&0x7ff)==iy-0x3c00
+		    && (((hx<<11)|(lx>>21))>(hy&0x7fffffff)
+			|| (((hx<<11)|(lx>>21))==(hy&0x7fffffff)
+			    && (lx<<11)>ly))))	{/* x < y, x -= ulp */
+		if(lx==0) hx -= 1;
+		lx -= 1;
+	    } else {				/* x > y, x += ulp */
+		lx += 1;
+		if(lx==0) hx += 1;
+	    }
+	}
+	hy = hx&0x7ff00000;
+	if(hy>=0x7ff00000) {
+	  x = x+x;	/* overflow  */
+	  return x;
+	}
+	if(hy<0x00100000) {
+	    volatile double u = x*x;		/* underflow */
+	    if(u==x) {
+		INSERT_WORDS(x,hx,lx);
+		return x;
+	    }
+	}
+	INSERT_WORDS(x,hx,lx);
+	return x;
+}
--- a/src/openlibm/ld80/s_nexttowardf.c
+++ b/src/openlibm/ld80/s_nexttowardf.c
@@ -0,0 +1,67 @@
+/* @(#)s_nextafter.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 <float.h>
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+float
+nexttowardf(float x, long double y)
+{
+	int32_t hx,ix,iy;
+	u_int32_t hy,ly,esy;
+
+	GET_FLOAT_WORD(hx,x);
+	GET_LDOUBLE_WORDS(esy,hy,ly,y);
+	ix = hx&0x7fffffff;		/* |x| */
+	iy = esy&0x7fff;		/* |y| */
+
+	if((ix>0x7f800000) ||			/* x is nan */
+	   (iy>=0x7fff&&((hy|ly)!=0)))		/* y is nan */
+	   return x+y;
+	if((long double) x==y) return y;	/* x=y, return y */
+	if(ix==0) {				/* x == 0 */
+	    volatile float u;
+	    SET_FLOAT_WORD(x,((esy&0x8000)<<16)|1);/* return +-minsub*/
+	    u = x;
+	    u = u * u;				/* raise underflow flag */
+	    return x;
+	}
+	if(hx>=0) {				/* x > 0 */
+	    if(esy>=0x8000||((ix>>23)&0xff)>iy-0x3f80
+	       || (((ix>>23)&0xff)==iy-0x3f80
+		   && ((ix&0x7fffff)<<8)>(hy&0x7fffffff))) {/* x > y, x -= ulp */
+		hx -= 1;
+	    } else {				/* x < y, x += ulp */
+		hx += 1;
+	    }
+	} else {				/* x < 0 */
+	    if(esy<0x8000||((ix>>23)&0xff)>iy-0x3f80
+	       || (((ix>>23)&0xff)==iy-0x3f80
+		   && ((ix&0x7fffff)<<8)>(hy&0x7fffffff))) {/* x < y, x -= ulp */
+		hx -= 1;
+	    } else {				/* x > y, x += ulp */
+		hx += 1;
+	    }
+	}
+	hy = hx&0x7f800000;
+	if(hy>=0x7f800000) {
+	  x = x+x;	/* overflow  */
+	  return x;
+	}
+	if(hy<0x00800000) {
+	    volatile float u = x*x;		/* underflow */
+	}
+	SET_FLOAT_WORD(x,hx);
+	return x;
+}
--- a/src/openlibm/ld80/s_remquol.c
+++ b/src/openlibm/ld80/s_remquol.c
@@ -0,0 +1,166 @@
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+#include <sys/types.h>
+#include <machine/ieee.h>
+
+#include <float.h>
+#include <openlibm_math.h>
+#include <stdint.h>
+
+#include "math_private.h"
+
+#define	BIAS (LDBL_MAX_EXP - 1)
+
+/*
+ * These macros add and remove an explicit integer bit in front of the
+ * fractional mantissa, if the architecture doesn't have such a bit by
+ * default already.
+ */
+#ifdef LDBL_IMPLICIT_NBIT
+#define	LDBL_NBIT	0
+#define	SET_NBIT(hx)	((hx) | (1ULL << LDBL_MANH_SIZE))
+#define	HFRAC_BITS	EXT_FRACHBITS
+#else
+#define	LDBL_NBIT	0x80000000
+#define	SET_NBIT(hx)	(hx)
+#define	HFRAC_BITS	(EXT_FRACHBITS - 1)
+#endif
+
+#define	MANL_SHIFT	(EXT_FRACLBITS - 1)
+
+static const long double Zero[] = {0.0L, -0.0L};
+
+/*
+ * Return the IEEE remainder and set *quo to the last n bits of the
+ * quotient, rounded to the nearest integer.  We choose n=31 because
+ * we wind up computing all the integer bits of the quotient anyway as
+ * a side-effect of computing the remainder by the shift and subtract
+ * method.  In practice, this is far more bits than are needed to use
+ * remquo in reduction algorithms.
+ *
+ * Assumptions:
+ * - The low part of the mantissa fits in a manl_t exactly.
+ * - The high part of the mantissa fits in an int64_t with enough room
+ *   for an explicit integer bit in front of the fractional bits.
+ */
+long double
+remquol(long double x, long double y, int *quo)
+{
+	int64_t hx,hz;	/* We need a carry bit even if LDBL_MANH_SIZE is 32. */
+	uint32_t hy;
+	uint32_t lx,ly,lz;
+	uint32_t esx, esy;
+	int ix,iy,n,q,sx,sxy;
+
+	GET_LDOUBLE_WORDS(esx,hx,lx,x);
+	GET_LDOUBLE_WORDS(esy,hy,ly,y);
+	sx = esx & 0x8000;
+	sxy = sx ^ (esy & 0x8000);
+	esx &= 0x7fff;				/* |x| */
+	esy &= 0x7fff;				/* |y| */
+	SET_LDOUBLE_EXP(x,esx);
+	SET_LDOUBLE_EXP(y,esy);
+
+    /* purge off exception values */
+	if((esy|hy|ly)==0 ||			/* y=0 */
+	   (esx == BIAS + LDBL_MAX_EXP) ||	/* or x not finite */
+	   (esy == BIAS + LDBL_MAX_EXP &&
+	    ((hy&~LDBL_NBIT)|ly)!=0))		/* or y is NaN */
+	    return (x*y)/(x*y);
+	if(esx<=esy) {
+	    if((esx<esy) ||
+	       (hx<=hy &&
+		(hx<hy ||
+		 lx<ly))) {
+		q = 0;
+		goto fixup;			/* |x|<|y| return x or x-y */
+	    }
+	    if(hx==hy && lx==ly) {
+		*quo = 1;
+		return Zero[sx!=0];		/* |x|=|y| return x*0*/
+	    }
+	}
+
+    /* determine ix = ilogb(x) */
+	if(esx == 0) {				/* subnormal x */
+	    x *= 0x1.0p512;
+	    GET_LDOUBLE_WORDS(esx,hx,lx,x);
+	    ix = esx - (BIAS + 512);
+	} else {
+	    ix = esx - BIAS;
+	}
+
+    /* determine iy = ilogb(y) */
+	if(esy == 0) {				/* subnormal y */
+	    y *= 0x1.0p512;
+	    GET_LDOUBLE_WORDS(esy,hy,ly,y);
+	    iy = esy - (BIAS + 512);
+	} else {
+	    iy = esy - BIAS;
+	}
+
+    /* set up {hx,lx}, {hy,ly} and align y to x */
+	hx = SET_NBIT(hx);
+	lx = SET_NBIT(lx);
+
+    /* fix point fmod */
+	n = ix - iy;
+	q = 0;
+
+	while(n--) {
+	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	    if(hz<0){hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;}
+	    else {hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz; q++;}
+	    q <<= 1;
+	}
+	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	if(hz>=0) {hx=hz;lx=lz;q++;}
+
+    /* convert back to floating value and restore the sign */
+	if((hx|lx)==0) {			/* return sign(x)*0 */
+	    *quo = (sxy ? -q : q);
+	    return Zero[sx!=0];
+	}
+	while(hx<(1ULL<<HFRAC_BITS)) {	/* normalize x */
+	    hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;
+	    iy -= 1;
+	}
+	if (iy < LDBL_MIN_EXP) {
+	    esx = (iy + BIAS + 512) & 0x7fff;
+	    SET_LDOUBLE_WORDS(x,esx,hx,lx);
+	    x *= 0x1p-512;
+	    GET_LDOUBLE_WORDS(esx,hx,lx,x);
+	} else {
+	    esx = (iy + BIAS) & 0x7fff;
+	}
+	SET_LDOUBLE_WORDS(x,esx,hx,lx);
+fixup:
+	y = fabsl(y);
+	if (y < LDBL_MIN * 2) {
+	    if (x+x>y || (x+x==y && (q & 1))) {
+		q++;
+		x-=y;
+	    }
+	} else if (x>0.5*y || (x==0.5*y && (q & 1))) {
+	    q++;
+	    x-=y;
+	}
+
+	GET_LDOUBLE_EXP(esx,x);
+	esx ^= sx;
+	SET_LDOUBLE_EXP(x,esx);
+
+	q &= 0x7fffffff;
+	*quo = (sxy ? -q : q);
+	return x;
+}
--- a/src/openlibm/ld80/s_tanhl.c
+++ b/src/openlibm/ld80/s_tanhl.c
@@ -0,0 +1,79 @@
+/* @(#)s_tanh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* tanhl(x)
+ * Return the Hyperbolic Tangent of x
+ *
+ * Method :
+ *				        x    -x
+ *				       e  - e
+ *	0. tanhl(x) is defined to be -----------
+ *				        x    -x
+ *				       e  + e
+ *	1. reduce x to non-negative by tanhl(-x) = -tanhl(x).
+ *	2.  0      <= x <= 2**-55 : tanhl(x) := x*(one+x)
+ *					         -t
+ *	    2**-55 <  x <=  1     : tanhl(x) := -----; t = expm1l(-2x)
+ *					        t + 2
+ *						      2
+ *	    1      <= x <=  23.0  : tanhl(x) := 1-  ----- ; t=expm1l(2x)
+ *						    t + 2
+ *	    23.0   <  x <= INF    : tanhl(x) := 1.
+ *
+ * Special cases:
+ *	tanhl(NaN) is NaN;
+ *	only tanhl(0)=0 is exact for finite argument.
+ */
+
+#include <openlibm_math.h>
+
+#include "math_private.h"
+
+static const long double one=1.0, two=2.0, tiny = 1.0e-4900L;
+
+long double
+tanhl(long double x)
+{
+	long double t,z;
+	int32_t se;
+	u_int32_t jj0,jj1,ix;
+
+    /* High word of |x|. */
+	GET_LDOUBLE_WORDS(se,jj0,jj1,x);
+	ix = se&0x7fff;
+
+    /* x is INF or NaN */
+	if(ix==0x7fff) {
+	    /* for NaN it's not important which branch: tanhl(NaN) = NaN */
+	    if (se&0x8000) return one/x-one;	/* tanhl(-inf)= -1; */
+	    else	   return one/x+one;	/* tanhl(+inf)=+1 */
+	}
+
+    /* |x| < 23 */
+	if (ix < 0x4003 || (ix == 0x4003 && jj0 < 0xb8000000u)) {/* |x|<23 */
+	    if ((ix|jj0|jj1) == 0)
+		return x;		/* x == +- 0 */
+	    if (ix<0x3fc8)		/* |x|<2**-55 */
+		return x*(one+tiny);	/* tanh(small) = small */
+	    if (ix>=0x3fff) {	/* |x|>=1  */
+		t = expm1l(two*fabsl(x));
+		z = one - two/(t+two);
+	    } else {
+		t = expm1l(-two*fabsl(x));
+		z= -t/(t+two);
+	    }
+    /* |x| > 23, return +-1 */
+	} else {
+	    z = one - tiny;		/* raised inexact flag */
+	}
+	return (se&0x8000)? -z: z;
+}
--- a/src/openlibm/ld80/s_tanhl.c.o
+++ b/src/openlibm/ld80/s_tanhl.c.o
--- a/src/openlibm/ld80/s_truncl.c
+++ b/src/openlibm/ld80/s_truncl.c
@@ -0,0 +1,72 @@
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ * From: @(#)s_floor.c 5.1 93/09/24
+ */
+
+/*
+ * truncl(x)
+ * Return x rounded toward 0 to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to truncl(x).
+ */
+
+#include <sys/types.h>
+//#include <machine/ieee.h>
+
+#include <float.h>
+#include <openlibm_math.h>
+#include <stdint.h>
+
+#include "math_private.h"
+
+#ifdef LDBL_IMPLICIT_NBIT
+#define	MANH_SIZE	(EXT_FRACHBITS + 1)
+#else
+#define	MANH_SIZE	EXT_FRACHBITS
+#endif
+
+static const long double huge = 1.0e300;
+static const float zero[] = { 0.0, -0.0 };
+
+long double
+truncl(long double x)
+{
+	int e, es;
+	uint32_t ix0, ix1;
+
+	GET_LDOUBLE_WORDS(es,ix0,ix1,x);
+	e = (es&0x7fff) - LDBL_MAX_EXP + 1;
+
+	if (e < MANH_SIZE - 1) {
+		if (e < 0) {			/* raise inexact if x != 0 */
+			if (huge + x > 0.0)
+				return (zero[(es&0x8000)!=0]);
+		} else {
+			uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
+			if (((ix0 & m) | ix1) == 0)
+				return (x);	/* x is integral */
+			if (huge + x > 0.0) {	/* raise inexact flag */
+				ix0 &= ~m;
+				ix1 = 0;
+			}
+		}
+	} else if (e < LDBL_MANT_DIG - 1) {
+		uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
+		if ((ix1 & m) == 0)
+			return (x);	/* x is integral */
+		if (huge + x > 0.0)		/* raise inexact flag */
+			ix1 &= ~m;
+	}
+	SET_LDOUBLE_WORDS(x,es,ix0,ix1);
+	return (x);
+}