Opened 9 years ago

Closed 4 years ago

# evaluate log gamma for complex input

Reported by: Owned by: kcrisman AlexGhitza critical sage-7.3 basic arithmetic lgamma log_gamma burcin, ktkohl, fstan, zimmerma, benjaminfjones, eviatarbach, paulmasson Eviatar Bach, Ralf Stephan Burcin Erocal, Ralf Stephan N/A c11c3ef (Commits) c11c3efddd8c6ae02df3fcf421616d5c665a0a62

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

### comment:1 follow-up: ↓ 2 Changed 9 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 9 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 8 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 8 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:6 Changed 8 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:8 Changed 7 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: ↓ 12 Changed 7 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 7 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 7 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 7 years ago by zimmerma

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

### comment:14 Changed 7 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 7 years ago by eviatarbach (previous) (diff)

### comment:15 Changed 7 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 7 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

### comment:17 Changed 7 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 7 years ago by eviatarbach

Let's try again:

Patchbot apply trac_12521_3.patch

### comment:19 Changed 7 years ago by jdemeyer

• Milestone changed from sage-5.11 to sage-5.12

### comment:20 Changed 7 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 7 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 7 years ago by vbraun_spam

• Milestone changed from sage-6.1 to sage-6.2

### comment:23 Changed 7 years ago by jdemeyer

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

### comment:24 Changed 7 years ago by vbraun_spam

• Milestone changed from sage-6.2 to sage-6.3

### comment:25 Changed 6 years ago by vbraun_spam

• Milestone changed from sage-6.3 to sage-6.4

### comment:26 Changed 6 years ago by rws

• Branch set to u/rws/evaluate_log_gamma_for_complex_input

### comment:27 Changed 6 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:

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

### comment:28 Changed 6 years ago by git

• Commit changed from d7731cb41c4b5c555e553996d132c725e73a3c5b to bacd2349b07fd6296737d56b6f430c08fe28c767

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

 ​bacd234 `Merge branch 'develop' into t/12521/evaluate_log_gamma_for_complex_input`

### comment:29 Changed 6 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 4 years ago by rws

• Description modified (diff)

### comment:32 Changed 4 years ago by git

• Commit changed from bacd2349b07fd6296737d56b6f430c08fe28c767 to c11c3efddd8c6ae02df3fcf421616d5c665a0a62

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`

### comment:33 Changed 4 years 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 4 years 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 4 years ago by paulmasson

• Status changed from needs_review to positive_review

### comment:36 Changed 4 years 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.