Mathematica code for Van Wijngaarden transformation?

Question

The Van Wijngaarden adaptation of the Euler transform is used to accelerate convergence of an alternating series. An efficient algorithm for this appears in Numerical Recipes.

Unfortunately, that numerical recipe appears to be written in C language. Does someone have an equivalent encoding in Mathematica language, or alternatively can someone translate from Numerical Recipes?

Answer 1

I believe there is a more concise way to express this in Mathematica.

AltSum[f_Function, n_Integer /; n > 1] := With[
  {row = Accumulate@Table[f[k], {k, 0, n}]},
  Last@Nest[ListConvolve[{1/2, 1/2}, #] &, row, Floor[2/3 n]]
]

The only two parameters this function gets is a function f that returns the single terms of the alternating sum and the parameter n that says how many terms we want to use to approximate the sum.

Referring to the wiki page of the Van Wijngaarden transformation, the two lines of code do the following:

• The row = ... part calculates the first row as given in the table on the wiki page. It uses a simple Table to calculate a list of each term and employs Accumulate to compute the list of partial sums {s0, s0+s1, s0+s1+s2, ...}
• In the second line, we use a convolution with the kernel {1/2, 1/2} to calculate the mean of neighboring elements. This call is done over and over again using Nest and per definition of Adriaan van Wijngaarden, is not to carried out through to the very end, but it stops two-thirds of the way (Floor[2/3*n]).

The final call to Last gives then the highlighted number of the 9th row in the example of the wiki page. To use this function with the example from the wikipedia page, you can do

AltSum[Function[k, (-1)^k/(2 k + 1)], 12] // N
(* 0.785398 *)

Note that if you use only exact numbers in your function, then Mathematica will compute the result exactly which leads to very large rational expressions. If you want to use approximate numbers from the start, then simply replace one of the numbers with either a machine precision real or with a number that has a specified accuracy. Here is the log(2) example given by Mats:

AltSum[Function[k, (-1``100)^(k + 2)/(k + 1)], 100]
(* 0.693147180559945309417232121458176568075500134360314076823448273125712684836516 *)

% - Log[2]
(* 5.8822702768263632319062866821*10^-50 *)

Answer 2

Here is (what I think is) my program for computing $\log(2)$ with Van Wijngaarden:

nn = 100;
S = N[Table[Sum[(-1)^(n + 1)/n, {n, 1, k}], {k, 1, nn}], 20];
S = Table[
   Last[S = 
     Table[(S[[k]] + S[[k + 1]])/2, {k, 1, Length[S] - 1}]], {i, 1, 
    nn - 1}];
S = Table[
   Last[S = 
     Table[(S[[k]] + S[[k + 1]])/2, {k, 1, Length[S] - 1}]], {i, 1, 
    nn - 2}];
S = Table[
   Last[S = 
     Table[(S[[k]] + S[[k + 1]])/2, {k, 1, Length[S] - 1}]], {i, 1, 
    nn - 3}];
S = Table[
   Last[S = 
     Table[(S[[k]] + S[[k + 1]])/2, {k, 1, Length[S] - 1}]], {i, 1, 
    nn - 4}];
S = Table[
   Last[S = 
     Table[(S[[k]] + S[[k + 1]])/2, {k, 1, Length[S] - 1}]], {i, 1, 
    nn - 5}];
S = Table[
   Last[S = 
     Table[(S[[k]] + S[[k + 1]])/2, {k, 1, Length[S] - 1}]], {i, 1, 
    nn - 6}];
S = Table[
   Last[S = 
     Table[(S[[k]] + S[[k + 1]])/2, {k, 1, Length[S] - 1}]], {i, 1, 
    nn - 7}];
S = Table[
   Last[S = 
     Table[(S[[k]] + S[[k + 1]])/2, {k, 1, Length[S] - 1}]], {i, 1, 
    nn - 8}];
S = Table[
   Last[S = 
     Table[(S[[k]] + S[[k + 1]])/2, {k, 1, Length[S] - 1}]], {i, 1, 
    nn - 9}];
S = Table[
   Last[S = 
     Table[(S[[k]] + S[[k + 1]])/2, {k, 1, Length[S] - 1}]], {i, 1, 
    nn - 10}];
Last[S]
N[Log[2], 20]

I leave it like this so that you can make the exercise and write it in compact form as a loop or something.

Answer 3

This doesn't answer your question, but you could consider using the function NumericalMath`NSequenceLimit (which used to be just SequenceLimit) that accelerates the convergence of sequences. For halirutan's Wikipedia example:

S = Accumulate @ Table[(-1)^k/(2 k + 1), {k, 0, 10}];

NumericalMath`NSequenceLimit[S]

0.7853982

and for Mats' example:

S = Accumulate @ Table[(-1)^(n + 1)/n, {n, 1, 10}];

NumericalMath`NSequenceLimit[S]

0.6931472