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.

[kcrisman]

comment:9

I'd like to have another reference for this (perhaps NIST or A&S) but this sounds plausible. So what is the relation between this and the mpfr versions?

[eviatarbach]

comment:10

Paul, could you please comment on this code for sage/rings/real_mpfr.pyx?

    def log_gamma(self):
        cdef RealNumber x = self._new()
        cdef int sign
        parent = (<RealField_class>self._parent)
        if parent.__prec > SIG_PREC_THRESHOLD: sig_on()
        mpfr_lgamma(x.value, &sign, self.value, parent.rnd)
        if parent.__prec > SIG_PREC_THRESHOLD: sig_off()
        if not mpfr_sgn(self.value) < 0:
            return x
        cdef RealNumber v = self._new()
        mpfr_div_si((<RealNumber>v).value, self.value, 2, parent.rnd)
        return parent.complex_field()(x, parent.pi() *
                                      (2 * (v - 1).ceil() + 1))

This correctly computes the principal branch. However, I'm concerned about precision; the result of lgamma is presumably guaranteed to be correct up to the precision of the input, but after doing some arithmetic it might not be. Is there any way to fix this?

[zimmermann6]

comment:11

Dear Eviatar,

I see no argument to tell the precision of the output in this function, thus I guess the result should be computed to the precision of the input.

Then the real part seems fine to me.

For the imaginary part, you create the variable v with the same precision than the input I guess.
In v-1 you might have a problem if v is very large, then v-1 rounded will equal v,
and (v-1).ceil() will equal v too. Thus the result will be too large by 2pi.

Paul

[zimmermann6]

comment:12

Replying to @kcrisman:

I'd like to have another reference for this (perhaps NIST or A&S) but this sounds plausible. So what is the relation between this and the mpfr versions?

there are some references on http://mathworld.wolfram.com/LogGammaFunction.html

Paul

[eviatarbach]

comment:13

Thank you. I suspected there might be a problem like this. I think maybe using MPFR for positive reals (since it's faster) and mpmath for negative reals might be best.

[eviatarbach]

Attachment: trac_12521.patch.gz

[eviatarbach]

comment:14

Patch uploaded. It uses mpmath to evaluate log_gamma in certain circumstances, while the MPFR implementation is left intact for positive reals. It also removes the lgamma and lngamma functions which were deprecated four years ago with #6992.

Patchbot apply trac_12521.2.patch

[burcin]

Reviewer: Burcin Erocal

[burcin]

comment:15

Attachment: trac_12521.2.patch.gz

The patchbot complains about old style doctest continuation.

Since you need to edit the patch anyway... In the py_lgamma function from pynac.pyx, you don't need to import call from mpmath or parent from sage.structure.coerce. These are already available as mpmath_utils.call and parent_c.

Thanks for the patch! Great work noticing that we don't actually compute the principal branch.

[burcin]

Author: Eviatar Bach

[eviatarbach]

comment:16

Thanks!

I had already fixed the doctest continuation and overwrote the previous patch, but the patchbot didn't test it again.

I fixed the import issue. I also noticed that the results for Python float and complex were inconsistent, so I fixed that and added some tests.

Patchbot apply trac_12521_3.patch

[jdemeyer]

comment:23

I am reverting the deprecation of lngamma() for PARI in #15767 (see comment [comment:20])

[rwst]

Branch: u/rws/evaluate_log_gamma_for_complex_input

[rwst]

Work Issues: speedup

[rwst]

Changed reviewer from Burcin Erocal to Burcin Erocal, Ralf Stephan

[rwst]

comment:27

Tests are good in functions, rings, symbolic. The speed issue in comment:21 still remains.

---

New commits:

f4b611d	12521: fixing numerical evaluation of log_gamma
d7731cb	12521: reviewer missed a part of the patch

[rwst]

Commit: d7731cb

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

bacd234

Merge branch 'develop' into t/12521/evaluate_log_gamma_for_complex_input

[sagetrac-git]

Changed commit from d7731cb to bacd234

[kcrisman]

comment:29

Has anyone checked the mathematics of this recently? Eviatar usually does a very good job on this stuff, but unfortunately Burcin didn't make it clear whether the rest was approved of other than the different imports.

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

d502104	Merge branch 'develop' into t/12521/evaluate_log_gamma_for_complex_input
c11c3ef	12521: simplifcations, doctest fixes

[sagetrac-git]

Changed commit from bacd234 to c11c3ef

[rwst]

Changed author from Eviatar Bach to Eviatar Bach, Ralf Stephan

[rwst]

Changed work issues from speedup to none

[paulmasson]

comment:34

With these changes the beta() function is still producing small imaginary parts around poles and not returning the appropriate signed Infinity at the poles. Presumably that should be on another ticket.

Otherwise looks fine to me.

[vbraun]

Changed branch from u/rws/evaluate_log_gamma_for_complex_input to c11c3ef