dead code in lgamma_r[f]?

[Bruce Evans]

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