%
% LaTeX source file
%
\documentstyle[A4,maple,12pt]{article}
%%
%% <SHAREFILE=numerics/trans/trans.tex >
%% <DESCRIBE>
%%       Package of routines for rational function approximations
%%       to rational points, Taylor and asymptotic series, and functions
%%       Includes special numerical versions of the
%%       approximation algorithms to use the performance of the
%%       hardware floating point arithmetic.
%%       AUTHOR: Johannes Grotendorst, j.grotendorst@kfa-juelich.de
%% </DESCRIBE>
%% <UPDATE=R4update >
\begin{document}
\font\mengft=cmss9 scaled \magstep1
\def \nzo{\hbox{\mengft{I \hskip -1.7mm N}$_0$}}
\newenvironment{listvar}[3]%
{\begin{list}{}
{\settowidth{\labelwidth}{#2 #1}
\setlength{\leftmargin}{\labelwidth}
\addtolength{\leftmargin}{\labelsep}
\parsep 0.2ex plus 0.1ex minus 0.1ex
\itemsep 0.0ex #3
\renewcommand{\makelabel}[1]{#2 ##1\hfill}}}%
{\end{list}}
\newcounter{saveeqn}
\newcommand{\alpheqn}{\setcounter{saveeqn}
{\value{equation}}
\addtocounter{saveeqn}{1}
\setcounter{equation}{0}
\renewcommand{\theequation}{\mbox{\arabic{saveeqn}\alph{equation}}}}
\newcommand {\reseteqn}{\setcounter{equation}{\value{saveeqn}}
\renewcommand{\theequation}{\arabic{equation}}}
%
% title
%
\title{The {\tt trans} package for computing approximations of functions,
series expansions and sequences}
\author{J. Grotendorst\\
        Zentralinstitut f\"ur Angewandte Mathematik\\
        Forschungszentrum J\"ulich GmbH\\
        D-52425 J\"ulich\\
        Fed. Rep. Germany\\
        e-mail: j.grotendorst@kfa.juelich.de}
\date{}
\maketitle
\begin{abstract}
The {\tt trans} package is a collection of Maple procedures for
computing rational approximations of functions and formal series
as well as for transforming series (sequences) numerically.  It
includes most of the currently known transformation techniques, mainly
nonlinear methods, for increasing the rate and expanding the
domain of the convergence of series (sequences). For numerical
purposes the generation of FORTRAN functions for the resulting rational
expressions is supported.
\end {abstract}
\thispagestyle{empty}
\pagebreak
\section{Program description}
The {\tt trans} package is loaded by executing the following Maple statements:
 
\begin{maple}
 
> read `trans.m`:
 
> with(trans);
 
  [aitken, epsilon, gb,  lev, ratgen, ratser, rho, rhoit, rich, sidi,
 
   theta, thetait]
 
 
\end{maple}
\noindent
The program {\bf ratser} is designed for computing rational approximations
of functions and formal series expansions with data type ``series''.
The program {\bf ratgen} transforms generalized series expansions, e.g.
asymptotic expansions or Puisieux series, with data type \`{}+\`{}.
{\bf ratser} has the following input parameters:
\begin{listvar}{Obligatory:}{\it}{\rm}
\item[Obligatory:]
\item[f]
an expression representing the function to be approximated or a series expansion
with data type ``series''
(has to be the first parameter)
\item[x=a]
the point about which to expand in a series (if a is not specified then the expansion
is about x=0)
\item[Optional:]
\item[type]
type of transformation
\item[norder]
order of the series expansion  used for the transformation
(default: 6)
\item[expts]
list with expansion points of data type ``real constant''
for segmented approximation
(default: [ ], empty list)
\item[ws]
flag for the generation of a FORTRAN function: fort or
nofort (default: nofort)
\item[precision]
type of the FORTRAN function: single or double (default:
single)
\item[evalorder]
computational form for the rational function:
horner or confrac (default: horner)
\item[fn]
filename for the FORTRAN function program. fn has to be of data type
``string'' (default: \`{}\`{}).
\end{listvar}
The formal parameters of procedure {\bf ratgen} are:
\begin{listvar}{Obligatory:}{\it}{\rm}
\item[Obligatory:]
\item[s]
generalized series with data type \`{}+\`{}
(has to be the first parameter)
\item[Optional:]
\item[type]
type of transformation
\item[norder]
order of the series expansion  used for the transformation
(default: number of terms of the series)
\item[ws]
flag for the generation of a FORTRAN function: fort or
nofort (default: nofort)
\item[precision]
type of the FORTRAN function: single or double (default:
single)
\item[evalorder]
computational form for the rational function:
horner (default: \`{}\`{})
\item[fn]
filename for the FORTRAN function program (default:
\`{}\`{}).
\end{listvar}
\pagebreak
The order of the optional parameters in {\bf ratser} and {\bf ratgen} is
arbitrary.
In both procedures the formal input parameter {\it type}
for the type of transformation can have one of the following values:
\begin{listvar}{thetit}{\tt}{\rm}
\item[u]      {\em u} transformation
\item[t]      {\em t} transformation
\item[d]      {\em d} transformation
\item[v]      {\em v} transformation
\item[su]     {\em su} transformation
\item[st]     {\em st} transformation
\item[sd]     {\em sd} transformation
\item[sv]     {\em sv} transformation
\item[sh]     Shanks transformation
\item[gb]     Germain-Bonne transformation
\item[rh]     $\varrho$ transformation
\item[rhit]   iterated ${\varrho}_{2}$ transformation
\item[r]      Richardson extrapolation process
\item[th]     $\vartheta$ algorithm
\item[thit]   iterated $\vartheta_{2}$ algorithm
\item[at]     iterated Aitken transformation
\end{listvar}
In case of the Richardson extrapolation (r), $\varrho$ transformation
(rh) or iterated $\varrho_{2}$ transformation (rhit)
 the user can specify a list of auxiliary parameters
$x_{n},\ n\in\nzo$
(default: $x_{n} = 1/(n+1),\ n\in\nzo$,
for the Richardson extrapolation process,
and $x_{n} = n+1,\ n\in\nzo$,
for the $\varrho$ transformation and the iterated $\varrho_{2}$
transformation).
Then the first element in the list characterizes the type of the
transformation, i.e. it must be equal to \`{}r\`{}, \`{}rh\`{}
or \`{}rhit\`{}.
The main structure of {\bf ratser}
and {\bf ratgen} is as follows: First it is checked whether forbidden
actual parameters were used as input parameters.
If an error occurs at this level the procedure terminates after an
error message. Then the sequence transformation specified via the
input parameter {\it type} is performed
by calling the appropriate procedure.
Finally, the resulting expression returned by this subprocedure
is normalized and, if specified by the user,
translated into a FORTRAN function program. The FORTRAN functions
can be generated as single (7 significant decimal digits)
or double precision  (16 significant decimal digits) functions.
In order to optimize the numerical evaluation of
rational functions two computational methods are supported:
(i) numerator and denominator polynomials are evaluated using nested
Horner rules, (ii) the associated continued fraction of the
rational function is evaluated.
Furthermore, procedure {\bf ratser} supports segmented rational
approximation
to improve the quality of the approximation over the whole interval.
With both procedures
{\bf ratser} and {\bf ratgen} there is associated a
remember table, i.e. whenever the value of
one of these procedures is computed, the Maple system makes an entry
in a table ``remembering'' the value of the procedure for the
given arguments. When the procedure is called again with the same
arguments, the Maple system retrieves the result from the procedure's
remember table. In this way, it is possible to avoid redundant
computations.
\pagebreak
 
\noindent
Examples (see also the help pages ?trans,ratser and ?trans,ratgen):
\begin{maple}
 
 
> ratser (GAMMA(x), x = x0, sh, horner, 2);
 
                                                         2
  - 2 Psi(x0) GAMMA(x0) + (Psi(1, x0) GAMMA(x0) - Psi(x0)  GAMMA(x0)) (x - x0)
  ----------------------------------------------------------------------------
                                                    2
                 - 2 Psi(x0) + (Psi(1, x0) + Psi(x0) ) (x - x0)
 
 
> ratser (GAMMA(x), x = x0, sh, fort, horner, double, [1,2,3], 4);
 
      function shtr (x)
      double precision shtr, x
      if (x .le. 150.D-2) then
      shtr = (-0.5556103712266652D2-0.9875532677166135D1*x-0.15983027224
     #83155D2*(x-1.D0)**2)/(-0.1779002397678055D2-0.476465458230521D2*x+
     #0.2123507182187591D2*(x-1.D0)**2)
      elseif (x .gt. 150.D-2 .and. x .le. 250.D-2) then
      shtr = (-0.3640020514529144D1-0.7908846765315307D1*x-0.44341599449
     #87097D1*(x-2.D0)**2)/(-0.2009285390476709D2+0.3175699298036668D0*x
     #+0.3445047845095364D1*(x-2.D0)**2)
      elseif (x .gt. 250.D-2) then
      shtr = (-0.3879262513353674D0-0.1169030154778331D2*x-0.58487629401
     #94525D1*(x-3.D0)**2)/(-0.4927524366144941D2+0.1051527607136892D2*x
     #-0.1578165631096409D1*(x-3.D0)**2)
      end if
      end
 
> series (sqrt(sin(x)), x, 8);
 
               1/2         5/2           9/2            13/2      15/2
              x    - 1/12 x    + 1/1440 x    - 1/24192 x     + O(x    )
 
> ratgen (", [rh, 1, 6, 9, 16]);
 
                        1/2          5/2         9/2       13/2
                846720 x    - 48960 x    - 1212 x    - 35 x
                -----------------------------------------------
                                               2
                               846720 + 21600 x
 
> asympt (Psi (2*exp(x))-x, x, 4);
 
                                  1           1            1
                      ln(2) - -------- - ---------- + O(-------)
                              4 exp(x)            2           4
                                         48 exp(x)      exp(x)
 
> ratgen (", th);
 
                                         2
                          48 ln(2) exp(x)  - 12 exp(x) - 1
                     1/48 --------------------------------
                                             2
 
 
 
\end{maple}
 
\pagebreak
 
\noindent
The programs
\begin{maple}
 
 aitken, epsilon, gb, lev, rho, rhoit, rich, sidi, theta, thetait
 
 
\end{maple}
represent the Maple implementations of the recursive algorithms
described in section 2.
These functions can be used for the transformation of numerical
series or sequences, i.e. for the purpose of accelerating the convergence of
slowly convergent series (sequences) or the summation of divergent series.
The following list describes the input parameters of the
procedures:
\begin{listvar}{aaaaaaa}{\bf}{\it}
\item[epsilon] (f, i=n..m, arith, tseq, ndig)
\item[aitken]  (f, i=n..m, arith, tseq, ndig)
\item[thetit]  (f, i=n..m, arith, tseq, ndig)
\item[theta]   (f, i=n..m, arith, tseq, ndig)
\item[rho]     (f, g, i=n..m, arith, tseq, ndig)
\item[rhoit]   (f, g, i=n..m, arith, tseq, ndig
\item[rich]    (f, g, i=n..m, arith, tseq, ndig)
\item[gb]      (f, i=n..m, arith, tseq, ndig)
\item[lev]     (f, i=n..m, type, arith, tseq, ndig)
\item[sidi]    (f, i=n..m, type, arith, tseq, ndig)
\end{listvar}
The meaning of the formal parameters is as follows:
\begin{listvar}{aaaaaaaaa}{\it}{\rm}
\item[Obligatory:]
\item[f] an expression, procedure or array representing the sequence under
consideration
\item [i=n..m] the sequence elements ${s}_{i}$ and transformed elements
$T({s}_{i})$ are computed from $0$ to $m$ but the resulting values
are displayed over the specified range $n..m$
\item [g] an expression, procedure or array representing an auxiliary sequence
\item [type]  type of Levin or Sidi transformation : $u, t, d, v$
\medskip
\item[Optional:]
\item [arith] a flag to specify the applied arithmetic system: ($h$) the hardware
floating-point system, ($s$) the Maple floating-point system or ($hs$) using
both arithmetic systems h and s
\item [tseq]  a flag to specify wether the sequence of partial sums ($sump$) or
the original sequence should be transformed ($seqp$)
\item [ndig] number of digits used
\end{listvar}
\pagebreak
\noindent
Examples (see also the help pages ?trans,sidi and ?trans,theatit):
\begin{maple}
 
# Summation of the divergent Euler series,
# alternating n! series
#
sidi ((-1)^i*i!, 'i' = 23..29, d, hs);
 
                                   23
          23, -.2477678962998852*10  , .5963473623419048, .5963473623231885
 
                                   24
          24,  .5956716121032509*10  , .5963473617912676, .5963473623231899
 
                                   26
          25, -.1491553843122774*10  , .5963473584425768, .5963473623231924
 
                                   27
          26,  .3883759226953779*10  , .5963473433582257, .5963473623231936
 
                                   29
          27, -.1050049352772297*10  , .5963472864218570, .5963473623231940
 
                                   30
          28,  .2943878510839909*10  , .5963470979147778, .5963473623231941
 
                                   31
          29, -.8547374142655711*10  , .5963465363580098, .5963473623231941
 
#Comparison:
 
> int(exp(-t)/(1+t), t=0 .. infinity);
 
                                  - exp(1) Ei(-1)
 
> evalf (");
                                 .5963473623231941
 
# Evaluation of the logarithmic convergent series for
# the lemniscate constant
#
> h0    := proc(i) option remember; h0(i-1)*(2*i-1)/(2*i) end:
> h1    := proc(i) h0(i)/(4*i+1) end:
> h0 (0):= 1:
 
> thetait (h1, 'i' = 10..16, hs);
 
            10, 1.225016204786238, 1.311028774265025, 1.311028771841582
 
            11, 1.228753718010522, 1.311028765393593, 1.311028775087538
 
            12, 1.232043111026792, 1.311028770521508, 1.311028777134897
 
            13, 1.234967281160994, 1.311028774467085, 1.311028777142468
 
            14, 1.237589140473052, 1.311028784224084, 1.311028777144742
 
            15, 1.239957410113943, 1.311028766572657, 1.311028777146056
 
            16, 1.242110486023042, 1.311028737698854, 1.311028777146059
 
# Comparison:
 
> evalf (GAMMA (1/4) ** 2 / (4*(2*Pi)^(1/2)));
 
                                 1.311028777146059
\end{maple}
\pagebreak
\section{Recursive algorithms}
In this section we give a short description of the algorithms for
computing the various sequence transformations included in the
package trans.
%
\subsection{Wynn's $\epsilon$ algorithm}
Wynn [1] showed that the Shanks transforms $e_{k}(s_{n})$ [2] can
be computed by an economic computational scheme usually being
called the $\epsilon$ algorithm. The $\epsilon$ algorithm
proceeds as follows. Starting from the boundary values
\alpheqn
\begin{equation}
{\epsilon}_{-1}^{(n)} =  0 \ {\rm and}\  {\epsilon}_{0}^{(n)} = s_{n},
 \ n\in\nzo,
\end{equation}
the nonlinear 4-term recursion
\begin{equation}
\epsilon_{k+1}^{(n)} = \epsilon_{k-1}^{(n+1)} + \frac{1}
                       {\epsilon_{k}^{(n+1)}-\epsilon_{k}^{(n)}},
                       \ k,n\in\nzo,
\end{equation}
\reseteqn
is used to construct the 2-dimensional array
$\{\epsilon_{k}^{(n)}\}_{k,n \geq 0}$,
called $\epsilon$ table, with
\begin{equation}
\epsilon_{2k}^{(n)} = e_{k}(s_{n}), \ k,n\in\nzo,
\end{equation}
and
\begin{equation}
\epsilon_{2k+1}^{(n)} = 1/e_{k}(\Delta s_{n}),\ k,n\in\nzo.
\end{equation}
\par\noindent
The computational tableau of algorithm (10)
can be displayed as follows:
\stepcounter{equation}
\leftline{}
$$\begin{array}{ccccccccccc}
& &        \epsilon_{0}^{(0)}   &  &      &          & & & & &\\
&                 &     & \searrow &  & &        &  & & &\\
\epsilon_{-1}^{(1)} &  & \longrightarrow &  & \epsilon_{1}^{(0)}
& & & & & &\\
&       &  & \nearrow  &        & \searrow  &       & & & &\\
& & \epsilon_{0}^{(1)} & & \longrightarrow &  & \epsilon_{2}^{(0)}
& & & &\\
&    &   & \searrow &  & \nearrow   &  &\searrow & & & \\
\epsilon_{-1}^{(2)} & & \longrightarrow & & \epsilon_{1}^{(1)}&
&\longrightarrow & &\epsilon_{3}^{(0)}        & & \\
& & & \nearrow & &   \searrow &    &   \nearrow &       & &\\
& &\epsilon_{0}^{(2)} & &  \longrightarrow & & \epsilon_{2}^{(1)}&
 & & &\epsilon_{4}^{(0)}  \\
&  & &          \searrow  & & \nearrow  &  &  & & & \vdots \\
\epsilon_{-1}^{(3)} & & \longrightarrow & & \epsilon_{1}^{(2)}&
& & &\epsilon_{3}^{(1)}     & &\\
&  & &     \nearrow  &  &  & &   & \vdots     &  & \\
& & \epsilon_{0}^{(3)} &  & & & \epsilon_{2}^{(2)}& & & & \\
&     & &            &     & &  \vdots     &  & & &\\
\epsilon_{-1}^{(4)} &  & & & \epsilon_{1}^{(3)}&  & & & & &\\
\vdots   &  & &     & \vdots &          &  & & & & \\
& & \epsilon_{0}^{(4)}&  & & & & &\\
& & \vdots &  &  & & & & \\
\end{array}$$
Because of the relations in (3) and (4)
only the elements $\epsilon_{2k}^{(n)}$ with even subscripts serve as
approximations to the limit whereas the elements $\epsilon_{2k+1}^
{(n)}$ with odd subscripts are only auxiliary quantities.
Therefore, if $n$, the index of the last sequence element $s_{n}$,
is even, $n=2\mu$, our approximation to the limit of the sequence
is the transformation
\begin{equation}
\{s_{0},s_{1}, ..., s_{2\mu}\} \longrightarrow
\epsilon_{2\mu}^{(0)}=e_{\mu}(s_{0}),
\end{equation}
and if $n$ is odd, $n=2\mu + 1$, we use the transformation
\begin{equation}
\{s_{1},s_{2}, ... ,s_{2\mu +1}\} \longrightarrow
\epsilon_{2\mu}^{(1)}=e_{\mu}(s_{1}).
\end{equation}
%
\subsection{Aitken's iterated $\Delta^{2}$ process}
The special Shanks transformation
\begin{equation}
e_{1}(s_{n})=s_{n}-\frac{{(\Delta s_{n})}^{2}}{{\Delta}^{2}s_{n}}
            ={\epsilon}_{2}^{(n)}, \ n\in\nzo,
\end{equation}
represents Aitken's famous $\Delta^{2}$ process [3].
Identifying the sequence elements $s_{n}$ with the
starting values ${\cal A}_{0}^{(n)}$ and then iterating
eq. (8) yields the following nonlinear recursive scheme:
\alpheqn
\begin{eqnarray}
{\cal A}_{0}^{(n)}   & = & s_{n}, \ n\in\nzo,\\
{\cal A}_{k+1}^{(n)} & = & {\cal A}_{k}^{(n)}-
                         \frac{{[\Delta {\cal A}_{k}^{(n)}]}^{2}}
                         {\Delta^{2}{\cal A}_{k}^{(n)}}, \ k,n\in\nzo.
\end{eqnarray}
\reseteqn
If we arrange the elements ${\cal A}_{k}^{(n)}$ of the 2-dimensional
array $\{{\cal A}_{k}^{(n)}\}_{k,n \geq 0}$
in a rectangular scheme,
\stepcounter{equation}
\leftline{}
$$\begin{array}{ccccc}
{\cal A}_{0}^{(0)} & {\cal A}_{1}^{(0)} & {\cal A}_{2}^{(0)} &
{\cal A}_{3}^{(0)} & \ldots \\
 &      &   &       &   \\
{\cal A}_{0}^{(1)} & {\cal A}_{1}^{(1)} & {\cal A}_{2}^{(1)} &
{\cal A}_{3}^{(1)} & \ldots \\
&     &           &     &    \\
{\cal A}_{0}^{(2)} & {\cal A}_{1}^{(2)} & {\cal A}_{2}^{(2)} &
{\cal A}_{3}^{(2)} & \ldots \\
&    &            &     &    \\
{\cal A}_{0}^{(3)} & {\cal A}_{1}^{(3)} & {\cal A}_{2}^{(3)} &
{\cal A}_{3}^{(3)} & \ldots \\
\vdots   &  \vdots & \vdots &   \vdots & \ddots
\end{array}\eqno(\theequation)$$
then the 4 elements, which are connected by the recurrence relation
(9b), form the following pattern:
\stepcounter{equation}
\leftline{}
$$\begin{array}{cc}
{\cal A}_{k}^{(n)}   & {\cal A}_{k+1}^{(n)} \\
 &           \\
{\cal A}_{k}^{(n+1)}   &  \\
&         \\
{\cal A}_{k}^{(n+2)}   & \\
\end{array}\eqno(\theequation)$$
It follows from the recursive algorithm (9) that the computation
of ${\cal A}_{k}^{(0)}$ requires the sequence elements $s_{0},s_{1},
...s_{2k}$. Therefore, as in Wynn's $\epsilon$ algorithm,
we distinguish two distinct approximations for the limit
of the sequence depending upon the index $n$ of the last
sequence element $s_{n}$ which was used by the algorithm. For
$n=2\mu$ we use
\begin{equation}
\{s_{0},s_{1}, ..., s_{2\mu}\} \longrightarrow {\cal A}_{2\mu}^{(0)}
\end{equation}
and for $n=2\mu+1$ we use
\begin{equation}
\{s_{1},s_{2}, ... ,s_{2\mu +1}\} \longrightarrow
{\cal A}_{2\mu}^{(1)}.
\end{equation}
%
\subsection{The $\varrho$ algorithm}
The $\varrho$ algorithm looks very much like
Wynn's $\epsilon$ algorithm:
With the initial values
\alpheqn
\begin{equation}
{\varrho}_{-1}^{(n)} =  0 \ {\rm and}\  {\varrho}_{0}^{(n)} = s_{n},
 \ n\in\nzo,
\end{equation}
the nonlinear 4-term recursion
\begin{equation}
\varrho_{k+1}^{(n)} = \varrho_{k-1}^{(n+1)} + \frac{x_{n+k+1}-x_{n}}
                       {\varrho_{k}^{(n+1)}-\varrho_{k}^{(n)}},
                       \ k,n\in\nzo,
\end{equation}
\reseteqn
is used to construct the $\varrho$ table, a 2-dimensional array
$\{\varrho_{k}^{(n)}\}_{k,n \geq 0}$ of transforms.
Algorithm (14) is a rhombus rule
and yields exactly the same computational
scheme as for the $\epsilon$ algorithm.
The only difference between the $\epsilon$ algorithm and the
$\varrho$ algorithm is that the $\varrho$ algorithm also involves
a sequence of interpolation points ${\{x_{n}\}}_{n \geq 0}$.
Here also, only the even order transforms $\varrho_{2k}^{(n)}$
are meaningful as approximations to the limit. Thus, we distinguish
the two transformations:
\begin{equation}
\{s_{0}, x_{0}; s_{1}, x_{1};...; s_{2\mu}, x_{2\mu}\} \longrightarrow
\varrho_{2\mu}^{(0)}
\end{equation}
for $n$ even, $n=2\mu$, and
\begin{equation}
\{s_{1},x_{1}; s_{2},x_{2};... ;s_{2\mu +1}, x_{2\mu +1}\}\longrightarrow
\varrho_{2\mu}^{(1)}
\end{equation}
for $n$ odd, $n=2\mu+1$.
With the choice $x_{n}=n+1,\ n\in\nzo$, the $\varrho$ transformation
assumes its original form  due to Wynn [4]:
\alpheqn
\begin{equation}
{\varrho}_{-1}^{(n)}  =  0, {\varrho}_{0}^{(n)} = s_{n},\ n\in\nzo,
\end{equation}
\begin{equation}
\varrho_{k+1}^{(n)}  =  \varrho_{k-1}^{(n+1)} + \frac{k+1}
                       {\varrho_{k}^{(n+1)}-\varrho_{k}^{(n)}},
                       \ k,n\in\nzo ,
\end{equation}
\reseteqn
which is very powerful
for logarithmically convergent sequences, i.e. sequences
${\{s_{n}\}}_{n=0}^{\infty}$ satisfying
\begin{equation}
\lim_{n \rightarrow \infty} \frac{s_{n+1}-s} {s_{n}-s} = 1.
\end{equation}
%
\subsection{The iterated $\varrho_{2}$ algorithm}
First we have to express the transform $\varrho_{2}^{(n)}$
in terms of some sequence elements $s_{n}$ and interpolation
points $x_{n}$. The resulting expression for $\varrho_{2}^{(n)}$
will then be iterated. For $k=2$ we get from the $\varrho$ algorithm
(14):
\begin{equation}
\varrho_{2}^{(n)}=s_{n+1} +
\frac{(x_{n+2}-x_{n})[\Delta s_{n+1}] [\Delta s_{n}]}
{[\Delta x_{n+1}][\Delta s_{n}]-[\Delta x_{n}][\Delta s_{n+1}]},
 \ n\in\nzo.
\end{equation}
Iterating this expression yields the following recursive algorithm
(see ref. [5], p. 234):
\alpheqn
\begin{eqnarray}
{\cal W}_{0}^{(n)} & = & s_{n},  \ n\in\nzo,\\
{\cal W}_{k+1}^{(n)} & = &
\frac{(x_{n+2k+2}-x_{n})[\Delta {\cal W}_{k}^{(n+1)}]
                        [\Delta {\cal W}_{k}^{(n)}]}
{(x_{n+2k+2}-x{n})[\Delta {\cal W}_{k}^{(n)}]-
(x_{n+2k+1}-x_{n})[\Delta {\cal W}_{k}^{(n+1)}]},\nonumber\\
 &   & \ k,n\in\nzo.
\end{eqnarray}
\reseteqn
With the choice $x_{n}=n+1,\ n\in\nzo$, the ${\cal W}$ transformation
assumes its standard form:
\alpheqn
\begin{eqnarray}
{\cal W}_{0}^{(n)} & = & s_{n},  \ n\in\nzo,\\
{\cal W}_{k+1}^{(n)} & = &
\frac{(2k+2)[\Delta {\cal W}_{k}^{(n+1)}] [\Delta {\cal W}_{k}^{(n)}]}
{(2k+1){\Delta}^{2} {\cal W}_{k}^{(n)}}, \ k,n\in\nzo.
\end{eqnarray}
\reseteqn
The recurrence formulas for Aitken's iterated ${\Delta}^{2}$ process
and the ${\cal W}_{k}^{(n)}$ transforms are structually identical.
Hence, the computational scheme for the transforms ${\cal W}_{k}^{(n)}$
is the same as for the transforms ${\cal A}_{k}^{(n)}$.
The approximation
to the limit depends upon the index $n$ of the last sequence element
$s_{n}$ which was used in the recursion. We distinguish
the two transformations:
\begin{equation}
\{s_{0}, x_{0}; s_{1}, x_{1};...; s_{2\mu}, x_{2\mu}\} \longrightarrow
{\cal W}_{2\mu}^{(0)}
\end{equation}
for $n$ even, $n=2\mu$, and
\begin{equation}
\{s_{1},x_{1}; s_{2},x_{2};... ;s_{2\mu +1}, x_{2\mu +1}\}\longrightarrow
{\cal W}_{2\mu}^{(1)}
\end{equation}
for $n$ odd, $n=2\mu+1$.
%
\subsection{Brezinski's $\vartheta$ algorithm}
It is well known that Wynn's $\epsilon$ algorithm accelerates linear
convergence quite efficiently and in addition is able to sum divergent
series. However, the $\epsilon$ algorithm is not able to accelerate
logarithmic convergence. In order to have a sequence transformation
capable to accelerate  linearly as well as logarithmically convergent
sequences (series) Brezinski [6] constructed a new algorithm, called
$\vartheta$ algorithm, by modifying the $\epsilon$ algorithm in a special
way. The defining equations of the recursive $\vartheta$ scheme are:
\alpheqn
\begin{eqnarray}
{\vartheta}_{-1}^{(n)} & = &  0, {\vartheta}_{0}^{(n)} = s_{n},\\
\vartheta_{2k+1}^{(n)} & = & \vartheta_{2k-1}^{(n+1)} + 1/
                         [\Delta \vartheta_{2k}^{(n)}],\\
\vartheta_{2k+2}^{(n)} & = & \vartheta_{2k}^{(n+1)} +
                             \frac{[\Delta \vartheta_{2k}^{(n+1)}]
                            [\Delta \vartheta_{2k+1}^{(n+1)}]}
                            {{\Delta}^{2} \vartheta_{2k+1}^{(n)}},
                            \ k,n\in\nzo.
\end{eqnarray}
\reseteqn
If we arrange the elements $\vartheta_{k}^{(n)}$ in a rectangular
scheme
\stepcounter{equation}
\leftline{}
$$\begin{array}{ccccc}
{\vartheta}_{0}^{(0)} & {\vartheta}_{1}^{(0)} &
{\vartheta}_{2}^{(0)} & {\vartheta}_{3}^{(0)} & \ldots \\
 &      &   &       &   \\
{\vartheta}_{0}^{(1)} & {\vartheta}_{1}^{(1)} &
{\vartheta}_{2}^{(1)} & {\vartheta}_{3}^{(1)} & \ldots \\
&     &           &     &    \\
{\vartheta}_{0}^{(2)} & {\vartheta}_{1}^{(2)} &
{\vartheta}_{2}^{(2)} & {\vartheta}_{3}^{(2)} & \ldots \\
&    &            &     &    \\
{\vartheta}_{0}^{(3)} & {\vartheta}_{1}^{(3)} &
{\vartheta}_{2}^{(3)} & {\vartheta}_{3}^{(3)} & \ldots \\
\vdots   &  \vdots & \vdots &   \vdots & \ddots
\end{array}\eqno(\theequation)$$
then the geometric pattern, which represents the recurrence
relation (24b), has the form
\leftline{}
\stepcounter{equation}
$$\begin{array}{ccc}
& {\vartheta}_{2k}^{(n)}   & {\vartheta}_{2k+1}^{(n)} \\
 & &           \\
{\vartheta}_{2k-1}^{(n+1)}   & {\vartheta}_{2k}^{(n+1)}  \\
\end{array}\eqno(\theequation)$$
and eq. (24c) is represented by
\stepcounter{equation}
\leftline{}
$$\begin{array}{ccc}
& {\vartheta}_{2k+1}^{(n)}   & {\vartheta}_{2k+2}^{(n)} \\
 & &           \\
{\vartheta}_{2k}^{(n+1)}   & {\vartheta}_{2k+1}^{(n+1)}  \\
 & &           \\
{\vartheta}_{2k}^{(n+2)}   & {\vartheta}_{2k+1}^{(n+2)}  \\
\end{array}\eqno(\theequation)$$
As with the $\epsilon$ algorithm, only transforms
$\vartheta_{2k}^{(n)}$ with even numbered subscripts serve as
approximations. The two geometric pattern (26) and (27) and the
corresponding eqs. (24b) and (24c) show that for the
computation of the transform $\vartheta_{2k}^{(0)}$ the sequence
elements $s_{0},s_{1}, ...s_{3k}$ have to be known. Hence,
our choice of approximants has to depend upon the index $n$ of the
last sequence element $s_{n}$ which was used in the recursion.
If $n$ is a multiple of 3, $n=3\mu$, we use the transformation
\begin{equation}
\{s_{0},s_{1}, ..., s_{3\mu}\} \longrightarrow
\vartheta_{2\mu}^{(0)},
\end{equation}
if $n=3\mu + 1$, we use the transformation
\begin{equation}
\{s_{1},s_{2}, ... ,s_{3\mu +1}\} \longrightarrow
\vartheta_{2\mu}^{(1)},
\end{equation}
and if $n=3\mu + 2$, we use
\begin{equation}
\{s_{2},s_{3}, ... ,s_{3\mu +2}\} \longrightarrow
\vartheta_{2\mu}^{(2)}.
\end{equation}
%
\subsection{The iterated $\vartheta_{2}$ algorithm}
>From the recursive scheme (24) we obtain
the following expression for $\vartheta_{2}^{(n)}$:
\begin{equation}
\vartheta_{2}^{(n)}=s_{n+1} -
\frac{[\Delta s_{n}] [\Delta s_{n+1}][{\Delta}^{2} s_{n+1}]}
{[\Delta s_{n+2}][{\Delta}^{2} s_{n}]-[\Delta s_{n}]
[{\Delta}^{2} s_{n+1}]},  \ n\in\nzo.
\end{equation}
Iterating eq. (31) we get the following recursive scheme (see ref. [5],
p. 282):
\alpheqn
\begin{eqnarray}
{\cal T}_{0}^{(n)}  & = &  s_{n},  \ n\in\nzo,\\
{\cal T}_{k+1}^{(n)} & = & {\cal T}_{k}^{(n+1)} -
\frac{[\Delta {\cal T}_{k}^{(n)}] [\Delta {\cal T}_{k}^{(n+1)}]
[{\Delta}^{2} {\cal T}_{k}^{(n+1)}]}
{[\Delta {\cal T}_{k}^{(n+2)}][{\Delta}^{2} {\cal T}_{k}^{(n)}]-
[\Delta {\cal T}_{k}^{(n)}][{\Delta}^{2} {\cal T}_{k}^{(n+1)}]},
\nonumber\\
 &   & \ k,n\in\nzo.
\end{eqnarray}
\reseteqn
The 5 elements, which are connected by the recurrence relation
(32b), form the following pattern in a rectangular scheme:
\stepcounter{equation}
\leftline{}
$$\begin{array}{cc}
{\cal T}_{k}^{(n)}   & {\cal T}_{k+1}^{(n)} \\
 &           \\
{\cal T}_{k}^{(n+1)}   &  \\
&         \\
{\cal T}_{k}^{(n+2)}   & \\
&         \\
{\cal T}_{k}^{(n+3)}   & \\
\end{array}\eqno(\theequation)$$
As in the case of the $\vartheta$ algorithm the computation
of the ${\cal T}_{k}^{(n)}$ transforms require the knowlegde of
the sequence elements
$s_{n},s_{n+1}, ...s_{n+3k}$. Therefore,
we distinguish the following approximations:
\begin{equation}
\{s_{0},s_{1}, ..., s_{3\mu}\} \longrightarrow
{\cal T}_{2\mu}^{(0)},
\end{equation}
for $n=3\mu$,
\begin{equation}
\{s_{1},s_{2}, ... ,s_{3\mu +1}\} \longrightarrow
{\cal T}_{2\mu}^{(1)},
\end{equation}
for $n=3\mu + 1$ and
\begin{equation}
\{s_{2},s_{3}, ... ,s_{3\mu +2}\} \longrightarrow
{\cal T}_{2\mu}^{(2)}.
\end{equation}
for $n=3\mu + 2$.
%
\subsection{The Richardson extrapolation process}
The linear Richardson extrapolation for sequences (series) [7]
can be computed by the following variant of Neville's recursive
scheme (see p. 6 of ref. [8]):
\alpheqn
\begin{eqnarray}
{\cal N}_{0}^{(n)} & = &  s_{n},  \ n\in\nzo, \\
{\cal N}_{k+1}^{(n)} & = & \frac{x_{n} {\cal N}_{k}^{(n+1)}-
                        x_{n+k+1} {\cal N}_{k}^{(n)}}
                       {x_{n} - x_{n+k+1}}, \ k,n\in\nzo
\end{eqnarray}
\reseteqn
For the special auxiliary sequence $x_{n} = \Delta s_{n}, \ n\in\nzo$,
we obtain the sequence transformation of Germain-Bonne [9].
In this case the recursive scheme (37) assumes the form
\alpheqn
\begin{eqnarray}
{\cal B}_{0}^{(n)} & = & s_{n},  \ n\in\nzo, \\
{\cal B}_{k+1}^{(n)}& = &\frac{\Delta s_{n} {\cal B}_{k}^{(n+1)}-
                        \Delta s_{n+k+1} {\cal B}_{k}^{(n)}}
                       {\Delta s_{n} - \Delta s_{n+k+1}},
                       \ k,n\in\nzo .
\end{eqnarray}
\reseteqn
The computational tableau of the deltoid scheme (38) can be
displayed as follows:
\stepcounter{equation}
\leftline{}
$$\begin{array}{ccccccccc}
{\cal N}_{0}^{(0)}   &  &   &  & & & & &\\
& \searrow &  & &  &  & & &\\
&  & {\cal N}_{1}^{(0)} & & & & & &\\
& \nearrow  &        & \searrow  &       & & & &\\
{\cal N}_{0}^{(1)} & &   &  & {\cal N}_{2}^{(0)} & & & &\\
& \searrow &  & \nearrow   &  &\searrow & & & \\
& & {\cal N}_{1}^{(1)}& &  & &{\cal N}_{3}^{(0)}& & \\
& \nearrow & &   \searrow &    &   \nearrow & \vdots   & &\\
{\cal N}_{0}^{(2)} & &    & & {\cal N}_{2}^{(1)} &  & & & \\
&          \searrow  & & \nearrow  & \vdots &  & & & \\
& & {\cal N}_{1}^{(2)}& & & &    & &\\
&     \nearrow  & \vdots &  & &   &  &  & \\
{\cal N}_{0}^{(3)} &  & & & & & & & \\
\vdots  &  &  & &   &  &   &  & \\
\end{array}\eqno(\theequation)$$
As approximation to the limit of the sequence we use
\begin{equation}
\{s_{1},x_{1}; s_{2},x_{2};... ;s_{n}, x_{n}\} \longrightarrow
{\cal N}_{n}^{(0)},
\end{equation}
and in case of Germain-Bonne's transformation we use
\begin{equation}
\{s_{1},\Delta s_{1}; s_{2},\Delta s_{2};...;
s_{n}, \Delta s_{n}\} \longrightarrow  {\cal B}_{n}^{(0)}.
\end{equation}
%
\subsection{Recursive evaluation of the Levin transformation}
The general sequence transformation of Levin [10] is given by
the following ratio of finite sums:
\begin{equation}
{\cal L}_{k}(s_{n}) =
\frac
{\sum\limits_{j=0}^{k} (-1)^{j} {k \choose j}
(n+j+1)^{k-1}
\frac{s_{n+j}}{R_{n+j}}}
{\sum\limits_{j=0}^{k} (-1)^{j} {k \choose j}
(n+j+1)^{k-1}
\frac{1}{R_{n+j}}}.
\end{equation}
Here, ${\{R_{n}\}}_{n \geq 0}$ is a sequence of remainder estimates
of a given sequence ${\{s_{n}\}}_{n \geq 0}$.
The numerator and denominator sums in eq. (42) can be computed
recursively with the help of a linear 3-term recurrence relation
of Fessler $et \ al$. [11]
\begin{equation}
L_{k+1}^{(n)} = L_{k}^{(n+1)}-
                \frac {(n+1) {(n+k+1)}^{k-1}}{(n+k+2)^{k}}
                L_{k}^{(n)}, \ k,n\in\nzo .
\end{equation}
The starting values for this recursive formula (43) are
\begin{equation}
L_{0}^{(n)} =  s_{n}/R_{n},  \ n\in\nzo,
\end{equation}
for the numerator in eq. (42) and
\begin{equation}
L_{0}^{(n)} =  1/R_{n},  \ n\in\nzo,
\end{equation}
for the denominator in eq. (42).
The computational tableau produced by the deltoid algorithm (43),
which either represent numerators or denominators of the general
Levin transformation (42), is exactly the same as in the case
of the Richardson extrapolation process. Our approximation
to the limit of the sequence ${\{s_{n}\}}_{n \geq 0}$
will be:
\begin{equation}
\{s_{0}, R_{0}; s_{1}, R_{1};...; s_{n}, R_{n}\} \longrightarrow
{\cal L}_{n}(s_{0}).
\end{equation}
%
\subsection{Recursive evaluation of the Sidi transformation}
Formally the transformation of Sidi can be obtained
by replacing the powers $(n+j+1)^{k-1}$ in the explicit expression
(42) for Levin's general sequence transformation by the Pochhammer
symbols $(n+j+1)_{k-1}$
\begin{equation}
{\cal S}_{k}(s_{n}) =
\frac
{\sum\limits_{j=0}^{k} (-1)^{j} {k \choose j}
(n+j+1)_{k-1}
\frac{s_{n+j}}{R_{n+j}}}
{\sum\limits_{j=0}^{k} (-1)^{j} {k \choose j}
(n+j+1)_{k-1}
\frac{1}{R_{n+j}}}.
\end{equation}
Here also, ${\{R_{n}\}}_{n \geq 0}$ is a sequence of remainder
estimates of a given sequence ${\{s_{n}\}}_{n \geq 0}$.
Transformation (47) was first used by Sidi [12] for deriving Pad\'{e}
approximants of some hypergeometric functions. Later Weniger [5]
introduced the transformation (47) as a nonlinear sequence transformation
for the acceleration of convergence and the summation of divergent
series.
The recursive computation of the numerator
and denominator sums in eq. (47) is performed
using the following linear 3-term recurrence relation
of Weniger (see ref. [5], p. 259):
\begin{equation}
S_{k+1}^{(n)} = S_{k}^{(n+1)}-
                \frac {(n+k)(n+k+1)}{(n+2k+1)(n+2k)}
                S_{k}^{(n)}, \ k,n\in\nzo .
\end{equation}
The starting values for the recursive formula in eq. (48) are
\begin{equation}
S_{0}^{(n)} =  s_{n}/R_{n},  \ n\in\nzo,
\end{equation}
in case of the numerator and
\begin{equation}
S_{0}^{(n)} =  1/R_{n},  \ n\in\nzo,
\end{equation}
in case of the denominator in eq. (47).
The computaional scheme for the deltoid algorithm (51)
is the same as for the general Levin transformation.
We use the following approximation to the limit of
the sequence ${\{s_{n}\}}_{n \geq 0}$:
\begin{equation}
\{s_{0}, R_{0}; s_{1}, R_{1};...; s_{n}, R_{n}\} \longrightarrow
{\cal S}_{n}(s_{0}).
\end{equation}
\pagebreak
\parskip=0pt plus 1pt
\begin{thebibliography}{wimp}
%\bibitem{jg1} J. Grotendorst, Comput. Phys. Commun. 59 (1990) 289.
%\bibitem{cb5} C. Brezinski, Numer. Math. 35 (1980) 176.
%\bibitem{th} T. H\aa vie, BIT 19 (1979) 204.
\bibitem{pw2} P. Wynn, Math. Tables Aids Comput. 10 (1956) 91.
\bibitem{ds} D. Shanks, J. Math. and Phys.
 (Cambridge, Mass.) 34 (1955) 1.
\bibitem{aca} A. C. Aitken, Proc. Roy. Soc. Edinburgh 46 (1926) 289.
\bibitem{pw1} P. Wynn, Proc. Camb. Phil. Soc. 52 (1956) 633.
\bibitem{ejw1} E. J. Weniger, Nonlinear Sequence Transformations for the
 Acceleration of Convergence and the Summation of Divergent Series,
 Comput. Phys. Rep. 10 (1989) 189.
\bibitem{cb4} C. Brezinski, C. R. Acad. Sc. Paris 273 (1971) 727.
\bibitem{lfr} L. F. Richardson, Phil. Trans. Roy. Soc. London A 226
 (1927) 229.
\bibitem{cb2} C. Brezinski, Algorithmes d'Acc\'{e}l\'{e}ration
 de la Convergence - \'{E}tude Num\'{e}rique (\'{E}ditions
 Technip, Paris, 1978).
\bibitem{bgb} B. Germain-Bonne, Estimation de la Limite des Suites et
 Formalisation de Pr\'{e}c\'{e}des d'Acc\'{e}l\'{e}ration
 de Convergence (Thes\`{e}, Universit\'{e} de Lille, 1978).
\bibitem{dl} D. Levin, Int. J. Comput. Math. B3 (1973) 371.
\bibitem{tf} T. Fessler, W. F. Ford and D. A. Smith, ACM Trans.
 Math. Software 9 (1983) 346.
\bibitem{as} A. Sidi, J. Comput. Appl. Math. 7 (1981) 37.
%\bibitem{bwc} B. W. Char, K. O. Geddes, G. H. Gonnet, M. B. Monagan
% and S. M. Watt, Maple Reference Manual 5th ed. (Watcom Pub.,
% Waterloo, Ontario, 1988).
\end{thebibliography}
\end{document}

