evaluate log gamma for complex input

[kcrisman]

Currently, we use MPFR or Ginac to evaluate log_gamma, but this returns NaN for negative input with even ceiling.

sage: log_gamma(-2.1)
NaN
sage: log_gamma(-3.1)
0.400311696703985
sage: log_gamma(-4.1)
NaN
sage: log_gamma(-5.1)
-2.63991581673655
sage: log_gamma(-21/10).n()
NaN
sage: get_systems('log_gamma(-21/10).n()')
['ginac']
sage: log_gamma(CC(-2.1))
1.53171380819509 + 3.14159265358979*I

We can use mpmath or something other trick to get this to work, now that #10075 has a nice symbolic function available. See #10072 for where we originally got better numerical evaluation.

sage: mpmath.loggamma(-2.1)
mpc(real='1.5317138081950856', imag='-9.4247779607693793')

Putting as defect because there is a log gamma for negative numbers, though we should talk about branches...

Apply: attachment: trac_12521_3.patch

There is now also the arb option:

sage: x=ComplexBallField(53)(-2.1)
sage: %timeit _ = x.log_gamma()
The slowest run took 8.11 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 8.15 µs per loop
sage: x.log_gamma()
[1.53171380819509 +/- 5.52e-15] + [-9.42477796076938 +/- 2.48e-15]*I

CC: @burcin @sagetrac-ktkohl @sagetrac-fstan @zimmermann6 @benjaminfjones @eviatarbach @paulmasson

Component: basic arithmetic

Keywords: lgamma log_gamma

Author: Eviatar Bach, Ralf Stephan

Branch/Commit: c11c3ef

Reviewer: Burcin Erocal, Ralf Stephan

Issue created by migration from https://trac.sagemath.org/ticket/12521

[zimmermann6]

comment:1

I cannot reproduce some examples from the description with Sage 4.8:

sage: from sage.misc.citation import get_systems
sage: get_systems('log_gamma(-21/10).n()')      
['MPFR']
sage: log_gamma(CC(-2.1))
NaN

Do you use a different version? Did you apply some patches?

Paul

[kcrisman]

comment:2

Do you use a different version? Did you apply some patches?

Yes, I guess this is after #10075. There is the same problem before that patch, but the systems and outputs will be different. In particular, inheriting from the parent of the input (which most generic Function classes do) if it has a method of the same name (which CC does) makes this change.

There is no patch, so this doesn't depend on #10075 in that sense, but I will put it anyway so that it's clear. Thanks!

[kcrisman]

Dependencies: #10075

[kcrisman]

Changed dependencies from #10075 to none

[kcrisman]

comment:3

Removing dependency so as not to confuse anyone, now that that is in a long-since released version of Sage.

[zimmermann6]

comment:4

by the way, MPFR provides two functions: mpfr_lngamma indeed returns NaN between
-2k-1 and -2k for k a non-negative integer.

However there is another function mpfr_lgamma which returns log(abs(gamma(x))),
and separately the sign of gamma(x). It suffices to add pi*I when the sign is
negative.

This would solve the problem at least for real input.

Paul

[eviatarbach]

comment:6

There is a related issue with the beta function:

sage: beta(-1.3,-0.4)
NaN
sage: (gamma(-1.3) * gamma(-0.4)) / gamma(-1.3 - 0.4) # expected functionality
-4.92909641669610

This is presumably because GiNaC uses exp(log_gamma(-1.3)+log_gamma(-0.4)-log_gamma(-1.3-0.4)) to evaluate, and gamma(-0.4) is negative. I'm confused as to why it works fine in the GiNaC interactive shell though.

[eviatarbach]

comment:8

Our current implementation is wrong. For the principal branch, log_gamma(z) should not be equal to log(gamma(z)) in general. See https://cs.uwaterloo.ca/research/tr/1994/23/CS-94-23.pdf:

ln_gamma(x) = ln(gamma(x)) + 2*pi*i*ceil(x/2 - 1)

Changing to critical because the documentation claims that it's the principal branch, which is mathematically wrong.