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
>>> 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).
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.5066284835166987585628174193828898874378396079338
>>> # 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