Levin's u-transformation Suppose I'm given a very slowly converging sequence $\sum_k a_k$. In the literature, the Levin u-transformation is mentioned as a good universal technique for convergence acceleration. 
I have difficulty in understanding how this method actually works and what the transformation actually is. Can anyone help? 
After understanding how the method works I would like to try it out on a computer. Does anyone know any software packages which have black-boxed this process? I have access to Mathematica, Maple, MATLAB and C++. 
Thank you very much in advance.
 A: As a reminder for everybody, the general Levin transformation for a sequence of partial sums $S_n=\sum\limits_{j=0}^n a_j$ looks like this:
$$\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}$, based on Weniger's short FORTRAN routine (comments are delimited by a backslash):
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

which uses the subroutine
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

I haven't bothered to implement the stopping rules described by Fessler, Ford, and Smith in that demo program, but it's doable. Translating that short routine to your language of choice should be straightforward.
As a note, the algorithm here looks very simple, due to the exploitation of the recursive identities satisfied by the forward differences.
