dead code in lgamma_r[f]?

[Steve Kargl]

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