Date: Tue, 29 Nov 2005 11:49:13 +1100 (EST) From: Bruce Evans <bde@zeta.org.au> To: Steve Kargl <sgk@troutmask.apl.washington.edu> Cc: cvs-src@freebsd.org, src-committers@freebsd.org, Andre Oppermann <andre@freebsd.org>, Bruce Evans <bde@freebsd.org>, cvs-all@freebsd.org Subject: Re: cvs commit: src/lib/msun/src e_lgammaf_r.c Message-ID: <20051129110058.T33820@delplex.bde.org> In-Reply-To: <20051128172718.GA59929@troutmask.apl.washington.edu> References: <200511280832.jAS8WGvs059057@repoman.freebsd.org> <438AD8FB.A8B96AB6@freebsd.org> <20051128172718.GA59929@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20051129110058.T33820>