Opened 7 years ago
Closed 6 years ago
#15033 closed defect (fixed)
Wrong limit value of expression involving gamma function
Reported by:  JGuzman  Owned by:  

Priority:  major  Milestone:  sage6.4 
Component:  calculus  Keywords:  maxima, gamma, limit 
Cc:  JGuzman  Merged in:  
Authors:  Peter Bruin  Reviewers:  KarlDieter Crisman 
Report Upstream:  Fixed upstream, but not in a stable release.  Work issues:  
Branch:  cab4585 (Commits)  Commit:  cab4585607523016a62f73e8e37462d05de4ced3 
Dependencies:  #13973, #16224  Stopgaps: 
Description
When computing the limit of the following function:
sage: f(x) = gamma(x+1/2)/gamma(x)/sqrt(x) sage: limit(f,x=oo)
Sage returns 0 in stead of one. See for example:
sage: plot(y,x,1,50)
or how the value approaches to one:
sage : for x in range(1, 20, 2): sage : print("f(%2d)=%f" %(i, f(i).n() ) )
f( 1)=0.886227 f( 3)=0.959369 f( 5)=0.975350 f( 7)=0.982316 f( 9)=0.986214 f(11)=0.988705 f(13)=0.990433 f(15)=0.991703 f(17)=0.992675 f(19)=0.993443
As described in sagesupport mailing list, this is probably a bug in Maxima, and it has been reported here
Tested in sage5.8
Change History (15)
comment:1 Changed 6 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:2 Changed 6 years ago by
 Milestone changed from sage6.2 to sage6.3
comment:3 Changed 6 years ago by
comment:4 Changed 6 years ago by
I didn't heard from the developers of Maxima so far...
comment:5 Changed 6 years ago by
After some Maxima debugging, it seems that the limit is evaluated via limit(g, x, inf)
, where g
can be defined by
g: 1/sqrt(x)*exp(x*log(2*x1)(x1/2)*log(x1)1/2)/2^x;
The fact that the two limits are equal follows from Stirling's approximation. And in fact this limit is also calculated incorrectly. In Maxima 5.33.0:
(%i1) g: 1/sqrt(x)*exp(x*log(2*x1)(x1/2)*log(x1)1/2)/2^x; x log(2 x  1)  (x  1/2) log(x  1)  1/2 %e (%o1)  x sqrt(x) 2 (%i2) limit(g, x, inf); (%o2) inf
The correct value is 1 (as it is for f
).
Another strange (but possibly related) thing is that if you simplify g
to
g1: 1/sqrt(x)*exp(x*log(x1/2)(x1/2)*log(x1)1/2);
then Maxima seems to be unable to compute the limit. [Edit: after a few minutes it actually returns the wrong result `inf`.]
comment:6 followup: ↓ 7 Changed 6 years ago by
I found the reason for the bug and uploaded a patch to the Maxima bug tracking page.
comment:7 in reply to: ↑ 6 ; followup: ↓ 8 Changed 6 years ago by
 Report Upstream changed from Reported upstream. No feedback yet. to Fixed upstream, but not in a stable release.
I found the reason for the bug and uploaded a patch to the Maxima bug tracking page.
Apparently it's just been committed!
comment:8 in reply to: ↑ 7 ; followup: ↓ 9 Changed 6 years ago by
Replying to kcrisman:
I found the reason for the bug and uploaded a patch to the Maxima bug tracking page.
Apparently it's just been committed!
Indeed; it's too late to add it to #13973, but should we maybe temporarily add this commit as a patch in our build/pkgs/maxima/patches
directory until the next time we upgrade Maxima?
comment:9 in reply to: ↑ 8 ; followup: ↓ 10 Changed 6 years ago by
I found the reason for the bug and uploaded a patch to the Maxima bug tracking page.
Apparently it's just been committed!
Indeed; it's too late to add it to #13973, but should we maybe temporarily add this commit as a patch in our
build/pkgs/maxima/patches
directory until the next time we upgrade Maxima?
I don't see why not, since it's already accepted upstream. Add a doctest and you'd be good to go.
comment:10 in reply to: ↑ 9 Changed 6 years ago by
Replying to kcrisman:
Indeed; it's too late to add it to #13973, but should we maybe temporarily add this commit as a patch in our
build/pkgs/maxima/patches
directory until the next time we upgrade Maxima?I don't see why not, since it's already accepted upstream. Add a doctest and you'd be good to go.
Hmm, it looks like the same domain: complex
slowness issue as in #13973 is going to hit us again. Trying this in Maxima 5.33.0 with the patch applied, the computation takes a few minutes with domain: complex
, but less than a second with domain: real
. (It does correctly return 1 in both cases.) Maybe we should try to find out if we can disable this domain: complex
setting in Sage for some things...
comment:11 Changed 6 years ago by
 Branch set to u/pbruin/15033maxima_limit_replace_logs
 Commit set to 0d266771007e39722defb3bab2200f36373ed834
 Dependencies set to #13973, #16224
 Status changed from new to needs_review
This branch adds the upstream patch and a doctest. The doctest temporarily switches to domain: real
to avoid the slowness caused by domain: complex
.
comment:12 Changed 6 years ago by
 Commit changed from 0d266771007e39722defb3bab2200f36373ed834 to cab4585607523016a62f73e8e37462d05de4ced3
Branch pushed to git repo; I updated commit sha1. New commits:
cab4585  Merge branch 'develop' into ticket/15033limit_replace_logs

comment:13 Changed 6 years ago by
 Milestone changed from sage6.3 to sage6.4
comment:14 Changed 6 years ago by
 Reviewers set to KarlDieter Crisman
 Status changed from needs_review to positive_review
Boy are you right about the domain:complex
slowness... Okay, this is good to go. Thanks!
comment:15 Changed 6 years ago by
 Branch changed from u/pbruin/15033maxima_limit_replace_logs to cab4585607523016a62f73e8e37462d05de4ced3
 Resolution set to fixed
 Status changed from positive_review to closed
This is bizarre. With Sage 6.2:
[Edit: this depends on the "domain: complex" setting. With "domain: real", Maxima 5.29.1 (in Sage 6.2) returns infinity, and only after asking whether 8 log(2)^2 and 8 log(2)^2 are integers.]
After #13973:
And indeed in Maxima 5.33.0:
So we go from one incorrect answer to another...