X-Recipient: archive-cygwin AT delorie DOT com
X-Spam-Check-By: sourceware.org
Message-ID: <ba40711f0710112101x51b7153bkd7d418aaa8deac84@mail.gmail.com>
Date: Fri, 12 Oct 2007 00:01:22 -0400
From: "Lev Bishop" <lev DOT bishop AT gmail DOT com>
To: cygwin AT cygwin DOT com
Subject: Re: A problem with tgamma function
In-Reply-To: <Pine.CYG.4.58.0710111825390.3032@homepc>
MIME-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: 7bit
Content-Disposition: inline
References: <Pine DOT CYG DOT 4 DOT 58 DOT 0710111825390 DOT 3032 AT homepc>
Mailing-List: contact cygwin-help AT cygwin DOT com; run by ezmlm
Precedence: bulk
List-Id: <cygwin.cygwin.com>
List-Unsubscribe: <mailto:cygwin-unsubscribe-archive-cygwin=delorie DOT com AT cygwin DOT com>
List-Subscribe: <mailto:cygwin-subscribe AT cygwin DOT com>
List-Archive: <http://sourceware.org/ml/cygwin/>
List-Post: <mailto:cygwin AT cygwin DOT com>
List-Help: <mailto:cygwin-help AT cygwin DOT com>, <http://sourceware.org/ml/#faqs>
Sender: cygwin-owner AT cygwin DOT com
Mail-Followup-To: cygwin AT cygwin DOT com
Delivered-To: mailing list cygwin AT cygwin DOT com
Note-from-DJ: This may be spam

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/