diff options
author | kargl <kargl@FreeBSD.org> | 2014-09-15 23:21:57 +0000 |
---|---|---|
committer | kargl <kargl@FreeBSD.org> | 2014-09-15 23:21:57 +0000 |
commit | 80b8071609d328fad184f145d51bb9d0fd790d9c (patch) | |
tree | 703e4d5f343e4dea078b2783f81a37b7d1b48b73 /lib/msun/src | |
parent | 971bc8fa84832ad2708b258da99f046bdc032055 (diff) | |
download | FreeBSD-src-80b8071609d328fad184f145d51bb9d0fd790d9c.zip FreeBSD-src-80b8071609d328fad184f145d51bb9d0fd790d9c.tar.gz |
* Makefile:
. Hook e_lgammal[_r].c to the build.
. Create man page links for lgammal[-r].3.
* Symbol.map:
. Sort lgammal to its rightful place.
. Add FBSD_1.4 section for the new lgamal_r symbol.
* ld128/e_lgammal_r.c:
. 128-bit implementataion of lgammal_r().
* ld80/e_lgammal_r.c:
. Intel 80-bit format implementation of lgammal_r().
* src/e_lgamma.c:
. Expose lgammal as a weak reference to lgamma for platforms
where long double is mapped to double.
* src/e_lgamma_r.c:
. Use integer literal constants instead of real literal constants.
Let compiler(s) do the job of conversion to the appropriate type.
. Expose lgammal_r as a weak reference to lgamma_r for platforms
where long double is mapped to double.
* src/e_lgammaf_r.c:
. Fixed the Cygnus Support conversion of e_lgamma_r.c to float.
This includes the generation of new polynomial and rational
approximations with fewer terms. For each approximation, include
a comment on an estimate of the accuracy over the relevant domain.
. Use integer literal constants instead of real literal constants.
Let compiler(s) do the job of conversion to the appropriate type.
This allows the removal of several explicit casts of double values
to float.
* src/e_lgammal.c:
. Wrapper for lgammal() about lgammal_r().
* src/imprecise.c:
. Remove the lgamma.
* src/math.h:
. Add a prototype for lgammal_r().
* man/lgamma.3:
. Document the new functions.
Reviewed by: bde
Diffstat (limited to 'lib/msun/src')
-rw-r--r-- | lib/msun/src/e_lgamma.c | 6 | ||||
-rw-r--r-- | lib/msun/src/e_lgamma_r.c | 19 | ||||
-rw-r--r-- | lib/msun/src/e_lgammaf_r.c | 171 | ||||
-rw-r--r-- | lib/msun/src/e_lgammal.c | 25 | ||||
-rw-r--r-- | lib/msun/src/imprecise.c | 1 | ||||
-rw-r--r-- | lib/msun/src/math.h | 6 |
6 files changed, 128 insertions, 100 deletions
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..ce2af94 100644 --- a/lib/msun/src/e_lgamma_r.c +++ b/lib/msun/src/e_lgamma_r.c @@ -83,6 +83,8 @@ __FBSDID("$FreeBSD$"); * */ +#include <float.h> + #include "math.h" #include "math_private.h" @@ -250,7 +252,7 @@ __ieee754_lgamma_r(double x, int *signgamp) p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); p = y*p1+p2; - r += (p-0.5*y); break; + r += (p-y/2); break; case 1: z = y*y; w = z*y; @@ -273,11 +275,11 @@ __ieee754_lgamma_r(double x, int *signgamp) r = half*y+p/q; z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ switch(i) { - case 7: z *= (y+6.0); /* FALLTHRU */ - case 6: z *= (y+5.0); /* FALLTHRU */ - case 5: z *= (y+4.0); /* FALLTHRU */ - case 4: z *= (y+3.0); /* FALLTHRU */ - case 3: z *= (y+2.0); /* FALLTHRU */ + case 7: z *= (y+6); /* FALLTHRU */ + case 6: z *= (y+5); /* FALLTHRU */ + case 5: z *= (y+4); /* FALLTHRU */ + case 4: z *= (y+3); /* FALLTHRU */ + case 3: z *= (y+2); /* FALLTHRU */ r += __ieee754_log(z); break; } /* 8.0 <= x < 2**58 */ @@ -293,3 +295,8 @@ __ieee754_lgamma_r(double x, int *signgamp) if(hx<0) r = nadj - r; return r; } + +#if (LDBL_MANT_DIG == 53) +__weak_reference(lgamma_r, lgammal_r); +#endif + diff --git a/lib/msun/src/e_lgammaf_r.c b/lib/msun/src/e_lgammaf_r.c index 9a7ab39..f940845 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) @@ -168,55 +160,50 @@ __ieee754_lgammaf_r(float x, int *signgamp) else {y = x; i=2;} } else { r = zero; - if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */ + if(ix>=0x3fdda618) {y=2-x;i=0;} /* [1.7316,2] */ else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */ else {y=x-one;i=2;} } switch(i) { case 0: z = y*y; - p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); - p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); + p1 = a0+z*(a2+z*a4); + p2 = z*(a1+z*(a3+z*a5)); p = y*p1+p2; - r += (p-(float)0.5*y); break; + r += (p-y/2); break; case 1: - z = y*y; - w = z*y; - p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ - p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); - p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); - p = z*p1-(tt-w*(p2+y*p3)); + p = t0+y*t1+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*t7))))); r += (tf + p); break; case 2: - p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); - p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); - r += (-(float)0.5*y + p1/p2); + p1 = y*(u0+y*(u1+y*u2)); + p2 = one+y*(v1+y*(v2+y*v3)); + r += (p1/p2-y/2); } } else if(ix<0x41000000) { /* x < 8.0 */ - i = (int)x; - y = x-(float)i; - p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); - q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); - r = half*y+p/q; + i = x; + y = x-i; + p = y*(s0+y*(s1+y*(s2+y*s3))); + q = one+y*(r1+y*(r2+y*r3)); + r = y/2+p/q; z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ switch(i) { - case 7: z *= (y+(float)6.0); /* FALLTHRU */ - case 6: z *= (y+(float)5.0); /* FALLTHRU */ - case 5: z *= (y+(float)4.0); /* FALLTHRU */ - case 4: z *= (y+(float)3.0); /* FALLTHRU */ - case 3: z *= (y+(float)2.0); /* FALLTHRU */ + case 7: z *= (y+6); /* FALLTHRU */ + case 6: z *= (y+5); /* FALLTHRU */ + case 5: z *= (y+4); /* FALLTHRU */ + case 4: z *= (y+3); /* FALLTHRU */ + case 3: z *= (y+2); /* FALLTHRU */ r += __ieee754_logf(z); break; } - /* 8.0 <= x < 2**58 */ - } else if (ix < 0x5c800000) { + /* 8.0 <= x < 2**24 */ + } else if (ix < 0x4b800000) { t = __ieee754_logf(x); z = one/x; y = z*z; - w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); + w = w0+z*(w1+y*w2); r = (x-half)*(t-one)+w; } else - /* 2**58 <= x <= inf */ + /* 2**24 <= x <= inf */ r = x*(__ieee754_logf(x)-one); if(hx<0) r = nadj - r; return r; 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/math.h b/lib/msun/src/math.h index e4e5682..2ff1731 100644 --- a/lib/msun/src/math.h +++ b/lib/msun/src/math.h @@ -496,8 +496,12 @@ long double tanhl(long double); long double tanl(long double); long double tgammal(long double); long double truncl(long double); - #endif /* __ISO_C_VISIBLE >= 1999 */ + +#if __BSD_VISIBLE +long double lgammal_r(long double, int *); +#endif + __END_DECLS #endif /* !_MATH_H_ */ |