AuroraMiddleware/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2025-01-19 07:00:08 +00:00
Code Issues Packages Projects Releases Wiki Activity

Merge with mainline.

Browse Source
This commit is contained in:
Andreas Jaeger 2002-08-21 08:00:55 +00:00
parent a4e8b66c28
commit 4d5dbeceac
1 changed files with 26 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
sysdeps/ieee754/ldbl-96/e_lgammal_r.c
View File

@ -18,9 +18,9 @@
*
* Method:
* 1. Argument Reduction for 0 < x <= 8
* Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
* reduce x to a number in [1.5,2.5] by
* lgamma(1+s) = log(s) + lgamma(s)
* Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
* reduce x to a number in [1.5,2.5] by
* lgamma(1+s) = log(s) + lgamma(s)
* for example,
* lgamma(7.3) = log(6.3) + lgamma(6.3)
* = log(6.3*5.3) + lgamma(5.3)
@ -50,13 +50,13 @@
* Let z = 1/x, then we approximation
* f(z) = lgamma(x) - (x-0.5)(log(x)-1)
* by
* 3 5 11
* 3 5 11
* w = w0 + w1*z + w2*z + w3*z + ... + w6*z
*
* 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))
* we have
* G(x) = pi/(sin(pi*x)*(-x)*G(-x))
* since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
* Hence, for x<0, signgam = sign(sin(pi*x)) and
* lgamma(x) = log(|Gamma(x)|)
@ -69,7 +69,7 @@
* lgamma(1)=lgamma(2)=0
* lgamma(x) ~ -log(x) for tiny x
* lgamma(0) = lgamma(inf) = inf
* lgamma(-integer) = +-inf
* lgamma(-integer) = +-inf
*
*/
@ -84,6 +84,7 @@ static long double
half = 0.5L,
one = 1.0L,
pi = 3.14159265358979323846264L,
two63 = 9.223372036854775808e18L,
/* lgam(1+x) = 0.5 x + x a(x)/b(x)
-0.268402099609375 <= x <= 0
@ -206,8 +207,7 @@ sin_pi (x)
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
i1 = (ix << 16) | (i0 >> 16);
ix = (ix << 16) | (i0 >> 16);
if (ix < 0x3ffd8000) /* 0.25 */
return __sinl (pi * x);
y = -x; /* x is assume negative */
@ -219,13 +219,25 @@ sin_pi (x)
z = __floorl (y);
if (z != y)
{ /* inexact anyway */
y *= half;
y = 2.0 * (y - __floorl (y)); /* y = |x| mod 2.0 */
n = (int) (y * 4.0);
y *= 0.5;
y = 2.0*(y - __floorl(y)); /* y = |x| mod 2.0 */
n = (int) (y*4.0);
}
else
{
return (zero + zero);
if (ix >= 0x403f8000) /* 2^64 */
{
y = zero; n = 0; /* y must be even */
}
else
{
if (ix < 0x403e8000) /* 2^63 */
z = y + two63; /* exact */
GET_LDOUBLE_WORDS (se, i0, i1, z);
n = i1 & 1;
y = n;
n <<= 2;
}
}
switch (n)
@ -267,6 +279,7 @@ __ieee754_lgammal_r (x, signgamp)
int i, ix;
u_int32_t se, i0, i1;
*signgamp = 1;
GET_LDOUBLE_WORDS (se, i0, i1, x);
ix = se & 0x7fff;
@ -276,7 +289,6 @@ __ieee754_lgammal_r (x, signgamp)
ix = (ix << 16) | (i0 >> 16);
/* purge off +-inf, NaN, +-0, and negative arguments */
*signgamp = 1;
if (ix >= 0x7fff0000)
return x * x;
@ -292,8 +304,6 @@ __ieee754_lgammal_r (x, signgamp)
}
if (se & 0x8000)
{
if (x == __floorl(x))
return x / zero;
t = sin_pi (x);
if (t == zero)
return one / fabsl (t); /* -integer */