Opened 8 years ago
Closed 3 years ago
#12521 closed defect (fixed)
evaluate log gamma for complex input
Reported by:  kcrisman  Owned by:  AlexGhitza 

Priority:  critical  Milestone:  sage7.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 )
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.52e15] + [9.42477796076938 +/ 2.48e15]*I
Attachments (3)
Change History (39)
comment:1 followup: ↓ 2 Changed 8 years ago by
comment:2 in reply to: ↑ 1 Changed 8 years ago by
 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 7 years ago by
 Dependencies #10075 deleted
Removing dependency so as not to confuse anyone, now that that is in a longsince released version of Sage.
comment:4 Changed 7 years ago by
by the way, MPFR provides two functions: mpfr_lngamma
indeed returns NaN
between
2k1
and 2k
for k
a nonnegative 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 7 years ago by
 Cc benjaminfjones added
comment:6 Changed 7 years ago by
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.30.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 7 years ago by
 Cc eviatarbach added
comment:8 Changed 6 years ago by
 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/CS9423.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 followup: ↓ 12 Changed 6 years ago by
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 6 years ago by
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 6 years ago by
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 v1
you might have a problem if v
is very large, then v1
rounded will equal v
,
and (v1).ceil()
will equal v
too. Thus the result will be too large by 2pi
.
Paul
comment:12 in reply to: ↑ 9 Changed 6 years ago by
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 6 years ago by
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 6 years ago by
comment:14 Changed 6 years ago by
 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
Changed 6 years ago by
comment:15 Changed 6 years ago by
 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 6 years ago by
 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 6 years ago by
comment:17 Changed 6 years ago by
 Description modified (diff)
Why is the patchbot trying to apply the last two patches instead of just trac_12521_3.patch?
comment:18 Changed 6 years ago by
Let's try again:
Patchbot apply trac_12521_3.patch
comment:19 Changed 6 years ago by
 Milestone changed from sage5.11 to sage5.12
comment:20 Changed 6 years ago by
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 6 years ago by
 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 6 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:23 Changed 6 years ago by
comment:24 Changed 6 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:25 Changed 5 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:26 Changed 5 years ago by
 Branch set to u/rws/evaluate_log_gamma_for_complex_input
comment:27 Changed 5 years ago by
 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 5 years ago by
 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 5 years ago by
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 3 years ago by
 Description modified (diff)
comment:31 Changed 3 years ago by
 Cc paulmasson added
comment:32 Changed 3 years ago by
 Commit changed from bacd2349b07fd6296737d56b6f430c08fe28c767 to c11c3efddd8c6ae02df3fcf421616d5c665a0a62
comment:33 Changed 3 years ago by
 Milestone changed from sage6.4 to sage7.3
 Status changed from needs_work to needs_review
 Work issues speedup deleted
comment:34 Changed 3 years ago by
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 3 years ago by
 Status changed from needs_review to positive_review
comment:36 Changed 3 years ago by
 Branch changed from u/rws/evaluate_log_gamma_for_complex_input to c11c3efddd8c6ae02df3fcf421616d5c665a0a62
 Resolution set to fixed
 Status changed from positive_review to closed
I cannot reproduce some examples from the description with Sage 4.8:
Do you use a different version? Did you apply some patches?
Paul