Opened 6 years ago

Closed 16 months ago

#12521 closed defect (fixed)

evaluate log gamma for complex input

Reported by: kcrisman Owned by: AlexGhitza
Priority: critical Milestone: sage-7.3
Component: basic arithmetic Keywords: lgamma log_gamma
Cc: burcin, ktkohl, fstan, zimmerma, benjaminfjones, eviatarbach, paulmasson Merged in:
Authors: Eviatar Bach, Ralf Stephan Reviewers: Burcin Erocal, Ralf Stephan
Report Upstream: N/A Work issues:
Branch: c11c3ef (Commits) Commit: c11c3efddd8c6ae02df3fcf421616d5c665a0a62
Dependencies: Stopgaps:

Description (last modified by rws)

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: 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

Attachments (3)

trac_12521.patch (10.0 KB) - added by eviatarbach 4 years ago.
trac_12521.2.patch (11.8 KB) - added by eviatarbach 4 years ago.
trac_12521_3.patch (13.0 KB) - added by eviatarbach 4 years ago.

Download all attachments as: .zip

Change History (39)

comment:1 follow-up: Changed 6 years ago by zimmerma

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

comment:2 in reply to: ↑ 1 Changed 6 years ago by kcrisman

  • Dependencies set to #10075

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!

comment:3 Changed 5 years ago by kcrisman

  • Dependencies #10075 deleted

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

comment:4 Changed 5 years ago by zimmerma

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

comment:5 Changed 5 years ago by benjaminfjones

  • Cc benjaminfjones added

comment:6 Changed 5 years ago by eviatarbach

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.

comment:7 Changed 5 years ago by eviatarbach

  • Cc eviatarbach added

comment:8 Changed 4 years ago by eviatarbach

  • Priority changed from major to critical

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.

comment:9 follow-up: Changed 4 years ago by 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?

comment:10 Changed 4 years ago by eviatarbach

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?

comment:11 Changed 4 years ago by zimmerma

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

comment:12 in reply to: ↑ 9 Changed 4 years ago by zimmerma

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

comment:13 Changed 4 years ago by eviatarbach

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.

Changed 4 years ago by eviatarbach

comment:14 Changed 4 years ago by eviatarbach

  • Status changed from new to needs_review

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

Last edited 4 years ago by eviatarbach (previous) (diff)

Changed 4 years ago by eviatarbach

comment:15 Changed 4 years ago by burcin

  • Authors set to Eviatar Bach
  • Description modified (diff)
  • Reviewers set to Burcin Erocal
  • Status changed from needs_review to needs_work

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.

comment:16 Changed 4 years ago by eviatarbach

  • Status changed from needs_work to needs_review

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

Changed 4 years ago by eviatarbach

comment:17 Changed 4 years ago by eviatarbach

  • Description modified (diff)

Why is the patchbot trying to apply the last two patches instead of just trac_12521_3.patch?

comment:18 Changed 4 years ago by eviatarbach

Let's try again:

Patchbot apply trac_12521_3.patch

comment:19 Changed 4 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:20 Changed 4 years ago by jdemeyer

I am against removing the lngamma() method and against the deprecation from #6992 for PARI. Inside PARI/GP, the name is lngamma() and that name should be available to Sage too.

comment:21 Changed 4 years ago by jdemeyer

  • Status changed from needs_review to needs_work

Also, it seems the algorithm for complex numbers is much faster than for real numbers, so we should probably use the complex algorithm for all inputs:

sage: a=RealField(53)(2.5); timeit('log_gamma(a)')
625 loops, best of 3: 118 µs per loop
sage: a=ComplexField(53)(2.5); timeit('log_gamma(a)')
625 loops, best of 3: 112 µs per loop
sage: a=RealField(10**4)(2.5); timeit('log_gamma(a)')
5 loops, best of 3: 506 ms per loop
sage: a=ComplexField(10**4)(2.5); timeit('log_gamma(a)')
625 loops, best of 3: 782 µs per loop

comment:22 Changed 4 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:23 Changed 4 years ago by jdemeyer

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

comment:24 Changed 3 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:25 Changed 3 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:26 Changed 3 years ago by rws

  • Branch set to u/rws/evaluate_log_gamma_for_complex_input

comment:27 Changed 3 years ago by rws

  • Commit set to d7731cb41c4b5c555e553996d132c725e73a3c5b
  • Reviewers changed from Burcin Erocal to Burcin Erocal, Ralf Stephan
  • Work issues set to speedup

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


New commits:

f4b611d12521: fixing numerical evaluation of log_gamma
d7731cb12521: reviewer missed a part of the patch

comment:28 Changed 3 years ago by git

  • Commit changed from d7731cb41c4b5c555e553996d132c725e73a3c5b to bacd2349b07fd6296737d56b6f430c08fe28c767

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

bacd234Merge branch 'develop' into t/12521/evaluate_log_gamma_for_complex_input

comment:29 Changed 3 years ago by kcrisman

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.

comment:30 Changed 16 months ago by rws

  • Description modified (diff)

comment:31 Changed 16 months ago by paulmasson

  • Cc paulmasson added

comment:32 Changed 16 months ago by git

  • Commit changed from bacd2349b07fd6296737d56b6f430c08fe28c767 to c11c3efddd8c6ae02df3fcf421616d5c665a0a62

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

d502104Merge branch 'develop' into t/12521/evaluate_log_gamma_for_complex_input
c11c3ef12521: simplifcations, doctest fixes

comment:33 Changed 16 months ago by rws

  • Authors changed from Eviatar Bach to Eviatar Bach, Ralf Stephan
  • Milestone changed from sage-6.4 to sage-7.3
  • Status changed from needs_work to needs_review
  • Work issues speedup deleted

comment:34 Changed 16 months ago by paulmasson

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.

comment:35 Changed 16 months ago by paulmasson

  • Status changed from needs_review to positive_review

comment:36 Changed 16 months ago by vbraun

  • Branch changed from u/rws/evaluate_log_gamma_for_complex_input to c11c3efddd8c6ae02df3fcf421616d5c665a0a62
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.