summaryrefslogtreecommitdiff
path: root/libm/e_lgamma_r.c
diff options
context:
space:
mode:
authorDenys Vlasenko <vda.linux@googlemail.com>2010-10-31 04:36:02 +0100
committerDenys Vlasenko <vda.linux@googlemail.com>2010-10-31 04:36:02 +0100
commit8b34cac8748067dbebc99f1b8450c5de91661db5 (patch)
tree645acac35767db45a53a68a561c6b307b32a16cf /libm/e_lgamma_r.c
parent7e308604ee0aa67bc37f37d6f69ed41481c427a0 (diff)
libm: fix tgamma to actually do return true gamma function
Signed-off-by: Denys Vlasenko <vda.linux@googlemail.com>
Diffstat (limited to 'libm/e_lgamma_r.c')
-rw-r--r--libm/e_lgamma_r.c42
1 files changed, 23 insertions, 19 deletions
diff --git a/libm/e_lgamma_r.c b/libm/e_lgamma_r.c
index 6d5d9c514..dbb0a0d5d 100644
--- a/libm/e_lgamma_r.c
+++ b/libm/e_lgamma_r.c
@@ -360,31 +360,35 @@ strong_alias(__ieee754_lgamma, gamma)
#endif
-/* FIXME! Looks like someone just used __ieee754_gamma_r,
- * believing it's a "true" gamma function, but it was not!
- * Our tgamma is WRONG.
- */
-
/* double tgamma(double x)
* Return the Gamma function of x.
*/
double tgamma(double x)
{
- double y;
- int local_signgam;
+ int sign_of_gamma;
+ int32_t hx;
+ u_int32_t lx;
- y = __ieee754_lgamma_r(x, &local_signgam); // was __ieee754_gamma_r
- if (local_signgam < 0)
- y = -y;
-#ifndef _IEEE_LIBM
- if (_LIB_VERSION == _IEEE_)
- return y;
- if (!isfinite(y) && isfinite(x)) {
- if (floor(x) == x && x <= 0.0)
- return __kernel_standard(x, x, 41); /* tgamma pole */
- return __kernel_standard(x, x, 40); /* tgamma overflow */
+ /* We don't have a real gamma implementation now. We'll use lgamma
+ and the exp function. But due to the required boundary
+ conditions we must check some values separately. */
+
+ EXTRACT_WORDS(hx, lx, x);
+
+ if (((hx & 0x7fffffff) | lx) == 0) {
+ /* Return value for x == 0 is Inf with divide by zero exception. */
+ return 1.0 / x;
}
-#endif
- return y;
+ if (hx < 0 && (u_int32_t)hx < 0xfff00000 && rint(x) == x) {
+ /* Return value for integer x < 0 is NaN with invalid exception. */
+ return (x - x) / (x - x);
+ }
+ if ((u_int32_t)hx == 0xfff00000 && lx == 0) {
+ /* x == -Inf. According to ISO this is NaN. */
+ return x - x;
+ }
+
+ x = exp(lgamma_r(x, &sign_of_gamma));
+ return sign_of_gamma >= 0 ? x : -x;
}
libm_hidden_def(tgamma)