Index: head/lib/msun/ld128/e_lgammal_r.c =================================================================== --- head/lib/msun/ld128/e_lgammal_r.c (revision 272844) +++ head/lib/msun/ld128/e_lgammal_r.c (revision 272845) @@ -1,333 +1,330 @@ /* @(#)e_lgamma_r.c 1.3 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __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; + y = y - z + n * 0.25; switch (n) { case 0: y = __kernel_sinl(pi*y,zero,0); break; case 1: case 2: y = __kernel_cosl(pi*(0.5-y),zero); break; - case 3: + case 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; + uint16_t hx,ix; - EXTRACT_LDBL128_WORDS(hx, lx, llx, x); + 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 */ + /* purge +-Inf and NaNs */ *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; - } + ix = hx&0x7fff; + if(ix==0x7fff) return x*x; - /* purge off tiny and negative arguments */ - if(fabsl(x)<0x1p-119L) { - if(hx&0x8000) { - *signgamp = -1; - return -logl(-x); - } else return -logl(x); + /* purge +-0 and tiny arguments */ + *signgamp = 1-2*(hx>>15); + if(ix<0x3fff-116) { /* |x|<2**-(p+3), return -log(|x|) */ + if((ix|lx|llx)==0) + return one/vzero; + return -logl(fabsl(x)); } + + /* purge negative integers and start evaluation for other x < 0 */ if(hx&0x8000) { - if(fabsl(x)>=0x1p112) + *signgamp = 1; + if(ix>=0x3fff+112) /* |x|>=2**(p-1), must be -integer */ return one/vzero; t = sin_pil(x); if(t==zero) return one/vzero; nadj = logl(pi/fabsl(t*x)); if(t=0.7315998077392578) {y = 1-x; i= 0;} - else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;} + if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;} + else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;} else {y = x; i=2;} } else { - r = 0; - if(x>=1.7316312789916992) {y=2-x;i=0;} - else if(x>=1.2316322326660156) {y=x-tc;i=1;} + r = 0; + if(x>=1.7316312789916992e+00) {y=2-x;i=0;} + else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;} else {y=x-1;i=2;} } switch(i) { case 0: z = y*y; p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*(a12+z*(a14+z*(a16+ z*(a18+z*(a20+z*a22)))))))))); p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*(a13+z*(a15+ z*(a17+z*(a19+z*(a21+z*a23))))))))))); p = y*p1+p2; - r += (p-y/2); break; + 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; + 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); + r += p1/p2-y/2; } } - else if(x<8) { + /* x < 8.0 */ + else if(ix<0x4002) { i = x; y = x-i; p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*(s6+y*(s7+y*(s8+ y*(s9+y*(s10+y*s11))))))))))); q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*(r7+y*(r8+ y*(r9+y*(r10+y*r11)))))))))); r = y/2+p/q; z = 1; /* lgamma(1+s) = log(s) + lgamma(s) */ switch(i) { case 7: z *= (y+6); /* FALLTHRU */ case 6: z *= (y+5); /* FALLTHRU */ case 5: z *= (y+4); /* FALLTHRU */ case 4: z *= (y+3); /* FALLTHRU */ case 3: z *= (y+2); /* FALLTHRU */ r += logl(z); break; } - } else if (x < 0x1p119L) { + /* 8.0 <= x < 2**(p+3) */ + } else if (ix<0x3fff+116) { t = logl(x); z = one/x; y = z*z; w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*(w6+y*(w7+y*(w8+ y*(w9+y*(w10+y*(w11+y*(w12+y*(w13+y*(w14+y*(w15+y*(w16+ y*(w17+y*w18))))))))))))))))); r = (x-half)*(t-one)+w; + /* 2**(p+3) <= x <= inf */ } else r = x*(logl(x)-1); if(hx&0x8000) r = nadj - r; return r; } Index: head/lib/msun/ld80/e_lgammal_r.c =================================================================== --- head/lib/msun/ld80/e_lgammal_r.c (revision 272844) +++ head/lib/msun/ld80/e_lgammal_r.c (revision 272845) @@ -1,349 +1,358 @@ /* @(#)e_lgamma_r.c 1.3 95/01/18 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); /* - * See s_lgamma_r.c for complete comments. + * See e_lgamma_r.c for complete comments. * * Converted to long double by Steven G. Kargl. */ -#include #ifdef __i386__ #include #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; + vz = y+0x1p63; + z = vz-0x1p63; if (z == y) return zero; vz = y+0x1p61; EXTRACT_LDBL80_WORDS(hx,n,vz); z = vz-0x1p61; if (z > y) { z -= 0.25; /* adjust to round down */ n--; } n &= 7; /* octant of y mod 2 */ y = y - z + n * 0.25; /* y mod 2 */ switch (n) { case 0: y = __kernel_sinl(pi*y,zero,0); break; case 1: case 2: y = __kernel_cosl(pi*(0.5-y),zero); break; case 3: case 4: y = __kernel_sinl(pi*(one-y),zero,0); break; case 5: case 6: y = -__kernel_cosl(pi*(y-1.5),zero); break; default: y = __kernel_sinl(pi*(y-2.0),zero,0); break; - } + } return -y; } long double lgammal_r(long double x, int *signgamp) { long double nadj,p,p1,p2,p3,q,r,t,w,y,z; uint64_t lx; int i; - uint16_t hx; + uint16_t hx,ix; EXTRACT_LDBL80_WORDS(hx,lx,x); - /* purge off +-inf, NaN, +-0 */ + /* purge +-Inf and NaNs */ *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; - } + ix = hx&0x7fff; + if(ix==0x7fff) return x*x; 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)); + /* purge +-0 and tiny arguments */ + *signgamp = 1-2*(hx>>15); + if(ix<0x3fff-67) { /* |x|<2**-(p+3), return -log(|x|) */ + if((ix|lx)==0) + RETURNI(one/vzero); + RETURNI(-logl(fabsl(x))); } + + /* purge negative integers and start evaluation for other x < 0 */ if(hx&0x8000) { - if(fabsl(x)>=0x1p63) /* |x|>=2**(p-1), must be -integer */ + *signgamp = 1; + if(ix>=0x3fff+63) /* |x|>=2**(p-1), must be -integer */ RETURNI(one/vzero); t = sin_pil(x); if(t==zero) RETURNI(one/vzero); /* -integer */ nadj = logl(pi/fabsl(t*x)); if(t=0.7315998077392578) {y = 1-x; i= 0;} - else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;} - else {y = x; i=2;} + else if(ix<0x4000) { + /* + * XXX Supposedly, one can use the following information to replace the + * XXX FP rational expressions. A similar approach is appropriate + * XXX for ld128, but one (may need?) needs to consider llx, too. + * + * 8.9999961853027344e-01 3ffe e666600000000000 + * 7.3159980773925781e-01 3ffe bb4a200000000000 + * 2.3163998126983643e-01 3ffc ed33080000000000 + * 1.7316312789916992e+00 3fff dda6180000000000 + * 1.2316322326660156e+00 3fff 9da6200000000000 + */ + if(x<8.9999961853027344e-01) { + r = -logl(x); + if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;} + else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;} + else {y = x; i=2;} } else { r = 0; - if(x>=1.7316312789916992) {y=2-x;i=0;} - else if(x>=1.2316322326660156) {y=x-tc;i=1;} + if(x>=1.7316312789916992e+00) {y=2-x;i=0;} + else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;} else {y=x-1;i=2;} } switch(i) { case 0: z = y*y; p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*a12))))); p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*a13)))))); p = y*p1+p2; - r += (p-y/2); break; + 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; + 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); + r += p1/p2-y/2; } } - else if(x<8) { + /* x < 8.0 */ + else if(ix<0x4002) { i = x; y = x-i; p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); q = 1+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*r7)))))); r = y/2+p/q; z = 1; /* lgamma(1+s) = log(s) + lgamma(s) */ switch(i) { case 7: z *= (y+6); /* FALLTHRU */ case 6: z *= (y+5); /* FALLTHRU */ case 5: z *= (y+4); /* FALLTHRU */ case 4: z *= (y+3); /* FALLTHRU */ case 3: z *= (y+2); /* FALLTHRU */ r += logl(z); break; } - /* 8.0 <= x < 2**70 */ - } else if (x < 0x1p70L) { + /* 8.0 <= x < 2**(p+3) */ + } else if (ix<0x3fff+67) { t = logl(x); z = one/x; y = z*z; w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*(w6+y*(w7+y*w8))))))); r = (x-half)*(t-one)+w; + /* 2**(p+3) <= x <= inf */ } else - /* 2**70 <= x <= inf */ r = x*(logl(x)-1); if(hx&0x8000) r = nadj - r; RETURNI(r); } Index: head/lib/msun/src/e_lgamma_r.c =================================================================== --- head/lib/msun/src/e_lgamma_r.c (revision 272844) +++ head/lib/msun/src/e_lgamma_r.c (revision 272845) @@ -1,302 +1,303 @@ /* @(#)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 __FBSDID("$FreeBSD$"); /* __ieee754_lgamma_r(x, signgamp) * Reentrant version of the logarithm of the Gamma function * with user provide pointer for the sign of Gamma(x). * * Method: * 1. Argument Reduction for 0 < x <= 8 * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may * reduce x to a number in [1.5,2.5] by * lgamma(1+s) = log(s) + lgamma(s) * for example, * lgamma(7.3) = log(6.3) + lgamma(6.3) * = log(6.3*5.3) + lgamma(5.3) * = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3) * 2. Polynomial approximation of lgamma around its * minimun ymin=1.461632144968362245 to maintain monotonicity. * On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use * Let z = x-ymin; * lgamma(x) = -1.214862905358496078218 + z^2*poly(z) * where * poly(z) is a 14 degree polynomial. * 2. Rational approximation in the primary interval [2,3] * We use the following approximation: * s = x-2.0; * lgamma(x) = 0.5*s + s*P(s)/Q(s) * with accuracy * |P/Q - (lgamma(x)-0.5s)| < 2**-61.71 * Our algorithms are based on the following observation * * zeta(2)-1 2 zeta(3)-1 3 * lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ... * 2 3 * * where Euler = 0.5771... is the Euler constant, which is very * close to 0.5. * * 3. For x>=8, we have * lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+.... * (better formula: * lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...) * Let z = 1/x, then we approximation * f(z) = lgamma(x) - (x-0.5)(log(x)-1) * by * 3 5 11 * w = w0 + w1*z + w2*z + w3*z + ... + w6*z * where * |w - f(z)| < 2**-58.74 * * 4. For negative x, since (G is gamma function) * -x*G(-x)*G(x) = pi/sin(pi*x), * we have * G(x) = pi/(sin(pi*x)*(-x)*G(-x)) * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 * Hence, for x<0, signgam = sign(sin(pi*x)) and * lgamma(x) = log(|Gamma(x)|) * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x); * Note: one should avoid compute pi*(-x) directly in the * computation of sin(pi*(-x)). * * 5. Special Cases * lgamma(2+s) ~ s*(1-Euler) for tiny s * lgamma(1) = lgamma(2) = 0 * lgamma(x) ~ -log(|x|) for tiny x * lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero * lgamma(inf) = inf * lgamma(-inf) = inf (bug for bug compatible with C99!?) */ #include #include "math.h" #include "math_private.h" static const volatile double vzero = 0; static const double zero= 0.00000000000000000000e+00, half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */ a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */ a2 = 6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */ a3 = 2.05808084325167332806e-02, /* 0x3F951322, 0xAC92547B */ a4 = 7.38555086081402883957e-03, /* 0x3F7E404F, 0xB68FEFE8 */ a5 = 2.89051383673415629091e-03, /* 0x3F67ADD8, 0xCCB7926B */ a6 = 1.19270763183362067845e-03, /* 0x3F538A94, 0x116F3F5D */ a7 = 5.10069792153511336608e-04, /* 0x3F40B6C6, 0x89B99C00 */ a8 = 2.20862790713908385557e-04, /* 0x3F2CF2EC, 0xED10E54D */ a9 = 1.08011567247583939954e-04, /* 0x3F1C5088, 0x987DFB07 */ a10 = 2.52144565451257326939e-05, /* 0x3EFA7074, 0x428CFA52 */ a11 = 4.48640949618915160150e-05, /* 0x3F07858E, 0x90A45837 */ tc = 1.46163214496836224576e+00, /* 0x3FF762D8, 0x6356BE3F */ tf = -1.21486290535849611461e-01, /* 0xBFBF19B9, 0xBCC38A42 */ /* tt = -(tail of tf) */ tt = -3.63867699703950536541e-18, /* 0xBC50C7CA, 0xA48A971F */ t0 = 4.83836122723810047042e-01, /* 0x3FDEF72B, 0xC8EE38A2 */ t1 = -1.47587722994593911752e-01, /* 0xBFC2E427, 0x8DC6C509 */ t2 = 6.46249402391333854778e-02, /* 0x3FB08B42, 0x94D5419B */ t3 = -3.27885410759859649565e-02, /* 0xBFA0C9A8, 0xDF35B713 */ t4 = 1.79706750811820387126e-02, /* 0x3F9266E7, 0x970AF9EC */ t5 = -1.03142241298341437450e-02, /* 0xBF851F9F, 0xBA91EC6A */ t6 = 6.10053870246291332635e-03, /* 0x3F78FCE0, 0xE370E344 */ t7 = -3.68452016781138256760e-03, /* 0xBF6E2EFF, 0xB3E914D7 */ t8 = 2.25964780900612472250e-03, /* 0x3F6282D3, 0x2E15C915 */ t9 = -1.40346469989232843813e-03, /* 0xBF56FE8E, 0xBF2D1AF1 */ t10 = 8.81081882437654011382e-04, /* 0x3F4CDF0C, 0xEF61A8E9 */ t11 = -5.38595305356740546715e-04, /* 0xBF41A610, 0x9C73E0EC */ t12 = 3.15632070903625950361e-04, /* 0x3F34AF6D, 0x6C0EBBF7 */ t13 = -3.12754168375120860518e-04, /* 0xBF347F24, 0xECC38C38 */ t14 = 3.35529192635519073543e-04, /* 0x3F35FD3E, 0xE8C2D3F4 */ u0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */ u1 = 6.32827064025093366517e-01, /* 0x3FE4401E, 0x8B005DFF */ u2 = 1.45492250137234768737e+00, /* 0x3FF7475C, 0xD119BD6F */ u3 = 9.77717527963372745603e-01, /* 0x3FEF4976, 0x44EA8450 */ u4 = 2.28963728064692451092e-01, /* 0x3FCD4EAE, 0xF6010924 */ u5 = 1.33810918536787660377e-02, /* 0x3F8B678B, 0xBF2BAB09 */ v1 = 2.45597793713041134822e+00, /* 0x4003A5D7, 0xC2BD619C */ v2 = 2.12848976379893395361e+00, /* 0x40010725, 0xA42B18F5 */ v3 = 7.69285150456672783825e-01, /* 0x3FE89DFB, 0xE45050AF */ v4 = 1.04222645593369134254e-01, /* 0x3FBAAE55, 0xD6537C88 */ v5 = 3.21709242282423911810e-03, /* 0x3F6A5ABB, 0x57D0CF61 */ s0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */ s1 = 2.14982415960608852501e-01, /* 0x3FCB848B, 0x36E20878 */ s2 = 3.25778796408930981787e-01, /* 0x3FD4D98F, 0x4F139F59 */ s3 = 1.46350472652464452805e-01, /* 0x3FC2BB9C, 0xBEE5F2F7 */ s4 = 2.66422703033638609560e-02, /* 0x3F9B481C, 0x7E939961 */ s5 = 1.84028451407337715652e-03, /* 0x3F5E26B6, 0x7368F239 */ s6 = 3.19475326584100867617e-05, /* 0x3F00BFEC, 0xDD17E945 */ r1 = 1.39200533467621045958e+00, /* 0x3FF645A7, 0x62C4AB74 */ r2 = 7.21935547567138069525e-01, /* 0x3FE71A18, 0x93D3DCDC */ r3 = 1.71933865632803078993e-01, /* 0x3FC601ED, 0xCCFBDF27 */ r4 = 1.86459191715652901344e-02, /* 0x3F9317EA, 0x742ED475 */ r5 = 7.77942496381893596434e-04, /* 0x3F497DDA, 0xCA41A95B */ r6 = 7.32668430744625636189e-06, /* 0x3EDEBAF7, 0xA5B38140 */ w0 = 4.18938533204672725052e-01, /* 0x3FDACFE3, 0x90C97D69 */ w1 = 8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */ w2 = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */ w3 = 7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */ w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */ w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */ w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */ /* * Compute sin(pi*x) without actually doing the pi*x multiplication. * sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is * the precision of x. */ static double sin_pi(double x) { volatile double vz; double y,z; int n; y = -x; vz = y+0x1p52; /* depend on 0 <= y < 0x1p52 */ z = vz-0x1p52; /* rint(y) for the above range */ if (z == y) return zero; vz = y+0x1p50; GET_LOW_WORD(n,vz); /* bits for rounded y (units 0.25) */ z = vz-0x1p50; /* y rounded to a multiple of 0.25 */ if (z > y) { z -= 0.25; /* adjust to round down */ n--; } n &= 7; /* octant of y mod 2 */ y = y - z + n * 0.25; /* y mod 2 */ switch (n) { case 0: y = __kernel_sin(pi*y,zero,0); break; case 1: case 2: y = __kernel_cos(pi*(0.5-y),zero); break; case 3: case 4: y = __kernel_sin(pi*(one-y),zero,0); break; case 5: case 6: y = -__kernel_cos(pi*(y-1.5),zero); break; default: y = __kernel_sin(pi*(y-2.0),zero,0); break; } return -y; } double __ieee754_lgamma_r(double x, int *signgamp) { - double t,y,z,nadj,p,p1,p2,p3,q,r,w; + double nadj,p,p1,p2,p3,q,r,t,w,y,z; int32_t hx; int i,ix,lx; EXTRACT_WORDS(hx,lx,x); - /* purge off +-inf, NaN, +-0, tiny and negative arguments */ + /* purge +-Inf and NaNs */ *signgamp = 1; ix = hx&0x7fffffff; if(ix>=0x7ff00000) return x*x; - if((ix|lx)==0) { - if(hx<0) - *signgamp = -1; - return one/vzero; + + /* purge +-0 and tiny arguments */ + *signgamp = 1-2*((uint32_t)hx>>31); + if(ix<0x3c700000) { /* |x|<2**-56, return -log(|x|) */ + if((ix|lx)==0) + return one/vzero; + return -__ieee754_log(fabs(x)); } - if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */ - if(hx<0) { - *signgamp = -1; - return -__ieee754_log(-x); - } else return -__ieee754_log(x); - } + + /* purge negative integers and start evaluation for other x < 0 */ if(hx<0) { + *signgamp = 1; if(ix>=0x43300000) /* |x|>=2**52, must be -integer */ return one/vzero; t = sin_pi(x); if(t==zero) return one/vzero; /* -integer */ nadj = __ieee754_log(pi/fabs(t*x)); if(t=0x3FE76944) {y = one-x; i= 0;} else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;} else {y = x; i=2;} } else { r = zero; if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */ else if(ix>=0x3FF3B4C4) {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))))); p = y*p1+p2; - r += (p-y/2); break; + r += p-y/2; break; case 1: z = y*y; w = z*y; p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); p = z*p1-(tt-w*(p2+y*p3)); - r += (tf + p); break; + r += tf + p; break; case 2: p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); - r += (-0.5*y + p1/p2); + r += p1/p2-y/2; } } - else if(ix<0x40200000) { /* x < 8.0 */ - i = (int)x; - y = x-(double)i; + /* x < 8.0 */ + else if(ix<0x40200000) { + i = x; + y = x-i; p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); - r = half*y+p/q; + r = y/2+p/q; z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ switch(i) { case 7: z *= (y+6); /* FALLTHRU */ case 6: z *= (y+5); /* FALLTHRU */ case 5: z *= (y+4); /* FALLTHRU */ case 4: z *= (y+3); /* FALLTHRU */ case 3: z *= (y+2); /* FALLTHRU */ r += __ieee754_log(z); break; } - /* 8.0 <= x < 2**58 */ - } else if (ix < 0x43900000) { + /* 8.0 <= x < 2**56 */ + } else if (ix < 0x43700000) { t = __ieee754_log(x); z = one/x; y = z*z; w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); r = (x-half)*(t-one)+w; } else - /* 2**58 <= x <= inf */ + /* 2**56 <= x <= inf */ r = x*(__ieee754_log(x)-one); if(hx<0) r = nadj - r; return r; } #if (LDBL_MANT_DIG == 53) __weak_reference(lgamma_r, lgammal_r); #endif Index: head/lib/msun/src/e_lgammaf_r.c =================================================================== --- head/lib/msun/src/e_lgammaf_r.c (revision 272844) +++ head/lib/msun/src/e_lgammaf_r.c (revision 272845) @@ -1,214 +1,215 @@ /* 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. */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ #include __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" static const volatile float vzero = 0; static const float zero= 0, half= 0.5, one = 1, pi = 3.1415927410e+00, /* 0x40490fdb */ /* * Domain y in [0x1p-27, 0.27], range ~[-3.4599e-10, 3.4590e-10]: * |(lgamma(2 - y) + 0.5 * y) / y - a(y)| < 2**-31.4 */ a0 = 7.72156641e-02, /* 0x3d9e233f */ a1 = 3.22467119e-01, /* 0x3ea51a69 */ a2 = 6.73484802e-02, /* 0x3d89ee00 */ a3 = 2.06395667e-02, /* 0x3ca9144f */ a4 = 6.98275631e-03, /* 0x3be4cf9b */ a5 = 4.11768444e-03, /* 0x3b86eda4 */ /* * Domain x in [tc-0.24, tc+0.28], range ~[-5.6577e-10, 5.5677e-10]: * |(lgamma(x) - tf) - t(x - tc)| < 2**-30.8. */ tc = 1.46163213e+00, /* 0x3fbb16c3 */ tf = -1.21486291e-01, /* 0xbdf8cdce */ t0 = -2.94064460e-11, /* 0xae0154b7 */ t1 = -2.35939837e-08, /* 0xb2caabb8 */ t2 = 4.83836412e-01, /* 0x3ef7b968 */ t3 = -1.47586212e-01, /* 0xbe1720d7 */ t4 = 6.46013096e-02, /* 0x3d844db1 */ t5 = -3.28450352e-02, /* 0xbd068884 */ t6 = 1.86483748e-02, /* 0x3c98c47a */ t7 = -9.89206228e-03, /* 0xbc221251 */ /* * Domain y in [-0.1, 0.232], range ~[-8.4931e-10, 8.7794e-10]: * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-31.2 */ u0 = -7.72156641e-02, /* 0xbd9e233f */ u1 = 7.36789703e-01, /* 0x3f3c9e40 */ u2 = 4.95649040e-01, /* 0x3efdc5b6 */ v1 = 1.10958421e+00, /* 0x3f8e06db */ v2 = 2.10598111e-01, /* 0x3e57a708 */ v3 = -1.02995494e-02, /* 0xbc28bf71 */ /* * Domain x in (2, 3], range ~[-5.5189e-11, 5.2317e-11]: * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-35.0 * with y = x - 2. */ s0 = -7.72156641e-02, /* 0xbd9e233f */ s1 = 2.69987404e-01, /* 0x3e8a3bca */ s2 = 1.42851010e-01, /* 0x3e124789 */ s3 = 1.19389519e-02, /* 0x3c439b98 */ r1 = 6.79650068e-01, /* 0x3f2dfd8c */ r2 = 1.16058730e-01, /* 0x3dedb033 */ r3 = 3.75673687e-03, /* 0x3b763396 */ /* * Domain z in [8, 0x1p24], range ~[-1.2640e-09, 1.2640e-09]: * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-29.6. */ w0 = 4.18938547e-01, /* 0x3ed67f1d */ w1 = 8.33332464e-02, /* 0x3daaaa9f */ w2 = -2.76129087e-03; /* 0xbb34f6c6 */ static float sin_pif(float x) { volatile float vz; float y,z; int n; y = -x; vz = y+0x1p23F; /* depend on 0 <= y < 0x1p23 */ z = vz-0x1p23F; /* rintf(y) for the above range */ if (z == y) return zero; vz = y+0x1p21F; GET_FLOAT_WORD(n,vz); /* bits for rounded y (units 0.25) */ z = vz-0x1p21F; /* y rounded to a multiple of 0.25 */ if (z > y) { z -= 0.25F; /* adjust to round down */ n--; } n &= 7; /* octant of y mod 2 */ y = y - z + n * 0.25F; /* y mod 2 */ switch (n) { case 0: y = __kernel_sindf(pi*y); break; case 1: case 2: y = __kernel_cosdf(pi*((float)0.5-y)); break; case 3: case 4: y = __kernel_sindf(pi*(one-y)); break; case 5: case 6: y = -__kernel_cosdf(pi*(y-(float)1.5)); break; default: y = __kernel_sindf(pi*(y-(float)2.0)); break; } return -y; } float __ieee754_lgammaf_r(float x, int *signgamp) { - float t,y,z,nadj,p,p1,p2,p3,q,r,w; + float nadj,p,p1,p2,p3,q,r,t,w,y,z; int32_t hx; int i,ix; GET_FLOAT_WORD(hx,x); - /* purge off +-inf, NaN, +-0, tiny and negative arguments */ + /* purge +-Inf and NaNs */ *signgamp = 1; ix = hx&0x7fffffff; if(ix>=0x7f800000) return x*x; - if(ix==0) { - if(hx<0) - *signgamp = -1; - return one/vzero; + + /* purge +-0 and tiny arguments */ + *signgamp = 1-2*((uint32_t)hx>>31); + if(ix<0x32000000) { /* |x|<2**-27, return -log(|x|) */ + if(ix==0) + return one/vzero; + return -__ieee754_logf(fabsf(x)); } - if(ix<0x35000000) { /* |x|<2**-21, return -log(|x|) */ - if(hx<0) { - *signgamp = -1; - return -__ieee754_logf(-x); - } else return -__ieee754_logf(x); - } + + /* purge negative integers and start evaluation for other x < 0 */ if(hx<0) { - if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ + *signgamp = 1; + if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ return one/vzero; t = sin_pif(x); if(t==zero) return one/vzero; /* -integer */ nadj = __ieee754_logf(pi/fabsf(t*x)); if(t=0x3f3b4a20) {y = one-x; i= 0;} else if(ix>=0x3e6d3308) {y= x-(tc-one); i=1;} else {y = x; i=2;} } else { r = zero; 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); p2 = z*(a1+z*(a3+z*a5)); p = y*p1+p2; - r += (p-y/2); break; + r += p-y/2; break; case 1: p = t0+y*t1+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*t7))))); - r += (tf + p); break; + r += tf + p; break; case 2: p1 = y*(u0+y*(u1+y*u2)); p2 = one+y*(v1+y*(v2+y*v3)); - r += (p1/p2-y/2); + r += p1/p2-y/2; } } - else if(ix<0x41000000) { /* x < 8.0 */ + /* x < 8.0 */ + else if(ix<0x41000000) { i = x; y = x-i; p = y*(s0+y*(s1+y*(s2+y*s3))); q = one+y*(r1+y*(r2+y*r3)); r = y/2+p/q; z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ switch(i) { case 7: z *= (y+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**24 */ - } else if (ix < 0x4b800000) { + /* 8.0 <= x < 2**27 */ + } else if (ix < 0x4d000000) { t = __ieee754_logf(x); z = one/x; y = z*z; w = w0+z*(w1+y*w2); r = (x-half)*(t-one)+w; } else - /* 2**24 <= x <= inf */ + /* 2**27 <= x <= inf */ r = x*(__ieee754_logf(x)-one); if(hx<0) r = nadj - r; return r; }