gammaln: why no negative arguments allowed?

[Fri Jan 13 19:33:45 GMT 2012 - OK] Removed unmeaningful tag

Based on MATLAB documentation for gammaln, it seems to me that the following two commands should generate the same answer:

>> log(gamma(-1.5))

ans = 0.8600

>> gammaln(-1.5)

ans = Inf

So they are actually not the same. It is clear from experimenting that gammaln does not like negative arguments, but the question is why, and also why does the documentation description of gammaln written so that it would seem that the above two answers should be the same?

Note: "doc gammaln" in v2011b returns the a description including the statement "gammaln(A) = log(gamma(A))" with no discussion of a restricted domain for A.

Is this a bug or is there a real mathematical issue here that I'm not seeing? The analytic continuation of the classic integral definition of the gamma function seems perfectly well behaved at -1.5, or any complex number barring non-positive integers. Why not have gammaln(-(abs(x)) just return log(gamma(-(abs(x)))) for negative numbers, and Inf for negative integers and zero?

PMM

[the cyclist]

Really interesting question. Nothing obvious comes to mind as to why MATLAB could not report the correct answer here. Out of curiosity, I tried the same calculation in Excel and in R, both of which have equivalent commands. Excel reports #NUM! for both gammaln(-1.5) and gammaln(-0.5). R reports lgamma(-1.5)=0.860047... and lgamma(-0.5)=1.265512... . I find the latter a bit dangerous, as it does not warn that there is an unreported imaginary part. I am not an expert in R, though, so there may be a better way to handle that case than my simple-minded approach. [lgamma(as.complex(-1.5)) gives an error, though.]

Answer 1

[Andrew Newell]

The functions gamma and gammaln, and their documentation, are restricted to the real numbers, but these functions can only be meaningfully discussed in the complex plane. Here is what Maple says about the two functions:

"The Gamma function is defined for Re(z)>0 by (the integral expression goes here) and is extended to the rest of the complex plane, less the non-positive integers, by analytic continuation. GAMMA has a simple pole at each of the points z=0,-1,-2, ..."

"For positive real arguments z, the lnGAMMA function is defined by:

lnGAMMA(z) = ln(GAMMA(z))

For complex z, Maple evaluates the principal branch of the log-Gamma function, which is defined by analytic continuation from the positive real axis. Each of the points z=0,-1,-2, ... is a singularity and a branch point, and the union of the branch cuts is the negative real axis. On the branch cuts, the values of lnGAMMA(z) are determined by continuity from above. (Note, therefore, that lnGAMMA href=""</a> ln@GAMMA in general)."

Note also that Maple does not get the same answers as MATLAB for lnGAMMA. For example, it returns 0.8600470154 - 6.283185307*i for lnGAMMA(-1.5). Clearly, MATLAB's answer of Inf is wrong.

EDIT: I have submitted a bug report on this issue.

[the cyclist]

Thanks for following up, Andrew. I'll be interested to hear the response from MW on this. If you could post a summary, that'd be great.

[Sean de Wolski]

We are aware of this bug.