diff options
| -rw-r--r-- | libm/Android.mk | 2 | ||||
| -rw-r--r-- | libm/include/math.h | 5 | ||||
| -rw-r--r-- | libm/upstream-freebsd/lib/msun/ld128/e_lgammal_r.c | 333 | ||||
| -rw-r--r-- | libm/upstream-freebsd/lib/msun/src/e_lgamma.c | 6 | ||||
| -rw-r--r-- | libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c | 25 | ||||
| -rw-r--r-- | libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c | 177 | ||||
| -rw-r--r-- | libm/upstream-freebsd/lib/msun/src/e_lgammal.c | 25 | ||||
| -rw-r--r-- | libm/upstream-freebsd/lib/msun/src/imprecise.c | 1 | ||||
| -rw-r--r-- | tests/math_test.cpp | 39 |
9 files changed, 511 insertions, 102 deletions
diff --git a/libm/Android.mk b/libm/Android.mk index 78223ca..60f4cb8 100644 --- a/libm/Android.mk +++ b/libm/Android.mk @@ -187,6 +187,7 @@ libm_ld_src_files = \ upstream-freebsd/lib/msun/src/e_atanhl.c \ upstream-freebsd/lib/msun/src/e_fmodl.c \ upstream-freebsd/lib/msun/src/e_hypotl.c \ + upstream-freebsd/lib/msun/src/e_lgammal.c \ upstream-freebsd/lib/msun/src/e_remainderl.c \ upstream-freebsd/lib/msun/src/e_sqrtl.c \ upstream-freebsd/lib/msun/src/s_asinhl.c \ @@ -224,6 +225,7 @@ libm_ld_src_files = \ libm_ld_src_files += \ upstream-freebsd/lib/msun/ld128/invtrig.c \ + upstream-freebsd/lib/msun/ld128/e_lgammal_r.c \ upstream-freebsd/lib/msun/ld128/k_cosl.c \ upstream-freebsd/lib/msun/ld128/k_sinl.c \ upstream-freebsd/lib/msun/ld128/k_tanl.c \ diff --git a/libm/include/math.h b/libm/include/math.h index 1fcc578..1542374 100644 --- a/libm/include/math.h +++ b/libm/include/math.h @@ -459,9 +459,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 + #if defined(__USE_GNU) void sincos(double, double*, double*); void sincosf(float, float*, float*); diff --git a/libm/upstream-freebsd/lib/msun/ld128/e_lgammal_r.c b/libm/upstream-freebsd/lib/msun/ld128/e_lgammal_r.c new file mode 100644 index 0000000..13f8f75 --- /dev/null +++ b/libm/upstream-freebsd/lib/msun/ld128/e_lgammal_r.c @@ -0,0 +1,333 @@ +/* @(#)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) { + if (hx&0x8000) + *signgamp = -1; + 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/libm/upstream-freebsd/lib/msun/src/e_lgamma.c b/libm/upstream-freebsd/lib/msun/src/e_lgamma.c index 4674d9b..43f5175 100644 --- a/libm/upstream-freebsd/lib/msun/src/e_lgamma.c +++ b/libm/upstream-freebsd/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/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c b/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c index 7a95ea4..a6fd691 100644 --- a/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c +++ b/libm/upstream-freebsd/lib/msun/src/e_lgamma_r.c @@ -83,6 +83,8 @@ __FBSDID("$FreeBSD$"); * */ +#include <float.h> + #include "math.h" #include "math_private.h" @@ -212,7 +214,11 @@ __ieee754_lgamma_r(double x, int *signgamp) *signgamp = 1; ix = hx&0x7fffffff; if(ix>=0x7ff00000) return x*x; - if((ix|lx)==0) return one/vzero; + if((ix|lx)==0) { + if(hx<0) + *signgamp = -1; + return one/vzero; + } if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */ if(hx<0) { *signgamp = -1; @@ -250,7 +256,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 +279,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 +299,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/libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c b/libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c index 9a7ab39..9d23053 100644 --- a/libm/upstream-freebsd/lib/msun/src/e_lgammaf_r.c +++ b/libm/upstream-freebsd/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) @@ -140,7 +132,11 @@ __ieee754_lgammaf_r(float x, int *signgamp) *signgamp = 1; ix = hx&0x7fffffff; if(ix>=0x7f800000) return x*x; - if(ix==0) return one/vzero; + if(ix==0) { + if(hx<0) + *signgamp = -1; + return one/vzero; + } if(ix<0x35000000) { /* |x|<2**-21, return -log(|x|) */ if(hx<0) { *signgamp = -1; @@ -168,55 +164,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/libm/upstream-freebsd/lib/msun/src/e_lgammal.c b/libm/upstream-freebsd/lib/msun/src/e_lgammal.c new file mode 100644 index 0000000..ebc2fc7 --- /dev/null +++ b/libm/upstream-freebsd/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/libm/upstream-freebsd/lib/msun/src/imprecise.c b/libm/upstream-freebsd/lib/msun/src/imprecise.c index 92fb2d0..08cd239 100644 --- a/libm/upstream-freebsd/lib/msun/src/imprecise.c +++ b/libm/upstream-freebsd/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/tests/math_test.cpp b/tests/math_test.cpp index b4604f8..538ebb1 100644 --- a/tests/math_test.cpp +++ b/tests/math_test.cpp @@ -1188,12 +1188,51 @@ TEST(math, lgamma_r) { ASSERT_EQ(1, sign); } +TEST(math, lgamma_r_17471883) { + int sign; + + sign = 0; + ASSERT_DOUBLE_EQ(HUGE_VAL, lgamma_r(0.0, &sign)); + ASSERT_EQ(1, sign); + sign = 0; + ASSERT_DOUBLE_EQ(HUGE_VAL, lgamma_r(-0.0, &sign)); + ASSERT_EQ(-1, sign); +} + TEST(math, lgammaf_r) { int sign; ASSERT_FLOAT_EQ(logf(24.0f), lgammaf_r(5.0f, &sign)); ASSERT_EQ(1, sign); } +TEST(math, lgammaf_r_17471883) { + int sign; + + sign = 0; + ASSERT_FLOAT_EQ(HUGE_VALF, lgammaf_r(0.0f, &sign)); + ASSERT_EQ(1, sign); + sign = 0; + ASSERT_FLOAT_EQ(HUGE_VALF, lgammaf_r(-0.0f, &sign)); + ASSERT_EQ(-1, sign); +} + +TEST(math, lgammal_r) { + int sign; + ASSERT_DOUBLE_EQ(log(24.0L), lgamma_r(5.0L, &sign)); + ASSERT_EQ(1, sign); +} + +TEST(math, lgammal_r_17471883) { + int sign; + + sign = 0; + ASSERT_DOUBLE_EQ(HUGE_VAL, lgammal_r(0.0L, &sign)); + ASSERT_EQ(1, sign); + sign = 0; + ASSERT_DOUBLE_EQ(HUGE_VAL, lgammal_r(-0.0L, &sign)); + ASSERT_EQ(-1, sign); +} + TEST(math, tgamma) { ASSERT_DOUBLE_EQ(24.0, tgamma(5.0)); } |