* Makefile:

  . Hook e_lgammal[_r].c to the build.
  . Create man page links for lgammal[-r].3.

* Symbol.map:
  . Sort lgammal to its rightful place.
  . Add FBSD_1.4 section for the new lgamal_r symbol.

* ld128/e_lgammal_r.c:
  . 128-bit implementataion of lgammal_r().

* ld80/e_lgammal_r.c:
  . Intel 80-bit format implementation of lgammal_r().

* src/e_lgamma.c:
  . Expose lgammal as a weak reference to lgamma for platforms
    where long double is mapped to double.

* src/e_lgamma_r.c:
  . Use integer literal constants instead of real literal constants.
    Let compiler(s) do the job of conversion to the appropriate type.
  . Expose lgammal_r as a weak reference to lgamma_r for platforms
    where long double is mapped to double.

* src/e_lgammaf_r.c:
  . Fixed the Cygnus Support conversion of e_lgamma_r.c to float.
    This includes the generation of new polynomial and rational
    approximations with fewer terms.  For each approximation, include
    a comment on an estimate of the accuracy over the relevant domain.
  . Use integer literal constants instead of real literal constants.
    Let compiler(s) do the job of conversion to the appropriate type.
    This allows the removal of several explicit casts of double values
    to float.

* src/e_lgammal.c:
  . Wrapper for lgammal() about lgammal_r().

* src/imprecise.c:
  . Remove the lgamma.

* src/math.h:
  . Add a prototype for lgammal_r().

* man/lgamma.3:
  . Document the new functions.

Reviewed by:	bde

11 files changed, 837 insertions, 113 deletions

diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index 654d0f5..e827767 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -98,6 +98,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 \
@@ -188,7 +189,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 \
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/ld128/e_lgammal_r.c b/lib/msun/ld128/e_lgammal_r.c
new file mode 100644
index 0000000..a612681
--- /dev/null
+++ b/lib/msun/ld128/e_lgammal_r.c
@@ -0,0 +1,329 @@
+/* @(#)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.25L;
+
+	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;
+
+	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
+
+	if((hx & 0x7fff) == 0x7fff) {	/* erfl(nan)=nan */
+		i = (hx>>15)<<1;
+		return (1-i)+one/x;	/* erfl(+-inf)=+-1 */
+	}
+
+    /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+	*signgamp = 1;
+	if((hx & 0x7fff) == 0x7fff)	/* x is +-Inf or NaN */
+		return x*x;
+	if((hx==0||hx==0x8000)&&lx==0) return one/vzero;
+
+   /* purge off tiny and negative arguments */
+	if(fabsl(x)<0x1p-119L) {
+	    if(hx&0x8000) {
+	        *signgamp = -1;
+	        return -logl(-x);
+	    } else return -logl(x);
+	}
+	if(hx&0x8000) {
+	    if(fabsl(x)>=0x1p112)
+		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;
+	}
+
+	if(x == 1 || x ==2) r = 0;
+	else if(x<2) {
+	    if(x<=0.8999996185302734) {
+		r = -logl(x);
+		if(x>=0.7315998077392578) {y = 1-x; i= 0;}
+		else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;}
+	  	else {y = x; i=2;}
+	    } else {
+	  	r = 0;
+	        if(x>=1.7316312789916992) {y=2-x;i=0;}
+	        else if(x>=1.2316322326660156) {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 += (-y/2 + p1/p2);
+	    }
+	}
+	else if(x<8) {
+	    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;
+	    }
+	} else if (x < 0x1p119L) {
+	    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;
+	} 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..bafa9a3
--- /dev/null
+++ b/lib/msun/ld80/e_lgammal_r.c
@@ -0,0 +1,345 @@
+/* @(#)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 s_lgamma_r.c for complete comments.
+ *
+ * Converted to long double by Steven G. Kargl.
+ */
+
+#include <float.h>
+#ifdef __i386__
+#include <ieeefp.h>
+#endif
+
+#include "fpmath.h"
+#include "math.h"
+#include "math_private.h"
+
+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+0x1p63L;
+	z = vz-0x1p63L;
+	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;
+
+	EXTRACT_LDBL80_WORDS(hx,lx,x);
+
+    /* purge off +-inf, NaN, +-0 */
+	*signgamp = 1;
+	if((hx & 0x7fff) == 0x7fff)	/* x is +-Inf or NaN */
+		return x*x;
+	if((hx==0||hx==0x8000)&&lx==0) return one/vzero;
+
+	ENTERI();
+
+    /* purge off tiny and negative arguments */
+	if(fabsl(x)<0x1p-70L) {		/* |x|<2**-70, return -log(|x|) */
+	    if(hx&0x8000) {
+	        *signgamp = -1;
+	        RETURNI(-logl(-x));
+	    } else RETURNI(-logl(x));
+	}
+	if(hx&0x8000) {
+	    if(fabsl(x)>=0x1p63) 	/* |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 off 1 and 2 */
+	if(x == 1 || x == 2) r = 0;
+    /* for x < 2.0 */
+	else if(x<2) {
+	    if(x<=0.8999996185302734) {	/* lgamma(x) = lgamma(x+1)-log(x) */
+		r = - logl(x);
+		if(x>=0.7315998077392578) {y = 1-x; i= 0;}
+		else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;}
+	  	else {y = x; i=2;}
+	    } else {
+		r = 0;
+	        if(x>=1.7316312789916992) {y=2-x;i=0;}
+	        else if(x>=1.2316322326660156) {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 += (-y/2 + p1/p2);
+	    }
+	}
+	else if(x<8) {
+	    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**70 */
+	} else if (x < 0x1p70L) {
+	    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;
+	} else 
+    /* 2**70 <= x <= inf */
+	    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 7a95ea4..ce2af94 100644
--- a/lib/msun/src/e_lgamma_r.c
+++ b/lib/msun/src/e_lgamma_r.c
@@ -83,6 +83,8 @@ __FBSDID("$FreeBSD$");
  *	
  */
 
+#include <float.h>
+
 #include "math.h"
 #include "math_private.h"
 
@@ -250,7 +252,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;
@@ -273,11 +275,11 @@ __ieee754_lgamma_r(double x, int *signgamp)
 	    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 */
+	    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 */
@@ -293,3 +295,8 @@ __ieee754_lgamma_r(double x, int *signgamp)
 	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 9a7ab39..f940845 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.
  */
 
 /*
@@ -22,72 +23,63 @@ __FBSDID("$FreeBSD$");
 static const volatile float vzero = 0;
 
 static const float
-zero=  0.0000000000e+00,
-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 */
+/*
+ * 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)
@@ -168,55 +160,50 @@ __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));
+		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;
+	    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**24 */
+	} else if (ix < 0x4b800000) {
 	    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**24 <= 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 e4e5682..2ff1731 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -496,8 +496,12 @@ 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 */
+
+#if __BSD_VISIBLE
+long double	lgammal_r(long double, int *);
+#endif
+
 __END_DECLS
 
 #endif /* !_MATH_H_ */