#===============================================================================
#                          Program TRANS
#===============================================================================
#
#             A Maple package for transforming sequences, series
#                                and functions
#
#
#    Copyright (C) 1990 by J Grotendorst, all rights reserved.
#
#    Error reports please to: J. Grotendorst
#                             Zentralinstitut fuer Angewandte Mathematik
#                             Forschungszentrum Juelich GmbH
#                             D-5170 Juelich
#                             Federal Republic of Germany
#                     E-mail: <ZDV071@DJUKFA11.BITNET>
#
#
#        modifications: 24 July 90: - generation of FORTRAN functions
#                                     with parameters
#                        8 Aug  90: - numerical evaluation of Maple
#                                     constants (Pi, E, ...) when
#                                     generating a FORTRAN function
#                                   - check the type `+` of the generalized
#                                     series
#                          Mar  91  - neccesary changes to use the
#                                     fortran function in MapleV
#                                   - rename and restructure of procedures
#                                     tay (now ratser) and asy (now ratgen)
#                                   - ratser applies to expansions with
#                                     data type series
#                                   - ratgen applies to generalized expansions
#                                     with data type `+`, e.g. asymptotic
#                                     expansions
#                          Oct 91   - changes due to the Maple package concept
#
#===============================================================================
#  MapleV is assumed
#===============================================================================
#
trans := 'trans':
trans := table ([ratser = proc() `trans/ratser` (args) end,
                ratgen  = proc() `trans/ratgen` (args) end,
                aitken  = proc() `trans/aitken` (args) end,
                eps     = proc() `trans/eps` (args) end,
                gb      = proc() `trans/gb` (args) end,
                lev     = proc() `trans/lev` (args) end,
                rho     = proc() `trans/rho` (args) end,
                rhoit   = proc() `trans/rhoit` (args) end,
                rich    = proc() `trans/rich` (args) end,
                sidi    = proc() `trans/sidi`(args) end,
                theta   = proc() `trans/theta` (args) end,
                thetait = proc() `trans/thetait` (args)end
                ]):
`trans/ratser` := proc (s)
#*******************************************************************************
#
#  procedure:   ratser
#
#  Description:
#
#  ratser transforms series expansions with data type "series" via
#  well-known linear as well as nonlinear sequence (series) transformations.
#  Generation of optimized FORTRAN code (functions) is supported.
#
#  Input parameters :
#
#  obligatory:
#
#  s           : series expansion of data type "series"
#
#  optional:
#
#  <typ>       : type of transformation: u, t, d, v, sh, gb, rh, r, at,
#                rhit, th, thit, su, st, sd, sv;
#                in case of Richardson extrapolation (r) or
#                rho transformation (rh) or
#                iterated rho_2 transformation (rhit)
#                a list of auxiliary parameters x_n can be specified.
#                (Default : x_n = 1/(n+1) for Richardson extrapolation
#                and x_n = n+1 for rho and rho_2 transformation).
#                Then the first element of the list characterizes
#                the type of transformation, i.e. it must be equal to
#                r, rh or rhit (Default: u)
#  <norder>    : order of transformation  (Default: number of terms of the
#                                          series expansion)
#  <expts>     : list of real expansion points for the rational approximations
#                                         (Default: [ ] )
#  <switch>    : switch for generation of a FORTRAN function: fort or nofort
#                                         (Default: nofort)
#  <precision> : type of the FORTRAN function: single or double
#                                         (Default : single)
#  <evalorder> : computational scheme: horner or confrac (Default : horner)
#
#  <fn>        : filename for the FORTRAN function (Default: `` =  terminal)
#
#*******************************************************************************
#
local norder ,typ, i, floats, x, s1, ndim1, ndim, app, resnum, resden,
      x0, x1, x2, expts, resarb, nexp, a, na, fn, exp1, expn, expk1,
      expk2, j, k, head, head1, head2, ndig, exp1m, exp1e, expnm, expne,
      expk1m, expk1e, expk2m, expk2e, echar, switch, evalorder, boolist, nd1,
      nlist, laux, typ1, a_k1, s_k1, indarg, indarg1, l, lcoeden,
      degnum, degden, appnum, appden, y, tra, old_precision;
#
#  initializations
#
expts     := [ ];
laux      := [ ];
norder    := 0;
typ       := 'u';
fn        := ``;
switch    := 'nofort';
evalorder := ``;
ndig      := 10;
echar     := `.E`;
a         := array ('sparse', 1 .. 8, []);
if assigned (precision) then
  old_precision := precision
else
  old_precision := 'precision'
fi;
precision := 'single';
#
#  check the input parameters
#
if not type (s, series) then
   ERROR (`expecting a series expansion of data type "series"`,
          `try ratgen() for generalized series with dataype ``+`` `)
elif nargs < 1 or 8 < nargs then
   ERROR (`wrong number of arguments`)
elif nargs > 1 then
   for i from 2 to nargs do
      if (args [ i ] = 'u') or (args [ i ]  = 't')
         or (args [ i ] = 'v') or (args [ i ] = 'd')
         or (args [ i ] = 'su') or (args [ i ] = 'st')
         or (args [ i ] = 'sd') or (args [ i ] = 'sv')
         or (args [ i ] = 'sh')
         or (args [ i ] = 'gb') or (args [ i ] = 'rh')
         or (args [ i ] = 'r') or (args [ i ] = 'at')
         or (args [ i ] = 'rhit') or (args [ i ] = 'th')
         or (args [ i ] = 'thit')
         and type (args [ i ], name) then
         if a[3] = 0 then
            a[1] := 1;
         fi;
         typ  := args [ i ];
      elif type (args [ i ], integer) and args [ i ] > 0 then
         a[2] := 1;
         norder := args [ i ];
      elif type (args [ i ], list) then
         nlist := nops (args [ i ]);
         if nlist > 0 then
            typ1 := op (1, args [ i ]);
            if typ1 = 'r' or typ1 = 'rh' or typ1 = 'rhit' then
               typ  := op (1, args [ i ]);
               if a[1] = 0 then
                  a[3] := 1;
               fi;
               laux := [seq(op (j, args [ i ]), j = 2 .. nlist)];
            else
               a[4]    := 1;
               boolist := map (type, args [ i ], realcons);
               x  := op (0, s);
               x1 := op (1, x);
               x2 := solve (x, x1);
               if type (x2, realcons) then
                  ERROR (`expansion point of the series must be indeterminate`)
               elif has (boolist, false) then
                  ERROR (`expansion points have to be of type real constant`)
               fi;
               expts := args [ i ];
            fi;
         fi;
      elif (args [ i ] = 'fort') or (args [ i ] = 'nofort')
         and type (args [ i ], name) then
         a[5]   := 1;
         switch := args [ i ];
      elif (args [ i ] = 'single') or (args [ i ] = 'double')
         and type (args [ i ], name) then
         a[6]      := 1;
         precision := args [ i ];
         if (precision = 'double') then
            ndig      := 20;
            echar     := `.D`;
         fi;
      elif (args [ i ] = `horner`) or (args [ i ] = `confrac`)
         and type (args [ i ], name) then
         a[7]      := 1;
         evalorder := args [ i ];
      elif type (args [ i ], name) then
         a[8] := 1;
         fn   := args [ i ];
      else
         ERROR (`invalid argument found`)
      fi;
   od;
#
#  check the number of input arguments
#
   na := sum (a ['i'], 'i' = 1 .. 8) + 1;
#
   if na <> nargs then
      ERROR (`invalid argument found`)
   fi;
fi;
#
# end of input check
#
# eliminate the order symbol
#
s1 := subs (O(1) = 0, s);
#
# determination of the expansion point
#
x := op (0, s1);
#
# check whether the series contains coefficients with data type "float"
#
if hastype (s1, float) or type (x, float) then
   x0 := x;
   x  := map (convert, x, rational);
   s1 := subs (x0 = x, s1);
   s1 := map (convert, s1, rational);
   floats := true;
else
   floats := false;
fi;
#
# determination of the independent variable
#
x1 := op (1, x);
x2 := solve (x, x1);
#
# determination of the series order
#
ndim1 := nops (s1) / 2 - 1;
#
if norder = 0 then
   ndim := ndim1;
else
   if norder <= ndim1 then
      ndim := norder
   else
      ERROR (`series order too small for specified transformation order`)
   fi;
fi;
#
# check the transformation order and parameters
#
if typ = 'gb' and ndim < ndim1 then
   nd1 := ndim + 1
else
   nd1 := ndim
fi;
#
if typ = 'rh' or typ = 'rhit' or typ = 'r' then
   if nops (laux) = 0 then
      if typ = 'rh' or typ = 'rhit' then
         laux := [seq (k + 1, k = 0 .. nd1)]
      elif typ = 'r' then
         laux := [seq (1 / (k + 1), k = 0 .. nd1)]
      fi
   elif nops (laux) > 0 and nops (laux) <= nd1 then
      if typ = 'rh' then
         ERROR (`too few parameters for the rho transformation`)
      elif typ = 'rhit' then
         ERROR (`too few parameters for the iterated rho_2 transformation`)
      elif typ = 'r' then
         ERROR (`too few parameters for the Richardson exptrapolation`)
      fi
   fi
fi;
#
# perform the specific series transformation
#
app := frontend (`trans/series_transform`, [s1, y, ndim, typ, laux],
                 [{`+`, `*`, list, series},{}]);
#
# normalization of the rational approximation
#
appnum  := expand (numer (app));
appden  := expand (denom (app));
degnum  := degree (appnum, y);
degden  := degree (appden, y);
lcoeden := tcoeff (appden, y);
#
if type (lcoeden, numeric) then
   resnum  := series (appnum / lcoeden, y, degnum + 1);
   resden  := series (appden / lcoeden, y, degden + 1)
else
   resnum  := series (appnum, y, degnum + 1);
   resden  := series (appden, y, degden + 1)
fi;
#
# indeterminants for the output expression
#
nexp := nops (expts);
if type (x2, name) and nexp > 0 then
   indarg  := indets (app, string) minus {constants, y, x2}
else
   indarg  := indets (app, string) minus {constants, y}
fi;
#
indarg1 := [x1, seq (l, l = indarg minus {x1})];
#
resnum := convert (subs (op (0, resnum) = x, resnum), polynom);
resden := convert (subs (op (0, resden) = x, resden), polynom);
#
if evalorder = `horner` then
   resnum  := convert (resnum, horner, indarg1);
   resden  := convert (resden, horner, indarg1);
   resarb  := resnum / resden
elif evalorder = `confrac` then
   resarb  := convert (resnum / resden, confrac, x1)
else
   resarb  := resnum / resden
fi;
#
# return the result
#
if switch = 'nofort' then
#
#  check whether the result contains a floating point number
#
   if floats then
      resarb := evalf (resarb)
   fi;
#  return a sequence of approximations
#
   if nexp > 0 then
      resarb := seq (subs (x2 = expts [ k ], resarb), k = 1 .. nexp)
   fi;
#
#  insert the actual arguments as table index and the result as
#  table value into the remember table
#
   `trans/ratser` (args) := resarb;
#
   precision := old_precision;
#
   RETURN (resarb)
else
#
#  determination of the fortran function header
#
   `trans/fortran_header` (typ, indarg1, precision, 'head', 'tra');
   head1 := head [1];
   head2 := head [2];
#
#  write the fortran function
#
   if nexp = 1 then
      if (fn <> ``) then writeto(fn) fi;
      lprint(`      `.head1);
      lprint(`      `.head2);
      `trans/my_fort`([tra = evalf (subs (x2 = expts [ 1 ], resarb), ndig)]);
      lprint(`      end`);
      if (fn <> ``) then writeto('terminal') fi;
   elif nexp > 1 then
      exp1  := evalf (expts [ 1 ] + (expts [ 2 ] -
               expts [ 1 ]) / 2, 3);
      expn  := evalf (expts [ nexp - 1 ] + (expts [ nexp ] -\
               expts [ nexp - 1 ]) / 2, 3);
      exp1m := op (1, exp1);
      expnm := op (1, expn);
#
      if exp1 = 0 then
         exp1e := 0
      else
         exp1e := op (2, exp1)
      fi;
#
      if expn = 0 then
         expne := 0
      else
         expne := op (2, expn)
      fi;
#
      if nexp = 2 then
         if (fn <> ``) then writeto(fn) fi;
         lprint(`      `.head1);
         lprint(`      `.head2);
         lprint(`      if (`.x1.` .le. `.exp1m.echar.exp1e.`) then `);
         `trans/my_fort`([tra = evalf (subs (x2 = expts [ 1 ], resarb), ndig)]);
         lprint(`      elseif (`.x1.` .gt. `.expnm.echar.expne.`) then `);
         `trans/my_fort`([tra = evalf (subs (x2 = expts [ nexp ], resarb),
                          ndig)]);
         lprint(`      end if`);
         lprint(`      end`);
         if (fn <> ``) then writeto('terminal') fi;
      elif nexp > 2 then
         if (fn <> ``) then writeto(fn) fi;
         lprint(`      `.head1);
         lprint(`      `.head2);
         lprint(`      if (`.x1.` .le. `.exp1m.echar.exp1e.`) then `);
         `trans/my_fort`([tra = evalf (subs (x2 = expts [ 1 ], resarb), ndig)]);
         for k from 2 to (nexp - 1) do
            expk1  := evalf (expts [ k - 1 ] + (expts [ k ] -\
                      expts [ k - 1 ]) / 2, 3):
            expk2  := evalf (expts [ k ] + (expts [ k + 1 ] -\
                      expts [ k ]) / 2, 3):
            expk1m := op (1, expk1):
            expk2m := op (1, expk2):
#
            if expk1 = 0 then
               expk1e := 0
            else
               expk1e := op (2, expk1)
            fi:
#
            if expk2 = 0 then
               expk2e := 0
            else
               expk2e := op (2, expk2)
            fi:
#
            lprint(`      elseif (`.x1.` .gt. `.expk1m.echar.
            expk1e.` .and. `.x1.` .le. `.expk2m.echar.expk2e.`) then`);
            `trans/my_fort`([tra = evalf (subs (x2 = expts [ k ], resarb),
                           ndig)])
         od;
         lprint(`      elseif (`.x1.` .gt. `.expnm.echar. expne.`) then `);
         `trans/my_fort`([tra = evalf (subs (x2 = expts [ nexp ], resarb),
                 ndig)]);
         lprint(`      end if`);
         lprint(`      end`);
         if (fn <> ``) then writeto('terminal') fi;
      fi
   else
      if (fn <> ``) then writeto(fn) fi;
      lprint(`      `.head1);
      lprint(`      `.head2);
      `trans/my_fort`([tra = evalf (resarb, ndig)]);
      lprint(`      end`);
      if (fn <> ``) then writeto('terminal') fi
   fi;
#
   precision := old_precision;
#
   RETURN ()
fi
end:
`trans/fortran_header` := proc (typ, parlist, prec, head12, nam)
#******************************************************************************
#
#  description:
#
#  This procedure evaluates the FORTRAN function header
#
#  input parameters:
#
#  typ     : type of transformation
#  parlist : list with the formal parameters for the fortran function
#  prec    : precision type of the fortran function
#  head12  : line 1 and 2 of the fortran function
#  nam     : name of the fortran function
#
#******************************************************************************
#
local parg, j, nam1, namtr;
#
# function name
#
namtr :='tr';
nam1  := typ.namtr;
#
parg   := ``;
for j from 1 to nops (parlist)
do
   b.j := op (j, parlist);
   if j = nops (parlist) then
      parg := cat (parg, b.j)
   else
      parg := cat (parg, b.j, `, `)
   fi;
od;
#
head12 [1] := `function `.nam1.` (`.parg.`)`;
#
# declaration for the parameters
#
if prec = 'single' then
   head12 [2] := `real `.nam1.`, `.parg
elif prec = 'double' then
   head12 [2] := `double precision `.nam1.`, `.parg
fi;
#
nam := nam1;
#
end:
`trans/series_transform` := proc (s, ep, nterm, typ, param)
#******************************************************************************
#
# description:
#
# This procedure performs the sequence transformation by calling
# the specific transformation procedure seccessively.
#
# input parameters:
#
# s    :  series expansion of type "series" or `+`
# ep   :  variable for "x - x0"
# nterm:  number of series terms
# type :  type of transformation
# param:  list with parameters for the rho, rhoit or r transformation
#
#******************************************************************************
#
local s_n, a_n, x_n, n, arr, xarr, ratapp;
#
# initialization of the first partial sum
#
s_n := 0;
#
#  loop over the sequence of partial sums
#
for n from 0 to nterm do
   if type (s, series) then
      a_n := op (2 * n + 1, s) * ep ** op (2 * n + 2, s)
   elif type (s, `+`) then
      a_n := op (n + 1, s)
   else
      ERROR (`series type is wrong`)
   fi;
#
   s_n := s_n + a_n;
#
#  call the algorithm to perform the specified transformation
#
   if typ = 'rh' then
      x_n := op (n + 1, param);
      ratapp  := `trans/rho` (s_n, x_n, n, arr, xarr)
   elif typ = 'sh' then
      ratapp  := `trans/eps` (s_n, n, arr)
   elif typ = 'u' or typ = 't' or typ = 'd' or typ = 'v' then
      ratapp  := `trans/lev` (s_n, a_n, n, arr, xarr, typ)
   elif typ = 'su' or typ = 'st' or typ = 'sd' or typ = 'sv' then
      ratapp  := `trans/sidi` (s_n, a_n, n, arr, xarr, typ)
   elif typ = 'r' then
      x_n := op (n + 1, param);
      ratapp  := `trans/rich` (s_n, x_n, n, arr, xarr)
   elif typ = 'gb' then
      ratapp := `trans/gb` (s_n, a_n, n, arr, xarr)
   elif typ = 'at' then
      ratapp  := `trans/aitken` (s_n, n, arr)
   elif typ = 'rhit' then
      x_n := op (n + 1, param);
      ratapp  := `trans/rhoit` (s_n, x_n, n, arr, xarr)
   elif typ = 'th' then
      ratapp  := `trans/theta` (s_n, n, arr, xarr)
   elif typ = 'thit' then
      ratapp  := `trans/thetait` (s_n, n, arr)
   fi;
od;
#
RETURN (ratapp);
end:
`trans/ratgen` := proc (s)
#*******************************************************************************
#
#  procedure:   ratgen
#
#  Description:
#
#  ratgen transforms generalized series expansions, e.g. asymptotic
#  expansions, with data type `+` via well-known linear as well as
#  nonlinear sequence (series) transformations. Generation of
#  optimized FORTRAN code (functions) is supported.
#
#  Input parameters :
#
#  obligatory:
#
#  s           : generalized series of dataype `+`
#
#  optional:
#
#  <typ>       : type of transformation: u, t, d, v, sh, gb, rh, r, at,
#                rhit, th ,thit, su, st, sd, sv;
#                in case of Richardson extrapolation (r) or
#                rho transformation (rh) or
#                iterated rho_2 transformation (rhit)
#                a list of auxiliary parameters x_n can be specified.
#                (Default : x_n = 1/(n+1) for Richardson extrapolation
#                and x_n = n+1 for rho and rho_2 transformation).
#                Then the first element of the list characterizes
#                the type of transformation, i.e. it must be equal to
#                r, rh or rhit (Default: u)
#  <norder>    : order of transformation  (Default: number of terms of the
#                                         generalized series)
#  <switch>    : switch for generation of a FORTRAN function: fort or nofort
#                                         (Default: nofort)
#  <precision> : type of the FORTRAN function: single or double
#                                         (Default: single)
#  <evalorder> : computational scheme: horner (Default: ``)
#
#  <fn>        : filename for FORTRAN function (Default: `` = terminal)
#
#*******************************************************************************
#
local s1, norder ,typ, typ1, i, j, k, floats, ndim1, ndim, app, resnum, resden,
      resarb, a, na, fn, ndig, switch, evalorder, lterm, nd1, nlist, laux,
      a_k1, s_k1, indarg, y, head, head1, head2, tra, old_precision;
#
# initializations
#
laux      := [ ];
norder    := 0;
typ       := 'u';
fn        := ``;
switch    := 'nofort';
evalorder := ``;
ndig      := 10;
a         := array ('sparse', 1 .. 7, []);
if assigned (precision) then
  old_precision := precision
else
  old_precision := 'precision'
fi;
precision := 'single';
#
#  check the number of the arguments and the type `+` (sum-of-products form)
#  of the generalized series
#
if not type (s, `+`) then
   ERROR (`expecting a series expansion of data type ``+`` `,
          `try ratser() for series expansions of data type series`)
elif nargs < 1 or 7 < nargs then
   ERROR (`wrong number of arguments`)
fi;
#
# eliminate the order symbol
#
lterm := op (nops(s),s);
if type(lterm, function) then
   if op(0, lterm) = O then
      s1 := subs (lterm = 0, s)
   fi;
else
   s1 := s
fi;
#
# check the other input parameters
#
if nargs > 1 then
   for i from 2 to nargs do
      if (args [ i ] = 'u') or (args [ i ]  = 't')
         or (args [ i ] = 'd') or (args [ i ] = 'v')
         or (args [ i ] = 'su') or (args [ i ] = 'st')
         or (args [ i ] = 'sd') or (args [ i ] = 'sv')
         or (args [ i ] = 'sh')
         or (args [ i ] = 'gb') or (args [ i ] = 'rh')
         or (args [ i ] = 'r') or (args [ i ] = 'at')
         or (args [ i ] = 'rhit') or (args [ i ] = 'th')
         or (args [ i ] = 'thit')
         and type (args [ i ], name) then
         if a[3] = 0 then
            a[1] := 1;
         fi;
         typ  := args [ i ];
      elif type (args [ i ], integer) and args [ i ] > 0 then
         a[2] := 1;
         norder := args [ i ];
      elif type (args [ i ], list) then
         nlist := nops (args [ i ]);
         if nlist > 0 then
            typ1 := op (1, args [ i ]);
            if typ1 = 'r' or typ1 = 'rh' or typ1 = 'rhit' then
               typ  := op (1, args [ i ]);
               if a[1] = 0 then
                  a[3] := 1;
               fi;
               laux := [seq (op (j, args [ i ]), j = 2 .. nlist)];
            fi;
         fi;
      elif (args [ i ] = 'fort') or (args [ i ] = 'nofort')
         and type (args [ i ], name) then
         a[4]   := 1;
         switch := args [ i ];
      elif (args [ i ] = 'single') or (args [ i ] = 'double')
         and type (args [ i ], name) then
         a[5]   := 1;
         precision := args [ i ];
         if (precision = 'double') then
            ndig      := 20;
         fi;
      elif (args [ i ] = `horner`)
         and type (args [ i ], name) then
         a[6]        := 1;
         evalorder   := args [ i ];
      elif type (args [ i ], name) then
         a[7] := 1;
         fn   := args [ i ];
      else
         ERROR (`invalid argument found`)
      fi;
   od;
#
   na := sum (a ['i'] ,'i' = 1 .. 7) + 1;
#
#  check the number of input parameters
#
   if na <> nargs then
      ERROR (`invalid argument found`)
   fi;
fi;
#
# check whether the series contains coefficients with data type float
#
if hastype (s1, float) then
   s1     := map (convert, s1, rational);
   floats := true;
else
   floats := false;
fi;
#
ndim1 := nops (s1) - 1;
#
if norder = 0 then
   ndim := ndim1;
else
   if norder <= ndim1 then
      ndim := norder;
   else
      ERROR (`series order too small for specified transformation order`)
   fi;
fi;
#
# check the transformation order and parameters
#
if typ = 'gb' and ndim < ndim1 then
   nd1 := ndim + 1
else
   nd1 := ndim
fi;
#
if typ = 'rh' or typ = 'rhit' or typ = 'r' then
   if nops (laux) = 0 then
      if typ = 'rh' or typ = 'rhit' then
         laux := [seq (k + 1, k = 0 .. nd1)]
      elif typ = 'r' then
         laux := [seq (1 / (k + 1), k = 0 .. nd1)]
      fi
   elif nops (laux) > 0 and nops (laux) <= nd1 then
      if typ = 'rh' then
         ERROR (`too few parameters for the rho transformation`)
      elif typ = 'rhit' then
         ERROR (`too few parameters for the iterated rho_2 transformation`)
      elif typ = 'r' then
         ERROR (`too few parameters for the Richardson exptrapolation`)
      fi
   fi
fi;
#
#  perform the specific series transformation
#
app := frontend (`trans/series_transform`, [s1, y, ndim, typ, laux],
                 [{`+`, `*`, list, series},{}]);
#
#  determination of the output form
#
resnum := expand (numer (app));
resden := expand (denom (app));
#
#  determination of the indeterminats for the output expresssion
#
indarg  := indets (app, string) minus {constants};
#
#  convert into Horner scheme
#
if evalorder = `horner` then
   resnum  := convert (resnum, horner, indarg);
   resden  := convert (resden, horner, indarg);
   resarb  := resnum / resden
else
   resarb := resnum / resden
fi;
#
# return the result
#
if switch = 'nofort' then
#
#  check whether the result contains a floating-point number
#
   if floats then
      resarb := evalf (resarb)
   fi;
#
#  insert the actual arguments as table index and the result as
#  table value into the remember table
#
   `trans/ratgen` (args) := resarb;
#
   precision := old_precision;
#
   RETURN (resarb)
else
#
#  determination of the FORTRAN function header
#
   `trans/fortran_header` (typ, indarg, precision, 'head', 'tra');
   head1 := head [1];
   head2 := head [2];
#
#  write the fortran function
#
   if (fn <> ``) then writeto(fn) fi;
   lprint(`      `.head1);
   lprint(`      `.head2);
   `trans/my_fort`([tra = evalf (resarb, ndig)]);
   lprint(`      end`);
   if (fn <> ``) then writeto('terminal') fi;
#
   precision := old_precision;
#
   RETURN()
fi
end:
`trans/eps` := proc (s_k, k, arr)
#***********************************************************************
#
#  Procedure  : eps
#
#  Programmer : Johannes Grotendorst
#
#  Date       : Juelich, 18 june 1990
#  Update     : 1 july 1990
#               normalize intermetiate expressions
#
#-----------------------------------------------------------------------
#
#  Description:
#
#  The epsilon algorithm due to Wynn is used to perform Shanks'
#  series transformation.
#
#  Input parameters :
#
#  s_k  : sequence element s_k of the sequence to be accelerated
#  k    : number of the last element
#  arr  : 1-dimensional array to store the actual counterdiagonal of
#         the epsilon table
#
#***********************************************************************
#
local j, aux1, aux2, diff1;
#
# initialization
#
arr [ k ]    := s_k;
#
# recursive computation of the epsilon table
#
if k = 0 then
   RETURN (s_k)
else
   aux2 := 0;
   for j from k by -1 to 1
   do
      aux1          := aux2;
      aux2          := arr [ j - 1 ];
      diff1         := normal (arr [ j ] - aux2);
      arr [ j - 1 ] := normal (aux1 + 1 / diff1)
   od
fi;
#
if modp (k, 2) = 0 then
   RETURN (arr [ 0 ])
else
   RETURN (arr [ 1 ])
fi
end:
`trans/rho` := proc (s_k, x_k, k, arr, xarr)
#***********************************************************************
#
#  Procedure  : rho
#
#  Programmer : Johannes Grotendorst
#
#  Date       : Juelich, 18 june 1990
#  Update     : 1 july 1990
#               normalize intermetiate expressions
#
#-----------------------------------------------------------------------
#
#  Description:
#
#  The rho algorithm due to Wynn is used to perform the rho
#  transformation.
#
#  Input parameters :
#
#  s_k  : sequence element s_k of the sequence to be accelerated
#  x_k  : element x_k of the auxiliary sequence
#  k    : number of the last element
#  arr  : 1-dimensional array to store the actual counterdiagonal of
#         the rho table
#  xarr : 1-dimensional array containing the auxiliary sequence x_k
#         required in the generalized version of the rho algorithm
#
#***********************************************************************
#
local j, aux1, aux2, diff1;
#
# initializations
#
arr  [ k ] := s_k;
xarr [ k ] := x_k;
#
# recursive computation of the rho table.
#
if k = 0 then
   RETURN (s_k)
else
   aux2 := 0;
   for j from k by -1 to 1
   do
      aux1          := aux2;
      aux2          := arr [ j - 1 ];
      diff1         := normal (arr [ j ] - aux2);
      arr [ j - 1 ] := normal (aux1 + (xarr [ k ] - xarr [ j - 1 ]) / diff1)
   od;
   if modp (k, 2) = 0 then
      RETURN (arr [ 0 ])
   else
      RETURN (arr [ 1 ])
   fi
fi
end:
`trans/rhoit` := proc (s_k, x_k, k, arr, xarr)
#***********************************************************************
#
#  Procedure  : rhoit
#
#  Programmer : Johannes Grotendorst
#
#  Date       : Juelich, 18 june 1990
#  Update     : 1 july 1990
#               normalize intermetiate expressions
#
#-----------------------------------------------------------------------
#
#  Description:
#
#  The iterated rho_2 algorithm is used to perform the series
#  transformation.
#
#  Input parameters :
#
#  s_k  : sequence element s_k of the sequence to be accelerated
#  x_k  : element x_k of the auxiliary sequence
#  k    : number of the last element
#  arr  : 1-dimensional array to store the actual counterdiagonal of
#         the table for the iterated rho_2 transformation
#  xar  : 1-dimensional array containing the auxiliary sequence x_k
#
#***********************************************************************
#
local j , diff1, diff2, m, dnom;
#
# initializations
#
arr  [ k ] := s_k;
xarr [ k ] := x_k;
#
# perform the iterated rho_2 algorithm
#
if k < 2 then
   RETURN (s_k)
else
   m := k;
   for j from 1 to iquo (k, 2)
   do
      m         := m - 2;
      diff1     := normal (arr [ m + 1 ] - arr [ m ]);
      diff2     := normal (arr [ m + 2 ] - arr [ m + 1 ]);
      dnom      := normal ((xarr [ k ] - xarr [ m + 1 ]) * diff1 -
                   (xarr [ k - 1 ] - xarr [ m ]) * diff2);
      arr [ m ] := normal (arr [ m + 1 ] + (xarr [ k ] - xarr [ m ]) *
                   diff1 * diff2 / dnom)
   od;
   if modp (k, 2) = 0 then
      RETURN (arr [ 0 ])
   else
      RETURN (arr [ 1 ])
   fi
fi
end:
`trans/theta` := proc (s_k, k, arr, xarr)
#***********************************************************************
#
#  Procedure  : theta
#
#  Programmer : Johannes Grotendorst
#
#  Date       : Juelich, 18 june 1990
#  Update     : 1 july 1990
#               normalize intermetiate expressions
#
#-----------------------------------------------------------------------
#
#  Description:
#
#  The theta algorithm due to Brezinski is used to perform the theta
#  transformation.
#
#  Input parameters :
#
#  s_k  : sequence element s_k of the sequence to be accelerated
#  k    : number of the last element
#  arr  : 1-dimensional array
#  xarr : 1-dimensional array
#
#***********************************************************************
#
local j, aux1, aux2, aux3, dnom, diff1, jmax;
#
# recursive computation of the theta table
#
if k = 0 then
   arr [ k ] := s_k;
   RETURN (s_k)
elif k > 0 then
   jmax := iquo (2 * k + 1, 3);
   if (modp (k, 2) = 0) then
      aux2      := 0;
      aux1      := arr [ 0 ];
      arr [ 0 ] := s_k;
      for j from 1 to jmax
      do
         aux3 := aux2;
         aux2 := aux1;
         if j < jmax then
            aux1 := arr [ j ]
         fi;
         if (modp (j, 2) = 0) then
            dnom        := normal (arr [ j - 1 ] - 2 *
                           xarr [ j - 1 ] + aux2);
            arr [ j ] := normal (aux3 + (xarr [ j - 2 ] - aux3) *
                         (arr [ j - 1 ] - xarr [ j - 1 ]) / dnom)
         else
            diff1     := normal (arr [ j - 1 ] - xarr [ j - 1 ]);
            arr [ j ] := normal (aux3 + 1 / diff1)
         fi
      od;
#
      if (modp (jmax, 2) = 0) then
         RETURN (arr [ jmax ])
      else
         RETURN (arr [ jmax - 1 ])
      fi
   else
      aux2       := 0;
      aux1       := xarr [ 0 ];
      xarr [ 0 ] := s_k;
      for j from 1 to jmax
      do
         aux3 := aux2;
         aux2 := aux1;
         if j < jmax then
            aux1 := xarr [ j ]
         fi;
         if (modp (j, 2) = 0) then
            dnom         := normal (xarr [ j - 1 ] - 2 *
                            arr [ j - 1 ] + aux2);
            xarr [ j ] := normal (aux3 + (arr [ j - 2 ] - aux3) *
                          (xarr [ j - 1 ] - arr [ j - 1 ]) / dnom)
         else
            diff1      := normal (xarr [ j - 1 ] - arr [ j - 1 ]);
            xarr [ j ] := normal (aux3 + 1 / diff1)
         fi
      od;
#
      if (modp (jmax, 2) = 0) then
         RETURN (xarr [ jmax ])
      else
         RETURN (xarr [ jmax - 1 ])
      fi
   fi
fi
end:
`trans/thetait` := proc (s_k, k, arr)
#***********************************************************************
#
#  Procedure  : thetait
#
#  Programmer : Johannes Grotendorst
#
#  Date       : Juelich, 18 june 1990
#  Update     : 1 july 1990
#               normalize intermetiate expressions
#
#-----------------------------------------------------------------------
#
#  Description:
#
#  The iterated theta_2 algorithm is used to perform the sequence
#  (series) transformation.
#
#  Input parameters :
#
#  s_k  : sequence element s_k of the sequence to be accelerated
#  k    : number of the last element
#  arr  : 1-dimensional array to store the actual counterdiagonal of
#         theta_2 table
#
#***********************************************************************
#
local j, m, dnom, diff0, diff1, diff2;
#
# initialization
#
arr [ k ] := s_k;
#
# recursive computation of the theta_2 table
#
if k < 3 then
   RETURN (s_k);
else
   m := k;
   for j from 1 to iquo (k, 3)
   do
      m         := m - 3;
      diff0     := normal (arr [ m + 1 ] - arr [ m ]);
      diff1     := normal (arr [ m + 2 ] - arr [ m + 1 ]);
      diff2     := normal (arr [ m + 3 ] - arr [ m + 2 ]);
      dnom      := normal (diff2 * (diff1 - diff0) -
                   diff0 * (diff2 - diff1));
      arr [ m ] := normal (arr [ m + 1 ] - diff0 * diff1 *
                   (diff2 - diff1) / dnom)
   od;
#
   RETURN (arr [modp (k, 3) ])
fi
end:
`trans/aitken` := proc (s_k, k, arr)
#***********************************************************************
#
#  Procedure  : aitken
#
#  Programmer : Johannes Grotendorst
#
#  Date       : Juelich, 18 june 1990
#  Update     : 1 july 1990
#               normalize intermetiate expressions
#
#-----------------------------------------------------------------------
#
#  Description:
#
#  A recursive algorithm is used to perform the iterated Aitken
#  transformation.
#
#  Input parameters :
#
#  s_k  : sequence element s_k of the sequence to be accelerated
#  k    : number of the last element
#  arr  : 1-dimensional array to store the actual counterdiagonal of
#         the Aitken table
#
#***********************************************************************
#
local j, dnom, m;
#
# initialization
#
arr [ k ] := s_k;
#
# perform the iterated Aitken transformation
#
if k < 2 then
   RETURN (s_k)
else
   m := k;
   for j from 1 to iquo (k, 2)
   do
      m         := m - 2;
      dnom      := normal (arr [ m + 2 ] - 2 * arr [ m + 1 ] + arr [ m ]);
      arr [ m ] := normal (arr [ m ] - (arr [ m ] - arr [ m + 1 ]) ** 2 / dnom)
   od;
#
   if modp (k, 2) = 0 then
      RETURN (arr [ 0 ])
   else
      RETURN (arr [ 1 ])
   fi
fi
end:
`trans/rich` := proc (s_k, x_k, k, arr, xarr)
#***********************************************************************
#
#  Procedure  : rich
#
#  Programmer : Johannes Grotendorst
#
#  Date       : Juelich, 18 june 1990
#  Update     : 1 july 1990
#               normalize intermetiate expressions
#
#-----------------------------------------------------------------------
#
#  Description:
#
#  Neville's recursive scheme is used to perform the Richardson
#  extrapolation process.
#
#  Input parameters :
#
#  s_k  : sequence element s_k of the sequence to be accelerated
#  x_k  : element x_k of the auxiliary sequence
#  k    : number of the last element
#  arr  : 1-dimensional array to store the actual counterdiagonal of
#         Neville's table
#  xarr : 1-dimensional array to store the auxiliary sequence
#
#***********************************************************************
#
local j, diff1;
#
# initializations
#
arr  [ k ] := s_k;
xarr [ k ] := x_k;
#
# recursive computation of Neville's table
#
for j from 1 to k
do
   diff1         := normal (xarr [ k - j ] - xarr [ k ]);
   arr [ k - j ] := normal ((xarr [ k - j ] * arr [ k - j + 1 ] -
                    xarr [ k ] * arr [ k - j ]) / diff1)
od;
#
RETURN (arr [ 0 ]);
#
end:
`trans/gb` := proc (s_k, x_k, k, arr, xarr)
#***********************************************************************
#
#  Procedure  : gb
#
#  Programmer : Johannes Grotendorst
#
#  Date       : Juelich, 18 june 1990
#  Update     : 1 july 1990
#               normalize intermetiate expressions
#
#-----------------------------------------------------------------------
#
#  Description:
#
#  Neville's recursive scheme is used to perform the transformation
#  of Germain-Bonne.
#
#  Input parameters :
#
#  s_k  : partial sum s_k of the series to be accelerated
#  x_k  : term x_k of the series
#  k    : number of the last partial sum
#  arr  : 1-dimensional array to store the actual counterdiagonal of
#         Neville's table
#  xarr : 1-dimensional array to store the terms of the series
#
#***********************************************************************
#
local j, diff1;
#
# initializations
#
arr  [ k ] := s_k;
xarr [ k ] := x_k;
#
# recursive computation of Neville's table
#
for j from 2 to k
do
   diff1         := normal (xarr [ k - j + 1] - xarr [ k ]);
   arr [ k - j ] := normal( (xarr [ k - j + 1 ] *
                    arr [ k - j + 1 ] - xarr [ k ] *
                    arr [ k - j ]) / diff1)
od;
RETURN (arr [ 0 ]);
end:
`trans/lev` := proc (s_k, a_k, k, arlo, arup, typ)
#***********************************************************************
#
#  Procedure  : lev
#
#  Programmer : Johannes Grotendorst
#
#  Date       : Juelich, 18 june 1990
#  Update     : 1 july 1990
#               normalize intermetiate expressions
#
#-----------------------------------------------------------------------
#
#  Description:
#
#  lev converts a series expansion to a rational function
#  using a series transformation of D. Levin (J. Comput. Math.,
#  Vol.B3, 1973, pp.371-388 ).
#
#  In Levin's transformations it is assumed that the partial sums s_n
#  of a series can be written as s_n = x + r_n (x limit or antilimit)
#  and that the first term of the Poincare-type asymptotic expansion
#  of the remainder r_n is proportional to:
#
#  (k+1) * a_k                            (u transformation)
#  a_k                                    (t transformation)
#  a_{k+1}                                (d transformation)
#  a_{k+1} * a_k / (a_k - a_{k+1})        (v transformation)
#
#  (a_k, k >= 0,  are the terms of the series)
#
#  Input parameters :
#
#  s_k   : partial sum s_k of the series to be accelerated
#  a_k   : term a_k of the series
#  k     : number of the last partial sum
#  arlo  : 1-dimensional array to store the counterdiagonal of the
#          numerator table
#  arup  : 1-dimensional array to store the counterdiagonal of the
#          denominator table
#  typ   : type of transformation: u, t, d or v
#
#***********************************************************************
#
local j, fact, k1;
#
#  determine the special part of the specified transformation
#
k1 := k;
#
if typ = 'u' then
   arlo [ k1 ] := 1 / ((k1 + 1) ** 2 * a_k)
elif typ = 't' then
   arlo [ k1 ] := 1 / ((k1 + 1) * a_k)
elif typ = 'd' then
   if k = 0 then
      arlo [ k1 ] := 1
   else
      k1          := k - 1;
      arlo [ k1 ] := 1 / ((k1 + 1) * a_k)
   fi
elif typ = 'v' then
   if k = 0 then
      a_k1        := a_k;
      arlo [ k1 ] := 1
   else
      k1          := k - 1;
      arlo [ k1 ] := (a_k1 - a_k) / ((k1 + 1) * a_k1 * a_k);
      a_k1        := a_k
   fi
fi;
#
if typ = 'd'  or typ = 'v' then
   if k = 0 then
      s_k1        := s_k;
      arup [ k1 ] := s_k
   else
      arup [ k1 ] := s_k1 * arlo [ k1 ];
      s_k1        := s_k
   fi;
else
   arup [ k1 ] := s_k * arlo [ k1 ]
fi;
#
for j from 1 to k1
do
   fact            := (k1 - j + 1) * k1 ** (j - 1) / (k1 + 1) ** j;
   arup [ k1 - j ] := normal (arup [ k1 - j + 1 ] - fact * arup [ k1 - j ]);
   arlo [ k1 - j ] := normal (arlo [ k1 - j + 1 ] - fact * arlo [ k1 - j ])
od;
RETURN (normal (arup [ 0 ] / arlo [ 0 ]));
end:
`trans/sidi` := proc (s_k, a_k, k, arlo, arup, typ)
#***********************************************************************
#
#  Procedure  : sidi
#
#  Programmer : Johannes Grotendorst
#
#  Date       : Juelich, 31 july 1990
#
#-----------------------------------------------------------------------
#
#  Description:
#
#  sidi converts a series expansion to a rational function
#  using a transformation due to Sidi
#  (A. Sidi J. Comput.Appl.Math. Vol. 7 (1981), 37-40)
#
#  As in Levin's transformations it is assumed that the partial sums s_n
#  of a series can be written as s_n = x + r_n (x limit or antilimit)
#  and that the first term of the asymptotic expansion
#  of the remainder r_n is proportional to:
#
#  (k+1) * a_k                            (su transformation)
#  a_k                                    (st transformation)
#  a_{k+1}                                (sd transformation)
#  a_{k+1} * a_k / (a_k - a_{k+1})        (sv transformation)
#
#  (a_k, k >= 0,  are the terms of the series)
#
#  Input parameters :
#
#  s_k   : partial sum s_k of the series to be accelerated
#  a_k   : term a_k of the series
#  k     : number of the last partial sum
#  arlo  : 1-dimensional array to store the counterdiagonal of the
#          numerator table
#  arup  : 1-dimensional array to store the counterdiagonal of the
#          denominator table
#  typ   : type of transformation: su, st, sd or sv
#
#***********************************************************************
#
local j, fact, k1;
#
#  determine the special part of the specified transformation
#
k1 := k;
#
if typ = 'su' then
   arlo [ k1 ] := 1 / ((k1 + 1) * a_k)
elif typ = 'st' then
   arlo [ k1 ] := 1 / a_k
elif typ = 'sd' then
   if k = 0 then
      arlo [ k1 ] := 1
   else
      k1          := k - 1;
      arlo [ k1 ] := 1 / a_k
   fi
elif typ = 'sv' then
   if k = 0 then
      arlo [ k1 ] := 1;
      a_k1        := a_k
   else
      arlo [ k1 ] := (a_k1 - a_k) / (a_k1 * a_k);
      a_k1        := a_k
   fi
fi;
#
if typ = 'sd'  or typ = 'sv' then
   if k = 0 then
      arup [ k ] := s_k;
      s_k1       := s_k
   else
      arup [ k1 ] := s_k1 * arlo [ k1 ];
      s_k1        := s_k
   fi;
else
   arup [ k1 ] := s_k * arlo [ k1 ]
fi;
#
if k1 > 0 then
   arup [ k1 - 1 ] := arup [ k1 ] - arup [ k1 - 1 ];
   arlo [ k1 - 1 ] := arlo [ k1 ] - arlo [ k1 - 1 ];
   for j from 2 to k1
   do
      fact            := (k1 - 1) * k1 / ((k1 + j - 2) * (k1 + j - 1));
      arup [ k1 - j ] := normal (arup [ k1 - j + 1 ] - fact * arup [ k1 - j ]);
      arlo [ k1 - j ] := normal (arlo [ k1 - j + 1 ] - fact * arlo [ k1 - j ])
   od;
fi;
RETURN (normal (arup [ 0 ] / arlo [ 0 ]));
end:
#
# difference between my_fort and the fortran function in MapleV:
# In my_fort an output_file is opended only if the output_filename
# in unequal to terminal.
#

macro(suffix=`fortran/suffix`); # Fix for Maple V release II
`trans/my_fort` := proc(x)
local k,error,old_Digits,old_Quiet,output_file,statseq,opt;
options `Copyright 1990 by the University of Waterloo`;
    fortran; # load fortran code to define `fortran/breakup` etc.
    old_Digits := Digits;
    old_Quiet := interface(quiet);
    if assigned(precision) and not member(precision,['single','double']) then
        ERROR(`bad setting for global precision`,precision)
    fi;
    if precision = 'double' then suffix := 'D'; Digits := 16
	else Digits := 7; suffix := 'E' fi;
    opt := false;
    output_file := 'terminal';
    for k in [args[2 .. nargs]] do
        if k = optimized then opt := true
        elif type(k,'identical(digits) = posint') then Digits := rhs(k)
        elif type(k,'identical(filename) = string') then
            output_file := op(2,k)
        else ERROR(`bad option`,k)
        fi
    od;
    if opt then statseq := [readlib(optimize)(x)]
    elif type(x,list(name = algebraic)) then statseq := x
    elif type(x,name) and type(x,'array(algebraic)') then
        statseq := [`fortran/arrayexpand`(x)]
    elif type(x,algebraic) then statseq := [t0 = x]
    else ERROR(`invalid argument`,x)
    fi;
    statseq := `fortran/breakup`(statseq);
    if 1 < printlevel then lprint(`fortran: generating output`) fi;
    interface(quiet = true);
#   writeto(output_file);
    if output_file <> `terminal` then writeto(output_file) fi;
    error := traperror(`fortran/statement`(op(statseq)));
#   writeto(terminal);
    if output_file <> `terminal` then writeto(terminal) fi;
    interface(quiet = old_Quiet);
    Digits := old_Digits;
    if error = lasterror then ERROR(error) else NULL fi
end:
#************************************************************************
#
# Help file for the package trans
#
#************************************************************************
`help/text/trans` := TEXT (
`HELP FOR: The rational approximation package`,
`   `,
`SYNOPSIS:`,
`To use a trans function either define that function alone using`,
`with(trans,<function>), or define all trans functions using with(trans)`,
`   `,
`The functions available are:`,
`   `,
`            aitken eps   gb   lev  ratser ratgen`,
`            rho    rhoit rich sidi theta  thetait`,
`   `,
`They can transform symbolic as well as numerical sequences (series).`,
`The elements s_n with n = 0, 1, 2, ... of the sequence to be`,
`transformed have be read in successively, starting with s_0.`,
`The obligatory input parameters are:`,
`   `,
`eps     (s_n, n, arr1)`,
`aitken  (s_n, n, arr1)`,
`thetait (s_n, n, arr1)`,
`theta   (s_n, n, arr1, arr2)`,
`rho     (s_n, x_n, n, arr1, arr2)`,
`rhoit   (s_n, x_n, n, arr1, arr2)`,
`rich    (s_n, x_n, n, arr1, arr2)`,
`gb      (s_n, a_n, n, arr1, arr2)`,
`lev     (s_n, a_n, n, arr1, arr2, typ)`,
`sidi    (s_n, a_n, n, arr1, arr2, typ)`,
`   `,
`The meaning of the formal parameters is as follows:`,
`   `,
`s_n    element of the sequence `,
`x_n    element of the auxiliary sequence`,
`a_n    term a_n = s_n - s_{n-1} of a series`,
`n      number of the last sequence element`,
`arr1   1-dimensional array to store intermediate`,
`       elements of the table of transforms`,
`arr2   1-dimensional array to store intermediate`,
`       elements of the table of transforms`,
`typ    type of transformation:`,
`       u, t, d, v in case of function lev`,
`       su, st, sd, sv in case of function sidi`,
`   `,
`The program ratser is designed to transform symbolically`,
`series expansions with  data type "series". The program ratgen`,
`transforms generalized series expansions, e.g. asymptotic`,
`expansions, with data type ``+``.`,
`   `,
`ratser has the following formal input parameters:`,
`   `,
`obligatory:`,
`s         series expansion with data type "series"`,
`          (has to be the first parameter)`,
`optional:`,
`typ       type of transformation`,
`norder    order of transformation (default: number of terms)`,
`expts     list with expansion points of data type real constant`,
`          for segmented approximation (default: [ ], empty list)`,
`switch    switch for the generation of a FORTRAN function:`,
`          fort or nofort (default: nofort)`,
`precision type of the FORTRAN function: single or double`,
`          (default: single)`,
`evalorder computational form for the rational function:`,
`          horner or confrac (default: horner)`,
`fn        filename for the FORTRAN function program. (default: ````)`,
`   `,
`Formal parameters of procedure ratgen:`,
`   `,
`obligatory:`,
`s         generalized series with data type ``+```,
`          (has to be the first parameter)`,
`optional:`,
`typ       type of transformation`,
`norder    order of transformation (default: number of terms)`,
`switch    switch for the generation of a FORTRAN function:`,
`          fort or nofort (default: nofort)`,
`precision type of the FORTRAN function: single or double`,
`          (default: single)`,
`evalorder computational form for the rational function:`,
`          horner (default: ````)`,
`fn        filename for the FORTRAN function program (default:````)`,
`   `,
`In both procedures the formal input parameter typ (type of transformation)`,
`can have one of the following values:`,
`   `,
`u      u transformation`,
`t      t transformation`,
`d      d transformation`,
`v      v transformation`,
`su     su transformation`,
`st     st transformation`,
`sd     sd transformation`,
`sv     sv transformation`,
`sh     Shanks transformation`,
`gb     Germain-Bonne transformation`,
`rh     rho transformation`,
`rhit   iterated rho_2 transformation`,
`r      Richardson extrapolation process`,
`th     theta algorithm`,
`thit   iterated theta_2 algorithm`,
`at     iterated Aitken transformation`,
`   `,
`In case of the Richardson extrapolation (r), rho transformation (rh)`,
`or iterated rho_2 transformation (rhit) the user can specify a list`,
`of auxiliary parameters x_n, (default: x_n = 1/(n+1) for the`,
`Richardson extrapolation process and x_n = n+1 for the rho`,
`transformation and the iterated rho_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.`,
`   `,
`EXAMPLES:`,
`> with(trans);`,
`> a1 := series(GAMMA(x), x = x0, 3);`,
`   `,
` a1 := GAMMA(x0) + Psi(x0) GAMMA(x0) (x - x0)`,
`   `,
`                                               2                    2`,
`      + (1/2 Psi(1, x0) GAMMA(x0) + 1/2 Psi(x0)  GAMMA(x0)) (x - x0)`,
`   `,
`                  3`,
`      + O((x - x0) )`,
`   `,
`> a2 := ratser (a1, sh, horner);`,
`   `,
`                                                                2`,
`a2 := (- 2 Psi(x0) GAMMA(x0) + (- Psi(1, x0) GAMMA(x0) + Psi(x0)  GAMMA(x0)) x0`,
`   `,
`                                      2`,
`     + (Psi(1, x0) GAMMA(x0) - Psi(x0)  GAMMA(x0)) x)`,
`   `,
`       /                                       2                            2`,
`      /  (- 2 Psi(x0) + (- Psi(1, x0) - Psi(x0) ) x0 + (Psi(1, x0) + Psi(x0) )`,
`     /`,
`   `,
`> ratser (a1, sh, fort, horner, double, [1,2,3]);`,
`   `,
`      function shtr (x)`,
`      double precision shtr, x`,
`      if (x .le. 150.D-2) then`,
`      shtr = (-0.157324813237442D0+0.1311756143040508D1*x)/(-0.823680660`,
`     #8528794D0+0.1978111990655945D1*x)`,
`      elseif (x .gt. 150.D-2 .and. x .le. 250.D-2) then`,
`      shtr = (-0.1777943615884081D1+0.4661874728435735D0*x)/(-0.24929299`,
`     #91902693D1+0.8236806608528794D0*x)`,
`      elseif (x .gt. 250.D-2) then`,
`      shtr = (-0.9515561668645066D0-0.9131937245097873D0*x)/(-0.55849636`,
`     #58050974D1+0.1246464995951347D1*x)`,
`      end if`,
`      end`,
`   `,
`> a3 := series (sqrt(sin(x)), x, 8);`,
`   `,
`               1/2         5/2           9/2            13/2      15/2`,
`        a3 := x    - 1/12 x    + 1/1440 x    - 1/24192 x     + O(x    )`,
`   `,
`> ratgen (a3, [rh, 1, 6, 9, 16]);`,
`   `,
`                       1/2          5/2         9/2       13/2`,
`               846720 x    - 34560 x    - 2412 x    - 35 x     `,
`               -----------------------------------------------`,
`                                              2 `,
`                              846720 + 36000 x `,
`   `,
`> a4 := asympt (Psi (2*exp(x))-x, x, 4);`,
`   `,
`                                  1           1            1`,
`                a4 := ln(2) - -------- - ---------- + O(-------)`,
`                              4 exp(x)            2           4`,
`                                         48 exp(x)      exp(x)`,
`   `,
`> ratgen (a4, th);`,
`   `,
`                                        2`,
`                         48 ln(2) exp(x)  - 12 exp(x) - 1`,
`                    1/48 --------------------------------`,
`                                            2`,
`                                      exp(x)`
):

macro(suffix=suffix);
save `trans.m`;
quit
