Date: Sat, 13 Sep 2014 12:51:53 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-numerics@freebsd.org Subject: Re: lgammal[_r] patch Message-ID: <20140913195153.GB2323@troutmask.apl.washington.edu> In-Reply-To: <20140913193634.GA2323@troutmask.apl.washington.edu> References: <20140913193634.GA2323@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20140913195153.GB2323>