As per the ISO standard, tgamma(x) should compute Gamma(x), while lgamma(x) computes ln(Gamma(x)). The current implementation aliases tgamma and lgamma, which is incorrect. This was identified in Sept 2009, as per the comment in e_lgamma_r.c: """ /* FIXME! Looks like someone just used __ieee754_gamma_r, * believing it's a "true" gamma function, but it was not! * Our tgamma is WRONG. */ """ (Denis Vlasenko, 2009-09-18) The function tgamma() thus needs to be replaced entirely. A "quick and dirty" fix is to rewrite as follows: """ double y; int local_signgam; y = __ieee754_lgamma_r(x, &local_signgam); return local_signgam * exp(y); """ This provides relatively accurate results for small inputs, but diverges for large inputs. An alternate implementation can be found in FreeBSD's source, in /scratch/usr/src/lib/msun/bsdsrc/b_tgamma.c.
Care to submit a tested, proper patch? See http://uclibc.org/developing.html#contrib for instructions. Thanks!