Limits
======

Richardson extrapolation (``limit``)
------------------------------------

Numerical limits of functions and sequences can be computed with ``limit``, which attempts to use Richardson extrapolation to accelerate their convergence. Here is one limit representation for e and one for pi:

    >>> from mpmath import *
    >>> mp.dps = 50
    >>> print limit(lambda n: (1+1/n)**n, inf)
    2.7182818284590452353602874713526624977572470937
    >>> fac = factorial
    >>> print limit(lambda n: 2**(4*n+1)*fac(n)**4/(2*n+1)/fac(2*n)**2, inf)
    3.1415926535897932384626433832795028841971693993751

Richardson extrapolation works well here, as both results are accurate to full precision. For comparison, plugging in a large n value directly gives only a few correct digits:

    >>> n = mpf(10**6)
    >>> nprint(e - (1+1/n)**n)
    1.35914e-6
    >>> nprint(pi - 2**(4*n+1)*fac(n)**4/(2*n+1)/fac(2*n)**2)
    7.85398e-7

Here is a simple limit of a function at a finite point (x = 1):

    >>> print limit(lambda x: (x**(mpf(1)/3)-1)/(x**(mpf(1)/4)-1), 1)
    1.3333333333333333333333333333333333333333333333333

As for ``diff``, the optional ``direction`` keyword argument can be used to specify the direction from which to approach x. By default (``direction = -1``) the point is approached from above.

Richardson extrapolation is not a universal solution. Here is one limit for which it does not work so well:

    >>> mp.dps = 15
    >>> print limit(lambda n: 1 - 2**(-n), inf)
    0.999906468619675

Indeed, this limit converges so quickly that acceleration is unnecessary; when applied, it almost does more harm than good.

The value returned by ``limit`` is also likely to be misleading if the real limit is divergent. Currently, ``limit`` provides no way to automatically estimate the error, so it must be used with care (trial and error often works).

Neat examples
-------------

Computing the constant factor in Stirling's formula, without and with acceleration:

    >>> mp.dps = 50
    >>> f = lambda n: factorial(n) / (sqrt(n)*(n/e)**n)
    >>> # only good to ~7 digits, despite computing (10**6)! which has
    >>> # over 5 million digits when written out exactly
    >>> print f(10**6)
    2.5066284835166987585628174193828898874378396032067
    >>> # fast and fully accurate
    >>> print limit(f, inf)
    2.5066282746310005024157652848110452530069867406099
    >>> print sqrt(2*pi)
    2.5066282746310005024157652848110452530069867406099

Calculating Euler's constant from its definition as the limiting difference between the harmonic series and the natural logarithm:

    >>> mp.dps = 50
    >>> f = lambda n: sum([mpf(1)/k for k in range(1,n)]) - log(n)
    >>> # this is agonizingly slow and only gives 5 digits
    >>> print f(10**5)
    0.57721066489319952727326209008239846278819149864013
    >>> # fast and fully accurate
    >>> print limit(f, inf)
    0.57721566490153286060651209008240243104215933593992
    >>> print euler
    0.57721566490153286060651209008240243104215933593992