ENH: implement the principle branch of the logarithm of Gamma.

[person142]

Closes gh-5840.

[person142]

Hm... the test failures are coming from test_basic.py and test_spherical_bessel.py. I'm seeing those failures on master too, so I don't think this code has anything to do with them...

Edit: forgot to rebuild before testing.

[pv]

The failures are real, and they are caused by this PR.

The mechanism is a bit unclear, but I believe it's because you are including "complex.h" which introduces C99 complex numbers --- it should use numpy complex number functions instead, see _complexstuff.pxd. (Whether the failing spherical bessel functions are completely stably implemented if they fail with C99 complex is then a different question.)

A different question is that gammaln(x+0j) does compute some imaginary part --- as it appears to differ by integer multiples, is there a need to reimplement the whole computation from scratch?

[person142]

Yes, including "complex.h" was the culprit; thanks for the hint.

It's true that one could compute loggamma from gammaln. This implementation of loggamma has a much larger range than gammaln, however. For example, gammaln(-250 + 250j) = (-inf+nan*j), whereas loggamma computes the correct result. (At least the result agrees with mpmath.) In fact, loggamma is correct even at numbers like -1e20 + 1e20j, far beyond where gammaln has conked out. So it might make more sense to compute gammaln(z) from loggamma.

[codecov-io]

@@            master   #5931   diff @@
======================================
  Files          238     238       
  Stmts        43744   43805    +61
  Branches      8213    8213       
  Methods          0       0       
======================================
+ Hit          34168   34230    +62
- Partial       2602    2604     +2
+ Missed        6974    6971     -3

Review entire Coverage Diff as of e49a067

Powered by Codecov. Updated on successful CI builds.

[pv]

Ok good. Then, perhaps we can just replace the complex version of gammaln (clngamma_wrap: D->D in generate_ufuncs.py) by loggamma. I think the intent of the code was to also implement the primary branch, but there's some implementation mistake for Re z < 0 which gets the branch cuts wrong.

[person142]

Right. Though maybe I'll make a proposal, and if you don't like it I'll merge things into gammaln.

Currently gammaln computes log(|Gamma(x)|) for real inputs and log(Gamma(z)) for complex inputs (or, as you pointed out, maybe a buggy implementation of the principle branch). That behavior is a little strange. There is a precedent of functions doing slightly different things for real vs. complex numbers, but it's things like np.sqrt(-1) returning nan versus np.sqrt(-1 + 0j) returning 1j, i.e. it's only a change of definition of the domain of the function. This is not what gammaln does--it actually computes a different function for real inputs.

This behavior was only documented in 0.18 (before that the documentation said that gammaln always computed log(|Gamma(z)|)), so if there were ever a time to fix the strange behavior now is it. The best fix I can see is to have two functions: loggamma, which computes the principle branch, and logabsgamma (maybe there's a better name), which computes log(|Gamma(x)|) and only takes real inputs. The original gammaln would have to be kept around for back-compat, but we could add a message to its documentation which discourages its use and redirects the user to loggamma or logabsgamma (along with a thoughtful discussion of when you would want to use each).

[pv]

As a data point, in a similar situation with the legendre P function lpmn where the real version computed something quite different from complex, we just dropped the complex implementation and added a new function, cf. 44e7c61. So, in principle, we could also deprecate complex argument in gammaln which leaves loggamma as the canonical choice.

[person142]

Updated in accordance with @pv's comment. I stripped out the complex part of gammln; since the complex valued behavior isn't documented in any release I wasn't sure if adding a deprecation warning was necessary.

[pv]

Perhaps gammaln should raise a deprecationwarning for complex inputs, rather than scrapping it right away. iow, "gammaln -> _gammaln" and add a Python wrapper function def gammaln in basic.py that checks if np.iscomplexobj and raises a warning if so. Then after a few releases we could remove the complex version.

[person142]

Updated.

[person142]

Updated again to explicitly mention the deprecation in the documentation of gammaln. Also fixed a bug I introduced in the gammaln wrapper and a mistake in the definition of loggamma in the documentation.

Also, @pv is correct that the current behavior of gammaln for complex arguments is a buggy implementation of loggamma and not an implementation of log(gamma(z)). (The source uses the Stirling series for loggamma, plus gammaln returns imaginary parts outside of the interval (-pi, pi]. So I guess that this pull request should really be a BUG instead of an ENH. Edit: for comparison, if you replace loggamma with gammaln in the TestSystematic.test_loggamma_complex test you get 1509/4150 failures.

A point worth noting is that for complex values gamma is computed by taking exp(gammaln(z)). (This happens inside the cgama function in specfun.f.) Some of the errors in gammaln for complex arguments propagate into gamma, which isn't caught by the current tests because gamma is only tested for real arguments (when it's using cephes anyway). If you change the mpmath tests to check ComplexArg() instead of Arg() there are 32 failures, but these can be fixed by using exp(loggamma(z)) instead. But maybe that change is best left for a separate pull request? (I think a good implementation would be a little more complicated than just returning exp(loggamma(z)); since the Gamma function blows up so quickly it would be good to do checks that return nan's right away instead of doing a bunch of unnecessary computations.)

[person142]

Updated yet again to improve accuracy around the pole at 0 and the zeros at 1 and 2.

[person142]

Updated with more improvements to accuracy in the z = 0, 1, 2 region and some simplifications to the code. Since there have been a number of changes I'll summarize where things stand:

• The complex behavior of gammaln has been moved into loggamma. Now gammaln computes log(|Gamma(x)|) and loggamma computes the principle branch of the logarithm of Gamma. A deprecation warning has been added for complex arguments to gammaln.
• The old specfun.f code computing the principle branch has been rewritten to improve accuracy. As mentioned above, the test TestSystematic.test_loggamma_complex goes from 1509/4150 failures with the old code to 0 failures with the new code.

I ran into a bug which (AFAICT) turned out to be caused by differences in accuracy of npy_clog around z = 1on different systems. (Compare the failures in 076d7a1 to the fix in acec387.) I fixed the issue by adding an expansion for log around z = 1 to the loggamma code, but perhaps there's a better way to do this?

[pv]

thanks, merged

[pv]

On the npy_clog issue --- may be a bug in the system libm, or in the numpy implementation of clog (depending of which one is used). So maybe report to numpy?