cvs commit: src/lib/msun/src e_lgammaf_r.c

[Bruce Evans]

On Mon, 28 Nov 2005, Steve Kargl wrote:

> On Mon, Nov 28, 2005 at 11:16:27AM +0100, Andre Oppermann wrote:
>> Bruce Evans wrote:
>>>   ...
>>>     lib/msun/src         e_lgammaf_r.c
>>>   Log:
>>>   Fixed about 50 million errors of infinity ulps and about 3 million errors
>>>   of between 1.0 and 1.8509 ulps for lgammaf(x) with x between -2**-21 and
>>>   -2**-70.
>>
>> What is an ULP and are you going to write a paper on how FreeBSD has
>> the best, fastest and most precise msun library of all OSs?
>
> Units in the Last Place.
> http://docs.sun.com/source/806-3568/ncg_goldberg.html

Yes, that is probably the best reference.  Knuth is also good.  See
also ieee(3) and <float.h>.  The expression b**(1-p) (FLT_EPSILON,
etc.) in <float.h> is 1 ulp according to 1 sentence in Knuth, but this
is only 1 ulp for values normalized to be >= 1.0 and < 2.0.  Generally,
an ulp is just like FLT_EPSILON, with the epsilon normalized to the
same precision as the value instead of the value normalized to the same
precision as the epsilon.

I don't like writing papers, and rarely read them these days.  ISTR
reading the Goldberg one when it was first published not long before
I almost stopped reading many technical paper papers.

Hopefully at least Sun still knows everything in docs.sun.com and has the
most precise if not the fastest libm.

The history of {t,l}gamma's precision is interesting.  The old BSD
libm one (actually not so old; it is by McIlroy in {Oct,Nov} 1992)
tries much harder than the (1993) fdlibm one to be precise.  We still
use it for tgamma.  For lgamma, I think it achieves more precision for
positive args.  On negative args, it comments that "all bets are off"
because cases like lgamma(-2.4...) have a result near 0.  Such cases
are especially hard to make precise and lgamma's implementation is
especially unsuitable for making them precise (it uses essentially
lgamma(x) = log(f(x)/g(x)) = log(f(x)) - log(g(x)), where f(x) and
g(x) are large, so almost all precision is lost to cancellation when
the difference is small).

Bruce