lib/glibc
lib/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2025-05-08 18:21:37 +03:00
Code

(__ieee754_lgammal_r): Remove test for negative integer arg; sin_pi does it correctly.

Browse Source
This commit is contained in:
Andreas Jaeger 2002-01-28 10:18:33 +00:00
parent fbee8a1eb1
commit 98dee2c20d
1 changed files with 8 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -18,9 +18,9 @@
* *
* Method: * Method:
* 1. Argument Reduction for 0 < x <= 8 * 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 * reduce x to a number in [1.5,2.5] by
* lgamma(1+s) = log(s) + lgamma(s) * lgamma(1+s) = log(s) + lgamma(s)
* for example, * for example,
* lgamma(7.3) = log(6.3) + lgamma(6.3) * lgamma(7.3) = log(6.3) + lgamma(6.3)
* = log(6.3*5.3) + lgamma(5.3) * = log(6.3*5.3) + lgamma(5.3)
@ -50,13 +50,13 @@
* Let z = 1/x, then we approximation * Let z = 1/x, then we approximation
* f(z) = lgamma(x) - (x-0.5)(log(x)-1) * f(z) = lgamma(x) - (x-0.5)(log(x)-1)
* by * by
* 3 5 11 * 3 5 11
* w = w0 + w1*z + w2*z + w3*z + ... + w6*z * w = w0 + w1*z + w2*z + w3*z + ... + w6*z
* *
* 4. For negative x, since (G is gamma function) * 4. For negative x, since (G is gamma function)
* -x*G(-x)*G(x) = pi/sin(pi*x), * -x*G(-x)*G(x) = pi/sin(pi*x),
* we have * we have
* G(x) = pi/(sin(pi*x)*(-x)*G(-x)) * G(x) = pi/(sin(pi*x)*(-x)*G(-x))
* since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 * 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)|) * lgamma(x) = log(|Gamma(x)|)
@ -69,7 +69,7 @@
* lgamma(1)=lgamma(2)=0 * lgamma(1)=lgamma(2)=0
* lgamma(x) ~ -log(x) for tiny x * lgamma(x) ~ -log(x) for tiny x
* lgamma(0) = lgamma(inf) = inf * lgamma(0) = lgamma(inf) = inf
* lgamma(-integer) = +-inf * lgamma(-integer) = +-inf
* *
*/ */
@ -84,7 +84,7 @@ static long double
half = 0.5L, half = 0.5L,
one = 1.0L, one = 1.0L,
pi = 3.14159265358979323846264L, pi = 3.14159265358979323846264L,
two63 = 9.223372036854775808e18L, two63 = 9.223372036854775808e18L,
/* lgam(1+x) = 0.5 x + x a(x)/b(x) /* lgam(1+x) = 0.5 x + x a(x)/b(x)
-0.268402099609375 <= x <= 0 -0.268402099609375 <= x <= 0
@ -304,8 +304,6 @@ __ieee754_lgammal_r (x, signgamp)
} }
if (se & 0x8000) if (se & 0x8000)
{ {
if (x == __floorl(x))
return x / zero;
t = sin_pi (x); t = sin_pi (x);
if (t == zero) if (t == zero)
return one / fabsl (t); /* -integer */ return one / fabsl (t); /* -integer */