Merge from head@274682

17 files changed, 1049 insertions, 249 deletions

diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index e8458d4..afb7067 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -100,6 +100,7 @@ COMMON_SRCS+=	s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c
 # If long double != double use these; otherwise, we alias the double versions.
 COMMON_SRCS+=	e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \
 	e_coshl.c e_fmodl.c e_hypotl.c \
+	e_lgammal.c e_lgammal_r.c \
 	e_remainderl.c e_sinhl.c e_sqrtl.c \
 	invtrig.c k_cosl.c k_sinl.c k_tanl.c \
 	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \
@@ -190,7 +191,8 @@ MLINKS+=ilogb.3 ilogbf.3 ilogb.3 ilogbl.3 \
 	ilogb.3 logb.3 ilogb.3 logbf.3 ilogb.3 logbl.3
 MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3
 MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3
-MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 lgamma.3 lgammaf.3 \
+MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 \
+	lgamma.3 lgammaf.3 lgamma.3 lgammal.3 \
 	lgamma.3 tgamma.3 lgamma.3 tgammaf.3
 MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log10l.3 \
 	log.3 log1p.3 log.3 log1pf.3 log.3 log1pl.3 \
@@ -219,4 +221,8 @@ MLINKS+=tan.3 tanf.3 tan.3 tanl.3
 MLINKS+=tanh.3 tanhf.3 tanh.3 tanhl.3
 MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3
 
+.include <src.opts.mk>
+
+.include <bsd.arch.inc.mk>
+
 .include <bsd.lib.mk>
diff --git a/lib/msun/Makefile.amd64 b/lib/msun/Makefile.amd64
new file mode 100644
index 0000000..dd0f5b0
--- /dev/null
+++ b/lib/msun/Makefile.amd64
@@ -0,0 +1,6 @@
+# $FreeBSD$
+
+.if ${MK_TESTS} != "no"
+SUBDIR+=	tests
+.endif
+
diff --git a/lib/msun/Makefile.i386 b/lib/msun/Makefile.i386
new file mode 100644
index 0000000..dd0f5b0
--- /dev/null
+++ b/lib/msun/Makefile.i386
@@ -0,0 +1,6 @@
+# $FreeBSD$
+
+.if ${MK_TESTS} != "no"
+SUBDIR+=	tests
+.endif
+
diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map
index e53ca07..d5b7f46 100644
--- a/lib/msun/Symbol.map
+++ b/lib/msun/Symbol.map
@@ -269,6 +269,7 @@ FBSD_1.3 {
 	erfl;
 	expl;
 	expm1l;
+	lgammal;
 	log10l;
 	log1pl;
 	log2l;
@@ -276,7 +277,11 @@ FBSD_1.3 {
 	sinhl;
 	tanhl;
 	/* Implemented as weak aliases for imprecise versions */
-	lgammal;
 	powl;
 	tgammal;
 };
+
+/* First added in 11.0-CURRENT */
+FBSD_1.4 {
+	lgammal_r;
+};
diff --git a/lib/msun/arm/fenv.c b/lib/msun/arm/fenv.c
index 2dd1933..16f1f48 100644
--- a/lib/msun/arm/fenv.c
+++ b/lib/msun/arm/fenv.c
@@ -30,7 +30,9 @@
 #define	__fenv_static
 #include "fenv.h"
 
-#if defined(__FreeBSD_ARCH_armv6__) || (defined(__ARM_ARCH) && __ARM_ARCH >= 6)
+#include <machine/acle-compat.h>
+
+#if __ARM_ARCH >= 6
 #define FENV_ARMv6
 #endif
 
diff --git a/lib/msun/ld128/e_lgammal_r.c b/lib/msun/ld128/e_lgammal_r.c
new file mode 100644
index 0000000..53d3af1
--- /dev/null
+++ b/lib/msun/ld128/e_lgammal_r.c
@@ -0,0 +1,330 @@
+/* @(#)e_lgamma_r.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/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_lgamma_r.c for complete comments.
+ *
+ * Converted to long double by Steven G. Kargl.
+ */
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+static const volatile double vzero = 0;
+
+static const double
+zero=  0,
+half=  0.5,
+one =  1;
+
+static const long double
+pi  =  3.14159265358979323846264338327950288e+00L;
+/*
+ * Domain y in [0x1p-119, 0.28], range ~[-1.4065e-36, 1.4065e-36]:
+ * |(lgamma(2 - y) + y / 2) / y - a(y)| < 2**-119.1
+ */
+static const long double
+a0  =  7.72156649015328606065120900824024296e-02L,
+a1  =  3.22467033424113218236207583323018498e-01L,
+a2  =  6.73523010531980951332460538330282217e-02L,
+a3  =  2.05808084277845478790009252803463129e-02L,
+a4  =  7.38555102867398526627292839296001626e-03L,
+a5  =  2.89051033074152328576829509522483468e-03L,
+a6  =  1.19275391170326097618357349881842913e-03L,
+a7  =  5.09669524743042462515256340206203019e-04L,
+a8  =  2.23154758453578096143609255559576017e-04L,
+a9  =  9.94575127818397632126978731542755129e-05L,
+a10 =  4.49262367375420471287545895027098145e-05L,
+a11 =  2.05072127845117995426519671481628849e-05L,
+a12 =  9.43948816959096748454087141447939513e-06L,
+a13 =  4.37486780697359330303852050718287419e-06L,
+a14 =  2.03920783892362558276037363847651809e-06L,
+a15 =  9.55191070057967287877923073200324649e-07L,
+a16 =  4.48993286185740853170657139487620560e-07L,
+a17 =  2.13107543597620911675316728179563522e-07L,
+a18 =  9.70745379855304499867546549551023473e-08L,
+a19 =  5.61889970390290257926487734695402075e-08L,
+a20 =  6.42739653024130071866684358960960951e-09L,
+a21 =  3.34491062143649291746195612991870119e-08L,
+a22 = -1.57068547394315223934653011440641472e-08L,
+a23 =  1.30812825422415841213733487745200632e-08L;
+/*
+ * Domain x in [tc-0.24, tc+0.28], range ~[-6.3201e-37, 6.3201e-37]:
+ * |(lgamma(x) - tf) - t(x - tc)| < 2**-120.3.
+ */
+static const long double
+tc  =  1.46163214496836234126265954232572133e+00L,
+tf  = -1.21486290535849608095514557177691584e-01L,
+tt  =  1.57061739945077675484237837992951704e-36L,
+t0  = -1.99238329499314692728655623767019240e-36L,
+t1  = -6.08453430711711404116887457663281416e-35L,
+t2  =  4.83836122723810585213722380854828904e-01L,
+t3  = -1.47587722994530702030955093950668275e-01L,
+t4  =  6.46249402389127526561003464202671923e-02L,
+t5  = -3.27885410884813055008502586863748063e-02L,
+t6  =  1.79706751152103942928638276067164935e-02L,
+t7  = -1.03142230366363872751602029672767978e-02L,
+t8  =  6.10053602051788840313573150785080958e-03L,
+t9  = -3.68456960831637325470641021892968954e-03L,
+t10 =  2.25976482322181046611440855340968560e-03L,
+t11 = -1.40225144590445082933490395950664961e-03L,
+t12 =  8.78232634717681264035014878172485575e-04L,
+t13 = -5.54194952796682301220684760591403899e-04L,
+t14 =  3.51912956837848209220421213975000298e-04L,
+t15 = -2.24653443695947456542669289367055542e-04L,
+t16 =  1.44070395420840737695611929680511823e-04L,
+t17 = -9.27609865550394140067059487518862512e-05L,
+t18 =  5.99347334438437081412945428365433073e-05L,
+t19 = -3.88458388854572825603964274134801009e-05L,
+t20 =  2.52476631610328129217896436186551043e-05L,
+t21 = -1.64508584981658692556994212457518536e-05L,
+t22 =  1.07434583475987007495523340296173839e-05L,
+t23 = -7.03070407519397260929482550448878399e-06L,
+t24 =  4.60968590693753579648385629003100469e-06L,
+t25 = -3.02765473778832036018438676945512661e-06L,
+t26 =  1.99238771545503819972741288511303401e-06L,
+t27 = -1.31281299822614084861868817951788579e-06L,
+t28 =  8.60844432267399655055574642052370223e-07L,
+t29 = -5.64535486432397413273248363550536374e-07L,
+t30 =  3.99357783676275660934903139592727737e-07L,
+t31 = -2.95849029193433121795495215869311610e-07L,
+t32 =  1.37790144435073124976696250804940384e-07L;
+/*
+ * Domain y in [-0.1, 0.232], range ~[-1.4046e-37, 1.4181e-37]:
+ * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-122.8
+ */
+static const long double
+u0  = -7.72156649015328606065120900824024311e-02L,
+u1  =  4.24082772271938167430983113242482656e-01L,
+u2  =  2.96194003481457101058321977413332171e+00L,
+u3  =  6.49503267711258043997790983071543710e+00L,
+u4  =  7.40090051288150177152835698948644483e+00L,
+u5  =  4.94698036296756044610805900340723464e+00L,
+u6  =  2.00194224610796294762469550684947768e+00L,
+u7  =  4.82073087750608895996915051568834949e-01L,
+u8  =  6.46694052280506568192333848437585427e-02L,
+u9  =  4.17685526755100259316625348933108810e-03L,
+u10 =  9.06361003550314327144119307810053410e-05L,
+v1  =  5.15937098592887275994320496999951947e+00L,
+v2  =  1.14068418766251486777604403304717558e+01L,
+v3  =  1.41164839437524744055723871839748489e+01L,
+v4  =  1.07170702656179582805791063277960532e+01L,
+v5  =  5.14448694179047879915042998453632434e+00L,
+v6  =  1.55210088094585540637493826431170289e+00L,
+v7  =  2.82975732849424562719893657416365673e-01L,
+v8  =  2.86424622754753198010525786005443539e-02L,
+v9  =  1.35364253570403771005922441442688978e-03L,
+v10 =  1.91514173702398375346658943749580666e-05L,
+v11 = -3.25364686890242327944584691466034268e-08L;
+/*
+ * Domain x in (2, 3], range ~[-1.3341e-36, 1.3536e-36]:
+ * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-120.1
+ * with y = x - 2.
+ */
+static const long double
+s0  = -7.72156649015328606065120900824024297e-02L,
+s1  =  1.23221687850916448903914170805852253e-01L,
+s2  =  5.43673188699937239808255378293820020e-01L,
+s3  =  6.31998137119005233383666791176301800e-01L,
+s4  =  3.75885340179479850993811501596213763e-01L,
+s5  =  1.31572908743275052623410195011261575e-01L,
+s6  =  2.82528453299138685507186287149699749e-02L,
+s7  =  3.70262021550340817867688714880797019e-03L,
+s8  =  2.83374000312371199625774129290973648e-04L,
+s9  =  1.15091830239148290758883505582343691e-05L,
+s10 =  2.04203474281493971326506384646692446e-07L,
+s11 =  9.79544198078992058548607407635645763e-10L,
+r1  =  2.58037466655605285937112832039537492e+00L,
+r2  =  2.86289413392776399262513849911531180e+00L,
+r3  =  1.78691044735267497452847829579514367e+00L,
+r4  =  6.89400381446725342846854215600008055e-01L,
+r5  =  1.70135865462567955867134197595365343e-01L,
+r6  =  2.68794816183964420375498986152766763e-02L,
+r7  =  2.64617234244861832870088893332006679e-03L,
+r8  =  1.52881761239180800640068128681725702e-04L,
+r9  =  4.63264813762296029824851351257638558e-06L,
+r10 =  5.89461519146957343083848967333671142e-08L,
+r11 =  1.79027678176582527798327441636552968e-10L;
+/*
+ * Domain z in [8, 0x1p70], range ~[-9.8214e-35, 9.8214e-35]:
+ * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-113.0
+ */
+static const long double
+w0  =  4.18938533204672741780329736405617738e-01L,
+w1  =  8.33333333333333333333333333332852026e-02L,
+w2  = -2.77777777777777777777777727810123528e-03L,
+w3  =  7.93650793650793650791708939493907380e-04L,
+w4  = -5.95238095238095234390450004444370959e-04L,
+w5  =  8.41750841750837633887817658848845695e-04L,
+w6  = -1.91752691752396849943172337347259743e-03L,
+w7  =  6.41025640880333069429106541459015557e-03L,
+w8  = -2.95506530801732133437990433080327074e-02L,
+w9  =  1.79644237328444101596766586979576927e-01L,
+w10 = -1.39240539108367641920172649259736394e+00L,
+w11 =  1.33987701479007233325288857758641761e+01L,
+w12 = -1.56363596431084279780966590116006255e+02L,
+w13 =  2.14830978044410267201172332952040777e+03L,
+w14 = -3.28636067474227378352761516589092334e+04L,
+w15 =  5.06201257747865138432663574251462485e+05L,
+w16 = -6.79720123352023636706247599728048344e+06L,
+w17 =  6.57556601705472106989497289465949255e+07L,
+w18 = -3.26229058141181783534257632389415580e+08L;
+
+static long double
+sin_pil(long double x)
+{
+	volatile long double vz;
+	long double y,z;
+	uint64_t lx, n;
+	uint16_t hx;
+
+	y = -x;
+
+	vz = y+0x1.p112;
+	z = vz-0x1.p112;
+	if (z == y)
+	    return zero;
+
+	vz = y+0x1.p110;
+	EXTRACT_LDBL128_WORDS(hx,lx,n,vz);
+	z = vz-0x1.p110;
+	if (z > y) {
+	    z -= 0.25;
+	    n--;
+	}
+	n &= 7;
+	y = y - z + n * 0.25;
+
+	switch (n) {
+	    case 0:   y =  __kernel_sinl(pi*y,zero,0); break;
+	    case 1:
+	    case 2:   y =  __kernel_cosl(pi*(0.5-y),zero); break;
+	    case 3:
+	    case 4:   y =  __kernel_sinl(pi*(one-y),zero,0); break;
+	    case 5:
+	    case 6:   y = -__kernel_cosl(pi*(y-1.5),zero); break;
+	    default:  y =  __kernel_sinl(pi*(y-2.0),zero,0); break;
+	    }
+	return -y;
+}
+
+long double
+lgammal_r(long double x, int *signgamp)
+{
+	long double nadj,p,p1,p2,p3,q,r,t,w,y,z;
+	uint64_t llx,lx;
+	int i;
+	uint16_t hx,ix;
+
+	EXTRACT_LDBL128_WORDS(hx,lx,llx,x);
+
+    /* purge +-Inf and NaNs */
+	*signgamp = 1;
+	ix = hx&0x7fff;
+	if(ix==0x7fff) return x*x;
+
+   /* purge +-0 and tiny arguments */
+	*signgamp = 1-2*(hx>>15);
+	if(ix<0x3fff-116) {		/* |x|<2**-(p+3), return -log(|x|) */
+	    if((ix|lx|llx)==0)
+		return one/vzero;
+	    return -logl(fabsl(x));
+	}
+
+    /* purge negative integers and start evaluation for other x < 0 */
+	if(hx&0x8000) {
+	    *signgamp = 1;
+	    if(ix>=0x3fff+112) 		/* |x|>=2**(p-1), must be -integer */
+		return one/vzero;
+	    t = sin_pil(x);
+	    if(t==zero) return one/vzero;
+	    nadj = logl(pi/fabsl(t*x));
+	    if(t<zero) *signgamp = -1;
+	    x = -x;
+	}
+
+    /* purge 1 and 2 */
+	if((ix==0x3fff || ix==0x4000) && (lx|llx)==0) r = 0;
+    /* for x < 2.0 */
+	else if(ix<0x4000) {
+	    if(x<=8.9999961853027344e-01) {
+		r = -logl(x);
+		if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;}
+		else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;}
+	  	else {y = x; i=2;}
+	    } else {
+		r = 0;
+	        if(x>=1.7316312789916992e+00) {y=2-x;i=0;}
+	        else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;}
+		else {y=x-1;i=2;}
+	    }
+	    switch(i) {
+	      case 0:
+		z = y*y;
+		p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*(a12+z*(a14+z*(a16+
+		    z*(a18+z*(a20+z*a22))))))))));
+		p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*(a13+z*(a15+
+		    z*(a17+z*(a19+z*(a21+z*a23)))))))))));
+		p  = y*p1+p2;
+		r  += p-y/2; break;
+	      case 1:
+		p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+
+		    y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+
+		    y*(t17+y*(t18+y*(t19+y*(t20+y*(t21+y*(t22+y*(t23+
+		    y*(t24+y*(t25+y*(t26+y*(t27+y*(t28+y*(t29+y*(t30+
+		    y*(t31+y*t32))))))))))))))))))))))))))))));
+		r += tf + p; break;
+	      case 2:
+		p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*(u7+
+		    y*(u8+y*(u9+y*u10))))))))));
+		p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7+
+		    y*(v8+y*(v9+y*(v10+y*v11))))))))));
+		r += p1/p2-y/2;
+	    }
+	}
+    /* x < 8.0 */
+	else if(ix<0x4002) {
+	    i = x;
+	    y = x-i;
+	    p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*(s6+y*(s7+y*(s8+
+		y*(s9+y*(s10+y*s11)))))))))));
+	    q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*(r7+y*(r8+
+		y*(r9+y*(r10+y*r11))))))))));
+	    r = y/2+p/q;
+	    z = 1;	/* lgamma(1+s) = log(s) + lgamma(s) */
+	    switch(i) {
+	    case 7: z *= (y+6);		/* FALLTHRU */
+	    case 6: z *= (y+5);		/* FALLTHRU */
+	    case 5: z *= (y+4);		/* FALLTHRU */
+	    case 4: z *= (y+3);		/* FALLTHRU */
+	    case 3: z *= (y+2);		/* FALLTHRU */
+		    r += logl(z); break;
+	    }
+    /* 8.0 <= x < 2**(p+3) */
+	} else if (ix<0x3fff+116) {
+	    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+y*(w8+
+		y*(w9+y*(w10+y*(w11+y*(w12+y*(w13+y*(w14+y*(w15+y*(w16+
+		y*(w17+y*w18)))))))))))))))));
+	    r = (x-half)*(t-one)+w;
+    /* 2**(p+3) <= x <= inf */
+	} else 
+	    r =  x*(logl(x)-1);
+	if(hx&0x8000) r = nadj - r;
+	return r;
+}
diff --git a/lib/msun/ld80/e_lgammal_r.c b/lib/msun/ld80/e_lgammal_r.c
new file mode 100644
index 0000000..1e65769
--- /dev/null
+++ b/lib/msun/ld80/e_lgammal_r.c
@@ -0,0 +1,358 @@
+/* @(#)e_lgamma_r.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/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+/*
+ * See e_lgamma_r.c for complete comments.
+ *
+ * Converted to long double by Steven G. Kargl.
+ */
+
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+static const volatile double vzero = 0;
+
+static const double
+zero=  0,
+half=  0.5,
+one =  1;
+
+static const union IEEEl2bits
+piu = LD80C(0xc90fdaa22168c235, 1,  3.14159265358979323851e+00L);
+#define	pi	(piu.e)
+/*
+ * Domain y in [0x1p-70, 0.27], range ~[-4.5264e-22, 4.5264e-22]:
+ * |(lgamma(2 - y) + y / 2) / y - a(y)| < 2**-70.9
+ */
+static const union IEEEl2bits
+a0u = LD80C(0x9e233f1bed863d26, -4,  7.72156649015328606028e-02L),
+a1u = LD80C(0xa51a6625307d3249, -2,  3.22467033424113218889e-01L),
+a2u = LD80C(0x89f000d2abafda8c, -4,  6.73523010531979398946e-02L),
+a3u = LD80C(0xa8991563eca75f26, -6,  2.05808084277991211934e-02L),
+a4u = LD80C(0xf2027e10634ce6b6, -8,  7.38555102796070454026e-03L),
+a5u = LD80C(0xbd6eb76dd22187f4, -9,  2.89051035162703932972e-03L),
+a6u = LD80C(0x9c562ab05e0458ed, -10,  1.19275351624639999297e-03L),
+a7u = LD80C(0x859baed93ee48e46, -11,  5.09674593842117925320e-04L),
+a8u = LD80C(0xe9f28a4432949af2, -13,  2.23109648015769155122e-04L),
+a9u = LD80C(0xd12ad0d9b93c6bb0, -14,  9.97387167479808509830e-05L),
+a10u= LD80C(0xb7522643c78a219b, -15,  4.37071076331030136818e-05L),
+a11u= LD80C(0xca024dcdece2cb79, -16,  2.40813493372040143061e-05L),
+a12u= LD80C(0xbb90fb6968ebdbf9, -19,  2.79495621083634031729e-06L),
+a13u= LD80C(0xba1c9ffeeae07b37, -17,  1.10931287015513924136e-05L);
+#define	a0	(a0u.e)
+#define	a1	(a1u.e)
+#define	a2	(a2u.e)
+#define	a3	(a3u.e)
+#define	a4	(a4u.e)
+#define	a5	(a5u.e)
+#define	a6	(a6u.e)
+#define	a7	(a7u.e)
+#define	a8	(a8u.e)
+#define	a9	(a9u.e)
+#define	a10	(a10u.e)
+#define	a11	(a11u.e)
+#define	a12	(a12u.e)
+#define	a13	(a13u.e)
+/*
+ * Domain x in [tc-0.24, tc+0.28], range ~[-6.1205e-22, 6.1205e-22]:
+ * |(lgamma(x) - tf) -  t(x - tc)| < 2**-70.5
+ */
+static const union IEEEl2bits
+tcu  = LD80C(0xbb16c31ab5f1fb71, 0,  1.46163214496836234128e+00L),
+tfu  = LD80C(0xf8cdcde61c520e0f, -4, -1.21486290535849608093e-01L),
+ttu  = LD80C(0xd46ee54b27d4de99, -69, -2.81152980996018785880e-21L),
+t0u  = LD80C(0x80b9406556a62a6b, -68,  3.40728634996055147231e-21L),
+t1u  = LD80C(0xc7e9c6f6df3f8c39, -67, -1.05833162742737073665e-20L),
+t2u  = LD80C(0xf7b95e4771c55d51, -2,  4.83836122723810583532e-01L),
+t3u  = LD80C(0x97213c6e35e119ff, -3, -1.47587722994530691476e-01L),
+t4u  = LD80C(0x845a14a6a81dc94b, -4,  6.46249402389135358063e-02L),
+t5u  = LD80C(0x864d46fa89997796, -5, -3.27885410884846056084e-02L),
+t6u  = LD80C(0x93373cbd00297438, -6,  1.79706751150707171293e-02L),
+t7u  = LD80C(0xa8fcfca7eddc8d1d, -7, -1.03142230361450732547e-02L),
+t8u  = LD80C(0xc7e7015ff4bc45af, -8,  6.10053603296546099193e-03L),
+t9u  = LD80C(0xf178d2247adc5093, -9, -3.68456964904901200152e-03L),
+t10u = LD80C(0x94188d58f12e5e9f, -9,  2.25976420273774583089e-03L),
+t11u = LD80C(0xb7cbaef14e1406f1, -10, -1.40224943666225639823e-03L),
+t12u = LD80C(0xe63a671e6704ea4d, -11,  8.78250640744776944887e-04L),
+t13u = LD80C(0x914b6c9cae61783e, -11, -5.54255012657716808811e-04L),
+t14u = LD80C(0xb858f5bdb79276fe, -12,  3.51614951536825927370e-04L),
+t15u = LD80C(0xea73e744c34b9591, -13, -2.23591563824520112236e-04L),
+t16u = LD80C(0x99aeabb0d67ba835, -13,  1.46562869351659194136e-04L),
+t17u = LD80C(0xd7c6938325db2024, -14, -1.02889866046435680588e-04L),
+t18u = LD80C(0xe24cb1e3b0474775, -15,  5.39540265505221957652e-05L);
+#define	tc	(tcu.e)
+#define	tf	(tfu.e)
+#define	tt	(ttu.e)
+#define	t0	(t0u.e)
+#define	t1	(t1u.e)
+#define	t2	(t2u.e)
+#define	t3	(t3u.e)
+#define	t4	(t4u.e)
+#define	t5	(t5u.e)
+#define	t6	(t6u.e)
+#define	t7	(t7u.e)
+#define	t8	(t8u.e)
+#define	t9	(t9u.e)
+#define	t10	(t10u.e)
+#define	t11	(t11u.e)
+#define	t12	(t12u.e)
+#define	t13	(t13u.e)
+#define	t14	(t14u.e)
+#define	t15	(t15u.e)
+#define	t16	(t16u.e)
+#define	t17	(t17u.e)
+#define	t18	(t18u.e)
+/*
+ * Domain y in [-0.1, 0.232], range ~[-8.1938e-22, 8.3815e-22]:
+ * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-71.2
+ */
+static const union IEEEl2bits
+u0u = LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L),
+u1u = LD80C(0x98280ee45e4ddd3d, -1,  5.94361239198682739769e-01L),
+u2u = LD80C(0xe330c8ead4130733, 0,  1.77492629495841234275e+00L),
+u3u = LD80C(0xd4a213f1a002ec52, 0,  1.66119622514818078064e+00L),
+u4u = LD80C(0xa5a9ca6f5bc62163, -1,  6.47122051417476492989e-01L),
+u5u = LD80C(0xc980e49cd5b019e6, -4,  9.83903751718671509455e-02L),
+u6u = LD80C(0xff636a8bdce7025b, -9,  3.89691687802305743450e-03L),
+v1u = LD80C(0xbd109c533a19fbf5, 1,  2.95413883330948556544e+00L),
+v2u = LD80C(0xd295cbf96f31f099, 1,  3.29039286955665403176e+00L),
+v3u = LD80C(0xdab8bcfee40496cb, 0,  1.70876276441416471410e+00L),
+v4u = LD80C(0xd2f2dc3638567e9f, -2,  4.12009126299534668571e-01L),
+v5u = LD80C(0xa07d9b0851070f41, -5,  3.91822868305682491442e-02L),
+v6u = LD80C(0xe3cd8318f7adb2c4, -11,  8.68998648222144351114e-04L);
+#define	u0	(u0u.e)
+#define	u1	(u1u.e)
+#define	u2	(u2u.e)
+#define	u3	(u3u.e)
+#define	u4	(u4u.e)
+#define	u5	(u5u.e)
+#define	u6	(u6u.e)
+#define	v1	(v1u.e)
+#define	v2	(v2u.e)
+#define	v3	(v3u.e)
+#define	v4	(v4u.e)
+#define	v5	(v5u.e)
+#define	v6	(v6u.e)
+/*
+ * Domain x in (2, 3], range ~[-3.3648e-22, 3.4416e-22]:
+ * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-72.3
+ * with y = x - 2.
+ */
+static const union IEEEl2bits
+s0u = LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L),
+s1u = LD80C(0xd3ff0dcc7fa91f94, -3,  2.07027640921219389860e-01L),
+s2u = LD80C(0xb2bb62782478ef31, -2,  3.49085881391362090549e-01L),
+s3u = LD80C(0xb49f7438c4611a74, -3,  1.76389518704213357954e-01L),
+s4u = LD80C(0x9a957008fa27ecf9, -5,  3.77401710862930008071e-02L),
+s5u = LD80C(0xda9b389a6ca7a7ac, -9,  3.33566791452943399399e-03L),
+s6u = LD80C(0xbc7a2263faf59c14, -14,  8.98728786745638844395e-05L),
+r1u = LD80C(0xbf5cff5b11477d4d, 0,  1.49502555796294337722e+00L),
+r2u = LD80C(0xd9aec89de08e3da6, -1,  8.50323236984473285866e-01L),
+r3u = LD80C(0xeab7ae5057c443f9, -3,  2.29216312078225806131e-01L),
+r4u = LD80C(0xf29707d9bd2b1e37, -6,  2.96130326586640089145e-02L),
+r5u = LD80C(0xd376c2f09736c5a3, -10,  1.61334161411590662495e-03L),
+r6u = LD80C(0xc985983d0cd34e3d, -16,  2.40232770710953450636e-05L),
+r7u = LD80C(0xe5c7a4f7fc2ef13d, -25, -5.34997929289167573510e-08L);
+#define	s0	(s0u.e)
+#define	s1	(s1u.e)
+#define	s2	(s2u.e)
+#define	s3	(s3u.e)
+#define	s4	(s4u.e)
+#define	s5	(s5u.e)
+#define	s6	(s6u.e)
+#define	r1	(r1u.e)
+#define	r2	(r2u.e)
+#define	r3	(r3u.e)
+#define	r4	(r4u.e)
+#define	r5	(r5u.e)
+#define	r6	(r6u.e)
+#define	r7	(r7u.e)
+/*
+ * Domain z in [8, 0x1p70], range ~[-3.0235e-22, 3.0563e-22]:
+ * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-71.7
+ */
+static const union IEEEl2bits
+w0u = LD80C(0xd67f1c864beb4a69, -2,  4.18938533204672741776e-01L),
+w1u = LD80C(0xaaaaaaaaaaaaaaa1, -4,  8.33333333333333332678e-02L),
+w2u = LD80C(0xb60b60b60b5491c9, -9, -2.77777777777760927870e-03L),
+w3u = LD80C(0xd00d00cf58aede4c, -11,  7.93650793490637233668e-04L),
+w4u = LD80C(0x9c09bf626783d4a5, -11, -5.95238023926039051268e-04L),
+w5u = LD80C(0xdca7cadc5baa517b, -11,  8.41733700408000822962e-04L),
+w6u = LD80C(0xfb060e361e1ffd07, -10, -1.91515849570245136604e-03L),
+w7u = LD80C(0xcbd5101bb58d1f2b, -8,  6.22046743903262649294e-03L),
+w8u = LD80C(0xad27a668d32c821b, -6, -2.11370706734662081843e-02L);
+#define	w0	(w0u.e)
+#define	w1	(w1u.e)
+#define	w2	(w2u.e)
+#define	w3	(w3u.e)
+#define	w4	(w4u.e)
+#define	w5	(w5u.e)
+#define	w6	(w6u.e)
+#define	w7	(w7u.e)
+#define	w8	(w8u.e)
+
+static long double
+sin_pil(long double x)
+{
+	volatile long double vz;
+	long double y,z;
+	uint64_t n;
+	uint16_t hx;
+
+	y = -x;
+
+	vz = y+0x1p63;
+	z = vz-0x1p63;
+	if (z == y)
+	    return zero;
+
+	vz = y+0x1p61;
+	EXTRACT_LDBL80_WORDS(hx,n,vz);
+	z = vz-0x1p61;
+	if (z > y) {
+	    z -= 0.25;			/* adjust to round down */
+	    n--;
+	}
+	n &= 7;				/* octant of y mod 2 */
+	y = y - z + n * 0.25;		/* y mod 2 */
+
+	switch (n) {
+	    case 0:   y =  __kernel_sinl(pi*y,zero,0); break;
+	    case 1:
+	    case 2:   y =  __kernel_cosl(pi*(0.5-y),zero); break;
+	    case 3:
+	    case 4:   y =  __kernel_sinl(pi*(one-y),zero,0); break;
+	    case 5:
+	    case 6:   y = -__kernel_cosl(pi*(y-1.5),zero); break;
+	    default:  y =  __kernel_sinl(pi*(y-2.0),zero,0); break;
+	    }
+	return -y;
+}
+
+long double
+lgammal_r(long double x, int *signgamp)
+{
+	long double nadj,p,p1,p2,p3,q,r,t,w,y,z;
+	uint64_t lx;
+	int i;
+	uint16_t hx,ix;
+
+	EXTRACT_LDBL80_WORDS(hx,lx,x);
+
+    /* purge +-Inf and NaNs */
+	*signgamp = 1;
+	ix = hx&0x7fff;
+	if(ix==0x7fff) return x*x;
+
+	ENTERI();
+
+    /* purge +-0 and tiny arguments */
+	*signgamp = 1-2*(hx>>15);
+	if(ix<0x3fff-67) {		/* |x|<2**-(p+3), return -log(|x|) */
+	    if((ix|lx)==0)
+		RETURNI(one/vzero);
+	    RETURNI(-logl(fabsl(x)));
+	}
+
+    /* purge negative integers and start evaluation for other x < 0 */
+	if(hx&0x8000) {
+	    *signgamp = 1;
+	    if(ix>=0x3fff+63) 		/* |x|>=2**(p-1), must be -integer */
+		RETURNI(one/vzero);
+	    t = sin_pil(x);
+	    if(t==zero) RETURNI(one/vzero); /* -integer */
+	    nadj = logl(pi/fabsl(t*x));
+	    if(t<zero) *signgamp = -1;
+	    x = -x;
+	}
+
+    /* purge 1 and 2 */
+	if((ix==0x3fff || ix==0x4000) && lx==0x8000000000000000ULL) r = 0;
+    /* for x < 2.0 */
+	else if(ix<0x4000) {
+    /*
+     * XXX Supposedly, one can use the following information to replace the
+     * XXX FP rational expressions.  A similar approach is appropriate
+     * XXX for ld128, but one (may need?) needs to consider llx, too.
+     *
+     * 8.9999961853027344e-01 3ffe e666600000000000
+     * 7.3159980773925781e-01 3ffe bb4a200000000000
+     * 2.3163998126983643e-01 3ffc ed33080000000000
+     * 1.7316312789916992e+00 3fff dda6180000000000
+     * 1.2316322326660156e+00 3fff 9da6200000000000
+     */
+	    if(x<8.9999961853027344e-01) {
+		r = -logl(x);
+		if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;}
+		else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;}
+		else {y = x; i=2;}
+	    } else {
+		r = 0;
+		if(x>=1.7316312789916992e+00) {y=2-x;i=0;}
+		else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;}
+		else {y=x-1;i=2;}
+	    }
+	    switch(i) {
+	      case 0:
+		z = y*y;
+		p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*a12)))));
+		p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*a13))))));
+		p  = y*p1+p2;
+		r  += p-y/2; break;
+	      case 1:
+		p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+
+		    y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+
+		    y*(t17+y*t18))))))))))))))));
+		r += tf + p; break;
+	      case 2:
+		p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*u6))))));
+		p2 = 1+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*v6)))));
+		r += p1/p2-y/2;
+	    }
+	}
+    /* x < 8.0 */
+	else if(ix<0x4002) {
+	    i = x;
+	    y = x-i;
+	    p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
+	    q = 1+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*r7))))));
+	    r = y/2+p/q;
+	    z = 1;	/* lgamma(1+s) = log(s) + lgamma(s) */
+	    switch(i) {
+	    case 7: z *= (y+6);		/* FALLTHRU */
+	    case 6: z *= (y+5);		/* FALLTHRU */
+	    case 5: z *= (y+4);		/* FALLTHRU */
+	    case 4: z *= (y+3);		/* FALLTHRU */
+	    case 3: z *= (y+2);		/* FALLTHRU */
+		    r += logl(z); break;
+	    }
+    /* 8.0 <= x < 2**(p+3) */
+	} else if (ix<0x3fff+67) {
+	    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+y*w8)))))));
+	    r = (x-half)*(t-one)+w;
+    /* 2**(p+3) <= x <= inf */
+	} else 
+	    r =  x*(logl(x)-1);
+	if(hx&0x8000) r = nadj - r;
+	RETURNI(r);
+}
diff --git a/lib/msun/man/lgamma.3 b/lib/msun/man/lgamma.3
index ea6eb6f..d244722 100644
--- a/lib/msun/man/lgamma.3
+++ b/lib/msun/man/lgamma.3
@@ -28,7 +28,7 @@
 .\"     from: @(#)lgamma.3	6.6 (Berkeley) 12/3/92
 .\" $FreeBSD$
 .\"
-.Dd January 14, 2005
+.Dd September 12, 2014
 .Dt LGAMMA 3
 .Os
 .Sh NAME
@@ -36,6 +36,8 @@
 .Nm lgamma_r ,
 .Nm lgammaf ,
 .Nm lgammaf_r ,
+.Nm lgammal ,
+.Nm lgammal_r ,
 .Nm gamma ,
 .Nm gamma_r ,
 .Nm gammaf ,
@@ -58,6 +60,10 @@
 .Fn lgammaf "float x"
 .Ft float
 .Fn lgammaf_r "float x" "int *signgamp"
+.Ft "long double"
+.Fn lgammal "long double x"
+.Ft "long double"
+.Fn lgammal_r "long double x" "int *signgamp"
 .Ft double
 .Fn gamma "double x"
 .Ft double
@@ -66,14 +72,15 @@
 .Fn gammaf "float x"
 .Ft float
 .Fn gammaf_r "float x" "int *signgamp"
-.Ft double
+.Ft "long double"
 .Fn tgamma "double x"
 .Ft float
 .Fn tgammaf "float x"
 .Sh DESCRIPTION
-.Fn lgamma x
+.Fn lgamma x ,
+.Fn lgammaf x ,
 and
-.Fn lgammaf x
+.Fn lgammal x
 .if t \{\
 return ln\||\(*G(x)| where
 .Bd -unfilled -offset indent
@@ -87,13 +94,15 @@ The external integer
 .Fa signgam
 returns the sign of \(*G(x).
 .Pp
-.Fn lgamma_r x signgamp
+.Fn lgamma_r x signgamp ,
+.Fn lgammaf_r x signgamp ,
 and
-.Fn lgammaf_r x signgamp
+.Fn lgammal_r x signgamp
 provide the same functionality as
-.Fn lgamma x
+.Fn lgamma x , 
+.Fn lgammaf x ,
 and
-.Fn lgammaf x
+.Fn lgammal x ,
 but the caller must provide an integer to store the sign of \(*G(x).
 .Pp
 The
@@ -115,6 +124,7 @@ are deprecated aliases for
 and
 .Fn lgammaf_r ,
 respectively.
+
 .Sh IDIOSYNCRASIES
 Do not use the expression
 .Dq Li signgam\(**exp(lgamma(x))
@@ -139,14 +149,18 @@ Exponentiation of
 will lose up to 10 significant bits.
 .Sh RETURN VALUES
 .Fn gamma ,
-.Fn gamma_r ,
 .Fn gammaf ,
+.Fn gammal ,
+.Fn gamma_r ,
 .Fn gammaf_r ,
+.Fn gammal_r ,
 .Fn lgamma ,
-.Fn lgamma_r ,
 .Fn lgammaf ,
+.Fn lgammal ,
+.Fn lgamma_r ,
+.Fn lgammaf_r ,
 and
-.Fn lgammaf_r
+.Fn lgammal_r
 return appropriate values unless an argument is out of range.
 Overflow will occur for sufficiently large positive values, and
 non-positive integers.
@@ -159,6 +173,7 @@ will underflow.
 The
 .Fn lgamma ,
 .Fn lgammaf ,
+.Fn lgammal ,
 .Fn tgamma ,
 and
 .Fn tgammaf
diff --git a/lib/msun/src/e_lgamma.c b/lib/msun/src/e_lgamma.c
index 4674d9b..43f5175 100644
--- a/lib/msun/src/e_lgamma.c
+++ b/lib/msun/src/e_lgamma.c
@@ -21,6 +21,8 @@ __FBSDID("$FreeBSD$");
  * Method: call __ieee754_lgamma_r
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -31,3 +33,7 @@ __ieee754_lgamma(double x)
 {
 	return __ieee754_lgamma_r(x,&signgam);
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(lgamma, lgammal);
+#endif
diff --git a/lib/msun/src/e_lgamma_r.c b/lib/msun/src/e_lgamma_r.c
index 1cff592..be70767 100644
--- a/lib/msun/src/e_lgamma_r.c
+++ b/lib/msun/src/e_lgamma_r.c
@@ -1,4 +1,3 @@
-
 /* @(#)e_lgamma_r.c 1.3 95/01/18 */
 /*
  * ====================================================
@@ -6,22 +5,21 @@
  *
  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice 
+ * software is freely granted, provided that this notice
  * is preserved.
  * ====================================================
- *
  */
 
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
 /* __ieee754_lgamma_r(x, signgamp)
- * Reentrant version of the logarithm of the Gamma function 
- * with user provide pointer for the sign of Gamma(x). 
+ * 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 
+ * 	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,
@@ -59,20 +57,20 @@ __FBSDID("$FreeBSD$");
  *	by
  *	  			    3       5             11
  *		w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
- *	where 
+ *	where
  *		|w - f(z)| < 2**-58.74
- *		
+ *
  *   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 
+ *	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 
+ *	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
@@ -80,14 +78,17 @@ __FBSDID("$FreeBSD$");
  *		lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero
  *		lgamma(inf) = inf
  *		lgamma(-inf) = inf (bug for bug compatible with C99!?)
- *	
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
+static const volatile double vzero = 0;
+
 static const double
-two52=  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
+zero=  0.00000000000000000000e+00,
 half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
 pi  =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
@@ -154,44 +155,40 @@ w4  = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */
 w5  =  8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */
 w6  = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
 
-static const double zero=  0.00000000000000000000e+00;
-
-	static double sin_pi(double x)
+/*
+ * Compute sin(pi*x) without actually doing the pi*x multiplication.
+ * sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is
+ * the precision of x.
+ */
+static double
+sin_pi(double x)
 {
+	volatile double vz;
 	double y,z;
-	int n,ix;
+	int n;
 
-	GET_HIGH_WORD(ix,x);
-	ix &= 0x7fffffff;
+	y = -x;
 
-	if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0);
-	y = -x;		/* x is assume negative */
+	vz = y+0x1p52;			/* depend on 0 <= y < 0x1p52 */
+	z = vz-0x1p52;			/* rint(y) for the above range */
+	if (z == y)
+	    return zero;
+
+	vz = y+0x1p50;
+	GET_LOW_WORD(n,vz);		/* bits for rounded y (units 0.25) */
+	z = vz-0x1p50;			/* y rounded to a multiple of 0.25 */
+	if (z > y) {
+	    z -= 0.25;			/* adjust to round down */
+	    n--;
+	}
+	n &= 7;				/* octant of y mod 2 */
+	y = y - z + n * 0.25;		/* y mod 2 */
 
-    /*
-     * argument reduction, make sure inexact flag not raised if input
-     * is an integer
-     */
-	z = floor(y);
-	if(z!=y) {				/* inexact anyway */
-	    y  *= 0.5;
-	    y   = 2.0*(y - floor(y));		/* y = |x| mod 2.0 */
-	    n   = (int) (y*4.0);
-	} else {
-            if(ix>=0x43400000) {
-                y = zero; n = 0;                 /* y must be even */
-            } else {
-                if(ix<0x43300000) z = y+two52;	/* exact */
-		GET_LOW_WORD(n,z);
-		n &= 1;
-                y  = n;
-                n<<= 2;
-            }
-        }
 	switch (n) {
 	    case 0:   y =  __kernel_sin(pi*y,zero,0); break;
-	    case 1:   
+	    case 1:
 	    case 2:   y =  __kernel_cos(pi*(0.5-y),zero); break;
-	    case 3:  
+	    case 3:
 	    case 4:   y =  __kernel_sin(pi*(one-y),zero,0); break;
 	    case 5:
 	    case 6:   y = -__kernel_cos(pi*(y-1.5),zero); break;
@@ -204,34 +201,38 @@ static const double zero=  0.00000000000000000000e+00;
 double
 __ieee754_lgamma_r(double x, int *signgamp)
 {
-	double t,y,z,nadj,p,p1,p2,p3,q,r,w;
+	double nadj,p,p1,p2,p3,q,r,t,w,y,z;
 	int32_t hx;
-	int i,lx,ix;
+	int i,ix,lx;
 
 	EXTRACT_WORDS(hx,lx,x);
 
-    /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+    /* purge +-Inf and NaNs */
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
 	if(ix>=0x7ff00000) return x*x;
-	if((ix|lx)==0) return one/zero;
-	if(ix<0x3b900000) {	/* |x|<2**-70, return -log(|x|) */
-	    if(hx<0) {
-	        *signgamp = -1;
-	        return -__ieee754_log(-x);
-	    } else return -__ieee754_log(x);
+
+    /* purge +-0 and tiny arguments */
+	*signgamp = 1-2*((uint32_t)hx>>31);
+	if(ix<0x3c700000) {	/* |x|<2**-56, return -log(|x|) */
+	    if((ix|lx)==0)
+	        return one/vzero;
+	    return -__ieee754_log(fabs(x));
 	}
+
+    /* purge negative integers and start evaluation for other x < 0 */
 	if(hx<0) {
+	    *signgamp = 1;
 	    if(ix>=0x43300000) 	/* |x|>=2**52, must be -integer */
-		return one/zero;
+		return one/vzero;
 	    t = sin_pi(x);
-	    if(t==zero) return one/zero; /* -integer */
+	    if(t==zero) return one/vzero; /* -integer */
 	    nadj = __ieee754_log(pi/fabs(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
 	}
 
-    /* purge off 1 and 2 */
+    /* purge 1 and 2 */
 	if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
     /* for x < 2.0 */
 	else if(ix<0x40000000) {
@@ -252,7 +253,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
 		p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
 		p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
 		p  = y*p1+p2;
-		r  += (p-0.5*y); break;
+		r  += p-y/2; break;
 	      case 1:
 		z = y*y;
 		w = z*y;
@@ -260,38 +261,43 @@ __ieee754_lgamma_r(double x, int *signgamp)
 		p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
 		p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
 		p  = z*p1-(tt-w*(p2+y*p3));
-		r += (tf + p); break;
-	      case 2:	
+		r += tf + p; break;
+	      case 2:
 		p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
 		p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
-		r += (-0.5*y + p1/p2);
+		r += p1/p2-y/2;
 	    }
 	}
-	else if(ix<0x40200000) { 			/* x < 8.0 */
-	    i = (int)x;
-	    y = x-(double)i;
+    /* x < 8.0 */
+	else if(ix<0x40200000) {
+	    i = x;
+	    y = x-i;
 	    p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
 	    q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
-	    r = half*y+p/q;
+	    r = y/2+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 */
+	    case 7: z *= (y+6);		/* FALLTHRU */
+	    case 6: z *= (y+5);		/* FALLTHRU */
+	    case 5: z *= (y+4);		/* FALLTHRU */
+	    case 4: z *= (y+3);		/* FALLTHRU */
+	    case 3: z *= (y+2);		/* FALLTHRU */
 		    r += __ieee754_log(z); break;
 	    }
-    /* 8.0 <= x < 2**58 */
-	} else if (ix < 0x43900000) {
+    /* 8.0 <= x < 2**56 */
+	} else if (ix < 0x43700000) {
 	    t = __ieee754_log(x);
 	    z = one/x;
 	    y = z*z;
 	    w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
 	    r = (x-half)*(t-one)+w;
-	} else 
-    /* 2**58 <= x <= inf */
+	} else
+    /* 2**56 <= x <= inf */
 	    r =  x*(__ieee754_log(x)-one);
 	if(hx<0) r = nadj - r;
 	return r;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(lgamma_r, lgammal_r);
+#endif
diff --git a/lib/msun/src/e_lgammaf_r.c b/lib/msun/src/e_lgammaf_r.c
index e2d90ef..9084e18 100644
--- a/lib/msun/src/e_lgammaf_r.c
+++ b/lib/msun/src/e_lgammaf_r.c
@@ -1,5 +1,6 @@
 /* e_lgammaf_r.c -- float version of e_lgamma_r.c.
  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Conversion to float fixed By Steven G. Kargl.
  */
 
 /*
@@ -19,107 +20,91 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+static const volatile float vzero = 0;
+
 static const float
-two23=  8.3886080000e+06, /* 0x4b000000 */
-half=  5.0000000000e-01, /* 0x3f000000 */
-one =  1.0000000000e+00, /* 0x3f800000 */
+zero=  0,
+half=  0.5,
+one =  1,
 pi  =  3.1415927410e+00, /* 0x40490fdb */
-a0  =  7.7215664089e-02, /* 0x3d9e233f */
-a1  =  3.2246702909e-01, /* 0x3ea51a66 */
-a2  =  6.7352302372e-02, /* 0x3d89f001 */
-a3  =  2.0580807701e-02, /* 0x3ca89915 */
-a4  =  7.3855509982e-03, /* 0x3bf2027e */
-a5  =  2.8905137442e-03, /* 0x3b3d6ec6 */
-a6  =  1.1927076848e-03, /* 0x3a9c54a1 */
-a7  =  5.1006977446e-04, /* 0x3a05b634 */
-a8  =  2.2086278477e-04, /* 0x39679767 */
-a9  =  1.0801156895e-04, /* 0x38e28445 */
-a10 =  2.5214456400e-05, /* 0x37d383a2 */
-a11 =  4.4864096708e-05, /* 0x383c2c75 */
-tc  =  1.4616321325e+00, /* 0x3fbb16c3 */
-tf  = -1.2148628384e-01, /* 0xbdf8cdcd */
-/* tt = -(tail of tf) */
-tt  =  6.6971006518e-09, /* 0x31e61c52 */
-t0  =  4.8383611441e-01, /* 0x3ef7b95e */
-t1  = -1.4758771658e-01, /* 0xbe17213c */
-t2  =  6.4624942839e-02, /* 0x3d845a15 */
-t3  = -3.2788541168e-02, /* 0xbd064d47 */
-t4  =  1.7970675603e-02, /* 0x3c93373d */
-t5  = -1.0314224288e-02, /* 0xbc28fcfe */
-t6  =  6.1005386524e-03, /* 0x3bc7e707 */
-t7  = -3.6845202558e-03, /* 0xbb7177fe */
-t8  =  2.2596477065e-03, /* 0x3b141699 */
-t9  = -1.4034647029e-03, /* 0xbab7f476 */
-t10 =  8.8108185446e-04, /* 0x3a66f867 */
-t11 = -5.3859531181e-04, /* 0xba0d3085 */
-t12 =  3.1563205994e-04, /* 0x39a57b6b */
-t13 = -3.1275415677e-04, /* 0xb9a3f927 */
-t14 =  3.3552918467e-04, /* 0x39afe9f7 */
-u0  = -7.7215664089e-02, /* 0xbd9e233f */
-u1  =  6.3282704353e-01, /* 0x3f2200f4 */
-u2  =  1.4549225569e+00, /* 0x3fba3ae7 */
-u3  =  9.7771751881e-01, /* 0x3f7a4bb2 */
-u4  =  2.2896373272e-01, /* 0x3e6a7578 */
-u5  =  1.3381091878e-02, /* 0x3c5b3c5e */
-v1  =  2.4559779167e+00, /* 0x401d2ebe */
-v2  =  2.1284897327e+00, /* 0x4008392d */
-v3  =  7.6928514242e-01, /* 0x3f44efdf */
-v4  =  1.0422264785e-01, /* 0x3dd572af */
-v5  =  3.2170924824e-03, /* 0x3b52d5db */
-s0  = -7.7215664089e-02, /* 0xbd9e233f */
-s1  =  2.1498242021e-01, /* 0x3e5c245a */
-s2  =  3.2577878237e-01, /* 0x3ea6cc7a */
-s3  =  1.4635047317e-01, /* 0x3e15dce6 */
-s4  =  2.6642270386e-02, /* 0x3cda40e4 */
-s5  =  1.8402845599e-03, /* 0x3af135b4 */
-s6  =  3.1947532989e-05, /* 0x3805ff67 */
-r1  =  1.3920053244e+00, /* 0x3fb22d3b */
-r2  =  7.2193557024e-01, /* 0x3f38d0c5 */
-r3  =  1.7193385959e-01, /* 0x3e300f6e */
-r4  =  1.8645919859e-02, /* 0x3c98bf54 */
-r5  =  7.7794247773e-04, /* 0x3a4beed6 */
-r6  =  7.3266842264e-06, /* 0x36f5d7bd */
-w0  =  4.1893854737e-01, /* 0x3ed67f1d */
-w1  =  8.3333335817e-02, /* 0x3daaaaab */
-w2  = -2.7777778450e-03, /* 0xbb360b61 */
-w3  =  7.9365057172e-04, /* 0x3a500cfd */
-w4  = -5.9518753551e-04, /* 0xba1c065c */
-w5  =  8.3633989561e-04, /* 0x3a5b3dd2 */
-w6  = -1.6309292987e-03; /* 0xbad5c4e8 */
-
-static const float zero=  0.0000000000e+00;
-
-	static float sin_pif(float x)
+/*
+ * Domain y in [0x1p-27, 0.27], range ~[-3.4599e-10, 3.4590e-10]:
+ * |(lgamma(2 - y) + 0.5 * y) / y - a(y)| < 2**-31.4
+ */
+a0  =  7.72156641e-02, /* 0x3d9e233f */
+a1  =  3.22467119e-01, /* 0x3ea51a69 */
+a2  =  6.73484802e-02, /* 0x3d89ee00 */
+a3  =  2.06395667e-02, /* 0x3ca9144f */
+a4  =  6.98275631e-03, /* 0x3be4cf9b */
+a5  =  4.11768444e-03, /* 0x3b86eda4 */
+/*
+ * Domain x in [tc-0.24, tc+0.28], range ~[-5.6577e-10, 5.5677e-10]:
+ * |(lgamma(x) - tf) - t(x - tc)| < 2**-30.8.
+ */
+tc  =  1.46163213e+00, /* 0x3fbb16c3 */
+tf  = -1.21486291e-01, /* 0xbdf8cdce */
+t0  = -2.94064460e-11, /* 0xae0154b7 */
+t1  = -2.35939837e-08, /* 0xb2caabb8 */
+t2  =  4.83836412e-01, /* 0x3ef7b968 */
+t3  = -1.47586212e-01, /* 0xbe1720d7 */
+t4  =  6.46013096e-02, /* 0x3d844db1 */
+t5  = -3.28450352e-02, /* 0xbd068884 */
+t6  =  1.86483748e-02, /* 0x3c98c47a */
+t7  = -9.89206228e-03, /* 0xbc221251 */
+/*
+ * Domain y in [-0.1, 0.232], range ~[-8.4931e-10, 8.7794e-10]:
+ * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-31.2
+ */
+u0  = -7.72156641e-02, /* 0xbd9e233f */
+u1  =  7.36789703e-01, /* 0x3f3c9e40 */
+u2  =  4.95649040e-01, /* 0x3efdc5b6 */
+v1  =  1.10958421e+00, /* 0x3f8e06db */
+v2  =  2.10598111e-01, /* 0x3e57a708 */
+v3  = -1.02995494e-02, /* 0xbc28bf71 */
+/*
+ * Domain x in (2, 3], range ~[-5.5189e-11, 5.2317e-11]:
+ * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-35.0
+ * with y = x - 2.
+ */
+s0 = -7.72156641e-02, /* 0xbd9e233f */
+s1 =  2.69987404e-01, /* 0x3e8a3bca */
+s2 =  1.42851010e-01, /* 0x3e124789 */
+s3 =  1.19389519e-02, /* 0x3c439b98 */
+r1 =  6.79650068e-01, /* 0x3f2dfd8c */
+r2 =  1.16058730e-01, /* 0x3dedb033 */
+r3 =  3.75673687e-03, /* 0x3b763396 */
+/*
+ * Domain z in [8, 0x1p24], range ~[-1.2640e-09, 1.2640e-09]:
+ * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-29.6.
+ */
+w0 =  4.18938547e-01, /* 0x3ed67f1d */
+w1 =  8.33332464e-02, /* 0x3daaaa9f */
+w2 = -2.76129087e-03; /* 0xbb34f6c6 */
+
+static float
+sin_pif(float x)
 {
+	volatile float vz;
 	float y,z;
-	int n,ix;
-
-	GET_FLOAT_WORD(ix,x);
-	ix &= 0x7fffffff;
-
-	if(ix<0x3e800000) return __kernel_sindf(pi*x);
-	y = -x;		/* x is assume negative */
-
-    /*
-     * argument reduction, make sure inexact flag not raised if input
-     * is an integer
-     */
-	z = floorf(y);
-	if(z!=y) {				/* inexact anyway */
-	    y  *= (float)0.5;
-	    y   = (float)2.0*(y - floorf(y));	/* y = |x| mod 2.0 */
-	    n   = (int) (y*(float)4.0);
-	} else {
-            if(ix>=0x4b800000) {
-                y = zero; n = 0;                 /* y must be even */
-            } else {
-                if(ix<0x4b000000) z = y+two23;	/* exact */
-		GET_FLOAT_WORD(n,z);
-		n &= 1;
-                y  = n;
-                n<<= 2;
-            }
-        }
+	int n;
+
+	y = -x;
+
+	vz = y+0x1p23F;			/* depend on 0 <= y < 0x1p23 */
+	z = vz-0x1p23F;			/* rintf(y) for the above range */
+	if (z == y)
+	    return zero;
+
+	vz = y+0x1p21F;
+	GET_FLOAT_WORD(n,vz);		/* bits for rounded y (units 0.25) */
+	z = vz-0x1p21F;			/* y rounded to a multiple of 0.25 */
+	if (z > y) {
+	    z -= 0.25F;			/* adjust to round down */
+	    n--;
+	}
+	n &= 7;				/* octant of y mod 2 */
+	y = y - z + n * 0.25F;		/* y mod 2 */
+
 	switch (n) {
 	    case 0:   y =  __kernel_sindf(pi*y); break;
 	    case 1:
@@ -137,34 +122,38 @@ static const float zero=  0.0000000000e+00;
 float
 __ieee754_lgammaf_r(float x, int *signgamp)
 {
-	float t,y,z,nadj,p,p1,p2,p3,q,r,w;
+	float nadj,p,p1,p2,p3,q,r,t,w,y,z;
 	int32_t hx;
 	int i,ix;
 
 	GET_FLOAT_WORD(hx,x);
 
-    /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+    /* purge +-Inf and NaNs */
 	*signgamp = 1;
 	ix = hx&0x7fffffff;
 	if(ix>=0x7f800000) return x*x;
-	if(ix==0) return one/zero;
-	if(ix<0x35000000) {	/* |x|<2**-21, return -log(|x|) */
-	    if(hx<0) {
-	        *signgamp = -1;
-	        return -__ieee754_logf(-x);
-	    } else return -__ieee754_logf(x);
+
+    /* purge +-0 and tiny arguments */
+	*signgamp = 1-2*((uint32_t)hx>>31);
+	if(ix<0x32000000) {		/* |x|<2**-27, return -log(|x|) */
+	    if(ix==0)
+	        return one/vzero;
+	    return -__ieee754_logf(fabsf(x));
 	}
+
+    /* purge negative integers and start evaluation for other x < 0 */
 	if(hx<0) {
-	    if(ix>=0x4b000000) 	/* |x|>=2**23, must be -integer */
-		return one/zero;
+	    *signgamp = 1;
+	    if(ix>=0x4b000000) 		/* |x|>=2**23, must be -integer */
+		return one/vzero;
 	    t = sin_pif(x);
-	    if(t==zero) return one/zero; /* -integer */
+	    if(t==zero) return one/vzero; /* -integer */
 	    nadj = __ieee754_logf(pi/fabsf(t*x));
 	    if(t<zero) *signgamp = -1;
 	    x = -x;
 	}
 
-    /* purge off 1 and 2 */
+    /* purge 1 and 2 */
 	if (ix==0x3f800000||ix==0x40000000) r = 0;
     /* for x < 2.0 */
 	else if(ix<0x40000000) {
@@ -175,55 +164,51 @@ __ieee754_lgammaf_r(float x, int *signgamp)
 	  	else {y = x; i=2;}
 	    } else {
 	  	r = zero;
-	        if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */
+	        if(ix>=0x3fdda618) {y=2-x;i=0;} /* [1.7316,2] */
 	        else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */
 		else {y=x-one;i=2;}
 	    }
 	    switch(i) {
 	      case 0:
 		z = y*y;
-		p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
-		p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
+		p1 = a0+z*(a2+z*a4);
+		p2 = z*(a1+z*(a3+z*a5));
 		p  = y*p1+p2;
-		r  += (p-(float)0.5*y); break;
+		r  += p-y/2; break;
 	      case 1:
-		z = y*y;
-		w = z*y;
-		p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12)));	/* parallel comp */
-		p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
-		p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
-		p  = z*p1-(tt-w*(p2+y*p3));
-		r += (tf + p); break;
+		p = t0+y*t1+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*t7)))));
+		r += tf + p; break;
 	      case 2:
-		p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
-		p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
-		r += (-(float)0.5*y + p1/p2);
+		p1 = y*(u0+y*(u1+y*u2));
+		p2 = one+y*(v1+y*(v2+y*v3));
+		r += p1/p2-y/2;
 	    }
 	}
-	else if(ix<0x41000000) { 			/* x < 8.0 */
-	    i = (int)x;
-	    y = x-(float)i;
-	    p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
-	    q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
-	    r = half*y+p/q;
+    /* x < 8.0 */
+	else if(ix<0x41000000) {
+	    i = x;
+	    y = x-i;
+	    p = y*(s0+y*(s1+y*(s2+y*s3)));
+	    q = one+y*(r1+y*(r2+y*r3));
+	    r = y/2+p/q;
 	    z = one;	/* lgamma(1+s) = log(s) + lgamma(s) */
 	    switch(i) {
-	    case 7: z *= (y+(float)6.0);	/* FALLTHRU */
-	    case 6: z *= (y+(float)5.0);	/* FALLTHRU */
-	    case 5: z *= (y+(float)4.0);	/* FALLTHRU */
-	    case 4: z *= (y+(float)3.0);	/* FALLTHRU */
-	    case 3: z *= (y+(float)2.0);	/* FALLTHRU */
+	    case 7: z *= (y+6);		/* FALLTHRU */
+	    case 6: z *= (y+5);		/* FALLTHRU */
+	    case 5: z *= (y+4);		/* FALLTHRU */
+	    case 4: z *= (y+3);		/* FALLTHRU */
+	    case 3: z *= (y+2);		/* FALLTHRU */
 		    r += __ieee754_logf(z); break;
 	    }
-    /* 8.0 <= x < 2**58 */
-	} else if (ix < 0x5c800000) {
+    /* 8.0 <= x < 2**27 */
+	} else if (ix < 0x4d000000) {
 	    t = __ieee754_logf(x);
 	    z = one/x;
 	    y = z*z;
-	    w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
+	    w = w0+z*(w1+y*w2);
 	    r = (x-half)*(t-one)+w;
 	} else
-    /* 2**58 <= x <= inf */
+    /* 2**27 <= x <= inf */
 	    r =  x*(__ieee754_logf(x)-one);
 	if(hx<0) r = nadj - r;
 	return r;
diff --git a/lib/msun/src/e_lgammal.c b/lib/msun/src/e_lgammal.c
new file mode 100644
index 0000000..ebc2fc7
--- /dev/null
+++ b/lib/msun/src/e_lgammal.c
@@ -0,0 +1,25 @@
+/* @(#)e_lgamma.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/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include "math.h"
+#include "math_private.h"
+
+extern int signgam;
+
+long double
+lgammal(long double x)
+{
+	return lgammal_r(x,&signgam);
+}
diff --git a/lib/msun/src/imprecise.c b/lib/msun/src/imprecise.c
index 92fb2d0..08cd239 100644
--- a/lib/msun/src/imprecise.c
+++ b/lib/msun/src/imprecise.c
@@ -60,5 +60,4 @@ DECLARE_WEAK(powl);
 	long double imprecise_ ## f ## l(long double v) { return f(v); }\
 	DECLARE_WEAK(f ## l)
 
-DECLARE_IMPRECISE(lgamma);
 DECLARE_IMPRECISE(tgamma);
diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h
index 32d01da..2ff1731 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -465,6 +465,7 @@ long double	frexpl(long double value, int *); /* fundamentally !__pure2 */
 long double	hypotl(long double, long double);
 int		ilogbl(long double) __pure2;
 long double	ldexpl(long double, int);
+long double	lgammal(long double);
 long long	llrintl(long double);
 long long	llroundl(long double);
 long double	log10l(long double);
@@ -481,6 +482,7 @@ long double	nextafterl(long double, long double);
 double		nexttoward(double, long double);
 float		nexttowardf(float, long double);
 long double	nexttowardl(long double, long double);
+long double	powl(long double, long double);
 long double	remainderl(long double, long double);
 long double	remquol(long double, long double, int *);
 long double	rintl(long double);
@@ -492,31 +494,14 @@ long double	sinl(long double);
 long double	sqrtl(long double);
 long double	tanhl(long double);
 long double	tanl(long double);
+long double	tgammal(long double);
 long double	truncl(long double);
-
 #endif /* __ISO_C_VISIBLE >= 1999 */
-__END_DECLS
 
-#endif /* !_MATH_H_ */
-
-/* separate header for cmath */
-#ifndef _MATH_EXTRA_H_
-#if __ISO_C_VISIBLE >= 1999
-#if _DECLARE_C99_LDBL_MATH
-
-#define _MATH_EXTRA_H_
-
-/*
- * extra long double versions of math functions for C99 and cmath
- */
-__BEGIN_DECLS
-
-long double	lgammal(long double);
-long double	powl(long double, long double);
-long double	tgammal(long double);
+#if __BSD_VISIBLE
+long double	lgammal_r(long double, int *);
+#endif
 
 __END_DECLS
 
-#endif /* !_DECLARE_C99_LDBL_MATH */
-#endif /* __ISO_C_VISIBLE >= 1999 */
-#endif /* !_MATH_EXTRA_H_ */
+#endif /* !_MATH_H_ */
diff --git a/lib/msun/src/s_tanh.c b/lib/msun/src/s_tanh.c
index f7b71c5..6d26c69 100644
--- a/lib/msun/src/s_tanh.c
+++ b/lib/msun/src/s_tanh.c
@@ -42,7 +42,8 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
-static const double one = 1.0, two = 2.0, tiny = 1.0e-300, huge = 1.0e300;
+static const volatile double tiny = 1.0e-300;
+static const double one = 1.0, two = 2.0, huge = 1.0e300;
 
 double
 tanh(double x)
diff --git a/lib/msun/src/s_tanhf.c b/lib/msun/src/s_tanhf.c
index 04f09c6..f537be4 100644
--- a/lib/msun/src/s_tanhf.c
+++ b/lib/msun/src/s_tanhf.c
@@ -19,7 +19,9 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
-static const float one=1.0, two=2.0, tiny = 1.0e-30, huge = 1.0e30;
+static const volatile float tiny = 1.0e-30;
+static const float one=1.0, two=2.0, huge = 1.0e30;
+
 float
 tanhf(float x)
 {
diff --git a/lib/msun/tests/Makefile b/lib/msun/tests/Makefile
new file mode 100644
index 0000000..4261e48
--- /dev/null
+++ b/lib/msun/tests/Makefile
@@ -0,0 +1,63 @@
+# $FreeBSD$
+
+OBJTOP=		${.OBJDIR:H:H:H}
+SRCTOP=		${.CURDIR:H:H:H}
+TESTSRC=	${SRCTOP}/contrib/netbsd-tests/lib/libm
+
+TESTSDIR=	${TESTSBASE}/lib/msun
+
+.if ${MACHINE} == "sparc" || ${MACHINE} == "i386" \
+	|| ${MACHINE} == "amd64" || ${MACHINE_CPU} == "arm"	\
+	|| ${MACHINE} == "sparc64"
+CFLAGS+=	-DHAVE_FENV_H
+.endif
+
+.if ${MACHINE} == "amd64" || ${MACHINE} == "i386"
+CFLAGS+=	-D__HAVE_LONG_DOUBLE
+.endif
+
+NETBSD_ATF_TESTS_C=	acos_test
+NETBSD_ATF_TESTS_C+=	asin_test
+NETBSD_ATF_TESTS_C+=	atan_test
+NETBSD_ATF_TESTS_C+=	cbrt_test
+NETBSD_ATF_TESTS_C+=	ceil_test
+NETBSD_ATF_TESTS_C+=	cos_test
+NETBSD_ATF_TESTS_C+=	cosh_test
+NETBSD_ATF_TESTS_C+=	erf_test
+NETBSD_ATF_TESTS_C+=	exp_test
+NETBSD_ATF_TESTS_C+=	fmod_test
+NETBSD_ATF_TESTS_C+=	infinity_test
+NETBSD_ATF_TESTS_C+=	ldexp_test
+NETBSD_ATF_TESTS_C+=	log_test
+NETBSD_ATF_TESTS_C+=	pow_test
+NETBSD_ATF_TESTS_C+=	precision_test
+NETBSD_ATF_TESTS_C+=	round_test
+NETBSD_ATF_TESTS_C+=	scalbn_test
+NETBSD_ATF_TESTS_C+=	sin_test
+NETBSD_ATF_TESTS_C+=	sinh_test
+NETBSD_ATF_TESTS_C+=	sqrt_test
+NETBSD_ATF_TESTS_C+=	tan_test
+NETBSD_ATF_TESTS_C+=	tanh_test
+
+CSTD=		c99
+
+LDADD+=		-lm
+DPADD+=		${LIBM}
+#COPTS+=	-Wfloat-equal
+
+# Copied from lib/msun/Makefile
+.if ${MACHINE_CPUARCH} == "i386"
+ARCH_SUBDIR= i387
+.else
+ARCH_SUBDIR= ${MACHINE_CPUARCH}
+.endif
+
+.include "../${ARCH_SUBDIR}/Makefile.inc"
+
+# XXX: for some odd reason float.h doesn't tell the full story about what the
+# precision is.
+CFLAGS+=	-DLDBL_PREC=${LDBL_PREC}
+
+.include <netbsd-tests.test.mk>
+
+.include <bsd.test.mk>