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

enter image description here

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?


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]])]];


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


  • $\begingroup$ See the Shanks implementation in this answer of the discussion "Numerical evaluation of a sum" or in this answer of "Limit of partial sums involving inverse squares" $\endgroup$ – Anton Antonov Apr 25 '17 at 17:48
  • $\begingroup$ @AntonAntonov Thank you very much, it helped! However I am not sure I managed to come up with the best recursive implementation, see edited question. $\endgroup$ – A. Florio Apr 26 '17 at 9:58
  • $\begingroup$ Wynn's algorithm does not seem to me an algorithm that readily lends itself to a (readable) functional implementation. You might be interested in the implementation here, which does not need a multidimensional table. $\endgroup$ – J. M.'s discontentment May 20 '17 at 17:27
  • $\begingroup$ @J.M. Thank you for the link! $\endgroup$ – A. Florio May 25 '17 at 9:04
  • $\begingroup$ You're welcome. If you want a slightly more flexible modification of Wynn $\varepsilon$, see this. $\endgroup$ – J. M.'s discontentment May 25 '17 at 9:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.