#10075 closed enhancement (fixed)
Make log gamma symbolic
Reported by: | kcrisman | Owned by: | burcin |
---|---|---|---|
Priority: | major | Milestone: | sage-5.0 |
Component: | symbolics | Keywords: | sd35.5 |
Cc: | ktkohl, benjaminfjones | Merged in: | sage-5.0.beta7 |
Authors: | Karen Kohl, Karl-Dieter Crisman | Reviewers: | Karl-Dieter Crisman, Benjamin Jones |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | #12507, #9130 | Stopgaps: |
Description (last modified by )
Currently, there is no way to send log_gamma
to Maxima, for instance. This can be fixed by following the models in the functions/ directory; it should be possible to make it a GinacFunction. Before doing so, though, one will have to resolve #10072, since the evaluation will be wrong (?) otherwise.
Apply only trac_10075.patch.
Attachments (3)
Change History (23)
comment:1 Changed 9 years ago by
- Cc ktkohl added
- Keywords sd35.5 added
comment:2 Changed 9 years ago by
- Cc benjaminfjones added
comment:3 Changed 9 years ago by
comment:4 Changed 9 years ago by
More precisely,
class Function_log_gamma(GinacFunction): def __init__(self): GinacFunction.__init__(self, "log_gamma", latex_name=r'\log\Gamma', ginac_name='lgamma', conversions={'mathematica':'LogGamma','maxima':'log_gamma'}) log_gamma = Function_log_gamma()
causes this failure.
comment:5 Changed 9 years ago by
- Milestone changed from sage-4.8 to sage-5.0
It works if you drop the ginac_name
argument. The function is named log_gamma
in pynac.
comment:6 Changed 9 years ago by
Weird. So what about things like
unsigned lgamma_serial "GiNaC::lgamma_SERIAL::serial" # logarithm of gamma function
?
Changed 9 years ago by
symbolic log_gamma (with modification of functions.rst in case merged before #9130)
Changed 9 years ago by
symbolic log_gamma (with modification of functions.rst in case merged after #9130)
comment:7 Changed 9 years ago by
- Status changed from new to needs_review
Load one of the above two patches depending on whether the functions.rst documentation file has been modified already (as in the combined patch for #9130).
The second patch file above (without functions.rst) was edited by hand from the first.
comment:8 Changed 9 years ago by
Burcin points out that
sage: log_gamma(-2.1) NaN
is not good. Sage itself does
sage: log(gamma(-2.1)) 1.53171380819509 + 3.14159265358979*I
but Wolfram Alpha says
1.53171... - 9.42478... i
so the branches seem to differ even there.
comment:9 follow-up: ↓ 10 Changed 9 years ago by
- Description modified (diff)
Since #9130 has positive review: Apply only trac_10075_log_gamma_without_functions.rst.patch.
One might think that Burcin's comment about -2.1
makes this 'needs work', but that is actually the current Sage behavior as well, so in principle that would be a different ticket, since making log_gamma
symbolic would not introduce a regression...
In fact,
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: get_systems('log_gamma(2.1)') ['MPFR']
Apparently this is how MPFR deals with this function. So maybe all is well?
comment:10 in reply to: ↑ 9 Changed 9 years ago by
Apparently this is how MPFR deals with this function. So maybe all is well?
I mean, for this ticket. Though we should not claim that it is evaluated by Ginac, because it isn't (all the above is in Sage with or without this patch).
Believe it or not:
Not any negative value, but in lngamma.c: /* if x < 0 and -2k-1 <= x <= -2k, then lngamma(x) = NaN */ probably because the gamma value is negative. This is because MPFR defines lngamma as log(gamma(x)) while the C standard defines it as log|gamma(x)|. I wonder if this should be regarded as a bug or if a new function (say, mpfr_lgamma) should be defined in MPFR (in which case, not before 2.3.0). Do other standards (other languages) define such a function, either as log(gamma(x)) or as log|gamma(x)|?
I'm cc:ing Paul Z. just to confirm that this is intended MPFR behavior. We should then open another ticket to make sure to use mpmath or ginac or something to get complex answers. We currently somehow use PARI to get the complex versions.
sage: log_gamma(CC(-2.1)) 1.53171380819509 + 3.14159265358979*I sage: from sage.misc.citation import get_systems sage: get_systems('log_gamma(CC(-2.1))') ['PARI', 'MPFR']
comment:11 Changed 9 years ago by
- Reviewers set to Karl-Dieter Crisman
Ok, here we go.
sage: log_gamma(-21/10).n() NaN sage: get_systems('log_gamma(-21/10).n()') ['ginac']
So both give NaN
, but we end up using RR.log_gamma()
as in the GinacFunction
code.
sage: log_gamma(-31/10).n() 0.400311696703985 sage: log_gamma(-3.1) 0.400311696703985 sage: a = RR(5) sage: a.log_gamma() 3.17805383034795
I don't see anything holding this up except cosmetics. I'll try to make a refreshed patch momentarily.
comment:12 Changed 9 years ago by
Okay, I've been messing with this for too long today.
sage: get_systems('log_gamma(SR(6))') ['ginac', 'GMP'] sage: get_systems('log_gamma(RR(6))') [] sage: get_systems('log_gamma(CC(6))') ['PARI', 'MPFR'] sage: get_systems('log_gamma(6.)') ['MPFR']
See also #10072, where a lot of the numerical evaluation was fixed. Anyway, updated patch with more explanation and other information coming up. It needs light review; no code was changed, only doctests.
I'm not sure I like the last doctest either
sage: conjugate(log_gamma(-2)) conjugate(+Infinity)
What is the conjugate of plus infinity? But I'll leave it for now, just to document it, unless someone has an objection, since we have in vanilla Sage
sage: conjugate(+Infinity) conjugate(+Infinity)
I've opened #12521 for the evaluation at negative input with even ceiling function issue (i.e., log_gamma(-2.1)
).
comment:13 Changed 9 years ago by
- Dependencies set to #12507, #9130
I'm marking this as 'needs review', because I did change a fair number of tests. This definitely depends on #9130 because of some doc fixes. Also, I am marking this as depending on #12507, because I don't want to bother fixing that doctest if no one else is either. However, I don't really care either way.
comment:14 Changed 9 years ago by
- Reviewers changed from Karl-Dieter Crisman to Karl-Dieter Crisman, Benjamin Jones
The latest patch trac_10075.patch failed to apply on top of 5.0.beta4 with this patch queue:
trac_12507_v2.patch trac_9130-beta_function.2.patch trac_9130-py_float_segfault.take2.patch trac_9130-reviewer.patch
the failure seems to be in all.py
~/sage/latest/devel/sage> hg qpush -v applying trac_10075.patch patching file sage/functions/all.py Hunk #1 FAILED at 15 1 out of 2 hunks FAILED -- saving rejects to file sage/functions/all.py.rej
Here is sage/functions/all.py.rej
--- all.py +++ all.py @@ -16,7 +16,7 @@ from other import ( ceil, floor, gamma, psi, factorial, - abs_symbolic, erf, sqrt, + abs_symbolic, erf, sqrt, log_gamma, gamma_inc, incomplete_gamma, arg, real_part, real, imag_part, imag, imaginary, conjugate)
comment:15 Changed 9 years ago by
Makes sense, since beta is in that list now. I was not careful enough about the dependencies, I guess. Coming up.
comment:16 Changed 9 years ago by
- Description modified (diff)
Okay, all should be well now? Sorry about that.
comment:17 Changed 9 years ago by
- Status changed from needs_review to positive_review
OK! All looks good now.
The patch now applies cleanly to 5.0.beta4 on top of the patch queue in my last comment. I've tested everything in sage/functions
and sage/symbolic
and running all tests now (I don't expect any problems). The docs look good too. Positive review.
comment:18 Changed 9 years ago by
Actually, one test did fail, but I don't think it's due to this patch (right?)
File "/home/jonesbe/sage/sage-5.0.beta4/devel/sage/sage/misc/trace.py", line 61: sage: print s.before[s.before.find('-'):] Expected: ---... ipdb> c 2 * 5 Got: ---------------------------------------------------------------------- | Sage Version 5.0.beta4, Release Date: 2012-02-14 | | Type notebook() for the GUI, and license() for information. | ---------------------------------------------------------------------- ********************************************************************** * * * Warning: this is a prerelease version, and it may be unstable. * * * ********************************************************************** trace('print factor(10)'); print 3+97 s c Loading Sage library. Current Mercurial branch is: ********************************************************************** 1 items had failures: 1 of 11 in __main__.example_1 ***Test Failed*** 1 failures. For whitespace errors, see the file /home/jonesbe/.sage//tmp/trace_30044.py [2.2 s]
I haven't seen that failure before.
comment:19 Changed 9 years ago by
- Merged in set to sage-5.0.beta7
- Resolution set to fixed
- Status changed from positive_review to closed
Sage gives this error message on startup:
with this change in functions/other.py: