Epsilon algorithm to accelerate sum's convergence in a functional approach

Question

I am trying to learn to use Mathematica in an efficient way. Thus, I decided to spend some time on functional programming. I would like to implement the so-called epsilon algorithm, that is used to accelerate convergence of sums (see http://www-m3.ma.tum.de/m3old/bornemann/challengebook/Chapter1/chall_int.pdf p.5 for example). The task is pretty simple. As an input you have $n$ partial sums $s_k$; the algorithm gives you back a $n\times(n+1)$ array, $\epsilon$.

You start by creating a $n\times(n+1)$ array of 0s and then you fill it (taken from the link above) as follow

I did not manage to compute this double loop using a functional approach.

Failed attempt

I tried to use MapIndexed but I could not make it work since one has to change the entries as one goes on and this affects the next value.

What seems to me a better approach

The last thing that occured to me is that this double loop can certainly be seen as a recursion. However I couldn't formulate it clear enough to try to actually code it...

Thank you in advance for your help!

Edit Following Anton's comment, I looked at the Shanks method implementation. It allowed me to code a recursive function for my epsilon algorithm (eal in the following code). However, compared to my procedural implementation (epsilonAlgo), it turns out to be about twice as slow. Would anyone have a suggestion to speed it up?

Definitions

eal[m_, k_, n_, tab_] := 0 /; m > n + 2 - k;
eal[m_, 1, n_, tab_] = 0;
eal[m_, 2, n_, tab_] := eal[m, 2, n, tab] = tab[[m]];
eal[m_, k_, n_, tab_] := eal[m, k, n, tab] = eal[m + 1, k - 2,n, tab] + 1/(eal[m + 1, k - 1, n, tab] - eal[m, k - 1, n, tab])

epsilonAlgo[tab_] := (Block[{int, l},
   l = Length[tab];
   int = Table[If[j == 2, Accumulate[tab][[i]], 0], {i, l}, {j, l + 1}];
   For[i = 3, i <= l + 1, ++i, 
    For[j = 1, j <= l + 2 - i, ++j, 
     int[[j]][[i]] = int[[j + 1]][[i - 2]] + 1/(int[[j + 1]][[i - 1]] - int[[j]][[i - 1]])]];
   int]
  )

Test

max=50
list = Table[(-1)^n/n^2, {n, 1., max}];
epsilonAlgo[list][[1]][[max + 1]] // AbsoluteTiming
eal[1, max + 1, Length[list], Accumulate[list]] // AbsoluteTiming

Adrien