Function: limitnum
Section: sums
C-Name: limitnum0
Prototype: GD0,L,DGp
Help: limitnum(expr,{k = 20},{alpha=1}): numerical limit of sequence expr
using Lagrange-Zagier extrapolation; k is a multiplier so that we extrapolate
from expr(k*n). Assume u(n) ~ sum a_i n^(-alpha*i).
Doc: Lagrange-Zagier numerical extrapolation of \var{expr}, corresponding to a
sequence
$u_n$, either given by a closure \kbd{n->u(n)} or by a vector of values
I.e., assuming that $u_n$ tends to a finite limit $\ell$, try to determine
$\ell$. This routine is purely numerical and heuristic, thus may or may not
work on your examples; $k$ is ignored if $u$ is given by a vector,
and otherwise is a multiplier such that we extrapolate from $u(kn)$.
Assume that $u_n$ has an asymptotic expansion in $n^{-\alpha}$ :
$$u_n = \ell + \sum_{i\geq 1} a_i n^{-i\alpha}$$
for some $a_i$.
\bprog
? limitnum(n -> n*sin(1/n))
%1 = 1.0000000000000000000000000000000000000
? limitnum(n -> (1+1/n)^n) - exp(1)
%2 = 0.E-37
? limitnum(n -> 2^(4*n+1)*(n!)^4 / (2*n)! /(2*n+1)! )
%3 = 3.1415926535897932384626433832795028842
? Pi
%4 = 3.1415926535897932384626433832795028842
@eprog\noindent
If $u_n$ is given by a vector, it must be long enough for the extrapolation
to make sense: at least $k$ times the current \kbd{realprecision}. The
preferred format is thus a closure, although it becomes inconvenient
when $u_n$ cannot be directly computed in time polynomial in $\log n$,
for instance if it is defined as a sum or by induction. In that case,
passing a vector of values is the best option. It usually pays off to
interpolate $u(kn)$ for some $k > 1$:
\bprog
? limitnum(vector(10,n,(1+1/n)^n))
*** ^--------------------
*** limitnum: non-existent component in limitnum: index < 20
\\ at this accuracy, we must have at least 20 values
? limitnum(vector(20,n,(1+1/n)^n)) - exp(1)
%5 = -2.05... E-20
? limitnum(vector(20,n, m=10*n;(1+1/m)^m)) - exp(1) \\ better accuracy
%6 = 0.E-37
? v = vector(20); s = 0;
? for(i=1,#v, s += 1/i; v[i]= s - log(i));
? limitnum(v) - Euler
%9 = -1.6... E-19
? V = vector(200); s = 0;
? for(i=1,#V, s += 1/i; V[i]= s);
? v = vector(#V \ 10, i, V[10*i] - log(10*i));
? limitnum(v) - Euler
%13 = 6.43... E-29
@eprog
\synt{limitnum}{void *E, GEN (*u)(void *,GEN,long), long muli, GEN alpha, long prec}, where \kbd{u(E, n, prec)} must return $u(n)$ in precision \kbd{prec}.
Also available is
\fun{GEN}{limitnum0}{GEN u, long muli, GEN alpha, long prec}, where $u$
must be a vector of sufficient length as above.