Date: Sat, 7 Dec 2013 17:29:39 -0800 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: freebsd-numerics@FreeBSD.org, Filipe Maia <filipe.c.maia@gmail.com> Subject: Re: dead code in lgamma_r[f]? Message-ID: <20131208012939.GA80145@troutmask.apl.washington.edu> In-Reply-To: <20131207064602.GA76042@troutmask.apl.washington.edu> References: <20131205173700.GA64575@troutmask.apl.washington.edu> <20131205182324.GA65135@troutmask.apl.washington.edu> <CAN5hRiV6N1arc4RDv=9JbRiXp-J9o3WzAbeZSOpAxks2ZeG%2B_w@mail.gmail.com> <20131205195225.GA65732@troutmask.apl.washington.edu> <20131206102724.W1033@besplex.bde.org> <20131206114426.I1329@besplex.bde.org> <20131207064602.GA76042@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Fri, Dec 06, 2013 at 10:46:02PM -0800, Steve Kargl wrote:
> On Fri, Dec 06, 2013 at 12:55:16PM +1100, Bruce Evans wrote:
> > On Fri, 6 Dec 2013, Bruce Evans wrote:
> >
> > > On Thu, 5 Dec 2013, Steve Kargl wrote:
> > > ...
> > >> If we again look at the code from __ieee754_lgamma_r(), we see
> > >> that sin_pi() is called if ix < 0x43300000, so by the time we
> > >> arrive at the 'if(ix<0x43300000)' statement we already know that
> > >> the condition is true.
> > >
> > > No, only for negative integers. hx<0 classifies negative values, and
> > > then ix>=0x43300000 classifies numbers that are so large negative that
> > > they must be integers, and the result is a sort of pole error. We
> > > are just filtering out this case, perhaps as an optimization.
> >
> > Oops, sin_pi() is only called for negative integers, so your change
> > seems to be correct. Just add a comment about the limited domain
> > of sin_pi() (it already has one saying that "x is assumed negative".
> >
>
> I wish to retract my earlier statement that after 2 additional
> years of reading fdlibm code that it was easier to work with.
> I spent the better part of Friday giving myself a headache trying
> to understand the algorithm for lgamma_r(). The code for x in
> the interval (0,2) does not match any comment in lgamma_r().
>
> I also think, but can't prove yet, that like erff() the
> polynomial and rational approximations in lgammaf_r() have
> too many terms.
>
I really need to stop looking at fdlibm code. The threshold(s) in
e_lgammaf.c could be chosen better (and actually documented). In
the float code
/* 8.0 <= x < 2**58 */
} else if (ix < 0x4e800000) { // (ix < 0x5c800000) {
t = __ieee754_logf(x);
z = one/x;
y = z*z;
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
r = (x-half)*(t-one)+w;
} else
/* 2**58 <= x <= inf */
r = x*(__ieee754_logf(x)-one);
the 2**58 is much too large. Asymptotically, we have
lgamma(x)~(x-0.5)log(x)-x = x*(log(x)-1) - 0.5*log(x)
where the 0.5*log(x) is insignificant when compared to 2**58.
I suspect that 2**58 is a sloppy threshold for double, which
can be reduced to 2**30 (or smaller) for float. I suppose I
could look for a tighter bound, but 2**(p+slop) with p, the
precision, probably is sufficient.
--
Steve
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20131208012939.GA80145>