\documentstyle[fullpage,maple,11pt]{article}
\title{Maple Package on Formal Power Series}%
\author{
             Dominik Gruntz\\
       Institute for Scientific Computing\\
       ETH Z\"urich\\
       CH-8092 Z\"urich \\
       {\tt gruntz@inf.ethz.ch}
\and
       Wolfram Koepf\\
       Konrad-Zuse-Zentrum f\"ur Informationstechnik\\
       Heilbronner Str.\ 10\\
       D-10711 Berlin\\
       {\tt koepf@zib-berlin.de}
}
\date{March 1, 1995}
\begin{document}
\maketitle
\newcommand{\heading}[1]{\section{#1}}

\newcommand{\C}{{\rm {\mbox{C{\llap{{\vrule height1.52ex}\kern.4em}}}}}}
\newcommand{\N} {{\rm {\mbox{\protect\makebox[.15em][l]{I}N}}}}
\newcommand{\Q} {{\rm {\mbox{Q{\llap{{\vrule height1.5ex}\kern.5em}}}}}}
\newcommand{\R} {{\rm {\mbox{\protect\makebox[.15em][l]{I}R}}}}
\newcommand{\D} {{\rm {\mbox{\protect\makebox[.15em][l]{I}D}}}}
\newcommand{\Z} {{\rm {\mbox{\protect\makebox[.2em][l]{\sf Z}\sf Z}}}}

\newcommand{\Maple}{{\sc Maple}}
\newtheorem{theorem}{\sc Theorem}
\newtheorem{lemma}[theorem]{\sc Lemma}
\newtheorem{definition}{\sc Definition}
\newcommand{\notexists}{\raisebox{1pt}{$/$}\!\!\!\exists\,\,}

\newcommand{\funkdeff}[4]{\left\{\begin{array}{ccc}
                                 #1 && \mbox{\rm{if} $#2$ } \\[4mm]
                                 #3 && \mbox{\rm{if} $#4$ } 
                                 \end{array}
                          \right.}
\def\erf{\mathop{\rm erf}\nolimits}
\def\arccosh{\mathop{\rm arccosh}\nolimits}
%\input{macros.tex}
\heading{Introduction}\noindent
Formal Laurent-Puiseux series of the form
\begin{equation}
f(x)=\sum\limits_{k=k_0}^{\infty}a_{k}x^{k/n}
\label{eq:formalLPS}
\end{equation}
with coefficients $a_{k}\in\C\;(k\in\Z)$ are important in many branches of
mathematics. \Maple\ supports the computation of {\em truncated\/} series
with its {\tt series} command, and through the {\tt powseries} package
infinite series are available. In the latter package the series is
represented as a table of coefficients that have already been determined
together with a function for computing additional coefficients. This
representation is known as {\em lazy evaluation\/}~\cite{Gruntz94}. But these
tools fail, if one is interested in an explicit formula for the coefficients
$a_k$.

In this article, which is a condensed version of~\cite{GruntzKoepf93}, we
describe the \Maple\ implementation of an algorithm presented
in~\cite{Koe92}--\cite{Koeortho} which computes an {\em exact\/} formal power
series~(FPS) of a given function. This procedure enables the user to
reproduce most of the results of the extensive bibliography on
series~\cite{Han}. We will give an overview of the algorithm and then present
some parts of it in more detail.

This package is available through the \Maple-share library with the name {\tt
FPS}. We will flavor this procedure with the following example.
\begin{maple}
> with(share):
> readshare(FPS, calculus):
> FormalPowerSeries(exp(x) * sin(x), x=0);

                   infinity
                    -----     k 1/2                k
                     \      (2 )    sin(1/4 k Pi) x
                      )     ------------------------
                     /                 k!
                    -----
                    k = 0
\end{maple}

The main method used works as follows: in a first step, a holonomic, i.e.\ a
homogeneous linear differential equation with polynomial coefficients~(DE) is
generated for a given function $f(x)$. Then this DE is converted to an
equivalent holonomic recurrence equation~(RE) for the corresponding
coefficients $a_k$. In a final step, this RE is solved if possible. Note,
that the class of functions for which a holonomic DE is found includes all
sums, products and compositions with algebraic functions of elementary
functions for which holonomic DEs are
known~\cite{Stanley80,Salvy92,Koediffgl}. In particular, the following
functions are supported: $\exp(x)$, $\sin(x)$, $\cos(x)$, $\arcsin(x)$,
$\arctan(x)$, $\erf(x)$, Bessel functions, the Hankel functions
$\mbox{Hankel1}(n,x)$ and $\mbox{Hankel2}(n,x)$~\cite[(9.1)]{AS}, the
Fibonacci polynomials $\mbox{Fibonacci}(n,x)$, the Whittaker functions
$\mbox{WhittakerM}(n,m,x)$ and $\mbox{WhittakerW}(n,m,x)$~\cite[(13.4)]{AS},
and families of orthogonal polynomials $\mbox{JacobiP}(n,a,b,x)$,
$\mbox{GegenbauerC}(n,a,x)$, $\mbox{ChebyshevT}(n,x)$,
$\mbox{ChebyshevU}(n,x)$, $\mbox{LegendreP(n,x)}, $\mbox{LaguerreL(n,a,x)},
and $\mbox{HermiteH}(n,x)$~\cite[(22.7) and (22.8)]{AS}. The {\tt FPS}
package just knows these functions under these names.

\heading{Search for the Holonomic DE}\noindent
In this section we present the algorithm that searches for a holonomic DE of
degree $k$ for a given function $f$. We set up the equation
\[
   f^{(k)}(x) + \sum_{j=0}^{k-1} A_j\,f^{(j)}(x) = 0
\]
and expand it. Then we collect the coefficients of all the rationally
dependent terms and equate them to zero. For testing whether two terms are
rationally dependent, we divide one by the other and test whether the
quotient is a rational function in $x$ or not. This is an easy and fast
approach. Of course, we could also use the Risch normalization
procedure~\cite{Risch,Bronstein} to generate the rationally independent
terms, but firstly this normalization is rather expensive and works only for
elementary functions and secondly, our simplified approach will not lead to
wrong results, it may at most happen that we miss a simpler solution, which,
in practice, rarely happens, however. This procedure by its own is available
under the name {\tt SimpleDE}.

It may happen, that for a given function $f(x)$ a DE of degree $k$ exists,
which has neither constant coefficients (which we call the {\em explike
case}), nor is the corresponding RE of hypergeometric type~(\ref{eq:RE1}) and
hence no closed form for the  formal Laurent-Puiseux series~(LPS) can be
computed. In this situation, the {\tt FPS} implementation looks for a DE of
higher degree which will have free parameters as from the existence of a DE
of degree $k$ follows the existence of families of DEs of higher degree.
These parameters can be set freely and we can try to choose them in such a
way, that we can use our tools to compute a formal LPS, i.e.\ we either need
a RE of hypergeometric type or a DE with constant coefficients.

In order to obtain a RE of hypergeometric type we convert the DE into the
corresponding RE and try to set all the coefficients but two to zero. For
example, let $f(x)=x^{-1}\,e^x\,\sin{x}$. We find DEs of degree 2 and of
degree 3, but neither of the corresponding REs can be solved. The RE which
corresponds to the DE of degree 4 has the form
 \[
   \sum_{j=0}^4 p_j(k)\,a_{k+j} = 0
\]
where
\begin{eqnarray*}
p_0(k) & = & 2 A_2+4 A_3+4\\
p_1(k) & = & - 4 A_2 - 10 A_3 - 16 - 2 A_3 k - 2 A_2 k\\
p_2(k) & = & 18 A_3 + 5 A_2 k + 6 A_2 + 48 + 8 k + 6 A_3 k + A_2 k^2\\
p_3(k) & = & (4 + k) (A_3 k^2  + 2 A_3 k - 24 - 3 A_3)\\
p_4(k) & = & (k + 5) (4 + k) (k^2 + k + 6)\\
\end{eqnarray*}
The solution of the equations $p_1(k)=0, p_2(k)=0, p_3(k)=0$ (forcing that
only two terms of the sum remain) with respect to $A_2$ and $A_3$ is
\[
A_{2} = -8 \frac{k^2+5k+12}{k^3+4 k^2+k-6},\quad
A_{3} =    \frac{24}{k^2+2k-3}.
% A[2] := -8*(k^2+5*k+12)/(k^3+4*k^2+k-6):
% A[3] := 24/(k^2+2*k-3):
% req := req * (k+2)*(k+3):
\]
This results (after multiplying by $\frac{(k+2)\,(k+3)}{k^2+k+6}$) in the
following RE of hypergeometric type for $f(x)$,
\[
\left(k+5\right)\left(k+4\right)\left(k+3\right)
\left(k+2\right)\,a_{k+4}
+4\,a_k = 0
\]
with symmetry number $m=4$~(\ref{eq:RE1}). Of course, we try to keep the
symmetry number, that is the number of resulting sums, as small as possible.
The final result for this example is
\begin{maple}
> FormalPowerSeries(exp(x)*sin(x)/x, x);

/infinity                          \   /infinity                              \
| -----       k   (- k)    k  (4 k)|   | -----       k   (- k)    k  (1 + 4 k)|
|  \      (-1)  64      256  x     |   |  \      (-1)  64      256  x         |
|   )     -------------------------| + |   )     -----------------------------|
|  /              (1 + 4 k)!       |   |  /           (1 + 4 k)! (2 k + 1)    |
| -----                            |   | -----                                |
\ k = 0                            /   \ k = 0                                /

       /infinity                               \
       | -----        k   (- k)    k  (4 k + 2)|
       |  \       (-1)  64      256  x         |
     + |   )     ------------------------------|
       |  /      (4 k + 3) (1 + 4 k)! (2 k + 1)|
       | -----                                 |
       \ k = 0                                 /
\end{maple}
Note that one of the four sums (with the powers $x^{4k+3}$) vanishes as a
result of a vanishing initial coefficient.

If this approach fails, we try to choose the parameters such that the DE gets
constant coefficients. Let us look at a rather similar example,
$f(x)=x\,e^{x}\,\sin(2x)$. We again find DEs of second and third order whose
corresponding REs cannot be solved. The DE of degree 4 has two free
parameters $A_2$ and $A_3$:
\begin{eqnarray*}
\lefteqn{x^2\,\frac{d^4}{dx^4}f(x) +
   A_{3}{x}^{2}\,\frac{d^3}{dx^3}f(x) +
   A_{2}{x}^{2}\,\frac{d^2}{dx^2}f(x) +}\\
\quad&&
   \left(
      (A_3-2\,A_2+12)\,x^2+(-6\,A_3-2\,A_2+4)\,x
   \right)\,\frac{d}{dx}f(x) +\\
\quad&&
   \left(
      (10\,A_3+5\,A_2-5)x^2+(14\,A_3+2\,A_2+28)x +(6\,A_3+2\,A_2-4)
   \right)f(x)=0
\:.
\end{eqnarray*}
The DE will have constant coefficients, if both $A_2$ and $A_3$ are constants
and if they meet the following equations 
\[
\begin{array}{rcc}
    6 A_3 + 2 A_2 -  4 & = & 0 \\
   14 A_3 + 2 A_2 + 28 & = & 0
\end{array}
\;.
\]
The solution of this system of equations is $A_2 = 14$ and $A_3 = -4$. If we
insert these values in the differential equation and divide by $x^2$, then we
get the DE
\[
   \frac{d^4}{dx^4}f(x) - 4\frac{d^3}{dx^3}f(x) + 14\frac{d^2}{dx^2}f(x)
   -20\frac{d}{dx}f(x) + 25f(x) = 0
\]
which has constant coefficients leading to a constant coefficient RE for
$b_k$ given by $f(x)=\sum\limits\frac{b_k}{k!}\,x^k$ that can be solved:
\begin{maple}
> FormalPowerSeries(x*exp(x)*sin(2*x),x);

infinity
 -----   /        k 1/2                            k 1/2                   \
  \      |      (5 )    cos(k arctan(2)) k       (5 )    sin(k arctan(2)) k|  k
   )     |- 2/5 -------------------------- + 1/5 --------------------------| x
  /      \                  k!                               k!            /
 -----
 k = 0
\end{maple}
In the {\tt gfun} package of the share library~\cite{Salvy92}, the procedure
${\tt `diffeq+diffeq`}(eq_1, eq_2, F(x))$ and ${\tt `diffeq*diffeq`}(eq_1,
eq_2, F(x))$ calculate holonomic DEs for the sum and product of solutions
$f_1$ and $f_2$ of $eq1$ and $eq2$, respectively. These procedures could also
be used to derive holonomic DEs for expressions given as sums and products of
elementary functions like the ones mentioned in the introduction. Note,
however, that {\tt SimpleDE} usually generates the holonomic DE of lowest
order, e.g.\ 
\begin{maple}
> SimpleDE(sqrt(1+x)+1/sqrt(1+x),x);

                                     /  d      \
                   2 (2 + x) (1 + x) |---- F(x)| - x F(x) = 0
                                     \ dx      /
\end{maple}
whereas {\tt `diffeq+diffeq`} calculates a holonomic DE valid for all linear
combinations of $f_1$ and $f_2$, e.g.\ 
% read `/home/rutishauser/ru4/monagan/share/5.3/share/gfun`;
\begin{maple}
> eq1 := SimpleDE(sqrt(1+x),x,F);
                                      /  d      \
                     eq1 := (2 + 2 x) |---- F(x)| - F(x) = 0
                                      \ dx      /

> eq2 := SimpleDE(1/sqrt(1+x),x,F);
                                       /  d      \
                      eq2 := (2 + 2 x) |---- F(x)| + F(x) = 0
                                       \ dx      /
> `diffeq+diffeq`(eq1, eq2, F(x));

                                                        2   (2)
           F(x) + (- 4 - 4 x) D(F)(x) + (- 4 - 8 x - 4 x ) D   (F)(x)
\end{maple}
An upper bound for the order of the differential equation can be computed a
priori~\cite{Stanley80}, however this bound is not sharp for all cases.
Therefore, and since the complexity of the algorithm increases rapidly with
the order, we stop the search if no DE of order five has been found. This
predefined upper bound can be overwritten by specifying an additional integer
argument when calling the {\tt FormalPowerSeries} function.

\heading{Conversion to the Holonomic Recurrence Equation}\noindent
It is well-known, that the conversion of a holonomic DE into the
corresponding holonomic RE is done by the substitution
\begin{equation}
x^l f^{(j)}(x) \mapsto (k\!+\!1\!-\!l)_j \cdot a_{k+j-l}
\label{eq:rule}
\end{equation}
into the DE.


The search of a DE and its conversion to the RE is directly available through
the command {\tt SimpleRE}. We see that {\tt AiryAi$(x^2)$}\footnote{{\tt
AiryAi} is the Airy wave function which is available in \Maple\ V Release 4.
In previous versions this function was called {\tt Ai}, but the derivative
was defined differently which led to a different DE and RE.} is of
hypergeometric type with symmetry number $m=6$, whereas the second RE is not
of hypergeometric type.
\begin{maple}
> SimpleRE(AiryAi(x^2), x);

                   (k - 1) (k + 1) a(k + 1) - 4 a(k - 5) = 0

> SimpleRE(x/(1-x-x^2), x);

             (1 - k) a(k) + (k - 1) a(k - 1) + (k - 1) a(k - 2) = 0
\end{maple}
The latter example is the generating function of the Fibonacci numbers, and
we get the expected RE. Note, that the common factor $(k-1)$ ensures, that
the RE holds $\forall\ k\in\Z$.

\heading{Solving a Recurrence Equation of Hypergeometric Type}\noindent
If the recurrence equation of a function $f(x)$ is of the form
\begin{equation}\label{eq:RE1} 
   Q(k)\, a_{k+m} = P(k)\, a_k
\end{equation}
($P,Q$ polynomials) then we call $f$ to be of hypergeometric type, and $m$
the symmetry number of the corresponding series representation. The explicit
formula for the coefficients can be found by the hypergeometric coefficient
formula and some initial conditions. Based on the analysis of the polynomials
$P(k)$ and $Q(k)$, we will convert the RE into one corresponding to a {\em
Taylor series\/} by applying a sequence of transformations stated in the
following lemma which preserve the hypergeometric type.

\begin{lemma}\label{Lemma:trafos}
Let $f$ be a formal Laurent series~(FLS) of hypergeometric type with
representation $\sum\limits_{k=k_{0}}^{\infty}a_{k}x^{k/n}$ and $a_{k_0}\neq
0$, whose coefficients $a_k$ satisfy a RE of the form
\begin{eqnarray*}
a_{k+m}&=&R(k)\,a_k\quad\:\mbox{for $k\geq k_0$}\\
a_{k}&=&A_{k}\quad\quad\quad\;\mbox{for $k=k_0,k_0+1,\ldots,k_0+m-1$},
\end{eqnarray*}
then the following functions are of hypergeometric type, too. Their
coefficients $b_k$ satisfy a RE whose relation to the RE of $f$ is also
given.
\[
\begin{array}{llllll}
\mbox{(a)} &\quad& x^n f(x)   &\quad& b_{k+m} = R(k-n)\,b_k & n \in \Z\\
\mbox{(b)} && f(A x)     && b_{k+m} = A^m\,R(k)\,b_k & A \in \C\\
\mbox{(c)} && f(x^n)     && b_{k+n\,m} = R(k/n)\,b_k & n \in \N\\
\mbox{(d)} && f(x^{1/m}) && b_{k+1} = R(k\,m)\,b_k &\\
\mbox{(e)} && f'(x)      && b_{k+m} = \frac{k+m+1}{k+1}\,R(k+1)\,b_k & \\
\end{array}
\]
\end{lemma}


First of all, we inspect the roots of $P(k)$ and $Q(k)$ of the
RE~(\ref{eq:RE1}). If there are any non-integer rational roots, then we know
that $f$ corresponds (possibly) to a Puiseux series. In this case we
transform $f$ to a function of Laurent type by an application of
transformation $(c)$, where $n$ is the least common multiple of the
denominators of all rational roots of $P(k)$ and $Q(k)$. The FLS we get when
we solve the transformed RE can be transformed back to the LPS of $f$ by
substituting $x$ by $x^{1/n}$.

Let us now assume that $f$ can be expanded in a FLS. We reduce this problem
to solving the RE of a function which has a FPS expansion. For that purpose
we remove the finite pole of $f$ at the origin by multiplying $f$ with a
suitable power of $x$ (transformation $(a)$). From the information of the RE
we can determine which power we have to use. Let us assume that $f$ may be
expanded in a Laurent series, i.e.\ that $\exists\ k_0\ :\ \forall\ k\leq
k_0\ :\ a_k=0$. From these known coefficients we can derive further ones
using the given RE in the form
\begin{equation}\label{eq:RE}
   a_{k+m} = \frac{P(k)}{Q(k)}\, a_k.
\end{equation}
If we know that $a_k=0$ then also $a_{k+m}=0$ provided that $Q(k) \neq 0$.
Let $k_{min}$ be the smallest integer root of $Q(k)$. Consequently $a_k=0\
\forall\ k < k_{min}+m$. From this it follows, that
\[
   g = x^{-(k_{min}+m)} \, f
\]
may be expanded into a FPS, i.e.\ $g$ has no pole at the origin, given the
assumption that $f$ may be expanded in a Laurent series. This latter
assumption can be tested by computing the limit of $g$ as $x$ tends to $0$.
This limit must be finite. (Since we also allow logarithmic singularities, we
test in fact whether the limit of $x\,g$ is 0.) If this is not the case, then
the assumption that $f$ may be expanded into a FLS is wrong and $f$ must have
an essential singularity, i.e.\ $\notexists k_0\ :\ \forall\ k\leq k_0\ :\
a_k=0$. Otherwise the FPS of $g$ exists and can again easily be transformed
back to the FLS of $f$.

What we finally have to show is how to solve a RE which corresponds to a
given function $f$ which has a FPS expansion. The RE~(\ref{eq:RE}) is valid
$\forall\ k\ :\ Q(k)\neq 0$, especially $\forall\ k>k_{max}$ where $k_{max}$
is the largest root of $Q(k)$. Consequently we must determine $a_k$ for
$k=0,1,\ldots,k_{max}+m$ and have to solve the hypergeometric RE for
$x^{-k}\,f(x^{1/m})$ for $k>k_{max}$. We investigate now how this condition
can be weakened.

The coefficient $a_k$ is given by the limit $\lim\limits_{x\to 0}
f^{(k)}(x)/k!$. Let us assume, that this limit is finite. Then the
hypergeometric RE can be solved in the case that $\forall\ j\geq 0\ :\
Q(k+j\,m)\neq 0$, otherwise simply $a_k\,x^k$ is added to the result.
Moreover note, that
\[
   \left(P(k)=0 \vee a_k=0 \right) \wedge 
   \forall\ j\geq 0\ :\ Q(k+j\,m)\neq0
   \wedge a_k \;\mbox{is finite}
   \Longrightarrow
   \forall\ j>0\ :\ a(k+j\,m) = 0
\]
and in this case the RE does not have to be solved and it is enough to add
$a_k\,x^k$ to the result. The indices $k+j\,m, j>0$ no longer need to be
considered.

If $a_k$ is infinite, then we found a logarithmic singularity which we can
remove by working with \[g=(x^{-k}\,f(x))'\] and by integrating and shifting
the resulting power series $S$. The RE of $g$ can be obtained from the RE of
$f$ by applying transformations $(a)$ and $(e)$. The constant term which we
lose by the differentiation can be determined by computing the limit
\[\lim_{x\to 0} f(x)/x^k-\int S(x)\,dx.\] Note that $g$ in general is a
function with a FLS expansion and we first must remove the pole to get a
function with a FPS expansion. This is the reason why we have chosen a
recursive implementation of the RE solver. Every application of a
transformation of Lemma~\ref{Lemma:trafos} is nothing but a recursive call of
the RE solver. One may inspect which steps the algorithm performs by
assigning the variable {\tt infolevel[FormalPowerSeries]}. As an example we
trace the computation for $f(x)=x^{-1}\,\sin{\sqrt{x}}$. A DE of degree 2 is
found whose corresponding RE is of hypergeometric type. From the root of
$2k+3$ it follows that the Puiseux number is 2 and so transformation $(c)$
with $n=2$ is applied. The resulting function is of Laurent type. The
smallest integer root of $Q(k)$ of the transformed RE is $-4$ and the
symmetry number is $m=2$, hence we multiply the function with $x^2$ and
adjust the RE accordingly. We end up with the function $\sin{x}$ whose FPS
can be computed directly. This result is then transformed back to the LPS of
$f(x)$ according to the two transformations we applied.
\begin{maple}
> infolevel[FormalPowerSeries] := 4:
> FormalPowerSeries(sin(sqrt(x))/x,x);
FPS/FPS:   looking for DE of degree   1
FPS/FPS:   looking for DE of degree   2
FPS/FPS:   DE of degree   2   found.
FPS/FPS:   DE = 
                     2
                  4 x  F''(x) + 10 x F'(x) + (2 + x) F(x) = 0

FPS/hypergeomRE:   RE is of hypergeometric type.
FPS/hypergeomRE:   Symmetry number m = 1
FPS/hypergeomRE:   RE: 2 (k + 2) (2 k + 3) a(k + 1) = - a(k)
FPS/hypergeomRE:   RE modified by k = 1/2*k
FPS/hypergeomRE:   => f := sin(x)/x^2

FPS/hypergeomRE:   RE is of hypergeometric type.
FPS/hypergeomRE:   Symmetry number m = 2
FPS/hypergeomRE:   RE: (k + 4) (k + 3) a(k + 2) = - a(k)
FPS/hypergeomRE:   working with x^2*f
FPS/hypergeomRE:   => f := sin(x)

FPS/hypergeomRE:   RE is of hypergeometric type.
FPS/hypergeomRE:   Symmetry number m = 2
FPS/hypergeomRE:   RE:  (k + 2) (k + 1) a(k + 2) = - a(k)
FPS/hypergeomRE:   RE valid for all k >=    0
FPS/hypergeomRE:   a(0) =   0
FPS/hypergeomRE:   a(2*j) = 0   for all j>0.
FPS/hypergeomRE:   a(1) =   1

                           infinity
                            -----       k  (k - 1/2)
                             \      (-1)  x
                              )     ----------------
                             /         (2 k + 1)!
                            -----
                            k = 0

> infolevel[FormalPowerSeries] := 1:
\end{maple}
\heading{Rational Algorithm}\noindent
If the given function is rational in $x$, first the complex partial fraction
decomposition of $f(x)$ is calculated. Each term of the form
$\frac{c}{(x-\alpha)^j}$ can be expanded by the binomial series whose
coefficients are
\[a_k = \frac{(-1)^j\,c}{\alpha^{j+k}}
          \left(\begin{array}{c}j+k-1\\k\end{array}\right).
\]
\begin{maple}
> FormalPowerSeries(1/((x-1)^2*(x-2)),x);

                infinity
                 -----                  k                  k   k
                  \            (k! - 2 2  k! + 2 (1 + k)! 2 ) x
                   )     - 1/2 ---------------------------------
                  /                           k
                 -----                       2  k!
                 k = 0

> FormalPowerSeries((1+x+x^2+x^3)/((x-1)*(x-2)),x);

                             /infinity                   \
                             | -----           k        k|
                             |  \          (8 2  - 15) x |
                     x + 4 + |   )     1/2 --------------|
                             |  /                 k      |
                             | -----             2       |
                             \ k = 0                     /

> FormalPowerSeries(C/(B*A - A*x - B*x+x^2),x);

                          infinity
                           -----         k      k   k
                            \      C (A A  - B B ) x
                             )     ------------------
                            /                  k    k
                           -----    (A - B) B B  A A
                           k = 0
\end{maple}
To get the complex partial fraction decomposition we must factor the
denominator which may be rather complicated, hence the rational algorithm to
compute the full partial fraction expansion presented in~\cite{Salvy} may be
used. The following example uses this code. This method can be forced to be
used, if the environment variable \verb|_EnvExplicit| is set to false.
\begin{maple}
> FormalPowerSeries(1/(x^4+x+1),x);

      infinity
       -----   /  -----                                 2           3\
        \      |   \            27 + 64 alpha - 48 alpha  + 36 alpha |  k
         )     |    )     1/229 -------------------------------------| x
        /      |   /                              (1 + k)            |
       -----   |  -----                      alpha                   |
       k = 0   \alpha = %1                                           /

                                 4           
%1 :=                   RootOf(_Z  + _Z + 1) 
\end{maple}
A closed form of the Fibonacci numbers can be derived by computing the FPS of
their generating function $x/(1-x-x^2)$.
\begin{maple}
> expand(FormalPowerSeries(x/(1-x-x^2),x));

             /infinity                                            \
             | -----   /          k                     k        \|
         1/2 |  \      |         x                     x         ||
    1/5 5    |   )     |------------------- - -------------------||
             |  /      |              1/2 k                 1/2 k||
             | -----   \(- 1/2 + 1/2 5   )    (- 1/2 - 1/2 5   ) /|
             \ k = 0                                              /
\end{maple}
The {\tt expand} command converts the coefficient in the usual notation. If
the factorization of the numerator is avoided, then the following result is
obtained:
\begin{maple}
> _EnvExplicit := false:
> FormalPowerSeries(x/(1-x-x^2),x);

        infinity
         -----   /            -----                                \
          \      |             \                         alpha - 2 |  k
           )     |              )                - 1/5 ------------| x
          /      |             /                            (1 + k)|
         -----   |            -----                    alpha       |
         k = 0   |                            2                    |
                 \alpha = RootOf(- 1 + _Z + _Z )                   /

> _EnvExplicit := '_EnvExplicit':
\end{maple}
 
\heading{Additional Examples}\noindent
In this section we present some more results generated with the procedure
{\tt FormalPowerSeries}.

First, we consider the FPS of the Airy wave function {\tt AiryAi}.
% `diff/AiryAi`      := (e,x) -> diff(e,x)*AiryAiPrime(e):
% `diff/AiryAiPrime` := (e,x) -> diff(e,x)*e*AiryAi(e):
% AiryAi(0)      :=  1/3^(2/3)/GAMMA(2/3):
% AiryAiPrime(0) := -1/3^(1/3)/GAMMA(1/3):
\begin{maple}
> FormalPowerSeries(AiryAi(x),x);

                      /infinity                      \
                      | -----        (- k)  (3 k)    |
                  1/3 |  \          9      x         |
                 3    |   )     ---------------------|
                      |  /      pochhammer(2/3, k) k!|
                      | -----                        |
                      \ k = 0                        /
             1/3 -------------------------------------
                               GAMMA(2/3)

                                     /infinity                      \
                                     | -----      (- k)  (3 k + 1)  |
                      1/6            |  \        9      x           |
                     3    GAMMA(2/3) |   )     ---------------------|
                                     |  /      pochhammer(4/3, k) k!|
                                     | -----                        |
                                     \ k = 0                        /
               - 1/2 ------------------------------------------------
                                            Pi
\end{maple}
For the next example the internal series command is not able to compute a
truncated power series as long a s the parameter $a$ is left as an unknown.
\begin{maple}
> FormalPowerSeries(x^a*sin(x^2),x);

                         infinity
                          -----       k  (4 k + 2 + a)
                           \      (-1)  x
                            )     --------------------
                           /           (1 + 2 k)!
                          -----
                          k = 0
\end{maple}


The following example has a logarithmic singularity at $x=0$ and hence the
FPS of $f'$ is computed.
\begin{maple}
> FormalPowerSeries(arcsech(x), x, right);

                              /infinity                             \
                              | -----               (- k)  (2 k + 2)|
                              |  \      (1 + 2 k)! 4      x         |
          ln(2) - ln(x) - 1/2 |   )     ----------------------------|
                              |  /             2                    |
                              | -----      (k!)  (k + 1) (2 k + 2)  |
                              \ k = 0                               /
\end{maple}

The next two examples are FPS expansions of two definite integrals, namely
the polylogarithm function and $\int\limits_0^x \erf(t)/t\,dt$. Note that for
the second integral we have used the inert function of {\tt Int} which is
only a placeholder and does not try to compute a closed form.
\begin{maple}
> FormalPowerSeries(dilog(1-x), x);

                               infinity
                                -----    (k + 1)
                                 \      x
                                  )     --------
                                 /             2
                                -----   (k + 1)
                                k = 0

> FormalPowerSeries(Int(erf(t)/t,t=0..x), x);

                            infinity
                             -----       k  (2 k + 1)
                              \      (-1)  x
                               )     ----------------
                              /                    2
                             -----     k! (2 k + 1)
                             k = 0
                          2 -------------------------
                                        1/2
                                      Pi
\end{maple}

The next two examples are FPS of functions which contain orthogonal
polynomials.
\begin{maple}
> FormalPowerSeries(exp(-x)*LaguerreL(n,a,x),x);
FPS/hypergeomRE:   provided that    -1    <=    min(-1,-a-1)

                       /infinity                                  \
                       | -----       k                           k|
                       |  \      (-1)  pochhammer(n + 1 + a, k) x |
    binomial(n + a, n) |   )     ---------------------------------|
                       |  /           pochhammer(1 + a, k) k!     |
                       | -----                                    |
                       \ k = 0                                    /

> FormalPowerSeries(exp(-x^2)*HermiteH(n,x),x);

                            /infinity                                \
                            | -----                          k  (2 k)|
                            |  \      (1/2 n + 1/2 + k)! (-4)  x     |
     cos(1/2 n Pi) (n + 1)! |   )     -------------------------------|
                            |  /            (n + 1 + 2 k) (2 k)!     |
                            | -----                                  |
                            \ k = 0                                  /
     -----------------------------------------------------------------
                          (1/2 n)! (1/2 n + 1/2)!

                                   /infinity                              \
                                   | -----       k               (2 k + 1)|
                                   |  \      (-4)  (1/2 n + k)! x         |
            sin(1/2 n Pi) (n + 1)! |   )     -----------------------------|
                                   |  /                (2 k + 1)!         |
                                   | -----                                |
                                   \ k = 0                                /
          + ---------------------------------------------------------------
                                (1/2 n + 1/2)! (1/2 n)!
\end{maple}

As we have seen above, we have three tools available to compute a formal
power series representation of a given function, and it may happen for some
examples, that several of these tools may be applied. Normally the solutions
of a RE of hypergeometric type leads to the simplest results, but this is not
true in general. Hence we added the possibility for the user to choose a
method by adding an additional argument which is either {\tt hypergeometric},
{\tt explike} or {\tt rational}.

\begin{maple}
> f := x*arctan(x)-1/2*ln(1+x^2):
> FormalPowerSeries(f, x, hypergeometric);

                            /infinity                  \
                            | -----        k  (2 k + 2)|
                            |  \       (-1)  x         |
                        1/2 |   )     -----------------|
                            |  /      (k + 1) (2 k + 1)|
                            | -----                    |
                            \ k = 0                    /

> FormalPowerSeries(f, x, rational);

                     infinity
                      -----               k    k   (k + 2)
                       \            ((- I)  + I ) x
                        )     1/2 -------------------------
                       /           k      k
                      -----       I  (- I)  (k + 1) (k + 2)
                      k = 0

> f := sin(x)*exp(x):
> FormalPowerSeries(f, x, explike);

                       infinity
                        -----     k 1/2                k
                         \      (2 )    sin(1/4 k Pi) x
                          )     ------------------------
                         /                 k!
                        -----
                        k = 0

> FormalPowerSeries(f, x, hypergeometric);

                /infinity                              \
                | -----       k   (- k)    k  (4 k + 1)|
                |  \      (-1)  64      256  x         |
                |   )     -----------------------------|
                |  /                (4 k + 1)!         |
                | -----                                |
                \ k = 0                                /

                       /infinity                              \
                       | -----       k   (- k)    k  (4 k + 2)|
                       |  \      (-1)  64      256  x         |
                     + |   )     -----------------------------|
                       |  /           (2 k + 1) (4 k + 1)!    |
                       | -----                                |
                       \ k = 0                                /

                       /infinity                               \
                       | -----        k   (- k)    k  (4 k + 3)|
                       |  \       (-1)  64      256  x         |
                     + |   )     ------------------------------|
                       |  /      (2 k + 1) (4 k + 1)! (4 k + 3)|
                       | -----                                 |
                       \ k = 0                                 /
\end{maple}

\heading{Asymptotic Series}\noindent
The same algorithm can also be used to compute asymptotic expansions. Note,
that we only look for asymptotic series of the form of a Laurent-Puiseux
series. 

One special thing of Laurent-Puiseux asymptotic expansions  is, that they are
only valid as long as the indeterminate approaches the expansion point from
one side. If the expansion point is not $\infty$, then one may specify with
an option {\tt right} or {\tt left} from which side one approaches the
expansion point.
\begin{maple}
> FormalPowerSeries(erf(x), x=infinity);

                                       1

> FormalPowerSeries(arctan(1/x), x=0, right);

                               /infinity                 \
                               | -----       k  (2 k + 1)|
                               |  \      (-1)  x         |
                      1/2 Pi - |   )     ----------------|
                               |  /           2 k + 1    |
                               | -----                   |
                               \ k = 0                   /

> FormalPowerSeries(exp(x), x=-infinity);

                                       0

> FormalPowerSeries(exp(x), x=infinity):
FPS/FPS:   ERROR:    essential singularity

> FormalPowerSeries(exp(x)*Ei(-x) + exp(-x)*Ei(x), x=infinity);

                       /infinity                          \
                       | -----                            |
                       |  \                      (2 k + 2)|
                     2 |   )     (1 + 2 k)! (1/x)         |
                       |  /                               |
                       | -----                            |
                       \ k = 0                            /

> FormalPowerSeries(exp(x)*(1-erf(sqrt(x))), x=infinity);

                infinity
                 -----       (- k)         (- k)      (1/2 + k)
                  \      (-1)      (2 k)! 4      (1/x)
                   )     --------------------------------------
                  /                        k!
                 -----
                 k = 0
                -----------------------------------------------
                                       1/2
                                     Pi
\end{maple}
By a plot of {\tt arctan(1/x)}, e.g, one may realize that, indeed, the
resulting series representation is one sided, only.

\heading{Conclusion}\noindent
We have presented an algorithm, and its implementation in \Maple, to compute
FPS. Since it is a goal of Computer Algebra to work with {\em formal\/}
objects, we think this is a very powerful and valuable addition.

The algorithm is beside of rational operations based on two basic tools:
solving systems of equations and computing symbolic limits. No truncated
series is ever computed. For the case of asymptotic series, one sided limits
are computed and for all other cases complex ones. The power of the procedure
for computing LPS stands and falls with the capabilities of the tool for
computing limits. The algorithms which are used in \Maple\ to compute limits
are described in~\cite{Geddes88,GonnetGruntz92}.

Further, in cases when the resulting RE cannot be solved explicitly, it can,
in principle be used to calculate the coefficients iteratively in a lazy
evaluation scheme. This is particularly efficient as the resulting RE always
is homogeneous and linear, so that each coefficient can be calculated by
finitely many of its predecessors. We note that the algorithm to find a
holonomic DE works, and so such a RE always exists, if the input function is
constructed by integration, differentiation, addition, multiplication, and
the composition with rational functions or rational powers, from functions
with the same property~\cite{Stanley80,Salvy92,Koediffgl}. This gives a huge
class of functions to which the method can be applied.

\begin{thebibliography}{99}
\bibitem{AS}
M. Abramowitz and I.A. Stegun, 
{\em Handbook of Mathematical Functions,\/}
Dover Publ., New York,
1964.

\bibitem{Bronstein}
M. Bronstein,
{\em Simplification of Real Elementary Functions,\/}
ISSAC'89, ed.\ G.H. Gonnet,
207 -- 211,
1989.

\bibitem{Salvy}
M. Bronstein and B. Salvy,
{\em Full Partial Fraction Decomposition of Rational Functions,\/}
ISSAC'93, ed.\ M. Bronstein, 
157 -- 160,
1993.

\bibitem{Geddes88}
K.O. Geddes and G.H. Gonnet,
{\em A New Algorithm for Computing Symbolic Limits using Hierarchical
Series},
Proceedings of the International Symposium ISSAC 1988,
Lecture Notes In Computer Science {\bf 358},
490 -- 495,
1988.

\bibitem{GonnetGruntz92}
G.H. Gonnet and D. Gruntz,
{\em Limit Computation in Computer Algebra\/},
Technical Report {\bf 187}, 
Department for Computer Science, ETH Z\"urich,
1992, (submitted to the JSC).

\bibitem{GruntzKoepf93}
D. Gruntz and W. Koepf,
{\em Formal Power Series,\/}
Preprint SC 93-31, ZIB, Konrad-Zuse-Zentrum f\"ur Informationstechnik Berlin,
1993.

\bibitem{Gruntz94}
D. Gruntz,
{\em Infinite Structures in Maple,\/}
MapleTech {\bf 1} (2),
pp.~19--30,
1994.

\bibitem{Han}
E.R. Hansen,
{\em A table of series and products,\/}
Prentice-Hall, Englewood Cliffs, NJ,
1975.

\bibitem{Koe92}
W. Koepf,
{\em Power series in Computer Algebra,\/}
J.\ Symb.\ Comp.\ {\bf 13},
581 -- 603,
1992.

\bibitem{Koe93a}
W. Koepf,
{\em Examples for the algorithmic calculation of formal Puiseux, Laurent and
power series,\/}
SIGSAM Bulletin {\bf 27}, 20 -- 32,
1993.

\bibitem{Koe93}
W. Koepf,
{\em Algorithmic development of power series\/},
In: Artificial intelligence and symbolic mathematical computing,
Proceedings of the International Conference AISMC-1, Karlsruhe, Germany,
August 1992, 
(J.~Calmet and J.A.~Campbell, ed.),
Lecture Notes in Computer Science {\bf 737}, Springer-Verlag
Berlin--Heidelberg,
195 -- 213,
1993.

\bibitem{Koediffgl}
W. Koepf and D. Schmersau,
{\em Spaces of functions satisfying simple differential equations,\/}
Konrad-Zuse-Zentrum Berlin (ZIB), Technical Report TR 94-2, 1994.

\bibitem{Koeortho}
W. Koepf,
{\em Algorithmic work with orthogonal polynomials and special functions,\/}
Konrad-Zuse-Zentrum Berlin (ZIB), Technical Report TR 94-5, 1994.

\bibitem{Risch}
R. Risch,
{\em Algebraic Properties of the Elementary Functions of Analysis,\/}
Amer.\ J.\ of Math.\ {\bf 101},
743 -- 759,
1979.

\bibitem{Salvy92}
B. Salvy and P. Zimmermann,
{\em GFUN: A package for the manipulation of generating and holonomic
functions in one variable\/},
AMS Transactions on Mathematical Software {\bf 20}, pp.~163--177, 1994.

\bibitem{Stanley80}
R.P. Stanley,
{\em Differentiably finite power series,\/}
Europ.~J.~Combinatorics {\bf 1}, pp.~175--188,
1980.

\end{thebibliography}
\end{document}

