lgammal[_r] patch
Steve Kargl
sgk at troutmask.apl.washington.edu
Sat Sep 13 19:51:54 UTC 2014
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
More information about the freebsd-numerics
mailing list