dead code in lgamma_r[f]?
Bruce Evans
brde at optusnet.com.au
Tue Dec 10 05:45:19 UTC 2013
On Sun, 8 Dec 2013, Steve Kargl wrote:
> On Sun, Dec 08, 2013 at 03:21:04PM +1100, Bruce Evans wrote:
>> ...
>> gamma() was the first part of libm that I looked closely at (even
>> before trig functions in 2004/5). I didn't understand it at first and
>> tried to write my own. My first idea was to use the asymptotic
>> expansion as a kernel for all cases, especially for complex args.
>> ...
> Interesting idea. For x in [2,8), fdlibm code (as you probably
> know) uses recursion to reduce the evaluation to a log of a
> product and lgamma(2+s) with s in [0,1). One can then show that
I remember a bit more about why I tried to use the asymptotic expansion.
The fdlibm method in [0,8] only works for real args.
> lgamma(2+s) = s * (1 - E) + s**2 * P(s) where E is Euler's constant
> and P(s) is a polynomial in s. But, instead of using a minimax
This works for complex s, but only for s near 0.
> polynomial approximation fdlibm code uses a rational approximation
> lgamma(2+s) = 0.5*s + p(s) / q(s). This is the origin of the
> s0, ..., s6, r1, ..., r6 coefficients.
It is indeed hard to see why the code matches the comment. The 0.5*s
expansion is presumably to calculate the leading term exactly. (All
ration approximations should do this, since p(s) / q(s) has an error
of a couple of ulps.) s * (1 - E) would need extra precision.
> A quick scan of n1256.pdf shows that C99 does not require
> lgamma of a complex argument, but it does reserve the
> names clgamma[fl].
Same in C11 (n1570.pdf).
>> % lg(z) =
>> % {
>> % local(bol, rnew);
>> %
>> % if (abs(z - 1) < 1.e-1, return (lngamma(z)););
>> % if (abs(z - 2) < 1.e-1, return (lngamma(z)););
>>
>> Near 1 and 2, I planned to use a poly approximation, but that is not
>> so easy to do in pari so I call the pari builtin lgamma() to fake it.
>
> fdlibm's lgammaf_r does not appear to have an issue near 1 or 2.
> Starting at the lower bound of each interval and using nextafterf
> to scan up to the upper bound, I'm seeing
I mean that these cases are easy (even for complex z) so I didn't bother
testing various implementations of them for accuracy.
>> % lgb(z) = (z - 0.5) * log(z) - z + 0.5 * log(2 * Pi) + \
>> % sum(n = 1, N, b[n] / z^(2 * n - 1))
>>
>> This is the asymptotic expansion. Pari syntax is only perfect for
>> numeric things like this (and probably for some number theory things
>> which I don't know well). b[n] is a table of coefficients.
> ...
> portion of your above code is approximated by a minimax polynomial.
>
> /* 8.0 <= x < 2**58 */
> } else if (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
This is essentially the same -- the approximation is mostly a minimax
poly in 1/z.
Another problem for complex args is that it would be surpising if
minimax polys worked for them. The polys use delicate cancelations
that probably only work for real args.
Bruce
More information about the freebsd-numerics
mailing list