Date: Thu, 25 Jun 2015 13:01:10 +0000 (UTC) From: Tijl Coosemans <tijl@FreeBSD.org> To: src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-stable@freebsd.org, svn-src-stable-10@freebsd.org Subject: svn commit: r284810 - in stable/10/lib/msun: . ld128 ld80 man src Message-ID: <201506251301.t5PD1AbQ066585@svn.freebsd.org>
next in thread | raw e-mail | index | archive | help
Author: tijl Date: Thu Jun 25 13:01:10 2015 New Revision: 284810 URL: https://svnweb.freebsd.org/changeset/base/284810 Log: MFC r271651, r271719, r272138, r272457, r272845, r275476, r275518, r275614, r275819, r276176, r278154, r278160, r278339, r279127, r279240, r279491, r279493, r279856, r283032, r284423, r284426, r284427, r284428 Merge changes to libm from the past 9 months. This includes improvements to the Bessel functions and adds the C99 function lgammal. Added: stable/10/lib/msun/ld128/e_lgammal_r.c - copied, changed from r271651, head/lib/msun/ld128/e_lgammal_r.c stable/10/lib/msun/ld80/e_lgammal_r.c - copied, changed from r271651, head/lib/msun/ld80/e_lgammal_r.c stable/10/lib/msun/src/e_lgammal.c - copied unchanged from r271651, head/lib/msun/src/e_lgammal.c Modified: stable/10/lib/msun/Makefile stable/10/lib/msun/Symbol.map stable/10/lib/msun/ld128/k_expl.h stable/10/lib/msun/ld80/k_expl.h stable/10/lib/msun/man/j0.3 stable/10/lib/msun/man/lgamma.3 stable/10/lib/msun/src/catrig.c stable/10/lib/msun/src/catrigf.c stable/10/lib/msun/src/e_j0.c stable/10/lib/msun/src/e_j0f.c stable/10/lib/msun/src/e_j1.c stable/10/lib/msun/src/e_j1f.c stable/10/lib/msun/src/e_jn.c stable/10/lib/msun/src/e_jnf.c stable/10/lib/msun/src/e_lgamma.c stable/10/lib/msun/src/e_lgamma_r.c stable/10/lib/msun/src/e_lgammaf_r.c stable/10/lib/msun/src/imprecise.c stable/10/lib/msun/src/k_exp.c stable/10/lib/msun/src/k_expf.c stable/10/lib/msun/src/math.h stable/10/lib/msun/src/math_private.h stable/10/lib/msun/src/s_ccosh.c stable/10/lib/msun/src/s_ccoshf.c stable/10/lib/msun/src/s_cexp.c stable/10/lib/msun/src/s_cexpf.c stable/10/lib/msun/src/s_conj.c stable/10/lib/msun/src/s_conjf.c stable/10/lib/msun/src/s_conjl.c stable/10/lib/msun/src/s_cproj.c stable/10/lib/msun/src/s_cprojf.c stable/10/lib/msun/src/s_cprojl.c stable/10/lib/msun/src/s_csinh.c stable/10/lib/msun/src/s_csinhf.c stable/10/lib/msun/src/s_csqrt.c stable/10/lib/msun/src/s_csqrtf.c stable/10/lib/msun/src/s_csqrtl.c stable/10/lib/msun/src/s_ctanh.c stable/10/lib/msun/src/s_ctanhf.c stable/10/lib/msun/src/s_scalbln.c Directory Properties: stable/10/ (props changed) Modified: stable/10/lib/msun/Makefile ============================================================================== --- stable/10/lib/msun/Makefile Thu Jun 25 11:32:41 2015 (r284809) +++ stable/10/lib/msun/Makefile Thu Jun 25 13:01:10 2015 (r284810) @@ -98,6 +98,7 @@ COMMON_SRCS+= s_copysignl.c s_fabsl.c s_ # If long double != double use these; otherwise, we alias the double versions. COMMON_SRCS+= e_acoshl.c e_acosl.c e_asinl.c e_atan2l.c e_atanhl.c \ e_coshl.c e_fmodl.c e_hypotl.c \ + e_lgammal.c e_lgammal_r.c \ e_remainderl.c e_sinhl.c e_sqrtl.c \ invtrig.c k_cosl.c k_sinl.c k_tanl.c \ s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \ @@ -188,7 +189,8 @@ MLINKS+=ilogb.3 ilogbf.3 ilogb.3 ilogbl. ilogb.3 logb.3 ilogb.3 logbf.3 ilogb.3 logbl.3 MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3 MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3 -MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 lgamma.3 lgammaf.3 \ +MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 \ + lgamma.3 lgammaf.3 lgamma.3 lgammal.3 \ lgamma.3 tgamma.3 lgamma.3 tgammaf.3 MLINKS+=log.3 log10.3 log.3 log10f.3 log.3 log10l.3 \ log.3 log1p.3 log.3 log1pf.3 log.3 log1pl.3 \ Modified: stable/10/lib/msun/Symbol.map ============================================================================== --- stable/10/lib/msun/Symbol.map Thu Jun 25 11:32:41 2015 (r284809) +++ stable/10/lib/msun/Symbol.map Thu Jun 25 13:01:10 2015 (r284810) @@ -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; +}; Copied and modified: stable/10/lib/msun/ld128/e_lgammal_r.c (from r271651, head/lib/msun/ld128/e_lgammal_r.c) ============================================================================== --- head/lib/msun/ld128/e_lgammal_r.c Mon Sep 15 23:21:57 2014 (r271651, copy source) +++ stable/10/lib/msun/ld128/e_lgammal_r.c Thu Jun 25 13:01:10 2015 (r284810) @@ -206,13 +206,13 @@ sin_pil(long double x) n--; } n &= 7; - y = y - z + n * 0.25L; + y = y - z + n * 0.25; switch (n) { case 0: y = __kernel_sinl(pi*y,zero,0); break; case 1: case 2: y = __kernel_cosl(pi*(0.5-y),zero); break; - case 3: + case 3: case 4: y = __kernel_sinl(pi*(one-y),zero,0); break; case 5: case 6: y = -__kernel_cosl(pi*(y-1.5),zero); break; @@ -221,37 +221,33 @@ sin_pil(long double x) return -y; } - long double lgammal_r(long double x, int *signgamp) { long double nadj,p,p1,p2,p3,q,r,t,w,y,z; uint64_t llx,lx; int i; - uint16_t hx; - - EXTRACT_LDBL128_WORDS(hx, lx, llx, x); + uint16_t hx,ix; - if((hx & 0x7fff) == 0x7fff) { /* erfl(nan)=nan */ - i = (hx>>15)<<1; - return (1-i)+one/x; /* erfl(+-inf)=+-1 */ - } + EXTRACT_LDBL128_WORDS(hx,lx,llx,x); - /* purge off +-inf, NaN, +-0, tiny and negative arguments */ + /* purge +-Inf and NaNs */ *signgamp = 1; - if((hx & 0x7fff) == 0x7fff) /* x is +-Inf or NaN */ - return x*x; - if((hx==0||hx==0x8000)&&lx==0) return one/vzero; - - /* purge off tiny and negative arguments */ - if(fabsl(x)<0x1p-119L) { - if(hx&0x8000) { - *signgamp = -1; - return -logl(-x); - } else return -logl(x); + 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) { - if(fabsl(x)>=0x1p112) + *signgamp = 1; + if(ix>=0x3fff+112) /* |x|>=2**(p-1), must be -integer */ return one/vzero; t = sin_pil(x); if(t==zero) return one/vzero; @@ -260,17 +256,19 @@ lgammal_r(long double x, int *signgamp) x = -x; } - if(x == 1 || x ==2) r = 0; - else if(x<2) { - if(x<=0.8999996185302734) { + /* 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>=0.7315998077392578) {y = 1-x; i= 0;} - else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;} + if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;} + else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;} else {y = x; i=2;} } else { - r = 0; - if(x>=1.7316312789916992) {y=2-x;i=0;} - else if(x>=1.2316322326660156) {y=x-tc;i=1;} + r = 0; + if(x>=1.7316312789916992e+00) {y=2-x;i=0;} + else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;} else {y=x-1;i=2;} } switch(i) { @@ -281,23 +279,24 @@ lgammal_r(long double x, int *signgamp) p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*(a13+z*(a15+ z*(a17+z*(a19+z*(a21+z*a23))))))))))); p = y*p1+p2; - r += (p-y/2); break; + r += p-y/2; break; case 1: p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+ y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+ y*(t17+y*(t18+y*(t19+y*(t20+y*(t21+y*(t22+y*(t23+ y*(t24+y*(t25+y*(t26+y*(t27+y*(t28+y*(t29+y*(t30+ y*(t31+y*t32)))))))))))))))))))))))))))))); - r += (tf + p); break; + r += tf + p; break; case 2: p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*(u7+ y*(u8+y*(u9+y*u10)))))))))); p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7+ y*(v8+y*(v9+y*(v10+y*v11)))))))))); - r += (-y/2 + p1/p2); + r += p1/p2-y/2; } } - else if(x<8) { + /* x < 8.0 */ + else if(ix<0x4002) { i = x; y = x-i; p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*(s6+y*(s7+y*(s8+ @@ -314,7 +313,8 @@ lgammal_r(long double x, int *signgamp) case 3: z *= (y+2); /* FALLTHRU */ r += logl(z); break; } - } else if (x < 0x1p119L) { + /* 8.0 <= x < 2**(p+3) */ + } else if (ix<0x3fff+116) { t = logl(x); z = one/x; y = z*z; @@ -322,6 +322,7 @@ lgammal_r(long double x, int *signgamp) 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; Modified: stable/10/lib/msun/ld128/k_expl.h ============================================================================== --- stable/10/lib/msun/ld128/k_expl.h Thu Jun 25 11:32:41 2015 (r284809) +++ stable/10/lib/msun/ld128/k_expl.h Thu Jun 25 13:01:10 2015 (r284810) @@ -322,7 +322,7 @@ __ldexp_cexpl(long double complex z, int 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 */ Copied and modified: stable/10/lib/msun/ld80/e_lgammal_r.c (from r271651, head/lib/msun/ld80/e_lgammal_r.c) ============================================================================== --- head/lib/msun/ld80/e_lgammal_r.c Mon Sep 15 23:21:57 2014 (r271651, copy source) +++ stable/10/lib/msun/ld80/e_lgammal_r.c Thu Jun 25 13:01:10 2015 (r284810) @@ -14,12 +14,11 @@ __FBSDID("$FreeBSD$"); /* - * See s_lgamma_r.c for complete comments. + * See e_lgamma_r.c for complete comments. * * Converted to long double by Steven G. Kargl. */ -#include <float.h> #ifdef __i386__ #include <ieeefp.h> #endif @@ -219,8 +218,8 @@ sin_pil(long double x) y = -x; - vz = y+0x1p63L; - z = vz-0x1p63L; + vz = y+0x1p63; + z = vz-0x1p63; if (z == y) return zero; @@ -243,7 +242,7 @@ sin_pil(long double x) 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; } @@ -253,27 +252,29 @@ lgammal_r(long double x, int *signgamp) long double nadj,p,p1,p2,p3,q,r,t,w,y,z; uint64_t lx; int i; - uint16_t hx; + uint16_t hx,ix; EXTRACT_LDBL80_WORDS(hx,lx,x); - /* purge off +-inf, NaN, +-0 */ + /* purge +-Inf and NaNs */ *signgamp = 1; - if((hx & 0x7fff) == 0x7fff) /* x is +-Inf or NaN */ - return x*x; - if((hx==0||hx==0x8000)&&lx==0) return one/vzero; + ix = hx&0x7fff; + if(ix==0x7fff) return x*x; ENTERI(); - /* purge off tiny and negative arguments */ - if(fabsl(x)<0x1p-70L) { /* |x|<2**-70, return -log(|x|) */ - if(hx&0x8000) { - *signgamp = -1; - RETURNI(-logl(-x)); - } else RETURNI(-logl(x)); + /* purge +-0 and tiny arguments */ + *signgamp = 1-2*(hx>>15); + if(ix<0x3fff-67) { /* |x|<2**-(p+3), return -log(|x|) */ + if((ix|lx)==0) + RETURNI(one/vzero); + RETURNI(-logl(fabsl(x))); } + + /* purge negative integers and start evaluation for other x < 0 */ if(hx&0x8000) { - if(fabsl(x)>=0x1p63) /* |x|>=2**(p-1), must be -integer */ + *signgamp = 1; + if(ix>=0x3fff+63) /* |x|>=2**(p-1), must be -integer */ RETURNI(one/vzero); t = sin_pil(x); if(t==zero) RETURNI(one/vzero); /* -integer */ @@ -282,19 +283,30 @@ lgammal_r(long double x, int *signgamp) x = -x; } - /* purge off 1 and 2 */ - if(x == 1 || x == 2) r = 0; + /* purge 1 and 2 */ + if((ix==0x3fff || ix==0x4000) && lx==0x8000000000000000ULL) r = 0; /* for x < 2.0 */ - else if(x<2) { - if(x<=0.8999996185302734) { /* lgamma(x) = lgamma(x+1)-log(x) */ - r = - logl(x); - if(x>=0.7315998077392578) {y = 1-x; i= 0;} - else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;} - else {y = x; i=2;} + else if(ix<0x4000) { + /* + * XXX Supposedly, one can use the following information to replace the + * XXX FP rational expressions. A similar approach is appropriate + * XXX for ld128, but one (may need?) needs to consider llx, too. + * + * 8.9999961853027344e-01 3ffe e666600000000000 + * 7.3159980773925781e-01 3ffe bb4a200000000000 + * 2.3163998126983643e-01 3ffc ed33080000000000 + * 1.7316312789916992e+00 3fff dda6180000000000 + * 1.2316322326660156e+00 3fff 9da6200000000000 + */ + if(x<8.9999961853027344e-01) { + r = -logl(x); + if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;} + else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;} + else {y = x; i=2;} } else { r = 0; - if(x>=1.7316312789916992) {y=2-x;i=0;} - else if(x>=1.2316322326660156) {y=x-tc;i=1;} + if(x>=1.7316312789916992e+00) {y=2-x;i=0;} + else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;} else {y=x-1;i=2;} } switch(i) { @@ -303,19 +315,20 @@ lgammal_r(long double x, int *signgamp) p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*a12))))); p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*a13)))))); p = y*p1+p2; - r += (p-y/2); break; + r += p-y/2; break; case 1: p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+ y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+ y*(t17+y*t18)))))))))))))))); - r += (tf + p); break; + r += tf + p; break; case 2: p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*u6)))))); p2 = 1+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*v6))))); - r += (-y/2 + p1/p2); + r += p1/p2-y/2; } } - else if(x<8) { + /* x < 8.0 */ + else if(ix<0x4002) { i = x; y = x-i; p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); @@ -330,15 +343,15 @@ lgammal_r(long double x, int *signgamp) case 3: z *= (y+2); /* FALLTHRU */ r += logl(z); break; } - /* 8.0 <= x < 2**70 */ - } else if (x < 0x1p70L) { + /* 8.0 <= x < 2**(p+3) */ + } else if (ix<0x3fff+67) { t = logl(x); z = one/x; y = z*z; w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*(w6+y*(w7+y*w8))))))); r = (x-half)*(t-one)+w; + /* 2**(p+3) <= x <= inf */ } else - /* 2**70 <= x <= inf */ r = x*(logl(x)-1); if(hx&0x8000) r = nadj - r; RETURNI(r); Modified: stable/10/lib/msun/ld80/k_expl.h ============================================================================== --- stable/10/lib/msun/ld80/k_expl.h Thu Jun 25 11:32:41 2015 (r284809) +++ stable/10/lib/msun/ld80/k_expl.h Thu Jun 25 13:01:10 2015 (r284810) @@ -299,7 +299,7 @@ __ldexp_cexpl(long double complex z, int 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 */ Modified: stable/10/lib/msun/man/j0.3 ============================================================================== --- stable/10/lib/msun/man/j0.3 Thu Jun 25 11:32:41 2015 (r284809) +++ stable/10/lib/msun/man/j0.3 Thu Jun 25 13:01:10 2015 (r284810) @@ -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 Modified: stable/10/lib/msun/man/lgamma.3 ============================================================================== --- stable/10/lib/msun/man/lgamma.3 Thu Jun 25 11:32:41 2015 (r284809) +++ stable/10/lib/msun/man/lgamma.3 Thu Jun 25 13:01:10 2015 (r284810) @@ -28,7 +28,7 @@ .\" from: @(#)lgamma.3 6.6 (Berkeley) 12/3/92 .\" $FreeBSD$ .\" -.Dd January 14, 2005 +.Dd September 12, 2014 .Dt LGAMMA 3 .Os .Sh NAME @@ -36,6 +36,8 @@ .Nm lgamma_r , .Nm lgammaf , .Nm lgammaf_r , +.Nm lgammal , +.Nm lgammal_r , .Nm gamma , .Nm gamma_r , .Nm gammaf , @@ -58,6 +60,10 @@ .Fn lgammaf "float x" .Ft float .Fn lgammaf_r "float x" "int *signgamp" +.Ft "long double" +.Fn lgammal "long double x" +.Ft "long double" +.Fn lgammal_r "long double x" "int *signgamp" .Ft double .Fn gamma "double x" .Ft double @@ -66,14 +72,15 @@ .Fn gammaf "float x" .Ft float .Fn gammaf_r "float x" "int *signgamp" -.Ft double +.Ft "long double" .Fn tgamma "double x" .Ft float .Fn tgammaf "float x" .Sh DESCRIPTION -.Fn lgamma x +.Fn lgamma x , +.Fn lgammaf x , and -.Fn lgammaf x +.Fn lgammal x .if t \{\ return ln\||\(*G(x)| where .Bd -unfilled -offset indent @@ -87,13 +94,15 @@ The external integer .Fa signgam returns the sign of \(*G(x). .Pp -.Fn lgamma_r x signgamp +.Fn lgamma_r x signgamp , +.Fn lgammaf_r x signgamp , and -.Fn lgammaf_r x signgamp +.Fn lgammal_r x signgamp provide the same functionality as -.Fn lgamma x +.Fn lgamma x , +.Fn lgammaf x , and -.Fn lgammaf x +.Fn lgammal x , but the caller must provide an integer to store the sign of \(*G(x). .Pp The @@ -115,6 +124,7 @@ are deprecated aliases for and .Fn lgammaf_r , respectively. + .Sh IDIOSYNCRASIES Do not use the expression .Dq Li signgam\(**exp(lgamma(x)) @@ -139,14 +149,18 @@ Exponentiation of will lose up to 10 significant bits. .Sh RETURN VALUES .Fn gamma , -.Fn gamma_r , .Fn gammaf , +.Fn gammal , +.Fn gamma_r , .Fn gammaf_r , +.Fn gammal_r , .Fn lgamma , -.Fn lgamma_r , .Fn lgammaf , +.Fn lgammal , +.Fn lgamma_r , +.Fn lgammaf_r , and -.Fn lgammaf_r +.Fn lgammal_r return appropriate values unless an argument is out of range. Overflow will occur for sufficiently large positive values, and non-positive integers. @@ -159,6 +173,7 @@ will underflow. The .Fn lgamma , .Fn lgammaf , +.Fn lgammal , .Fn tgamma , and .Fn tgammaf Modified: stable/10/lib/msun/src/catrig.c ============================================================================== --- stable/10/lib/msun/src/catrig.c Thu Jun 25 11:32:41 2015 (r284809) +++ stable/10/lib/msun/src/catrig.c Thu Jun 25 13:01:10 2015 (r284810) @@ -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))); } Modified: stable/10/lib/msun/src/catrigf.c ============================================================================== --- stable/10/lib/msun/src/catrigf.c Thu Jun 25 11:32:41 2015 (r284809) +++ stable/10/lib/msun/src/catrigf.c Thu Jun 25 13:01:10 2015 (r284810) @@ -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))); } Modified: stable/10/lib/msun/src/e_j0.c ============================================================================== --- stable/10/lib/msun/src/e_j0.c Thu Jun 25 11:32:41 2015 (r284809) +++ stable/10/lib/msun/src/e_j0.c Thu Jun 25 13:01:10 2015 (r284810) @@ -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 */ }; *** DIFF OUTPUT TRUNCATED AT 1000 LINES ***
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201506251301.t5PD1AbQ066585>