* 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
This commit is contained in:
Steve Kargl
parent
e2cc4003e2
commit
f7efd14df1
@ -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 \
|
||||
|
||||
@ -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;
|
||||
};
|
||||
|
||||
329
lib/msun/ld128/e_lgammal_r.c
Normal file
329
lib/msun/ld128/e_lgammal_r.c
Normal file
@ -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;
|
||||
}
|
||||
345
lib/msun/ld80/e_lgammal_r.c
Normal file
345
lib/msun/ld80/e_lgammal_r.c
Normal file
@ -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);
|
||||
}
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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;
|
||||
|
||||
25
lib/msun/src/e_lgammal.c
Normal file
25
lib/msun/src/e_lgammal.c
Normal file
@ -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);
|
||||
}
|
||||
@ -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);
|
||||
|
||||
@ -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_ */
|
||||
|
||||
Loading…
Reference in New Issue
Block a user