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: sage-6.4
Component: calculus Keywords: maxima, gamma, limit
Cc: JGuzman Merged in:
Authors: Peter Bruin Reviewers: Karl-Dieter 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 sage-support mailing list, this is probably a bug in Maxima, and it has been reported here

Tested in sage-5.8

Change History (15)

comment:1 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:2 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:3 Changed 6 years ago by pbruin

This is bizarre. With Sage 6.2:

sage: f(x) = gamma(x+1/2)/gamma(x)/sqrt(x)
sage: limit(f,x=oo)
x |--> 0

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

sage: f(x) = gamma(x+1/2)/gamma(x)/sqrt(x)
sage: sage.calculus.calculus.maxima('domain: real')  # takes forever otherwise
real
sage: limit(f,x=oo)
x |--> +Infinity

And indeed in Maxima 5.33.0:

(%i1) f: gamma(x+1/2)/gamma(x)/sqrt(x);
                                           1
                                 gamma(x + -)
                                           2
(%o1)                          ----------------
                               sqrt(x) gamma(x)
(%i2) limit(f, x, inf);
(%o2)                                 inf

So we go from one incorrect answer to another...

Last edited 6 years ago by pbruin (previous) (diff)

comment:4 Changed 6 years ago by JGuzman

I didn't heard from the developers of Maxima so far...

comment:5 Changed 6 years ago by pbruin

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*x-1)-(x-1/2)*log(x-1)-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*x-1)-(x-1/2)*log(x-1)-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(x-1/2)-(x-1/2)*log(x-1)-1/2);

then Maxima seems to be unable to compute the limit. [Edit: after a few minutes it actually returns the wrong result `inf`.]

Last edited 6 years ago by pbruin (previous) (diff)

comment:6 follow-up: Changed 6 years ago by pbruin

I found the reason for the bug and uploaded a patch to the Maxima bug tracking page.

comment:7 in reply to: ↑ 6 ; follow-up: Changed 6 years ago by kcrisman

  • 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 ; follow-up: Changed 6 years ago by pbruin

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 ; follow-up: Changed 6 years ago by 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?

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 pbruin

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 pbruin

  • Authors set to Peter Bruin
  • Branch set to u/pbruin/15033-maxima_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 git

  • Commit changed from 0d266771007e39722defb3bab2200f36373ed834 to cab4585607523016a62f73e8e37462d05de4ced3

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

cab4585Merge branch 'develop' into ticket/15033-limit_replace_logs

comment:13 Changed 6 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:14 Changed 6 years ago by kcrisman

  • Reviewers set to Karl-Dieter 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 vbraun

  • Branch changed from u/pbruin/15033-maxima_limit_replace_logs to cab4585607523016a62f73e8e37462d05de4ced3
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.