A problem with tgamma function

Lev Bishop lev.bishop@gmail.com
Fri Oct 12 04:01:00 GMT 2007


On 10/11/07, Angelo Graziosi wrote:
> > That is a bug report against glibc.  You cannot draw any inferences
> > between that and newlib/Cygwin because they are totally separate and
> > unrelated code bases.

Technically, when it comes to maths, both lean very heavily on fdlibm,
so they really aren't unrelated codebases.

> > That would be expected given Lev's indication of a fix in newlib
>
> What I wished to stress was that the 'tgamma' function of newlib gives the
> same results of the glibc-tgamma, for which a bug report has been sent.
>
> And this regardless the "Lev's indication of a fix in newlib".

There are two completely different problems. The newlib segfault
problem was due to a missing include file, causing a missing prototype
for __iee754_gammaf_r(), causing C's default type rules to infer a
type of double for the first parameter to __ieee754_gammaf_r() instead
of float, in turn causing the write intended for &signgam to instead
hit a random location and segfault/coredump. This problem is fixed in
CVS.

As a second issue, the glibc implementation of tgammaf() is accurate
to worse than 10^-6 although the glibc documentation claims it should
be 1ulp (about 10^-7). This is because __ieee754_gammaf_r(x) is
implemented as
__ieee754_expf (__ieee754_lgammaf_r (x))
(with an XXX FIXME comment). Newlib does exactly the same thing (minus
the coment). The trouble with this is that you lose some accuracy. The
error can be as large as log(tgamma(x)) ulps, which in the worst case
(log(FLT_MAX)) is about 89 ulps. Actually, in practice the error
doesn't seem to be quite as bad -- the worst I've seen is about 64
ulps -- but it's still not as accurate as you'd like.

I'm sure glibc and newlib would both appreciate a good algorithm for
tgamma(), if you felt like submitting one...

Lev

--
Unsubscribe info:      http://cygwin.com/ml/#unsubscribe-simple
Problem reports:       http://cygwin.com/problems.html
Documentation:         http://cygwin.com/docs.html
FAQ:                   http://cygwin.com/faq/



More information about the Cygwin mailing list