Skip site navigation (1)Skip section navigation (2)
Date:      Mon, 15 Sep 2014 23:21:57 +0000 (UTC)
From:      Steve Kargl <kargl@FreeBSD.org>
To:        src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-head@freebsd.org
Subject:   svn commit: r271651 - in head/lib/msun: . ld128 ld80 man src
Message-ID:  <201409152321.s8FNLvcX073887@svn.freebsd.org>

next in thread | raw e-mail | index | archive | help
Author: kargl
Date: Mon Sep 15 23:21:57 2014
New Revision: 271651
URL: http://svnweb.freebsd.org/changeset/base/271651

Log:
  * 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

Added:
  head/lib/msun/ld128/e_lgammal_r.c   (contents, props changed)
  head/lib/msun/ld80/e_lgammal_r.c   (contents, props changed)
  head/lib/msun/src/e_lgammal.c   (contents, props changed)
Modified:
  head/lib/msun/Makefile
  head/lib/msun/Symbol.map
  head/lib/msun/man/lgamma.3
  head/lib/msun/src/e_lgamma.c
  head/lib/msun/src/e_lgamma_r.c
  head/lib/msun/src/e_lgammaf_r.c
  head/lib/msun/src/imprecise.c
  head/lib/msun/src/math.h

Modified: head/lib/msun/Makefile
==============================================================================
--- head/lib/msun/Makefile	Mon Sep 15 22:32:35 2014	(r271650)
+++ head/lib/msun/Makefile	Mon Sep 15 23:21:57 2014	(r271651)
@@ -98,6 +98,7 @@ COMMON_SRCS+=	s_copysignl.c s_fabsl.c s_
 # 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.
 	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 \

Modified: head/lib/msun/Symbol.map
==============================================================================
--- head/lib/msun/Symbol.map	Mon Sep 15 22:32:35 2014	(r271650)
+++ head/lib/msun/Symbol.map	Mon Sep 15 23:21:57 2014	(r271651)
@@ -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;
+};

Added: head/lib/msun/ld128/e_lgammal_r.c
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ head/lib/msun/ld128/e_lgammal_r.c	Mon Sep 15 23:21:57 2014	(r271651)
@@ -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;
+}

Added: head/lib/msun/ld80/e_lgammal_r.c
==============================================================================
--- /dev/null	00:00:00 1970	(empty, because file is newly added)
+++ head/lib/msun/ld80/e_lgammal_r.c	Mon Sep 15 23:21:57 2014	(r271651)
@@ -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);
+}

Modified: head/lib/msun/man/lgamma.3
==============================================================================
--- head/lib/msun/man/lgamma.3	Mon Sep 15 22:32:35 2014	(r271650)
+++ head/lib/msun/man/lgamma.3	Mon Sep 15 23:21:57 2014	(r271651)
@@ -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

Modified: head/lib/msun/src/e_lgamma.c
==============================================================================
--- head/lib/msun/src/e_lgamma.c	Mon Sep 15 22:32:35 2014	(r271650)
+++ head/lib/msun/src/e_lgamma.c	Mon Sep 15 23:21:57 2014	(r271651)
@@ -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

Modified: head/lib/msun/src/e_lgamma_r.c
==============================================================================
--- head/lib/msun/src/e_lgamma_r.c	Mon Sep 15 22:32:35 2014	(r271650)
+++ head/lib/msun/src/e_lgamma_r.c	Mon Sep 15 23:21:57 2014	(r271651)
@@ -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 *signga
 		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 *signga
 	    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 *signga
 	if(hx<0) r = nadj - r;
 	return r;
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(lgamma_r, lgammal_r);
+#endif
+

Modified: head/lib/msun/src/e_lgammaf_r.c
==============================================================================
--- head/lib/msun/src/e_lgammaf_r.c	Mon Sep 15 22:32:35 2014	(r271650)
+++ head/lib/msun/src/e_lgammaf_r.c	Mon Sep 15 23:21:57 2014	(r271651)
@@ -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 */

*** DIFF OUTPUT TRUNCATED AT 1000 LINES ***



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201409152321.s8FNLvcX073887>