diff options
43 files changed, 1297 insertions, 511 deletions
diff --git a/lib/msun/Makefile b/lib/msun/Makefile index 4473623..f34255f 100644 --- a/lib/msun/Makefile +++ b/lib/msun/Makefile @@ -98,6 +98,7 @@ COMMON_SRCS+= s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c # If long double != double use these; otherwise, we alias the double versions. COMMON_SRCS+= e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \ e_coshl.c e_fmodl.c e_hypotl.c \ + e_lgammal.c e_lgammal_r.c \ e_remainderl.c e_sinhl.c e_sqrtl.c \ invtrig.c k_cosl.c k_sinl.c k_tanl.c \ s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \ @@ -188,7 +189,8 @@ MLINKS+=ilogb.3 ilogbf.3 ilogb.3 ilogbl.3 \ ilogb.3 logb.3 ilogb.3 logbf.3 ilogb.3 logbl.3 MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3 MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3 -MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 lgamma.3 lgammaf.3 \ +MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 \ + lgamma.3 lgammaf.3 lgamma.3 lgammal.3 \ lgamma.3 tgamma.3 lgamma.3 tgammaf.3 MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log10l.3 \ log.3 log1p.3 log.3 log1pf.3 log.3 log1pl.3 \ diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map index e53ca07..d5b7f46 100644 --- a/lib/msun/Symbol.map +++ b/lib/msun/Symbol.map @@ -269,6 +269,7 @@ FBSD_1.3 { erfl; expl; expm1l; + lgammal; log10l; log1pl; log2l; @@ -276,7 +277,11 @@ FBSD_1.3 { sinhl; tanhl; /* Implemented as weak aliases for imprecise versions */ - lgammal; powl; tgammal; }; + +/* First added in 11.0-CURRENT */ +FBSD_1.4 { + lgammal_r; +}; diff --git a/lib/msun/ld128/e_lgammal_r.c b/lib/msun/ld128/e_lgammal_r.c new file mode 100644 index 0000000..53d3af1 --- /dev/null +++ b/lib/msun/ld128/e_lgammal_r.c @@ -0,0 +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 <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.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 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,ix; + + EXTRACT_LDBL128_WORDS(hx,lx,llx,x); + + /* purge +-Inf and NaNs */ + *signgamp = 1; + ix = hx&0x7fff; + if(ix==0x7fff) return x*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) { + *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<zero) *signgamp = -1; + x = -x; + } + + /* purge 1 and 2 */ + if((ix==0x3fff || ix==0x4000) && (lx|llx)==0) r = 0; + /* for x < 2.0 */ + else if(ix<0x4000) { + 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.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; + 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 += p1/p2-y/2; + } + } + /* 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; + } + /* 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; +} diff --git a/lib/msun/ld128/k_expl.h b/lib/msun/ld128/k_expl.h index a5668fd..e0a48fc 100644 --- a/lib/msun/ld128/k_expl.h +++ b/lib/msun/ld128/k_expl.h @@ -322,7 +322,7 @@ __ldexp_cexpl(long double complex z, int expt) scale2 = 1; SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt); - return (cpackl(cos(y) * exp_x * scale1 * scale2, + return (CMPLXL(cos(y) * exp_x * scale1 * scale2, sinl(y) * exp_x * scale1 * scale2)); } #endif /* _COMPLEX_H */ diff --git a/lib/msun/ld80/e_lgammal_r.c b/lib/msun/ld80/e_lgammal_r.c new file mode 100644 index 0000000..1e65769 --- /dev/null +++ b/lib/msun/ld80/e_lgammal_r.c @@ -0,0 +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 <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +/* + * See e_lgamma_r.c for complete comments. + * + * Converted to long double by Steven G. Kargl. + */ + +#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+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,ix; + + EXTRACT_LDBL80_WORDS(hx,lx,x); + + /* purge +-Inf and NaNs */ + *signgamp = 1; + ix = hx&0x7fff; + if(ix==0x7fff) return x*x; + + ENTERI(); + + /* 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) { + *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<zero) *signgamp = -1; + x = -x; + } + + /* purge 1 and 2 */ + if((ix==0x3fff || ix==0x4000) && lx==0x8000000000000000ULL) r = 0; + /* for x < 2.0 */ + 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.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; + 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 += p1/p2-y/2; + } + } + /* 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**(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 + r = x*(logl(x)-1); + if(hx&0x8000) r = nadj - r; + RETURNI(r); +} diff --git a/lib/msun/ld80/k_expl.h b/lib/msun/ld80/k_expl.h index ebfb9a8..9b081fa 100644 --- a/lib/msun/ld80/k_expl.h +++ b/lib/msun/ld80/k_expl.h @@ -299,7 +299,7 @@ __ldexp_cexpl(long double complex z, int expt) scale2 = 1; SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt); - return (cpackl(cos(y) * exp_x * scale1 * scale2, + return (CMPLXL(cos(y) * exp_x * scale1 * scale2, sinl(y) * exp_x * scale1 * scale2)); } #endif /* _COMPLEX_H */ diff --git a/lib/msun/man/j0.3 b/lib/msun/man/j0.3 index 91849da..7e1b790 100644 --- a/lib/msun/man/j0.3 +++ b/lib/msun/man/j0.3 @@ -28,7 +28,7 @@ .\" from: @(#)j0.3 6.7 (Berkeley) 4/19/91 .\" $FreeBSD$ .\" -.Dd February 18, 2008 +.Dd March 10, 2015 .Dt J0 3 .Os .Sh NAME @@ -77,24 +77,17 @@ The functions .Fn j0 , .Fn j0f , -.Fn j1 +.Fn j1 , and .Fn j1f -compute the -.Em Bessel function of the first kind of the order -0 and the -.Em order -1, respectively, -for the -real value +compute the Bessel function of the first kind of orders +0 and 1 for the real value .Fa x ; the functions .Fn jn and .Fn jnf -compute the -.Em Bessel function of the first kind of the integer -.Em order +compute the Bessel function of the first kind of the integer order .Fa n for the real value .Fa x . @@ -105,13 +98,8 @@ The functions .Fn y1 , and .Fn y1f -compute the linearly independent -.Em Bessel function of the second kind of the order -0 and the -.Em order -1, respectively, -for the -positive +compute the linearly independent Bessel function of the second kind +of orders 0 and 1 for the positive .Em real value .Fa x ; @@ -119,9 +107,7 @@ the functions .Fn yn and .Fn ynf -compute the -.Em Bessel function of the second kind for the integer -.Em order +compute the Bessel function of the second kind for the integer order .Fa n for the positive .Em real @@ -141,11 +127,20 @@ and .Fn ynf . If .Fa x -is negative, these routines will generate an invalid exception and -return \*(Na. +is negative, including -\*(If, these routines will generate an invalid +exception and return \*(Na. +If +.Fa x +is \*(Pm0, these routines +will generate a divide-by-zero exception and return -\*(If. If .Fa x -is 0 or a sufficiently small positive number, these routines +is a sufficiently small positive number, then +.Fn y1 , +.Fn y1f , +.Fn yn , +and +.Fn ynf will generate an overflow exception and return -\*(If. .Sh SEE ALSO .Xr math 3 diff --git a/lib/msun/man/lgamma.3 b/lib/msun/man/lgamma.3 index ea6eb6f..d244722 100644 --- a/lib/msun/man/lgamma.3 +++ b/lib/msun/man/lgamma.3 @@ -28,7 +28,7 @@ .\" from: @(#)lgamma.3 6.6 (Berkeley) 12/3/92 .\" $FreeBSD$ .\" -.Dd January 14, 2005 +.Dd September 12, 2014 .Dt LGAMMA 3 .Os .Sh NAME @@ -36,6 +36,8 @@ .Nm lgamma_r , .Nm lgammaf , .Nm lgammaf_r , +.Nm lgammal , +.Nm lgammal_r , .Nm gamma , .Nm gamma_r , .Nm gammaf , @@ -58,6 +60,10 @@ .Fn lgammaf "float x" .Ft float .Fn lgammaf_r "float x" "int *signgamp" +.Ft "long double" +.Fn lgammal "long double x" +.Ft "long double" +.Fn lgammal_r "long double x" "int *signgamp" .Ft double .Fn gamma "double x" .Ft double @@ -66,14 +72,15 @@ .Fn gammaf "float x" .Ft float .Fn gammaf_r "float x" "int *signgamp" -.Ft double +.Ft "long double" .Fn tgamma "double x" .Ft float .Fn tgammaf "float x" .Sh DESCRIPTION -.Fn lgamma x +.Fn lgamma x , +.Fn lgammaf x , and -.Fn lgammaf x +.Fn lgammal x .if t \{\ return ln\||\(*G(x)| where .Bd -unfilled -offset indent @@ -87,13 +94,15 @@ The external integer .Fa signgam returns the sign of \(*G(x). .Pp -.Fn lgamma_r x signgamp +.Fn lgamma_r x signgamp , +.Fn lgammaf_r x signgamp , and -.Fn lgammaf_r x signgamp +.Fn lgammal_r x signgamp provide the same functionality as -.Fn lgamma x +.Fn lgamma x , +.Fn lgammaf x , and -.Fn lgammaf x +.Fn lgammal x , but the caller must provide an integer to store the sign of \(*G(x). .Pp The @@ -115,6 +124,7 @@ are deprecated aliases for and .Fn lgammaf_r , respectively. + .Sh IDIOSYNCRASIES Do not use the expression .Dq Li signgam\(**exp(lgamma(x)) @@ -139,14 +149,18 @@ Exponentiation of will lose up to 10 significant bits. .Sh RETURN VALUES .Fn gamma , -.Fn gamma_r , .Fn gammaf , +.Fn gammal , +.Fn gamma_r , .Fn gammaf_r , +.Fn gammal_r , .Fn lgamma , -.Fn lgamma_r , .Fn lgammaf , +.Fn lgammal , +.Fn lgamma_r , +.Fn lgammaf_r , and -.Fn lgammaf_r +.Fn lgammal_r return appropriate values unless an argument is out of range. Overflow will occur for sufficiently large positive values, and non-positive integers. @@ -159,6 +173,7 @@ will underflow. The .Fn lgamma , .Fn lgammaf , +.Fn lgammal , .Fn tgamma , and .Fn tgammaf diff --git a/lib/msun/src/catrig.c b/lib/msun/src/catrig.c index 200977c..c0f5f55 100644 --- a/lib/msun/src/catrig.c +++ b/lib/msun/src/catrig.c @@ -286,19 +286,19 @@ casinh(double complex z) if (isnan(x) || isnan(y)) { /* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */ if (isinf(x)) - return (cpack(x, y + y)); + return (CMPLX(x, y + y)); /* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */ if (isinf(y)) - return (cpack(y, x + x)); + return (CMPLX(y, x + x)); /* casinh(NaN + I*0) = NaN + I*0 */ if (y == 0) - return (cpack(x + x, y)); + return (CMPLX(x + x, y)); /* * All other cases involving NaN return NaN + I*NaN. * C99 leaves it optional whether to raise invalid if one of * the arguments is not NaN, so we opt not to raise it. */ - return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { @@ -307,7 +307,7 @@ casinh(double complex z) w = clog_for_large_values(z) + m_ln2; else w = clog_for_large_values(-z) + m_ln2; - return (cpack(copysign(creal(w), x), copysign(cimag(w), y))); + return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y))); } /* Avoid spuriously raising inexact for z = 0. */ @@ -325,7 +325,7 @@ casinh(double complex z) ry = asin(B); else ry = atan2(new_y, sqrt_A2my2); - return (cpack(copysign(rx, x), copysign(ry, y))); + return (CMPLX(copysign(rx, x), copysign(ry, y))); } /* @@ -335,9 +335,9 @@ casinh(double complex z) double complex casin(double complex z) { - double complex w = casinh(cpack(cimag(z), creal(z))); + double complex w = casinh(CMPLX(cimag(z), creal(z))); - return (cpack(cimag(w), creal(w))); + return (CMPLX(cimag(w), creal(w))); } /* @@ -370,19 +370,19 @@ cacos(double complex z) if (isnan(x) || isnan(y)) { /* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */ if (isinf(x)) - return (cpack(y + y, -INFINITY)); + return (CMPLX(y + y, -INFINITY)); /* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */ if (isinf(y)) - return (cpack(x + x, -y)); + return (CMPLX(x + x, -y)); /* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */ if (x == 0) - return (cpack(pio2_hi + pio2_lo, y + y)); + return (CMPLX(pio2_hi + pio2_lo, y + y)); /* * All other cases involving NaN return NaN + I*NaN. * C99 leaves it optional whether to raise invalid if one of * the arguments is not NaN, so we opt not to raise it. */ - return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { @@ -392,18 +392,18 @@ cacos(double complex z) ry = creal(w) + m_ln2; if (sy == 0) ry = -ry; - return (cpack(rx, ry)); + return (CMPLX(rx, ry)); } /* Avoid spuriously raising inexact for z = 1. */ if (x == 1 && y == 0) - return (cpack(0, -y)); + return (CMPLX(0, -y)); /* All remaining cases are inexact. */ raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) - return (cpack(pio2_hi - (x - pio2_lo), -y)); + return (CMPLX(pio2_hi - (x - pio2_lo), -y)); do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); if (B_is_usable) { @@ -419,7 +419,7 @@ cacos(double complex z) } if (sy == 0) ry = -ry; - return (cpack(rx, ry)); + return (CMPLX(rx, ry)); } /* @@ -437,15 +437,15 @@ cacosh(double complex z) ry = cimag(w); /* cacosh(NaN + I*NaN) = NaN + I*NaN */ if (isnan(rx) && isnan(ry)) - return (cpack(ry, rx)); + return (CMPLX(ry, rx)); /* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */ /* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */ if (isnan(rx)) - return (cpack(fabs(ry), rx)); + return (CMPLX(fabs(ry), rx)); /* cacosh(0 + I*NaN) = NaN + I*NaN */ if (isnan(ry)) - return (cpack(ry, ry)); - return (cpack(fabs(ry), copysign(rx, cimag(z)))); + return (CMPLX(ry, ry)); + return (CMPLX(fabs(ry), copysign(rx, cimag(z)))); } /* @@ -475,16 +475,16 @@ clog_for_large_values(double complex z) * this method is still poor since it is uneccessarily slow. */ if (ax > DBL_MAX / 2) - return (cpack(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x))); + return (CMPLX(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x))); /* * Avoid overflow when x or y is large. Avoid underflow when x or * y is small. */ if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) - return (cpack(log(hypot(x, y)), atan2(y, x))); + return (CMPLX(log(hypot(x, y)), atan2(y, x))); - return (cpack(log(ax * ax + ay * ay) / 2, atan2(y, x))); + return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x))); } /* @@ -575,30 +575,30 @@ catanh(double complex z) /* This helps handle many cases. */ if (y == 0 && ax <= 1) - return (cpack(atanh(x), y)); + return (CMPLX(atanh(x), y)); /* To ensure the same accuracy as atan(), and to filter out z = 0. */ if (x == 0) - return (cpack(x, atan(y))); + return (CMPLX(x, atan(y))); if (isnan(x) || isnan(y)) { /* catanh(+-Inf + I*NaN) = +-0 + I*NaN */ if (isinf(x)) - return (cpack(copysign(0, x), y + y)); + return (CMPLX(copysign(0, x), y + y)); /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */ if (isinf(y)) - return (cpack(copysign(0, x), + return (CMPLX(copysign(0, x), copysign(pio2_hi + pio2_lo, y))); /* * All other cases involving NaN return NaN + I*NaN. * C99 leaves it optional whether to raise invalid if one of * the arguments is not NaN, so we opt not to raise it. */ - return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) - return (cpack(real_part_reciprocal(x, y), + return (CMPLX(real_part_reciprocal(x, y), copysign(pio2_hi + pio2_lo, y))); if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { @@ -623,7 +623,7 @@ catanh(double complex z) else ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; - return (cpack(copysign(rx, x), copysign(ry, y))); + return (CMPLX(copysign(rx, x), copysign(ry, y))); } /* @@ -633,7 +633,7 @@ catanh(double complex z) double complex catan(double complex z) { - double complex w = catanh(cpack(cimag(z), creal(z))); + double complex w = catanh(CMPLX(cimag(z), creal(z))); - return (cpack(cimag(w), creal(w))); + return (CMPLX(cimag(w), creal(w))); } diff --git a/lib/msun/src/catrigf.c b/lib/msun/src/catrigf.c index 08ebef7..ca84ce2 100644 --- a/lib/msun/src/catrigf.c +++ b/lib/msun/src/catrigf.c @@ -156,12 +156,12 @@ casinhf(float complex z) if (isnan(x) || isnan(y)) { if (isinf(x)) - return (cpackf(x, y + y)); + return (CMPLXF(x, y + y)); if (isinf(y)) - return (cpackf(y, x + x)); + return (CMPLXF(y, x + x)); if (y == 0) - return (cpackf(x + x, y)); - return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXF(x + x, y)); + return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { @@ -169,7 +169,7 @@ casinhf(float complex z) w = clog_for_large_values(z) + m_ln2; else w = clog_for_large_values(-z) + m_ln2; - return (cpackf(copysignf(crealf(w), x), + return (CMPLXF(copysignf(crealf(w), x), copysignf(cimagf(w), y))); } @@ -186,15 +186,15 @@ casinhf(float complex z) ry = asinf(B); else ry = atan2f(new_y, sqrt_A2my2); - return (cpackf(copysignf(rx, x), copysignf(ry, y))); + return (CMPLXF(copysignf(rx, x), copysignf(ry, y))); } float complex casinf(float complex z) { - float complex w = casinhf(cpackf(cimagf(z), crealf(z))); + float complex w = casinhf(CMPLXF(cimagf(z), crealf(z))); - return (cpackf(cimagf(w), crealf(w))); + return (CMPLXF(cimagf(w), crealf(w))); } float complex @@ -214,12 +214,12 @@ cacosf(float complex z) if (isnan(x) || isnan(y)) { if (isinf(x)) - return (cpackf(y + y, -INFINITY)); + return (CMPLXF(y + y, -INFINITY)); if (isinf(y)) - return (cpackf(x + x, -y)); + return (CMPLXF(x + x, -y)); if (x == 0) - return (cpackf(pio2_hi + pio2_lo, y + y)); - return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXF(pio2_hi + pio2_lo, y + y)); + return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { @@ -228,16 +228,16 @@ cacosf(float complex z) ry = crealf(w) + m_ln2; if (sy == 0) ry = -ry; - return (cpackf(rx, ry)); + return (CMPLXF(rx, ry)); } if (x == 1 && y == 0) - return (cpackf(0, -y)); + return (CMPLXF(0, -y)); raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) - return (cpackf(pio2_hi - (x - pio2_lo), -y)); + return (CMPLXF(pio2_hi - (x - pio2_lo), -y)); do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); if (B_is_usable) { @@ -253,7 +253,7 @@ cacosf(float complex z) } if (sy == 0) ry = -ry; - return (cpackf(rx, ry)); + return (CMPLXF(rx, ry)); } float complex @@ -266,12 +266,12 @@ cacoshf(float complex z) rx = crealf(w); ry = cimagf(w); if (isnan(rx) && isnan(ry)) - return (cpackf(ry, rx)); + return (CMPLXF(ry, rx)); if (isnan(rx)) - return (cpackf(fabsf(ry), rx)); + return (CMPLXF(fabsf(ry), rx)); if (isnan(ry)) - return (cpackf(ry, ry)); - return (cpackf(fabsf(ry), copysignf(rx, cimagf(z)))); + return (CMPLXF(ry, ry)); + return (CMPLXF(fabsf(ry), copysignf(rx, cimagf(z)))); } static float complex @@ -291,13 +291,13 @@ clog_for_large_values(float complex z) } if (ax > FLT_MAX / 2) - return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1, + return (CMPLXF(logf(hypotf(x / m_e, y / m_e)) + 1, atan2f(y, x))); if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) - return (cpackf(logf(hypotf(x, y)), atan2f(y, x))); + return (CMPLXF(logf(hypotf(x, y)), atan2f(y, x))); - return (cpackf(logf(ax * ax + ay * ay) / 2, atan2f(y, x))); + return (CMPLXF(logf(ax * ax + ay * ay) / 2, atan2f(y, x))); } static inline float @@ -346,22 +346,22 @@ catanhf(float complex z) ay = fabsf(y); if (y == 0 && ax <= 1) - return (cpackf(atanhf(x), y)); + return (CMPLXF(atanhf(x), y)); if (x == 0) - return (cpackf(x, atanf(y))); + return (CMPLXF(x, atanf(y))); if (isnan(x) || isnan(y)) { if (isinf(x)) - return (cpackf(copysignf(0, x), y + y)); + return (CMPLXF(copysignf(0, x), y + y)); if (isinf(y)) - return (cpackf(copysignf(0, x), + return (CMPLXF(copysignf(0, x), copysignf(pio2_hi + pio2_lo, y))); - return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); + return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) - return (cpackf(real_part_reciprocal(x, y), + return (CMPLXF(real_part_reciprocal(x, y), copysignf(pio2_hi + pio2_lo, y))); if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { @@ -381,13 +381,13 @@ catanhf(float complex z) else ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; - return (cpackf(copysignf(rx, x), copysignf(ry, y))); + return (CMPLXF(copysignf(rx, x), copysignf(ry, y))); } float complex catanf(float complex z) { - float complex w = catanhf(cpackf(cimagf(z), crealf(z))); + float complex w = catanhf(CMPLXF(cimagf(z), crealf(z))); - return (cpackf(cimagf(w), crealf(w))); + return (CMPLXF(cimagf(w), crealf(w))); } diff --git a/lib/msun/src/e_j0.c b/lib/msun/src/e_j0.c index 8320f25..a253ed1 100644 --- a/lib/msun/src/e_j0.c +++ b/lib/msun/src/e_j0.c @@ -62,7 +62,9 @@ __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" -static double pzero(double), qzero(double); +static __inline double pzero(double), qzero(double); + +static const volatile double vone = 1, vzero = 0; static const double huge = 1e300, @@ -115,7 +117,7 @@ __ieee754_j0(double x) if(ix<0x3f200000) { /* |x| < 2**-13 */ if(huge+x>one) { /* raise inexact if x != 0 */ if(ix<0x3e400000) return one; /* |x|<2**-27 */ - else return one - 0.25*x*x; + else return one - x*x/4; } } z = x*x; @@ -150,10 +152,16 @@ __ieee754_y0(double x) EXTRACT_WORDS(hx,lx,x); ix = 0x7fffffff&hx; - /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ - if(ix>=0x7ff00000) return one/(x+x*x); - if((ix|lx)==0) return -one/zero; - if(hx<0) return zero/zero; + /* + * y0(NaN) = NaN. + * y0(Inf) = 0. + * y0(-Inf) = NaN and raise invalid exception. + */ + if(ix>=0x7ff00000) return vone/(x+x*x); + /* y0(+-0) = -inf and raise divide-by-zero exception. */ + if((ix|lx)==0) return -one/vzero; + /* y0(x<0) = NaN and raise invalid exception. */ + if(hx<0) return vzero/vzero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) * where x0 = x-pi/4 @@ -268,7 +276,8 @@ static const double pS2[5] = { 1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */ }; - static double pzero(double x) +static __inline double +pzero(double x) { const double *p,*q; double z,r,s; @@ -278,7 +287,7 @@ static const double pS2[5] = { if(ix>=0x40200000) {p = pR8; q= pS8;} else if(ix>=0x40122E8B){p = pR5; q= pS5;} else if(ix>=0x4006DB6D){p = pR3; q= pS3;} - else if(ix>=0x40000000){p = pR2; q= pS2;} + else {p = pR2; q= pS2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); @@ -363,7 +372,8 @@ static const double qS2[6] = { -5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */ }; - static double qzero(double x) +static __inline double +qzero(double x) { const double *p,*q; double s,r,z; @@ -373,7 +383,7 @@ static const double qS2[6] = { if(ix>=0x40200000) {p = qR8; q= qS8;} else if(ix>=0x40122E8B){p = qR5; q= qS5;} else if(ix>=0x4006DB6D){p = qR3; q= qS3;} - else if(ix>=0x40000000){p = qR2; q= qS2;} + else {p = qR2; q= qS2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); diff --git a/lib/msun/src/e_j0f.c b/lib/msun/src/e_j0f.c index c45faf3..3e2f7e8 100644 --- a/lib/msun/src/e_j0f.c +++ b/lib/msun/src/e_j0f.c @@ -16,10 +16,16 @@ #include <sys/cdefs.h> __FBSDID("$FreeBSD$"); +/* + * See e_j0.c for complete comments. + */ + #include "math.h" #include "math_private.h" -static float pzerof(float), qzerof(float); +static __inline float pzerof(float), qzerof(float); + +static const volatile float vone = 1, vzero = 0; static const float huge = 1e30, @@ -62,17 +68,17 @@ __ieee754_j0f(float x) * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) */ - if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x); + if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(x); /* |x|>2**49 */ else { u = pzerof(x); v = qzerof(x); z = invsqrtpi*(u*cc-v*ss)/sqrtf(x); } return z; } - if(ix<0x39000000) { /* |x| < 2**-13 */ + if(ix<0x3b000000) { /* |x| < 2**-9 */ if(huge+x>one) { /* raise inexact if x != 0 */ - if(ix<0x32000000) return one; /* |x|<2**-27 */ - else return one - (float)0.25*x*x; + if(ix<0x39800000) return one; /* |x|<2**-12 */ + else return one - x*x/4; } } z = x*x; @@ -107,10 +113,9 @@ __ieee754_y0f(float x) GET_FLOAT_WORD(hx,x); ix = 0x7fffffff&hx; - /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ - if(ix>=0x7f800000) return one/(x+x*x); - if(ix==0) return -one/zero; - if(hx<0) return zero/zero; + if(ix>=0x7f800000) return vone/(x+x*x); + if(ix==0) return -one/vzero; + if(hx<0) return vzero/vzero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) * where x0 = x-pi/4 @@ -136,14 +141,14 @@ __ieee754_y0f(float x) if ((s*c)<zero) cc = z/ss; else ss = z/cc; } - if(ix>0x80000000) z = (invsqrtpi*ss)/sqrtf(x); + if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */ else { u = pzerof(x); v = qzerof(x); z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); } return z; } - if(ix<=0x32000000) { /* x < 2**-27 */ + if(ix<=0x39000000) { /* x < 2**-13 */ return(u00 + tpi*__ieee754_logf(x)); } z = x*x; @@ -224,7 +229,8 @@ static const float pS2[5] = { 1.4657617569e+01, /* 0x416a859a */ }; - static float pzerof(float x) +static __inline float +pzerof(float x) { const float *p,*q; float z,r,s; @@ -232,9 +238,9 @@ static const float pS2[5] = { GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; if(ix>=0x41000000) {p = pR8; q= pS8;} - else if(ix>=0x40f71c58){p = pR5; q= pS5;} - else if(ix>=0x4036db68){p = pR3; q= pS3;} - else if(ix>=0x40000000){p = pR2; q= pS2;} + else if(ix>=0x409173eb){p = pR5; q= pS5;} + else if(ix>=0x4036d917){p = pR3; q= pS3;} + else {p = pR2; q= pS2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); @@ -319,7 +325,8 @@ static const float qS2[6] = { -5.3109550476e+00, /* 0xc0a9f358 */ }; - static float qzerof(float x) +static __inline float +qzerof(float x) { const float *p,*q; float s,r,z; @@ -327,9 +334,9 @@ static const float qS2[6] = { GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; if(ix>=0x41000000) {p = qR8; q= qS8;} - else if(ix>=0x40f71c58){p = qR5; q= qS5;} - else if(ix>=0x4036db68){p = qR3; q= qS3;} - else if(ix>=0x40000000){p = qR2; q= qS2;} + else if(ix>=0x409173eb){p = qR5; q= qS5;} + else if(ix>=0x4036d917){p = qR3; q= qS3;} + else {p = qR2; q= qS2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); diff --git a/lib/msun/src/e_j1.c b/lib/msun/src/e_j1.c index 63800ad..74a7c69 100644 --- a/lib/msun/src/e_j1.c +++ b/lib/msun/src/e_j1.c @@ -62,7 +62,9 @@ __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" -static double pone(double), qone(double); +static __inline double pone(double), qone(double); + +static const volatile double vone = 1, vzero = 0; static const double huge = 1e300, @@ -147,10 +149,16 @@ __ieee754_y1(double x) EXTRACT_WORDS(hx,lx,x); ix = 0x7fffffff&hx; - /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if(ix>=0x7ff00000) return one/(x+x*x); - if((ix|lx)==0) return -one/zero; - if(hx<0) return zero/zero; + /* + * y1(NaN) = NaN. + * y1(Inf) = 0. + * y1(-Inf) = NaN and raise invalid exception. + */ + if(ix>=0x7ff00000) return vone/(x+x*x); + /* y1(+-0) = -inf and raise divide-by-zero exception. */ + if((ix|lx)==0) return -one/vzero; + /* y1(x<0) = NaN and raise invalid exception. */ + if(hx<0) return vzero/vzero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ s = sin(x); c = cos(x); @@ -262,7 +270,8 @@ static const double ps2[5] = { 8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */ }; - static double pone(double x) +static __inline double +pone(double x) { const double *p,*q; double z,r,s; @@ -272,7 +281,7 @@ static const double ps2[5] = { if(ix>=0x40200000) {p = pr8; q= ps8;} else if(ix>=0x40122E8B){p = pr5; q= ps5;} else if(ix>=0x4006DB6D){p = pr3; q= ps3;} - else if(ix>=0x40000000){p = pr2; q= ps2;} + else {p = pr2; q= ps2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); @@ -358,7 +367,8 @@ static const double qs2[6] = { -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */ }; - static double qone(double x) +static __inline double +qone(double x) { const double *p,*q; double s,r,z; @@ -368,7 +378,7 @@ static const double qs2[6] = { if(ix>=0x40200000) {p = qr8; q= qs8;} else if(ix>=0x40122E8B){p = qr5; q= qs5;} else if(ix>=0x4006DB6D){p = qr3; q= qs3;} - else if(ix>=0x40000000){p = qr2; q= qs2;} + else {p = qr2; q= qs2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); diff --git a/lib/msun/src/e_j1f.c b/lib/msun/src/e_j1f.c index 88e2d83..ec7f381 100644 --- a/lib/msun/src/e_j1f.c +++ b/lib/msun/src/e_j1f.c @@ -16,10 +16,16 @@ #include <sys/cdefs.h> __FBSDID("$FreeBSD$"); +/* + * See e_j1.c for complete comments. + */ + #include "math.h" #include "math_private.h" -static float ponef(float), qonef(float); +static __inline float ponef(float), qonef(float); + +static const volatile float vone = 1, vzero = 0; static const float huge = 1e30, @@ -63,7 +69,7 @@ __ieee754_j1f(float x) * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) */ - if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(y); + if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(y); /* |x|>2**49 */ else { u = ponef(y); v = qonef(y); z = invsqrtpi*(u*cc-v*ss)/sqrtf(y); @@ -71,7 +77,7 @@ __ieee754_j1f(float x) if(hx<0) return -z; else return z; } - if(ix<0x32000000) { /* |x|<2**-27 */ + if(ix<0x39000000) { /* |x|<2**-13 */ if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */ } z = x*x; @@ -104,10 +110,9 @@ __ieee754_y1f(float x) GET_FLOAT_WORD(hx,x); ix = 0x7fffffff&hx; - /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if(ix>=0x7f800000) return one/(x+x*x); - if(ix==0) return -one/zero; - if(hx<0) return zero/zero; + if(ix>=0x7f800000) return vone/(x+x*x); + if(ix==0) return -one/vzero; + if(hx<0) return vzero/vzero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ s = sinf(x); c = cosf(x); @@ -129,14 +134,14 @@ __ieee754_y1f(float x) * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one. */ - if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x); + if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */ else { u = ponef(x); v = qonef(x); z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); } return z; } - if(ix<=0x24800000) { /* x < 2**-54 */ + if(ix<=0x33000000) { /* x < 2**-25 */ return(-tpi/x); } z = x*x; @@ -219,7 +224,8 @@ static const float ps2[5] = { 8.3646392822e+00, /* 0x4105d590 */ }; - static float ponef(float x) +static __inline float +ponef(float x) { const float *p,*q; float z,r,s; @@ -227,9 +233,9 @@ static const float ps2[5] = { GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; if(ix>=0x41000000) {p = pr8; q= ps8;} - else if(ix>=0x40f71c58){p = pr5; q= ps5;} - else if(ix>=0x4036db68){p = pr3; q= ps3;} - else if(ix>=0x40000000){p = pr2; q= ps2;} + else if(ix>=0x409173eb){p = pr5; q= ps5;} + else if(ix>=0x4036d917){p = pr3; q= ps3;} + else {p = pr2; q= ps2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); @@ -315,17 +321,18 @@ static const float qs2[6] = { -4.9594988823e+00, /* 0xc09eb437 */ }; - static float qonef(float x) +static __inline float +qonef(float x) { const float *p,*q; float s,r,z; int32_t ix; GET_FLOAT_WORD(ix,x); ix &= 0x7fffffff; - if(ix>=0x40200000) {p = qr8; q= qs8;} - else if(ix>=0x40f71c58){p = qr5; q= qs5;} - else if(ix>=0x4036db68){p = qr3; q= qs3;} - else if(ix>=0x40000000){p = qr2; q= qs2;} + if(ix>=0x41000000) {p = qr8; q= qs8;} + else if(ix>=0x409173eb){p = qr5; q= qs5;} + else if(ix>=0x4036d917){p = qr3; q= qs3;} + else {p = qr2; q= qs2;} /* ix>=0x40000000 */ z = one/(x*x); r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); diff --git a/lib/msun/src/e_jn.c b/lib/msun/src/e_jn.c index 8b0bc62..eefd4ff 100644 --- a/lib/msun/src/e_jn.c +++ b/lib/msun/src/e_jn.c @@ -43,6 +43,8 @@ __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" +static const volatile double vone = 1, vzero = 0; + static const double invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */ @@ -220,10 +222,12 @@ __ieee754_yn(int n, double x) EXTRACT_WORDS(hx,lx,x); ix = 0x7fffffff&hx; - /* if Y(n,NaN) is NaN */ + /* yn(n,NaN) = NaN */ if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x; - if((ix|lx)==0) return -one/zero; - if(hx<0) return zero/zero; + /* yn(n,+-0) = -inf and raise divide-by-zero exception. */ + if((ix|lx)==0) return -one/vzero; + /* yn(n,x<0) = NaN and raise invalid exception. */ + if(hx<0) return vzero/vzero; sign = 1; if(n<0){ n = -n; diff --git a/lib/msun/src/e_jnf.c b/lib/msun/src/e_jnf.c index f564aec..965feeb 100644 --- a/lib/msun/src/e_jnf.c +++ b/lib/msun/src/e_jnf.c @@ -16,9 +16,15 @@ #include <sys/cdefs.h> __FBSDID("$FreeBSD$"); +/* + * See e_jn.c for complete comments. + */ + #include "math.h" #include "math_private.h" +static const volatile float vone = 1, vzero = 0; + static const float two = 2.0000000000e+00, /* 0x40000000 */ one = 1.0000000000e+00; /* 0x3F800000 */ @@ -172,10 +178,9 @@ __ieee754_ynf(int n, float x) GET_FLOAT_WORD(hx,x); ix = 0x7fffffff&hx; - /* if Y(n,NaN) is NaN */ if(ix>0x7f800000) return x+x; - if(ix==0) return -one/zero; - if(hx<0) return zero/zero; + if(ix==0) return -one/vzero; + if(hx<0) return vzero/vzero; sign = 1; if(n<0){ n = -n; diff --git a/lib/msun/src/e_lgamma.c b/lib/msun/src/e_lgamma.c index 4674d9b..43f5175 100644 --- a/lib/msun/src/e_lgamma.c +++ b/lib/msun/src/e_lgamma.c @@ -21,6 +21,8 @@ __FBSDID("$FreeBSD$"); * Method: call __ieee754_lgamma_r */ +#include <float.h> + #include "math.h" #include "math_private.h" @@ -31,3 +33,7 @@ __ieee754_lgamma(double x) { return __ieee754_lgamma_r(x,&signgam); } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(lgamma, lgammal); +#endif diff --git a/lib/msun/src/e_lgamma_r.c b/lib/msun/src/e_lgamma_r.c index 7a95ea4..be70767 100644 --- a/lib/msun/src/e_lgamma_r.c +++ b/lib/msun/src/e_lgamma_r.c @@ -1,4 +1,3 @@ - /* @(#)e_lgamma_r.c 1.3 95/01/18 */ /* * ==================================================== @@ -6,22 +5,21 @@ * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== - * */ #include <sys/cdefs.h> __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). + * 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 + * 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, @@ -59,20 +57,20 @@ __FBSDID("$FreeBSD$"); * by * 3 5 11 * w = w0 + w1*z + w2*z + w3*z + ... + w6*z - * where + * 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 + * 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 + * 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 @@ -80,9 +78,10 @@ __FBSDID("$FreeBSD$"); * lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero * lgamma(inf) = inf * lgamma(-inf) = inf (bug for bug compatible with C99!?) - * */ +#include <float.h> + #include "math.h" #include "math_private.h" @@ -187,9 +186,9 @@ sin_pi(double x) switch (n) { case 0: y = __kernel_sin(pi*y,zero,0); break; - case 1: + case 1: case 2: y = __kernel_cos(pi*(0.5-y),zero); break; - case 3: + 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; @@ -202,24 +201,28 @@ sin_pi(double x) 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) return one/vzero; - if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */ - if(hx<0) { - *signgamp = -1; - return -__ieee754_log(-x); - } else return -__ieee754_log(x); + + /* 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)); } + + /* 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); @@ -229,7 +232,7 @@ __ieee754_lgamma_r(double x, int *signgamp) x = -x; } - /* purge off 1 and 2 */ + /* purge 1 and 2 */ if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0; /* for x < 2.0 */ else if(ix<0x40000000) { @@ -250,7 +253,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; @@ -258,38 +261,43 @@ __ieee754_lgamma_r(double x, int *signgamp) 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; - case 2: + 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.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 */ - } 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 */ + } else + /* 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 diff --git a/lib/msun/src/e_lgammaf_r.c b/lib/msun/src/e_lgammaf_r.c index 9a7ab39..9084e18 100644 --- a/lib/msun/src/e_lgammaf_r.c +++ b/lib/msun/src/e_lgammaf_r.c @@ -1,5 +1,6 @@ /* e_lgammaf_r.c -- float version of e_lgamma_r.c. * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + * Conversion to float fixed By Steven G. Kargl. */ /* @@ -22,72 +23,63 @@ __FBSDID("$FreeBSD$"); static const volatile float vzero = 0; static const float -zero= 0.0000000000e+00, -half= 5.0000000000e-01, /* 0x3f000000 */ -one = 1.0000000000e+00, /* 0x3f800000 */ +zero= 0, +half= 0.5, +one = 1, pi = 3.1415927410e+00, /* 0x40490fdb */ -a0 = 7.7215664089e-02, /* 0x3d9e233f */ -a1 = 3.2246702909e-01, /* 0x3ea51a66 */ -a2 = 6.7352302372e-02, /* 0x3d89f001 */ -a3 = 2.0580807701e-02, /* 0x3ca89915 */ -a4 = 7.3855509982e-03, /* 0x3bf2027e */ -a5 = 2.8905137442e-03, /* 0x3b3d6ec6 */ -a6 = 1.1927076848e-03, /* 0x3a9c54a1 */ -a7 = 5.1006977446e-04, /* 0x3a05b634 */ -a8 = 2.2086278477e-04, /* 0x39679767 */ -a9 = 1.0801156895e-04, /* 0x38e28445 */ -a10 = 2.5214456400e-05, /* 0x37d383a2 */ -a11 = 4.4864096708e-05, /* 0x383c2c75 */ -tc = 1.4616321325e+00, /* 0x3fbb16c3 */ -tf = -1.2148628384e-01, /* 0xbdf8cdcd */ -/* tt = -(tail of tf) */ -tt = 6.6971006518e-09, /* 0x31e61c52 */ -t0 = 4.8383611441e-01, /* 0x3ef7b95e */ -t1 = -1.4758771658e-01, /* 0xbe17213c */ -t2 = 6.4624942839e-02, /* 0x3d845a15 */ -t3 = -3.2788541168e-02, /* 0xbd064d47 */ -t4 = 1.7970675603e-02, /* 0x3c93373d */ -t5 = -1.0314224288e-02, /* 0xbc28fcfe */ -t6 = 6.1005386524e-03, /* 0x3bc7e707 */ -t7 = -3.6845202558e-03, /* 0xbb7177fe */ -t8 = 2.2596477065e-03, /* 0x3b141699 */ -t9 = -1.4034647029e-03, /* 0xbab7f476 */ -t10 = 8.8108185446e-04, /* 0x3a66f867 */ -t11 = -5.3859531181e-04, /* 0xba0d3085 */ -t12 = 3.1563205994e-04, /* 0x39a57b6b */ -t13 = -3.1275415677e-04, /* 0xb9a3f927 */ -t14 = 3.3552918467e-04, /* 0x39afe9f7 */ -u0 = -7.7215664089e-02, /* 0xbd9e233f */ -u1 = 6.3282704353e-01, /* 0x3f2200f4 */ -u2 = 1.4549225569e+00, /* 0x3fba3ae7 */ -u3 = 9.7771751881e-01, /* 0x3f7a4bb2 */ -u4 = 2.2896373272e-01, /* 0x3e6a7578 */ -u5 = 1.3381091878e-02, /* 0x3c5b3c5e */ -v1 = 2.4559779167e+00, /* 0x401d2ebe */ -v2 = 2.1284897327e+00, /* 0x4008392d */ -v3 = 7.6928514242e-01, /* 0x3f44efdf */ -v4 = 1.0422264785e-01, /* 0x3dd572af */ -v5 = 3.2170924824e-03, /* 0x3b52d5db */ -s0 = -7.7215664089e-02, /* 0xbd9e233f */ -s1 = 2.1498242021e-01, /* 0x3e5c245a */ -s2 = 3.2577878237e-01, /* 0x3ea6cc7a */ -s3 = 1.4635047317e-01, /* 0x3e15dce6 */ -s4 = 2.6642270386e-02, /* 0x3cda40e4 */ -s5 = 1.8402845599e-03, /* 0x3af135b4 */ -s6 = 3.1947532989e-05, /* 0x3805ff67 */ -r1 = 1.3920053244e+00, /* 0x3fb22d3b */ -r2 = 7.2193557024e-01, /* 0x3f38d0c5 */ -r3 = 1.7193385959e-01, /* 0x3e300f6e */ -r4 = 1.8645919859e-02, /* 0x3c98bf54 */ -r5 = 7.7794247773e-04, /* 0x3a4beed6 */ -r6 = 7.3266842264e-06, /* 0x36f5d7bd */ -w0 = 4.1893854737e-01, /* 0x3ed67f1d */ -w1 = 8.3333335817e-02, /* 0x3daaaaab */ -w2 = -2.7777778450e-03, /* 0xbb360b61 */ -w3 = 7.9365057172e-04, /* 0x3a500cfd */ -w4 = -5.9518753551e-04, /* 0xba1c065c */ -w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ -w6 = -1.6309292987e-03; /* 0xbad5c4e8 */ +/* + * Domain y in [0x1p-27, 0.27], range ~[-3.4599e-10, 3.4590e-10]: + * |(lgamma(2 - y) + 0.5 * y) / y - a(y)| < 2**-31.4 + */ +a0 = 7.72156641e-02, /* 0x3d9e233f */ +a1 = 3.22467119e-01, /* 0x3ea51a69 */ +a2 = 6.73484802e-02, /* 0x3d89ee00 */ +a3 = 2.06395667e-02, /* 0x3ca9144f */ +a4 = 6.98275631e-03, /* 0x3be4cf9b */ +a5 = 4.11768444e-03, /* 0x3b86eda4 */ +/* + * Domain x in [tc-0.24, tc+0.28], range ~[-5.6577e-10, 5.5677e-10]: + * |(lgamma(x) - tf) - t(x - tc)| < 2**-30.8. + */ +tc = 1.46163213e+00, /* 0x3fbb16c3 */ +tf = -1.21486291e-01, /* 0xbdf8cdce */ +t0 = -2.94064460e-11, /* 0xae0154b7 */ +t1 = -2.35939837e-08, /* 0xb2caabb8 */ +t2 = 4.83836412e-01, /* 0x3ef7b968 */ +t3 = -1.47586212e-01, /* 0xbe1720d7 */ +t4 = 6.46013096e-02, /* 0x3d844db1 */ +t5 = -3.28450352e-02, /* 0xbd068884 */ +t6 = 1.86483748e-02, /* 0x3c98c47a */ +t7 = -9.89206228e-03, /* 0xbc221251 */ +/* + * Domain y in [-0.1, 0.232], range ~[-8.4931e-10, 8.7794e-10]: + * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-31.2 + */ +u0 = -7.72156641e-02, /* 0xbd9e233f */ +u1 = 7.36789703e-01, /* 0x3f3c9e40 */ +u2 = 4.95649040e-01, /* 0x3efdc5b6 */ +v1 = 1.10958421e+00, /* 0x3f8e06db */ +v2 = 2.10598111e-01, /* 0x3e57a708 */ +v3 = -1.02995494e-02, /* 0xbc28bf71 */ +/* + * Domain x in (2, 3], range ~[-5.5189e-11, 5.2317e-11]: + * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-35.0 + * with y = x - 2. + */ +s0 = -7.72156641e-02, /* 0xbd9e233f */ +s1 = 2.69987404e-01, /* 0x3e8a3bca */ +s2 = 1.42851010e-01, /* 0x3e124789 */ +s3 = 1.19389519e-02, /* 0x3c439b98 */ +r1 = 6.79650068e-01, /* 0x3f2dfd8c */ +r2 = 1.16058730e-01, /* 0x3dedb033 */ +r3 = 3.75673687e-03, /* 0x3b763396 */ +/* + * Domain z in [8, 0x1p24], range ~[-1.2640e-09, 1.2640e-09]: + * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-29.6. + */ +w0 = 4.18938547e-01, /* 0x3ed67f1d */ +w1 = 8.33332464e-02, /* 0x3daaaa9f */ +w2 = -2.76129087e-03; /* 0xbb34f6c6 */ static float sin_pif(float x) @@ -130,25 +122,29 @@ sin_pif(float x) 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) return one/vzero; - if(ix<0x35000000) { /* |x|<2**-21, return -log(|x|) */ - if(hx<0) { - *signgamp = -1; - return -__ieee754_logf(-x); - } else return -__ieee754_logf(x); + + /* 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)); } + + /* 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 */ @@ -157,7 +153,7 @@ __ieee754_lgammaf_r(float x, int *signgamp) x = -x; } - /* purge off 1 and 2 */ + /* purge 1 and 2 */ if (ix==0x3f800000||ix==0x40000000) r = 0; /* for x < 2.0 */ else if(ix<0x40000000) { @@ -168,55 +164,51 @@ __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)); - r += (tf + p); break; + 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; + /* 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+(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**27 */ + } else if (ix < 0x4d000000) { 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**27 <= x <= inf */ r = x*(__ieee754_logf(x)-one); if(hx<0) r = nadj - r; return r; diff --git a/lib/msun/src/e_lgammal.c b/lib/msun/src/e_lgammal.c new file mode 100644 index 0000000..ebc2fc7 --- /dev/null +++ b/lib/msun/src/e_lgammal.c @@ -0,0 +1,25 @@ +/* @(#)e_lgamma.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#include "math.h" +#include "math_private.h" + +extern int signgam; + +long double +lgammal(long double x) +{ + return lgammal_r(x,&signgam); +} diff --git a/lib/msun/src/imprecise.c b/lib/msun/src/imprecise.c index 92fb2d0..08cd239 100644 --- a/lib/msun/src/imprecise.c +++ b/lib/msun/src/imprecise.c @@ -60,5 +60,4 @@ DECLARE_WEAK(powl); long double imprecise_ ## f ## l(long double v) { return f(v); }\ DECLARE_WEAK(f ## l) -DECLARE_IMPRECISE(lgamma); DECLARE_IMPRECISE(tgamma); diff --git a/lib/msun/src/k_exp.c b/lib/msun/src/k_exp.c index f592f69..ecef54c 100644 --- a/lib/msun/src/k_exp.c +++ b/lib/msun/src/k_exp.c @@ -103,6 +103,6 @@ __ldexp_cexp(double complex z, int expt) half_expt = expt - half_expt; INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); - return (cpack(cos(y) * exp_x * scale1 * scale2, + return (CMPLX(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2)); } diff --git a/lib/msun/src/k_expf.c b/lib/msun/src/k_expf.c index 548a008..f8c254d 100644 --- a/lib/msun/src/k_expf.c +++ b/lib/msun/src/k_expf.c @@ -82,6 +82,6 @@ __ldexp_cexpf(float complex z, int expt) half_expt = expt - half_expt; SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23); - return (cpackf(cosf(y) * exp_x * scale1 * scale2, + return (CMPLXF(cosf(y) * exp_x * scale1 * scale2, sinf(y) * exp_x * scale1 * scale2)); } diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h index 3ab76f8..be8240a 100644 --- a/lib/msun/src/math.h +++ b/lib/msun/src/math.h @@ -500,8 +500,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_ */ diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h index 8af2c65..afaf201 100644 --- a/lib/msun/src/math_private.h +++ b/lib/msun/src/math_private.h @@ -454,9 +454,15 @@ typedef union { * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product. * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted * to -0.0+I*0.0. + * + * The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL() + * to construct complex values. Compilers that conform to the C99 + * standard require the following functions to avoid the above issues. */ + +#ifndef CMPLXF static __inline float complex -cpackf(float x, float y) +CMPLXF(float x, float y) { float_complex z; @@ -464,9 +470,11 @@ cpackf(float x, float y) IMAGPART(z) = y; return (z.f); } +#endif +#ifndef CMPLX static __inline double complex -cpack(double x, double y) +CMPLX(double x, double y) { double_complex z; @@ -474,9 +482,11 @@ cpack(double x, double y) IMAGPART(z) = y; return (z.f); } +#endif +#ifndef CMPLXL static __inline long double complex -cpackl(long double x, long double y) +CMPLXL(long double x, long double y) { long_double_complex z; @@ -484,6 +494,8 @@ cpackl(long double x, long double y) IMAGPART(z) = y; return (z.f); } +#endif + #endif /* _COMPLEX_H */ #ifdef __GNUCLIKE_ASM diff --git a/lib/msun/src/s_ccosh.c b/lib/msun/src/s_ccosh.c index 9ea962b..2843198 100644 --- a/lib/msun/src/s_ccosh.c +++ b/lib/msun/src/s_ccosh.c @@ -32,6 +32,8 @@ * * Exceptional values are noted in the comments within the source code. * These values and the return value were taken from n1124.pdf. + * The sign of the result for some exceptional values is unspecified but + * must satisfy both cosh(conj(z)) == conj(cosh(z)) and cosh(-z) == cosh(z). */ #include <sys/cdefs.h> @@ -62,49 +64,48 @@ ccosh(double complex z) /* Handle the nearly-non-exceptional cases where x and y are finite. */ if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) - return (cpack(cosh(x), x * y)); - if (ix < 0x40360000) /* small x: normal case */ - return (cpack(cosh(x) * cos(y), sinh(x) * sin(y))); + return (CMPLX(cosh(x), x * y)); + if (ix < 0x40360000) /* |x| < 22: normal case */ + return (CMPLX(cosh(x) * cos(y), sinh(x) * sin(y))); /* |x| >= 22, so cosh(x) ~= exp(|x|) */ if (ix < 0x40862e42) { /* x < 710: exp(|x|) won't overflow */ h = exp(fabs(x)) * 0.5; - return (cpack(h * cos(y), copysign(h, x) * sin(y))); + return (CMPLX(h * cos(y), copysign(h, x) * sin(y))); } else if (ix < 0x4096bbaa) { /* x < 1455: scale to avoid overflow */ - z = __ldexp_cexp(cpack(fabs(x), y), -1); - return (cpack(creal(z), cimag(z) * copysign(1, x))); + z = __ldexp_cexp(CMPLX(fabs(x), y), -1); + return (CMPLX(creal(z), cimag(z) * copysign(1, x))); } else { /* x >= 1455: the result always overflows */ h = huge * x; - return (cpack(h * h * cos(y), h * sin(y))); + return (CMPLX(h * h * cos(y), h * sin(y))); } } /* - * cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0. - * The sign of 0 in the result is unspecified. Choice = normally - * the same as dNaN. Raise the invalid floating-point exception. + * cosh(+-0 +- I Inf) = dNaN + I (+-)(+-)0. + * The sign of 0 in the result is unspecified. Choice = product + * of the signs of the argument. Raise the invalid floating-point + * exception. * - * cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0. - * The sign of 0 in the result is unspecified. Choice = normally - * the same as d(NaN). + * cosh(+-0 +- I NaN) = d(NaN) + I (+-)(+-)0. + * The sign of 0 in the result is unspecified. Choice = product + * of the signs of the argument. */ - if ((ix | lx) == 0 && iy >= 0x7ff00000) - return (cpack(y - y, copysign(0, x * (y - y)))); + if ((ix | lx) == 0) /* && iy >= 0x7ff00000 */ + return (CMPLX(y - y, x * copysign(0, y))); /* * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0. * - * cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0. - * The sign of 0 in the result is unspecified. + * cosh(NaN +- I 0) = d(NaN) + I (+-)(+-)0. + * The sign of 0 in the result is unspecified. Choice = product + * of the signs of the argument. */ - if ((iy | ly) == 0 && ix >= 0x7ff00000) { - if (((hx & 0xfffff) | lx) == 0) - return (cpack(x * x, copysign(0, x) * y)); - return (cpack(x * x, copysign(0, (x + x) * y))); - } + if ((iy | ly) == 0) /* && ix >= 0x7ff00000 */ + return (CMPLX(x * x, copysign(0, x) * y)); /* * cosh(x +- I Inf) = dNaN + I dNaN. @@ -114,8 +115,8 @@ ccosh(double complex z) * Optionally raises the invalid floating-point exception for finite * nonzero x. Choice = don't raise (except for signaling NaNs). */ - if (ix < 0x7ff00000 && iy >= 0x7ff00000) - return (cpack(y - y, x * (y - y))); + if (ix < 0x7ff00000) /* && iy >= 0x7ff00000 */ + return (CMPLX(y - y, x * (y - y))); /* * cosh(+-Inf + I NaN) = +Inf + I d(NaN). @@ -126,10 +127,10 @@ ccosh(double complex z) * * cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y) */ - if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { + if (ix == 0x7ff00000 && lx == 0) { if (iy >= 0x7ff00000) - return (cpack(x * x, x * (y - y))); - return (cpack((x * x) * cos(y), x * sin(y))); + return (CMPLX(INFINITY, x * (y - y))); + return (CMPLX(INFINITY * cos(y), x * sin(y))); } /* @@ -143,7 +144,7 @@ ccosh(double complex z) * Optionally raises the invalid floating-point exception for finite * nonzero y. Choice = don't raise (except for signaling NaNs). */ - return (cpack((x * x) * (y - y), (x + x) * (y - y))); + return (CMPLX((x * x) * (y - y), (x + x) * (y - y))); } double complex @@ -151,5 +152,5 @@ ccos(double complex z) { /* ccos(z) = ccosh(I * z) */ - return (ccosh(cpack(-cimag(z), creal(z)))); + return (ccosh(CMPLX(-cimag(z), creal(z)))); } diff --git a/lib/msun/src/s_ccoshf.c b/lib/msun/src/s_ccoshf.c index 1de9ad4..eeed92f 100644 --- a/lib/msun/src/s_ccoshf.c +++ b/lib/msun/src/s_ccoshf.c @@ -25,7 +25,7 @@ */ /* - * Hyperbolic cosine of a complex argument. See s_ccosh.c for details. + * Float version of ccosh(). See s_ccosh.c for details. */ #include <sys/cdefs.h> @@ -55,50 +55,47 @@ ccoshf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) - return (cpackf(coshf(x), x * y)); - if (ix < 0x41100000) /* small x: normal case */ - return (cpackf(coshf(x) * cosf(y), sinhf(x) * sinf(y))); + return (CMPLXF(coshf(x), x * y)); + if (ix < 0x41100000) /* |x| < 9: normal case */ + return (CMPLXF(coshf(x) * cosf(y), sinhf(x) * sinf(y))); /* |x| >= 9, so cosh(x) ~= exp(|x|) */ if (ix < 0x42b17218) { /* x < 88.7: expf(|x|) won't overflow */ - h = expf(fabsf(x)) * 0.5f; - return (cpackf(h * cosf(y), copysignf(h, x) * sinf(y))); + h = expf(fabsf(x)) * 0.5F; + return (CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y))); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ - z = __ldexp_cexpf(cpackf(fabsf(x), y), -1); - return (cpackf(crealf(z), cimagf(z) * copysignf(1, x))); + z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1); + return (CMPLXF(crealf(z), cimagf(z) * copysignf(1, x))); } else { /* x >= 192.7: the result always overflows */ h = huge * x; - return (cpackf(h * h * cosf(y), h * sinf(y))); + return (CMPLXF(h * h * cosf(y), h * sinf(y))); } } - if (ix == 0 && iy >= 0x7f800000) - return (cpackf(y - y, copysignf(0, x * (y - y)))); + if (ix == 0) /* && iy >= 0x7f800000 */ + return (CMPLXF(y - y, x * copysignf(0, y))); - if (iy == 0 && ix >= 0x7f800000) { - if ((hx & 0x7fffff) == 0) - return (cpackf(x * x, copysignf(0, x) * y)); - return (cpackf(x * x, copysignf(0, (x + x) * y))); - } + if (iy == 0) /* && ix >= 0x7f800000 */ + return (CMPLXF(x * x, copysignf(0, x) * y)); - if (ix < 0x7f800000 && iy >= 0x7f800000) - return (cpackf(y - y, x * (y - y))); + if (ix < 0x7f800000) /* && iy >= 0x7f800000 */ + return (CMPLXF(y - y, x * (y - y))); - if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) { + if (ix == 0x7f800000) { if (iy >= 0x7f800000) - return (cpackf(x * x, x * (y - y))); - return (cpackf((x * x) * cosf(y), x * sinf(y))); + return (CMPLXF(INFINITY, x * (y - y))); + return (CMPLXF(INFINITY * cosf(y), x * sinf(y))); } - return (cpackf((x * x) * (y - y), (x + x) * (y - y))); + return (CMPLXF((x * x) * (y - y), (x + x) * (y - y))); } float complex ccosf(float complex z) { - return (ccoshf(cpackf(-cimagf(z), crealf(z)))); + return (ccoshf(CMPLXF(-cimagf(z), crealf(z)))); } diff --git a/lib/msun/src/s_cexp.c b/lib/msun/src/s_cexp.c index abe178f..9e11d51 100644 --- a/lib/msun/src/s_cexp.c +++ b/lib/msun/src/s_cexp.c @@ -50,22 +50,22 @@ cexp(double complex z) /* cexp(x + I 0) = exp(x) + I 0 */ if ((hy | ly) == 0) - return (cpack(exp(x), y)); + return (CMPLX(exp(x), y)); EXTRACT_WORDS(hx, lx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ if (((hx & 0x7fffffff) | lx) == 0) - return (cpack(cos(y), sin(y))); + return (CMPLX(cos(y), sin(y))); if (hy >= 0x7ff00000) { if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ - return (cpack(y - y, y - y)); + return (CMPLX(y - y, y - y)); } else if (hx & 0x80000000) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ - return (cpack(0.0, 0.0)); + return (CMPLX(0.0, 0.0)); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ - return (cpack(x, y - y)); + return (CMPLX(x, y - y)); } } @@ -84,6 +84,6 @@ cexp(double complex z) * - x = NaN (spurious inexact exception from y) */ exp_x = exp(x); - return (cpack(exp_x * cos(y), exp_x * sin(y))); + return (CMPLX(exp_x * cos(y), exp_x * sin(y))); } } diff --git a/lib/msun/src/s_cexpf.c b/lib/msun/src/s_cexpf.c index 0e30d08..0138cd1 100644 --- a/lib/msun/src/s_cexpf.c +++ b/lib/msun/src/s_cexpf.c @@ -50,22 +50,22 @@ cexpf(float complex z) /* cexp(x + I 0) = exp(x) + I 0 */ if (hy == 0) - return (cpackf(expf(x), y)); + return (CMPLXF(expf(x), y)); GET_FLOAT_WORD(hx, x); /* cexp(0 + I y) = cos(y) + I sin(y) */ if ((hx & 0x7fffffff) == 0) - return (cpackf(cosf(y), sinf(y))); + return (CMPLXF(cosf(y), sinf(y))); if (hy >= 0x7f800000) { if ((hx & 0x7fffffff) != 0x7f800000) { /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ - return (cpackf(y - y, y - y)); + return (CMPLXF(y - y, y - y)); } else if (hx & 0x80000000) { /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ - return (cpackf(0.0, 0.0)); + return (CMPLXF(0.0, 0.0)); } else { /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ - return (cpackf(x, y - y)); + return (CMPLXF(x, y - y)); } } @@ -84,6 +84,6 @@ cexpf(float complex z) * - x = NaN (spurious inexact exception from y) */ exp_x = expf(x); - return (cpackf(exp_x * cosf(y), exp_x * sinf(y))); + return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y))); } } diff --git a/lib/msun/src/s_conj.c b/lib/msun/src/s_conj.c index 5770c29..d0dd051 100644 --- a/lib/msun/src/s_conj.c +++ b/lib/msun/src/s_conj.c @@ -34,5 +34,5 @@ double complex conj(double complex z) { - return (cpack(creal(z), -cimag(z))); + return (CMPLX(creal(z), -cimag(z))); } diff --git a/lib/msun/src/s_conjf.c b/lib/msun/src/s_conjf.c index b090760..c4cf127 100644 --- a/lib/msun/src/s_conjf.c +++ b/lib/msun/src/s_conjf.c @@ -34,5 +34,5 @@ float complex conjf(float complex z) { - return (cpackf(crealf(z), -cimagf(z))); + return (CMPLXF(crealf(z), -cimagf(z))); } diff --git a/lib/msun/src/s_conjl.c b/lib/msun/src/s_conjl.c index 0e431ef..a4604b6 100644 --- a/lib/msun/src/s_conjl.c +++ b/lib/msun/src/s_conjl.c @@ -34,5 +34,5 @@ long double complex conjl(long double complex z) { - return (cpackl(creall(z), -cimagl(z))); + return (CMPLXL(creall(z), -cimagl(z))); } diff --git a/lib/msun/src/s_cproj.c b/lib/msun/src/s_cproj.c index 8e9404c..2f0db6e 100644 --- a/lib/msun/src/s_cproj.c +++ b/lib/msun/src/s_cproj.c @@ -39,7 +39,7 @@ cproj(double complex z) if (!isinf(creal(z)) && !isinf(cimag(z))) return (z); else - return (cpack(INFINITY, copysign(0.0, cimag(z)))); + return (CMPLX(INFINITY, copysign(0.0, cimag(z)))); } #if LDBL_MANT_DIG == 53 diff --git a/lib/msun/src/s_cprojf.c b/lib/msun/src/s_cprojf.c index 68ea77b..bd27bdc 100644 --- a/lib/msun/src/s_cprojf.c +++ b/lib/msun/src/s_cprojf.c @@ -39,5 +39,5 @@ cprojf(float complex z) if (!isinf(crealf(z)) && !isinf(cimagf(z))) return (z); else - return (cpackf(INFINITY, copysignf(0.0, cimagf(z)))); + return (CMPLXF(INFINITY, copysignf(0.0, cimagf(z)))); } diff --git a/lib/msun/src/s_cprojl.c b/lib/msun/src/s_cprojl.c index 07385bc..5cccc91 100644 --- a/lib/msun/src/s_cprojl.c +++ b/lib/msun/src/s_cprojl.c @@ -39,5 +39,5 @@ cprojl(long double complex z) if (!isinf(creall(z)) && !isinf(cimagl(z))) return (z); else - return (cpackl(INFINITY, copysignl(0.0, cimagl(z)))); + return (CMPLXL(INFINITY, copysignl(0.0, cimagl(z)))); } diff --git a/lib/msun/src/s_csinh.c b/lib/msun/src/s_csinh.c index c192f30..7af3634 100644 --- a/lib/msun/src/s_csinh.c +++ b/lib/msun/src/s_csinh.c @@ -32,6 +32,8 @@ * * Exceptional values are noted in the comments within the source code. * These values and the return value were taken from n1124.pdf. + * The sign of the result for some exceptional values is unspecified but + * must satisfy both sinh(conj(z)) == conj(sinh(z)) and sinh(-z) == -sinh(z). */ #include <sys/cdefs.h> @@ -62,48 +64,45 @@ csinh(double complex z) /* Handle the nearly-non-exceptional cases where x and y are finite. */ if (ix < 0x7ff00000 && iy < 0x7ff00000) { if ((iy | ly) == 0) - return (cpack(sinh(x), y)); - if (ix < 0x40360000) /* small x: normal case */ - return (cpack(sinh(x) * cos(y), cosh(x) * sin(y))); + return (CMPLX(sinh(x), y)); + if (ix < 0x40360000) /* |x| < 22: normal case */ + return (CMPLX(sinh(x) * cos(y), cosh(x) * sin(y))); /* |x| >= 22, so cosh(x) ~= exp(|x|) */ if (ix < 0x40862e42) { /* x < 710: exp(|x|) won't overflow */ h = exp(fabs(x)) * 0.5; - return (cpack(copysign(h, x) * cos(y), h * sin(y))); + return (CMPLX(copysign(h, x) * cos(y), h * sin(y))); } else if (ix < 0x4096bbaa) { /* x < 1455: scale to avoid overflow */ - z = __ldexp_cexp(cpack(fabs(x), y), -1); - return (cpack(creal(z) * copysign(1, x), cimag(z))); + z = __ldexp_cexp(CMPLX(fabs(x), y), -1); + return (CMPLX(creal(z) * copysign(1, x), cimag(z))); } else { /* x >= 1455: the result always overflows */ h = huge * x; - return (cpack(h * cos(y), h * h * sin(y))); + return (CMPLX(h * cos(y), h * h * sin(y))); } } /* - * sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN. - * The sign of 0 in the result is unspecified. Choice = normally - * the same as dNaN. Raise the invalid floating-point exception. + * sinh(+-0 +- I Inf) = +-0 + I dNaN. + * The sign of 0 in the result is unspecified. Choice = same sign + * as the argument. Raise the invalid floating-point exception. * - * sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN). - * The sign of 0 in the result is unspecified. Choice = normally - * the same as d(NaN). + * sinh(+-0 +- I NaN) = +-0 + I d(NaN). + * The sign of 0 in the result is unspecified. Choice = same sign + * as the argument. */ - if ((ix | lx) == 0 && iy >= 0x7ff00000) - return (cpack(copysign(0, x * (y - y)), y - y)); + if ((ix | lx) == 0) /* && iy >= 0x7ff00000 */ + return (CMPLX(x, y - y)); /* * sinh(+-Inf +- I 0) = +-Inf + I +-0. * * sinh(NaN +- I 0) = d(NaN) + I +-0. */ - if ((iy | ly) == 0 && ix >= 0x7ff00000) { - if (((hx & 0xfffff) | lx) == 0) - return (cpack(x, y)); - return (cpack(x, copysign(0, y))); - } + if ((iy | ly) == 0) /* && ix >= 0x7ff00000 */ + return (CMPLX(x + x, y)); /* * sinh(x +- I Inf) = dNaN + I dNaN. @@ -113,45 +112,45 @@ csinh(double complex z) * Optionally raises the invalid floating-point exception for finite * nonzero x. Choice = don't raise (except for signaling NaNs). */ - if (ix < 0x7ff00000 && iy >= 0x7ff00000) - return (cpack(y - y, x * (y - y))); + if (ix < 0x7ff00000) /* && iy >= 0x7ff00000 */ + return (CMPLX(y - y, y - y)); /* * sinh(+-Inf + I NaN) = +-Inf + I d(NaN). - * The sign of Inf in the result is unspecified. Choice = normally - * the same as d(NaN). + * The sign of Inf in the result is unspecified. Choice = same sign + * as the argument. * - * sinh(+-Inf +- I Inf) = +Inf + I dNaN. - * The sign of Inf in the result is unspecified. Choice = always +. - * Raise the invalid floating-point exception. + * sinh(+-Inf +- I Inf) = +-Inf + I dNaN. + * The sign of Inf in the result is unspecified. Choice = same sign + * as the argument. Raise the invalid floating-point exception. * * sinh(+-Inf + I y) = +-Inf cos(y) + I Inf sin(y) */ - if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { + if (ix == 0x7ff00000 && lx == 0) { if (iy >= 0x7ff00000) - return (cpack(x * x, x * (y - y))); - return (cpack(x * cos(y), INFINITY * sin(y))); + return (CMPLX(x, y - y)); + return (CMPLX(x * cos(y), INFINITY * sin(y))); } /* - * sinh(NaN + I NaN) = d(NaN) + I d(NaN). + * sinh(NaN1 + I NaN2) = d(NaN1, NaN2) + I d(NaN1, NaN2). * - * sinh(NaN +- I Inf) = d(NaN) + I d(NaN). + * sinh(NaN +- I Inf) = d(NaN, dNaN) + I d(NaN, dNaN). * Optionally raises the invalid floating-point exception. * Choice = raise. * - * sinh(NaN + I y) = d(NaN) + I d(NaN). + * sinh(NaN + I y) = d(NaN) + I d(NaN). * Optionally raises the invalid floating-point exception for finite * nonzero y. Choice = don't raise (except for signaling NaNs). */ - return (cpack((x * x) * (y - y), (x + x) * (y - y))); + return (CMPLX((x + x) * (y - y), (x * x) * (y - y))); } double complex csin(double complex z) { - /* csin(z) = -I * csinh(I * z) */ - z = csinh(cpack(-cimag(z), creal(z))); - return (cpack(cimag(z), -creal(z))); + /* csin(z) = -I * csinh(I * z) = I * conj(csinh(I * conj(z))). */ + z = csinh(CMPLX(cimag(z), creal(z))); + return (CMPLX(cimag(z), creal(z))); } diff --git a/lib/msun/src/s_csinhf.c b/lib/msun/src/s_csinhf.c index c523125..c9d9d3c 100644 --- a/lib/msun/src/s_csinhf.c +++ b/lib/msun/src/s_csinhf.c @@ -25,7 +25,7 @@ */ /* - * Hyperbolic sine of a complex argument z. See s_csinh.c for details. + * Float version of csinh(). See s_csinh.c for details. */ #include <sys/cdefs.h> @@ -55,51 +55,48 @@ csinhf(float complex z) if (ix < 0x7f800000 && iy < 0x7f800000) { if (iy == 0) - return (cpackf(sinhf(x), y)); - if (ix < 0x41100000) /* small x: normal case */ - return (cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y))); + return (CMPLXF(sinhf(x), y)); + if (ix < 0x41100000) /* |x| < 9: normal case */ + return (CMPLXF(sinhf(x) * cosf(y), coshf(x) * sinf(y))); /* |x| >= 9, so cosh(x) ~= exp(|x|) */ if (ix < 0x42b17218) { /* x < 88.7: expf(|x|) won't overflow */ - h = expf(fabsf(x)) * 0.5f; - return (cpackf(copysignf(h, x) * cosf(y), h * sinf(y))); + h = expf(fabsf(x)) * 0.5F; + return (CMPLXF(copysignf(h, x) * cosf(y), h * sinf(y))); } else if (ix < 0x4340b1e7) { /* x < 192.7: scale to avoid overflow */ - z = __ldexp_cexpf(cpackf(fabsf(x), y), -1); - return (cpackf(crealf(z) * copysignf(1, x), cimagf(z))); + z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1); + return (CMPLXF(crealf(z) * copysignf(1, x), cimagf(z))); } else { /* x >= 192.7: the result always overflows */ h = huge * x; - return (cpackf(h * cosf(y), h * h * sinf(y))); + return (CMPLXF(h * cosf(y), h * h * sinf(y))); } } - if (ix == 0 && iy >= 0x7f800000) - return (cpackf(copysignf(0, x * (y - y)), y - y)); + if (ix == 0) /* && iy >= 0x7f800000 */ + return (CMPLXF(x, y - y)); - if (iy == 0 && ix >= 0x7f800000) { - if ((hx & 0x7fffff) == 0) - return (cpackf(x, y)); - return (cpackf(x, copysignf(0, y))); - } + if (iy == 0) /* && ix >= 0x7f800000 */ + return (CMPLXF(x + x, y)); - if (ix < 0x7f800000 && iy >= 0x7f800000) - return (cpackf(y - y, x * (y - y))); + if (ix < 0x7f800000) /* && iy >= 0x7f800000 */ + return (CMPLXF(y - y, y - y)); - if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) { + if (ix == 0x7f800000) { if (iy >= 0x7f800000) - return (cpackf(x * x, x * (y - y))); - return (cpackf(x * cosf(y), INFINITY * sinf(y))); + return (CMPLXF(x, y - y)); + return (CMPLXF(x * cosf(y), INFINITY * sinf(y))); } - return (cpackf((x * x) * (y - y), (x + x) * (y - y))); + return (CMPLXF((x + x) * (y - y), (x * x) * (y - y))); } float complex csinf(float complex z) { - z = csinhf(cpackf(-cimagf(z), crealf(z))); - return (cpackf(cimagf(z), -crealf(z))); + z = csinhf(CMPLXF(cimagf(z), crealf(z))); + return (CMPLXF(cimagf(z), crealf(z))); } diff --git a/lib/msun/src/s_csqrt.c b/lib/msun/src/s_csqrt.c index 18a7ae3..75e12d8 100644 --- a/lib/msun/src/s_csqrt.c +++ b/lib/msun/src/s_csqrt.c @@ -58,12 +58,12 @@ csqrt(double complex z) /* Handle special cases. */ if (z == 0) - return (cpack(0, b)); + return (CMPLX(0, b)); if (isinf(b)) - return (cpack(INFINITY, b)); + return (CMPLX(INFINITY, b)); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ - return (cpack(a, t)); /* return NaN + NaN i */ + return (CMPLX(a, t)); /* return NaN + NaN i */ } if (isinf(a)) { /* @@ -73,9 +73,9 @@ csqrt(double complex z) * csqrt(-inf + y i) = 0 + inf i */ if (signbit(a)) - return (cpack(fabs(b - b), copysign(a, b))); + return (CMPLX(fabs(b - b), copysign(a, b))); else - return (cpack(a, copysign(b - b, b))); + return (CMPLX(a, copysign(b - b, b))); } /* * The remaining special case (b is NaN) is handled just fine by @@ -94,10 +94,10 @@ csqrt(double complex z) /* Algorithm 312, CACM vol 10, Oct 1967. */ if (a >= 0) { t = sqrt((a + hypot(a, b)) * 0.5); - result = cpack(t, b / (2 * t)); + result = CMPLX(t, b / (2 * t)); } else { t = sqrt((-a + hypot(a, b)) * 0.5); - result = cpack(fabs(b) / (2 * t), copysign(t, b)); + result = CMPLX(fabs(b) / (2 * t), copysign(t, b)); } /* Rescale. */ diff --git a/lib/msun/src/s_csqrtf.c b/lib/msun/src/s_csqrtf.c index da7fe18..7ee513f 100644 --- a/lib/msun/src/s_csqrtf.c +++ b/lib/msun/src/s_csqrtf.c @@ -49,12 +49,12 @@ csqrtf(float complex z) /* Handle special cases. */ if (z == 0) - return (cpackf(0, b)); + return (CMPLXF(0, b)); if (isinf(b)) - return (cpackf(INFINITY, b)); + return (CMPLXF(INFINITY, b)); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ - return (cpackf(a, t)); /* return NaN + NaN i */ + return (CMPLXF(a, t)); /* return NaN + NaN i */ } if (isinf(a)) { /* @@ -64,9 +64,9 @@ csqrtf(float complex z) * csqrtf(-inf + y i) = 0 + inf i */ if (signbit(a)) - return (cpackf(fabsf(b - b), copysignf(a, b))); + return (CMPLXF(fabsf(b - b), copysignf(a, b))); else - return (cpackf(a, copysignf(b - b, b))); + return (CMPLXF(a, copysignf(b - b, b))); } /* * The remaining special case (b is NaN) is handled just fine by @@ -80,9 +80,9 @@ csqrtf(float complex z) */ if (a >= 0) { t = sqrt((a + hypot(a, b)) * 0.5); - return (cpackf(t, b / (2.0 * t))); + return (CMPLXF(t, b / (2.0 * t))); } else { t = sqrt((-a + hypot(a, b)) * 0.5); - return (cpackf(fabsf(b) / (2.0 * t), copysignf(t, b))); + return (CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b))); } } diff --git a/lib/msun/src/s_csqrtl.c b/lib/msun/src/s_csqrtl.c index dd18e1e..ea2ef27 100644 --- a/lib/msun/src/s_csqrtl.c +++ b/lib/msun/src/s_csqrtl.c @@ -58,12 +58,12 @@ csqrtl(long double complex z) /* Handle special cases. */ if (z == 0) - return (cpackl(0, b)); + return (CMPLXL(0, b)); if (isinf(b)) - return (cpackl(INFINITY, b)); + return (CMPLXL(INFINITY, b)); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ - return (cpackl(a, t)); /* return NaN + NaN i */ + return (CMPLXL(a, t)); /* return NaN + NaN i */ } if (isinf(a)) { /* @@ -73,9 +73,9 @@ csqrtl(long double complex z) * csqrt(-inf + y i) = 0 + inf i */ if (signbit(a)) - return (cpackl(fabsl(b - b), copysignl(a, b))); + return (CMPLXL(fabsl(b - b), copysignl(a, b))); else - return (cpackl(a, copysignl(b - b, b))); + return (CMPLXL(a, copysignl(b - b, b))); } /* * The remaining special case (b is NaN) is handled just fine by @@ -94,10 +94,10 @@ csqrtl(long double complex z) /* Algorithm 312, CACM vol 10, Oct 1967. */ if (a >= 0) { t = sqrtl((a + hypotl(a, b)) * 0.5); - result = cpackl(t, b / (2 * t)); + result = CMPLXL(t, b / (2 * t)); } else { t = sqrtl((-a + hypotl(a, b)) * 0.5); - result = cpackl(fabsl(b) / (2 * t), copysignl(t, b)); + result = CMPLXL(fabsl(b) / (2 * t), copysignl(t, b)); } /* Rescale. */ diff --git a/lib/msun/src/s_ctanh.c b/lib/msun/src/s_ctanh.c index d427e28..f5b9bdd 100644 --- a/lib/msun/src/s_ctanh.c +++ b/lib/msun/src/s_ctanh.c @@ -25,7 +25,7 @@ */ /* - * Hyperbolic tangent of a complex argument z = x + i y. + * Hyperbolic tangent of a complex argument z = x + I y. * * The algorithm is from: * @@ -44,15 +44,15 @@ * * tanh(z) = sinh(z) / cosh(z) * - * sinh(x) cos(y) + i cosh(x) sin(y) + * sinh(x) cos(y) + I cosh(x) sin(y) * = --------------------------------- - * cosh(x) cos(y) + i sinh(x) sin(y) + * cosh(x) cos(y) + I sinh(x) sin(y) * - * cosh(x) sinh(x) / cos^2(y) + i tan(y) + * cosh(x) sinh(x) / cos^2(y) + I tan(y) * = ------------------------------------- * 1 + sinh^2(x) / cos^2(y) * - * beta rho s + i t + * beta rho s + I t * = ---------------- * 1 + beta s^2 * @@ -85,16 +85,16 @@ ctanh(double complex z) ix = hx & 0x7fffffff; /* - * ctanh(NaN + i 0) = NaN + i 0 + * ctanh(NaN +- I 0) = d(NaN) +- I 0 * - * ctanh(NaN + i y) = NaN + i NaN for y != 0 + * ctanh(NaN + I y) = d(NaN,y) + I d(NaN,y) for y != 0 * * The imaginary part has the sign of x*sin(2*y), but there's no * special effort to get this right. * - * ctanh(+-Inf +- i Inf) = +-1 +- 0 + * ctanh(+-Inf +- I Inf) = +-1 +- I 0 * - * ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite + * ctanh(+-Inf + I y) = +-1 + I 0 sin(2y) for y finite * * The imaginary part of the sign is unspecified. This special * case is only needed to avoid a spurious invalid exception when @@ -102,26 +102,27 @@ ctanh(double complex z) */ if (ix >= 0x7ff00000) { if ((ix & 0xfffff) | lx) /* x is NaN */ - return (cpack(x, (y == 0 ? y : x * y))); + return (CMPLX((x + 0) * (y + 0), + y == 0 ? y : (x + 0) * (y + 0))); SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */ - return (cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); + return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); } /* - * ctanh(x + i NAN) = NaN + i NaN - * ctanh(x +- i Inf) = NaN + i NaN + * ctanh(x + I NaN) = d(NaN) + I d(NaN) + * ctanh(x +- I Inf) = dNaN + I dNaN */ if (!isfinite(y)) - return (cpack(y - y, y - y)); + return (CMPLX(y - y, y - y)); /* - * ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the + * ctanh(+-huge +- I y) ~= +-1 +- I 2sin(2y)/exp(2x), using the * approximation sinh^2(huge) ~= exp(2*huge) / 4. * We use a modified formula to avoid spurious overflow. */ - if (ix >= 0x40360000) { /* x >= 22 */ + if (ix >= 0x40360000) { /* |x| >= 22 */ double exp_mx = exp(-fabs(x)); - return (cpack(copysign(1, x), + return (CMPLX(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx)); } @@ -131,14 +132,14 @@ ctanh(double complex z) s = sinh(x); rho = sqrt(1 + s * s); /* = cosh(x) */ denom = 1 + beta * s * s; - return (cpack((beta * rho * s) / denom, t / denom)); + return (CMPLX((beta * rho * s) / denom, t / denom)); } double complex ctan(double complex z) { - /* ctan(z) = -I * ctanh(I * z) */ - z = ctanh(cpack(-cimag(z), creal(z))); - return (cpack(cimag(z), -creal(z))); + /* ctan(z) = -I * ctanh(I * z) = I * conj(ctanh(I * conj(z))) */ + z = ctanh(CMPLX(cimag(z), creal(z))); + return (CMPLX(cimag(z), creal(z))); } diff --git a/lib/msun/src/s_ctanhf.c b/lib/msun/src/s_ctanhf.c index 4be28d8..520bf77 100644 --- a/lib/msun/src/s_ctanhf.c +++ b/lib/msun/src/s_ctanhf.c @@ -51,18 +51,19 @@ ctanhf(float complex z) if (ix >= 0x7f800000) { if (ix & 0x7fffff) - return (cpackf(x, (y == 0 ? y : x * y))); + return (CMPLXF((x + 0) * (y + 0), + y == 0 ? y : (x + 0) * (y + 0))); SET_FLOAT_WORD(x, hx - 0x40000000); - return (cpackf(x, + return (CMPLXF(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y)))); } if (!isfinite(y)) - return (cpackf(y - y, y - y)); + return (CMPLXF(y - y, y - y)); - if (ix >= 0x41300000) { /* x >= 11 */ + if (ix >= 0x41300000) { /* |x| >= 11 */ float exp_mx = expf(-fabsf(x)); - return (cpackf(copysignf(1, x), + return (CMPLXF(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx)); } @@ -71,14 +72,14 @@ ctanhf(float complex z) s = sinhf(x); rho = sqrtf(1 + s * s); denom = 1 + beta * s * s; - return (cpackf((beta * rho * s) / denom, t / denom)); + return (CMPLXF((beta * rho * s) / denom, t / denom)); } float complex ctanf(float complex z) { - z = ctanhf(cpackf(-cimagf(z), crealf(z))); - return (cpackf(cimagf(z), -crealf(z))); + z = ctanhf(CMPLXF(cimagf(z), crealf(z))); + return (CMPLXF(cimagf(z), crealf(z))); } diff --git a/lib/msun/src/s_scalbln.c b/lib/msun/src/s_scalbln.c index 2aa6900..5e4e42e 100644 --- a/lib/msun/src/s_scalbln.c +++ b/lib/msun/src/s_scalbln.c @@ -27,32 +27,28 @@ #include <sys/cdefs.h> __FBSDID("$FreeBSD$"); -#include <limits.h> #include <math.h> +#define NMAX 65536 +#define NMIN -65536 + double -scalbln (double x, long n) +scalbln(double x, long n) { - int in; - in = (n > INT_MAX) ? INT_MAX : (n < INT_MIN) ? INT_MIN : n; - return (scalbn(x, in)); + return (scalbn(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n)); } float -scalblnf (float x, long n) +scalblnf(float x, long n) { - int in; - in = (n > INT_MAX) ? INT_MAX : (n < INT_MIN) ? INT_MIN : n; - return (scalbnf(x, in)); + return (scalbnf(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n)); } long double -scalblnl (long double x, long n) +scalblnl(long double x, long n) { - int in; - in = (n > INT_MAX) ? INT_MAX : (n < INT_MIN) ? INT_MIN : n; - return (scalbnl(x, in)); + return (scalbnl(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n)); } |
