#2621 gamma limit error

limit (28)

Maxima 5.29.1 in ECL. Apparently the limit is actually 1, see https://groups.google.com/forum/#!topic/sage-support/y9UnyWbUagY
(%i2) y:gamma(x+1/2)/(sqrt(x)*gamma(x));
gamma(x + -)
(%o2) ----------------
sqrt(x) gamma(x)
(%i3) limit(y,x,0);
(%o3) 0


  • Peter Bruin

    Peter Bruin - 2014-05-18

    It seems to me that the limit as x -> 0 is correctly computed as 0.
    The question originally asked in the discussion linked to, however,
    was about the limit as x -> infinity.

    For any real number a, Stirling's approximation shows that the

    f: gamma(x + a)/(x^a*gamma(x));

    tends to 1 as x -> infinity. This is also what

    limit(x, a, inf);

    gives you (after asking a few questions). However, for specific
    rational but non-integral values of a the result returned by Maxima
    seems to be either 0 or inf. Besides the above example (using
    limit(y, x, inf)) one also gets, for example,

    (%i23) f: gamma(x - 2/5)/(x^(-2/5)*gamma(x));
                                    2/5           2
                                   x    gamma(x - -)
    (%o23)                         -----------------
    (%i24) limit(f,x,inf);
    (%o24)                                 0

    I tried to do a bit of debugging. It seems that the limit (in the
    case a = 1/2, say) is computed via limit(g, x, inf), where

    g: 1/sqrt(x)*exp(x*log(2*x-1)-(x-1/2)*log(x-1)-1/2)/2^x;

    as obtained by Stirling's approximation. Indeed, trying to compute
    this limit also yields infinity instead of the correct value 1.

    When we instead simplify the expression for g to

    g1: 1/sqrt(x)*exp(x*log(x-1/2)-(x-1/2)*log(x-1)-1/2);

    the limit is still computed as infinity, but this time it takes
    several minutes. I don't know if the slowness is related to the
    incorrect answer. Finally, further simplifying g1 to

    g2: exp(x*log(x-1/2)-(x-1/2)*log(x-1)-1/2-log(x)/2);

    and computing limit(g2, x, inf) does instantaneously return the
    correct answer 1.

    (All of the above is in Maxima 5.33.0 on GCL.)


    Peter Bruin

  • Peter Bruin

    Peter Bruin - 2014-05-20

    I finally discovered the reason for this bug: the series expansion
    step in Gruntz's algorithm (see his thesis [1]) is not implemented
    correctly. On page 49, it is stated that when expanding a logarithm,
    log w must be replaced by h. In Gruntz's Maple implementation
    (appendix of [1]), he does this by means of a magical declaration
    ln(W) := e3; that is somehow picked up by Series().

    In the Maxima implementation, the substitution log w -> h is only done
    after expanding the whole expression. The problem is that by time,
    the log w may have been transformed, e.g. because it appears inside an

    The patch limit-rewrite-logs.patch [see message below] fixes this,
    at least for "pure exp-log functions" (the ones for which Gruntz's
    algorithm was designed) by doing a suitable substitution on relevant
    subexpressions of the form log(f(x)) before the series expansion.

    The patch currently keeps the existing substition after the series
    expansion because (1) if the expression involves other transcendental
    functions whose Taylor expansion involves logarithms, we still want to
    transform those if we see them, and (2) not doing this substitution
    led to errors even in cases where it is apparently unnecessary,
    e.g. in the test

    integrate(x^6/(1 + x + x^2)^(15/2), x, 0, inf);

    I hope this is an acceptable fix for this bug.

    [1] D. Gruntz, On computing limits in a symbolic manipulation
    system. Ph.D. thesis, ETH Zürich, 1996,

    Last edit: Peter Bruin 2014-05-20
  • Peter Bruin

    Peter Bruin - 2014-05-20

    Here is a slightly simplified version of the patch. [Edit: simplified once more]

    Last edit: Peter Bruin 2014-05-20
  • Robert Dodier

    Robert Dodier - 2014-05-27
    • labels: --> limit
    • status: open --> closed
  • Robert Dodier

    Robert Dodier - 2014-05-27

    Fixed by commit [b0579c08a] (applied patch).


Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

No, thanks