lgammal[_r] patch

[Steve Kargl]

On Sat, Sep 13, 2014 at 12:36:34PM -0700, Steve Kargl wrote:
> Following my .sig is a patch that contains the ld80 and ld128
> implementations of lgamma[_r]().  It also contains an update
> to lgammaf_r() and lgamma_r().  Specifically, it has
> 

Yes, a follow-up to my own post.

The polynomial approximation for the "t" coefficients has been changed
from whatever occurs in lgamma_r() to one that actually works in
lgamma[fl]_r().  To explain, I tried to match the approximation in
lgamma_r() in the other functions in tc-0.23 <= x <= tc+0.27, but my
attempts were fraught with instabilities.  Specifically, I tried
approximations of the form

1. f(x) = (lgamma(x) - lgamma(tc)) / (x - tc)**2 = t(x - tc)
2. f(x) = (lgamma(x) - lgamma(tc)) = (x - tc)**2 * t(x - tc)
3. f(x) = (lgamma(x) - lgamma(tc)) = t(x - tc)

For 1., there were stability issue for x -> tc.

For 2., the special Remes algorithm I developed would loose an
extrema.  I speculate, but haven't pursued, that multiple zeros
occurring in the rhs at x = tc is the root of the problem.

For 3., I had sensitive issues near the endpoints of the domain.
So, I adjusted the domain to [tc-0.24, tc+0.28].  This yielded the
coefficients found in the code.

-- 
steve