Como un recordatorio para todos, el general Levin transformación de una secuencia de sumas parciales $S_n=\sum\limits_{j=0}^n a_j$ tiene este aspecto:
$$\mathcal{L}_k^{(n)}=\frac{\Delta^k\left((n+b)^{k-1}\frac{S_n}{g(n)}\right)}{\Delta^k\left((n+b)^{k-1}\frac1{g(n)}\right)}$$
where $\Delta^k$ is the usual $k$-th forward difference operator, $g(n)$ is an auxiliary sequence, and $b$ is an adjustable real parameter that is not a nonpositive integer; more explicitly, we have
$$\mathcal{L}_k^{(n)}=\frac{\sum\limits_{j=0}^k (-1)^j\binom{k}{j}(n+j+b)^{k-1}\frac{S_{n+j}}{g(n+j)}}{\sum\limits_{j=0}^k (-1)^j\binom{k}{j}(n+j+b)^{k-1}\frac1{g(n+j)}}$$
Often, $b$ is taken to be $1$, and $n$ is taken to be $0$. The various Levin transformations correspond to different choices of the auxiliary sequence $g(n)$; to give two of the simplest cases, the $u$-transformation, for instance, takes $g(n)=(n+b)a_n$, while the $t$-transformation uses $g(n)=a_n$.
The article you should be looking at, apart from David Levin's original paper, is E.J. Weniger's Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series; in there, he gives a short FORTRAN routine for implementing Levin's transformations. A more elaborate implementation by Fessler, Ford, and Smith is available at Netlib.
Just to show how easy it is to implement the Levin transformation, here's a demonstration program I wrote for the TI-83 Plus calculator for summing the alternating harmonic series $\sum\limits_{k=1}^\infty \frac{(-1)^{k+1}}{k}$, basado en Weniger corta rutina de FORTRAN (los comentarios están delimitados por una barra diagonal inversa):
PROGRAM:XTRPOL
Prompt N \\ number of terms of the series to use
1→P:0→S \\ P: alternating sign S: stores partial sums
For(K,1,N)
P/K→U \\ K-th term of the series, change to sum a different series
S+U→S
1→B \\ adjustable parameter for the Levin transformation
(B+K-1)U→T \\ Levin u-transform; for t version, remove the (B+K-1)
prgmLEVINT
Pause Y
-P→P
End
Y
que utiliza la subrutina
PROGRAM:LEVINT
(B+K-1)ֿ¹→W
W/T→U
If K=1
ClrList ∟DL,∟NL
U→∟DL(K)
SU→∟NL(K)
1-W→V
For(J,K-1,1,-1)
(B+J-1)W→U
∟NL(J+1)-U∟NL(J)→∟NL(J)
∟DL(J+1)-U∟DL(J)→∟DL(J)
WV→W
End
10^(99)→Y
If abs(∟DL(1))≥10^-99
∟NL(1)/∟DL(1)→Y
Y
Yo no he molestado a aplicar las reglas de detención descritas por Fessler, Ford, y Smith, en que programa de demostración, pero es factible. La traducción que corta rutina de su idioma de elección debe ser sencillo.
Como nota, el algoritmo aquí se ve muy simple, debido a la explotación de la recursivo identidades satisfecho por el avance de las diferencias.
