reduce-historical/r36/doc/SPECFN.TEX

\documentstyle[11pt,reduce]{article}
\title{{\tt SPECFN}: Special Functions Package for REDUCE}
\date{}
\author{Chris Cannam\\[0.05in]
Konrad--Zuse--Zentrum f\"ur Informationstechnik Berlin \\
Heilbronner Strasse 10 \\
D--10711 Berlin -- Wilmersdorf \\
Federal Republic of Germany \\[0.05in]
E--mail: neun@sc.zib--berlin.de \\[0.05in]
Version 2.0, November 1994}
\begin{document}
\maketitle
\index{SPECFN package}
\section{Introduction}

This package provides the 'common' special functions for REDUCE.
The names of the operators and implementation details can be
found in this document.

Due to the enormous number of special functions
a package for special functions is never complete.
Several users pointed out that important classes
of special functions were missing in the first version.
These comments and other hints from a number of contributors
and users were very helpful.

The first version of this package was developed while the author
worked as a student exchange grantee at ZIB Berlin in 1992/93.
The package is maintained by ZIB Berlin after the author left the ZIB.
Therefore, please direct comments, hints and bug reports etc. to
{\tt neun@sc.ZIB-Berlin.de}.

This package is designed to provide algebraic and numeric manipulations of
several common special functions, namely:

\begin{itemize}
\item Bernoulli numbers and Euler numbers;
\item Stirling numbers;
\item Binomial Coefficients;
\item Pochhammer notation;
\item The Gamma function;
\item The psi function and its derivatives;
\item The Riemann Zeta function;
\item The Bessel functions J and Y of the first and second kinds;
\item The modified Bessel functions I and K;
\item The Hankel functions H1 and H2;
\item The Kummer hypergeometric functions M and U;
\item The Beta function, and Struve, Lommel and Whittaker functions;
\item The Airy funcions;
\item The Exponential Integral, the Sine and Cosine Integrals;
\item The Hyperbolic Sine and Cosine Integrals;
\item The Fresnel Integrals and the Error function;
\item The Dilog function;
\item Hermite Polynomials;
\item Jacobi Polynomials;
\item Legendre Polynomials;
\item Associated Legendre Functions (Spherical and Solid Harmonics)
\item Laguerre Polynomials;
\item Chebyshev Polynomials;
\item Gegenbauer Polynomials;
\item Euler Polynomials;
\item Bernoulli Polynomials;
\item (Jacobi's) Elliptic Functions;
\item Elliptic Integrals;
\item 3j and 6j symbols , Clebsch-Gordan coefficients;
\item and some well-known constants.
\end{itemize}

All algorithms whose sources are uncredited are culled from series or
expressions found in the Dover Handbook of Mathematical
Functions\cite{Abramowitz:72}.

There is a nice collection of plot calls for special functions
in the file \$reduce/plot/specplot.tst.
These examples will reproduce a number of well-known pictures from
\cite{Abramowitz:72}.

\section{Compatibility with earlier \REDUCE\ versions}

For {PSL} versions,
this package is intended to be used with the new \REDUCE\ bigfloat
mechanisms which is distributed together with REDUCE 3.5 and later versions.
The package does work with the earlier bigfloat implementations,
but in order to
ensure that it works efficiently with the new versions, it has not been
optimized for the old.

\section{Simplification and Approximation}

All of the operators supported by this package have certain algebraic
simplification rules to handle special cases, poles, derivatives and so
on.  Such rules are applied whenever they are appropriate.  However, if
the {\tt ROUNDED} switch is on, numeric evaluation is also carried out.
Unless otherwise stated below, the result of an application of a special
function operator to real or complex numeric arguments in rounded mode
will be approximated numerically whenever it is possible to do so.  All
approximations are to the current precision.

Most algebraic simplifications within the special function package
are defined in the form of a \REDUCE\ ruleset. Therefore, in order to
get a quick insight into the simplification rules one can use the
ShowRules operator, e.g.\\

\begin{verbatim}

ShowRules BesselI;

                          1          ~z     - ~z
{besseli(~n,~z) => ---------------*(e   - e     )
                    sqrt(pi*2*~z)

                           1
  when numberp(~n) and ~n=---,
                           2

                          1          ~z     - ~z
 besseli(~n,~z) => ---------------*(e   + e     )
                    sqrt(pi*2*~z)

                              1
  when numberp(~n) and ~n= - ---,
                              2

 besseli(~n,~z) => 0

  when numberp(~z) and ~z=0 and numberp(~n) and ~n neq 0,

 besseli(~n,~z) => besseli( - ~n,~z) when numberp(~n)

  and impart(~n)=0 and ~n=floor(~n) and ~n<0,

 besseli(~n,~z) => do*i(~n,~z)

  when numberp(~n) and numberp(~z) and *rounded,

 df(besseli(~n,~z),~z)

      besseli(~n - 1,~z) + besseli(~n + 1,~z)
  => -----------------------------------------,
                         2

 df(besseli(~n,~z),~z)

  => besseli(1,~z) when numberp(~n) and ~n=0}

\end{verbatim}

Several \REDUCE\ packages (such as Sum or Limits) obtain different
(hopefully better)
results for the algebraic simplifications when the SPECFN package
is loaded, because the later package contains some information which
may be useful and directly applicable for other packages, e.g.:

\begin{verbatim}
sum(1/k^s,k,1,infinity); % will be evaluated to

zeta(s)
\end{verbatim}

\ttindex{savesfs}
A record is kept of all values previously approximated, so that should a
value be required which has already been computed to the current
precision or greater, it can be simply looked up.  This can result in
some storage overheads, particularly if many values are computed which
will not be needed again.  In this case, the switch {\tt savesfs} may be
turned off in order to inhibit the storage of approximated values.  The
switch is on by default.


\section{Constants}
\ttindex{Euler\_Gamma}\ttindex{Khinchin}\ttindex{Golden\_Ratio}
\ttindex{Catalan}

Some well-known constants are defined in the special function package.
Important properties of these constants which can be used to define them
are also known. Numerical values are computed  at arbitrary precision
if the switch ROUNDED is on.

\begin{itemize}
\item Euler\_Gamma : Euler's constants, also available as -$\psi(1)$;
\item Catalan : Catalan's constant;
\item Khinchin : Khinchin's constant , defined in  \cite{Khinchin:64}.
 (takes a lot of time to compute);
\item Golden\_Ratio : $\frac{1 + \sqrt{5}}{2}$
\end{itemize}


\section{Bernoulli Numbers and Euler Numbers}
\ttindex{Bernoulli}\index{Bernoulli numbers}
\ttindex{Euler}\index{Euler numbers}

The unary operator {\tt Bernoulli} provides notation and computation for
Bernoulli numbers.  {\tt Bernoulli(n)} evaluates to the $n$th Bernoulli
number; all of the odd Bernoulli numbers, except {\tt Bernoulli(1)}, are
zero.

The algorithms are based upon those by Herbert Wilf, presented by Sandra
Fillebrown \cite{Fillebrown:92}.  If the {\tt ROUNDED} switch is off,
the algorithms are exactly those; if it is on, some further rounding may
be done to prevent computation of redundant digits.  Hence, these
functions are particularly fast when used to approximate the Bernoulli
numbers in rounded mode.

Euler numbers are computed by the unary operator Euler, which
return the nth Euler number. The computation is derived
directly from Pascal's triangle of binomial coefficients.


\section{Stirling Numbers}
\index{Stirling Numbers}\ttindex{Stirling1}\ttindex{Stirling2}
Stirling numbers of the first and second kind are computed
by the binary operators {\tt Stirling1} and {\tt Stirling2}
using explicit formulae.


\section{The $\Gamma$ Function, and Related Functions}
\ttindex{Gamma}\index{$\Gamma$ function}\index{Gamma function}

\subsection{The $\Gamma$ Function}

This is represented by the unary operator {\tt Gamma}.

Initial transformations applied with {\tt ROUNDED} off are: $\Gamma(n)$ for
integral $n$ is computed, $\Gamma(n+1/2)$ for integral $n$ is rewritten to
an expression in $\sqrt\pi$, $\Gamma(n+1/m)$ for natural $n$ and $m$ a
positive integral power of 2 less than or equal to 64 is rewritten to an
expression in $\Gamma(1/m)$, expressions with arguments at which there is a
pole are replaced by {\tt INFINITY}, and those with a negative (real)
argument are rewritten so as to have positive arguments.

The algorithm used for numerical approximation is an implementation of an
asymptotic series for $ln(\Gamma)$, with a scaling factor obtained from
the Pochhammer functions.

An expression for $\Gamma'(z)$ in terms of $\Gamma$ and $\psi$ is
included.



\subsection{The Pochhammer Notation}
\ttindex{Pochhammer}\index{Pochhammer notation}

The Pochhammer notation $(a)_k$ is supported by the binary operator {\tt
Pochhammer}.  With {\tt ROUNDED} off, this expression is evaluated
numerically if $a$ and $k$ are both integral, and otherwise may be
simplified where appropriate.  The simplification rules are based upon
algorithms supplied by Wolfram Koepf \cite{Koepf:92}.



\subsection{The Digamma Function, $\psi$}
\ttindex{PSI}\index{$\psi$ function}\index{psi function}\index{Digamma
function}

This is represented by the unary operator {\tt PSI}.

Initial transformations for $\psi$ are applied on a similar basis to
those for $\Gamma$; where possible, $\psi(x)$ is rewritten in
terms of $\psi(1)$ and $\psi(\frac{1}{2})$, and expressions with negative
arguments are rewritten to have positive ones.

Numerical evaluation of $\psi$ is only carried out if the argument is
real. The algorithm used is based upon an asymptotic series, with a
suitable scaler.

Relations for the derivative and integral of $\psi$ are included.


\subsection{The Polygamma Functions, $\psi^{(n)}$}
\ttindex{Polygamma}\index{$\psi^{(n)}$ functions}\index{Polygamma
functions}

The $n$th derivative of the $\psi$ function is represented by the
binary operator {\tt Polygamma}, whose first argument is $n$.

Initial manipulations on $\psi^{(n)}$ are few; where the second argument
is $1$ or $3/2$, the expression is rewritten to one involving the
Riemann $\zeta$ function, and when the first is zero it is rewritten to
$\psi$; poles are also handled.

Numerical evaluation is only carried out with real arguments. The
algorithm used is again an asymptotic series with a scaling factor; for
negative (second) arguments, a Reflection Formula is used, introducing a
term in the $n$th derivative of $\cot(z\pi)$.

Simple relations for derivatives and integrals are provided.


\subsection{The Riemann $\zeta$ Function}
\ttindex{Zeta}\index{Riemann Zeta function}\index{Zeta
function}\index{$\zeta$ function}

This is represented by the unary operator {\tt Zeta}.

With {\tt ROUNDED} off, $\zeta(z)$ is evaluated numerically for even
integral arguments in the range $-31 < z < 31$, and for odd integral
arguments in the range $-30 < z < 16$.  Outside this range the values
become a little unwieldy.

Numerical evaluation of $\zeta$ is only carried out if the argument is real.
The algorithms used for $\zeta$ are: for odd integral arguments, an
expression relating $\zeta(n)$ with $\psi^{n-1}(3)$; for even arguments, a
trivial relationship with the Bernoulli numbers; and for other arguments the
approach is either (for larger arguments) to take the first few primes in
the standard over-all-primes expansion, and then continue with the defining
series with natural numbers not divisible by these primes, or (for smaller
arguments) to use a fast-converging series obtained from \cite{Bender:78}.

There are no rules for differentiation or integration of $\zeta$.


\section{Bessel Functions}
\ttindex{BesselJ}\ttindex{BesselY}\ttindex{BesselI}\ttindex{BesselK}\ttindex{Hankel1}\ttindex{Hankel2}\index{Bessel
functions}\index{Hankel functions}

Support is provided for the Bessel functions J and Y, the modified
Bessel functions I and K, and the Hankel functions of the first and
second kinds.  The relevant operators are, respectively, {\tt BesselJ},
{\tt BesselY}, {\tt BesselI}, {\tt BesselK}, {\tt Hankel1} and {\tt
Hankel2}, which are all binary operators.

The following initial transformations are performed:

\begin{itemize}
\item trivial cases or poles of J, Y, I and K are handled;
\item J, Y, I and K with negative first argument are transformed to have
positive first argument;
\item J with negative second argument is transformed for positive second
argument;
\item Y or K with non-integral or complex second argument is transformed
into an expression in J or I respectively;
\item derivatives of J, Y and I are carried out;
\item derivatives of K with zero first argument are carried out;
\item derivatives of Hankel functions are carried out.
\end{itemize}
Also, if the {\tt COMPLEX} switch is on and {\tt ROUNDED} is off,
expressions in Hankel functions are rewritten in terms of Bessels.

No numerical approximation is provided for the Bessel K function, or for
the Hankel functions for anything other than special cases.  The
algorithms used for the other Bessels are generally implementations of
standard ascending series for J, Y and I, together with asymptotic
series for J and Y; usually, the asymptotic series are tried first, and
if the argument is too small for them to attain the current precision,
the standard series are applied.  An obvious optimization prevents an
attempt with the asymptotic series if it is clear from the outset that
it will fail.

There are no rules for the integration of Bessel and Hankel functions.


\section{Hypergeometric and Other Functions}
\ttindex{Beta}\ttindex{KummerM}\ttindex{KummerU}\ttindex{StruveH}
\ttindex{StruveL}\ttindex{Lommel1}\ttindex{Lommel2}\ttindex{WhittakerM}
\ttindex{WhittakerW}\index{Beta function}\index{$B$ function}
\index{Kummer functions}\index{Struve functions}\index{Lommel functions}
\index{Whittaker functions}\index{Hypergeometric functions}

This package also provides some support for other functions, in the form
of algebraic simplifications:
\begin{itemize}
\item The Beta function, a variation upon the $\Gamma$
function\cite{Abramowitz:72}, with the binary operator {\tt Beta};
\item The Struve {\bf H} and {\bf L} functions, through the binary
operators {\tt StruveH} and {\tt StruveL}, for which manipulations are
provided to handle special cases, simplify to more readily handled
functions where appropriate, and differentiate with respect to the second
argument;
\item The Lommel functions of the first and second kinds, through the
ternary operators {\tt Lommel1} and {\tt Lommel2}, for which manipulations
are provided to handle special cases and simplify where appropriate;
\item The Kummer confluent hypergeometric functions M and U (the
hypergeometric ${_1F_1}$ or $\Phi$, and $z^{-a}{_2F_0}$ or $\Psi$,
respectively),
with the ternary operators {\tt KummerM} and {\tt KummerU}, for which
there are manipulations for special cases and simplifications, derivatives
and, for the M function, numerical approximations for real arguments;
\item The Whittaker M and W functions, variations upon the Kummer
functions, which, with the ternary operators {\tt WhittakerM} and {\tt
WhittakerW}, simplify to expressions in the Kummer functions.
\end{itemize}

\section{Integral Functions}

The SPECFN package includes manipulation and a limited numerical
evaluation for some Integral functions, namely

erf, erfc, Si, Shi, si, Ci, Chi, Ei, Fresnel\_C and Fresnel\_S.

The definitions from integral, the derviatives and some limits are
known together with some simple properties such as symmetry
conditions.

The numerical approximation for the Integral functions suffer
from the fact that the precision is not set correctly for
values of the argument above 10.0 (approx.) and from the
usage of summations even for large arguments.

\section{Airy Functions}
\ttindex{Airy\_Ai}\ttindex{Airy\_Bi}\ttindex{Airy\_Aiprime}
\ttindex{Airy\_Biprime}\index{Airy Functions}

Support is provided for the Airy Functions Ai and Bi and for the
Airyprime Functions Aiprime and Biprime. The relevant operators are
respectively {\tt Airy\_Ai}, {\tt Airy\_Bi}, {\tt Airy\_Aiprime}, and
{\tt Airy\_Biprime}, which are all unary operators with one argument.

The following cases can be performed:

\begin{itemize}
\item Trivial cases of Airy\_Ai and Airy\_Bi and their primes are handled.
\item All cases can handle both complex and real arguments.
\item The Airy Functions can also be represented in terms of Bessel
Functions by activating an inactive rule set.
\end{itemize}

In order to activate the Airy Function to Bessel Rules one should type: \\
{\tt let Airy2Bessel\_rules;}. As a result the Airy\_Ai function,
for example will be calculated using the formula :- \\
\\
{\tt Ai(z) = } $\frac{1}{3}$\( \sqrt{z} \)[{\bf {\sl I}}$_{-1/3}$($\zeta$)
- {\bf {\sl I}}$_{1/3}$({$\zeta$})] , where
 $\zeta$ =  $\frac{2}{3} z^{\frac{2}{3}}$\\
\\
\underline{{\tt Note}}:- In order to obtain satisfactory approximations
to results both the {\tt COMPLEX} and {\tt ROUNDED} switches must be on.

The algorithms used for the Airy Functions are implementations of
standard ascending series, together with asymptotic series. At some
point it is better to use the asymptotic approach, rather than the
series. This value is calculated by the program and depends on the given
precision.

There are no rules for the integration of Airy Functions.

\section{Polynomial Functions}

Two groups are defined, some well-known orthogonal Polynomials
(Hermite, Jacobi, Legendre, Laguerre, Chebyshev, Gegenbauer)
and Euler and Bernoulli Polynomials.
The names of the \REDUCE\  operator are build by adding a P
to the name
of the polynomials, e.g. EulerP implements the Euler polynomials.
Most definitions are equivalent to \cite{Abramowitz:72}, except
for the ternary (associated) Legendre Polynomials.

\begin{verbatim}
P(n,m,x) = (-1)^m *(1-x^2)^(m/2)*df(legendreP (n,x),x,m)
\end{verbatim}

\section{Spherical and Solid Harmonics}
\ttindex{SolidHarmonicY} \ttindex{SphericalHarmonicY}

The relevant operators are, respectively,\\
{\tt SolidHarmonicY} and {\tt SphericalHarmonicY}.

The SolidHarmonicY operator implements the Solid Harmonics
described below. It expects 6 parameter, namely n,m,x,y,z and r2
and returns a polynomial in x,y,z and r2.

The operator
SphericalHarmonicY is a special case of SolidHarmonicY
with the usual definition:

\begin{verbatim}
algebraic procedure SphericalHarmonicY(n,m,theta,phi);
        SolidHarmonicY(n,m,sin(theta)*cos(phi),
                sin(theta)*sin(phi),cos(theta),1)$
\end{verbatim}


Solid Harmonics of order n (Laplace polynomials)
are homogeneous polynomials of degree n in x,y,z
which are solutions of Laplace equation:-

\begin{verbatim}
       df(P,x,2) + df(P,y,2) + df(P,z,2) = 0.
\end{verbatim}

There are 2*n+1 independent such polynomials for any given $n >=0$
and with:-

\begin{verbatim}
       w!0 = z, w!+ = i*(x-i*y)/2, w!- = i*(x+i*y)/2,
\end{verbatim}

they are given by the Fourier integral:-

\begin{verbatim}
 S(n,m,w!-,w!0,w!+) =

       (1/(2*pi)) *
       for u:=-pi:pi integrate (w!0 + w!+ * exp(i*u) + w!- *
           exp(-i*u))^n * exp(i*m*u) * du;
\end{verbatim}

which is obviously zero if $|m| > n$ since then all terms in the
expanded integrand contain the factor exp(i*k*u) with k neq 0,

S(n,m,x,y,z) is proportional to
\begin{verbatim}
     r^n * Legendre(n,m,cos theta) * exp(i*phi)
\end{verbatim}

Let r2 = $x^2 + y^2 + z^2$ and r = sqrt(r2).

The spherical harmonics are simply the restriction of the solid
harmonics to the surface of the unit sphere and the set of all
spherical harmonics {$n >=0; -n <= m =< n$} form a complete orthogonal
basis on it, i.e. $<n,m|n',m'>$ = Kronecker\_delta(n,n') *
Kronecker\_delta(m,m') using
$<...|...>$ to designate the scalar product
of functions over the spherical surface.

The coefficients of the solid harmonics are normalised in what
follows to yield an ortho-normal system of spherical harmonics.

Given their polynomial nature, there are many recursions formulae
for the solid harmonics and any recursion valid for Legendre functions
can be 'translated' into solid harmonics. However the direct proof is
usually far simpler using Laplace's definition.

It is also clear that all differentiations of solid harmonics are
trivial, qua polynomials.

Some substantial reduction in the symbolic form would occur if one
maintained throughout the recursions the symbol r2 (r cannot occur
as it is not rational in x,y,z). Formally the solid harmonics appear
in this guise as more compact polynomials in (x,y,z,r2).

Only two recursions are needed:-

(i) along the diagonal (n,n);

(ii) along a line of constant n: (m,m),(m+1,m),...,(n,m).

Numerically these recursions are stable.

For $m < 0$ one has:-

\begin{verbatim}
      S(n,m,x,y,z) = (-1)^m * S(n,-m,x,-y,z);
\end{verbatim}

\section{Jacobi's Elliptic Functions}

The following functions have been implemented:

\begin{itemize}
\item The Twelve Jacobi Functions
\item Arithmetic Geometric Mean
\item Descending Landen Transformation
\end{itemize}

\subsection{Jacobi Functions}

The following Jacobi functions are available:-

\begin{itemize}
\item Jacobisn(u,m)
\item Jacobidn(u,m)
\item Jacobicn(u,m)
\item Jacobicd(u,m)
\item Jacobisd(u,m)
\item Jacobind(u,m)
\item Jacobidc(u,m)
\item Jacobinc(u,m)
\item Jacobisc(u,m)
\item Jacobins(u,m)
\item Jacobids(u,m)
\item Jacobics(u,m)
\end{itemize}

They will be evaluated numerically if the {\tt rounded} switch is used.

\subsection{Amplitude}

The amplitude of u can be evaluated using the
{\tt JacobiAmplitude(u,m)} command.

\subsection{Arithmetic Geometric Mean}

A procedure to evaluate the AGM of initial values \(a_0,b_0,c_0\)
exists as \\
{\tt AGM\_function(\(a_0,b_0,c_0\))} and will return \\
$\{ N, AGM, \{ a_N, \ldots ,a_0\}, \{ b_N, \ldots ,b_0\},
\{c_N, \ldots ,c_0\}\}$,
where N is the number of steps to compute the AGM to the
desired acuracy. \\

To determine the Elliptic Integrals K($m$), E($m$) we use initial values
\(a_0 = 1\); \(b_0 = \sqrt{1-m}\) ; \(c_0 = \sqrt{m}\).


\subsection{Descending Landen Transformation}

The procedure to evaluate the Descending Landen Transformation of
phi and alpha uses the following equations:\\
\indent \ \ \ \ \( (1+sin \alpha_{n+1})(1+cos \alpha_n)=2 \)
\ \ \ \ where \(\alpha_{n+1}<\alpha_n\) \\
\indent \ \ \ \ \(tan(\phi_{n+1}-\phi_n)=cos \alpha_n tan \phi_n \)
\ \ \ where \(\phi_{n+1}>\phi_n\) \\

It can be called using {\tt landentrans($\phi_0$,$\alpha_0$)}
and will return \\
$\{\{\phi_0, \ldots ,\phi_n\},\{\alpha_0, \ldots ,\alpha_n\}\}$.

\section{Elliptic Integrals}
The following functions have been implemented:

\begin{itemize}
\item Elliptic Integrals of the First Kind
\item Elliptic Integrals of the Second Kind
%\item Ellpitic Integrals of the Third Kind
\item Jacobi $\theta$ Functions
\item Jacobi $\zeta$ Function
\end{itemize}

\subsection{Elliptic F}

The Elliptic F function can be used as {\tt EllipticF($\phi$,m)} and
will return the value of the {\underline {Elliptic Integral of the
First Kind}}.

\subsection{Elliptic K}

The Elliptic K function can be used as {\tt EllipticK(m)} and will
return the value of the {\underline {Complete Elliptic Integral of the
First Kind}}, K. It is often used in the calculation of other elliptic
functions

\subsection{Elliptic K$'$}

The Elliptic K$'$ function can be used as {\tt EllipticK!$'$(m)} and will
return the value K($1-m$).

\subsection{Elliptic E}

The Elliptic E function comes with two different numbers of arguments:

It can be used with two arguments as {\tt EllipticE($\phi$,m)}
and will return the value
of the {\underline {Elliptic Integral of the Second Kind}}.

The Elliptic E function can also be used as {\tt EllipticE(m)} and
will return the value of the {\underline {Complete Elliptic Integral
of the Second Kind}}, E.

%\section{Ellpitic $\Pi$}
%
%The Elliptic $\pi$ function can be used as {\tt EllipticPi( )} and
%will return the value of the {\underline {Elliptic Integral of the
%Third Kind}}.
%
\subsection{Elliptic $\Theta$ Functions}

This can be used as {\tt EllipticTheta(a,u,m)}, where $a$ is the index
for the theta functions ($a = 1,2,3$ or $4$) and will return $H$;
$H_1$; $\Theta_1$; $\Theta$. (Also denoted in some texts as
$\vartheta_1$; $\vartheta_2$; $\vartheta_3$; $\vartheta_4$.)

\subsection{Jacobi's Zeta Function Z }

This can be used as {\tt JacobiZeta(u,m)} and will return
Jacobi Zeta. Note: the operator {\tt Zeta} will invoke Riemann's
$\zeta$ function.

\section{Lambert's W function}

Lambert's W function is the inverse of the function  $w*e^w$.
Therefore it is an important contribution for the solve package.
The function is studied extensively in \cite{Corless:92}.
The current implementation will compute the principal branch in
ROUNDED mode only.

\section{3j symbols and Clebsch-Gordan Coefficients}
The operators {\tt ThreeJSymbol}, {\tt Clebsch\_Gordan} are defined like
in \cite{Landolt:68} or \cite{Edmonds:57}. ThreeJSymbol expects as arguments
three lists of values \{$j_i,m_i$\}, e.g.

\begin{verbatim}
ThreeJSymbol({J+1,M},{J,-M},{1,0});
Clebsch_Gordan({2,0},{2,0},{2,0});
\end{verbatim}

\section{6j symbols }
The operator {\tt SixJSymbol} is defined like
in \cite{Landolt:68} or \cite{Edmonds:57}.
SixJSymbol expects two lists of values \{$j_1,j_2,j_3$\} and
 \{$l_1,l_2,l_3$\} as arguments, e.g.

\begin{verbatim}
SixJSymbol({7,6,3},{2,4,6});
\end{verbatim}

In the current implementation of the SixJSymbol, there is only a limited
reasoning about the minima and maxima of the summation using
the INEQ package, such that in most
cases the special 6j-symbols (see e.g. \cite{Landolt:68})
will not be found.



\section{Acknowledgements}

The contributions of Kerry Gaskell, Matthew Rebbeck, Lisa Temme
and Stephen Scowcroft (all students from the University of Bath
on placement in ZIB Berlin for one year) were very helpful
to augment the package. The advise of Ren\'e Grognard (CSIRO , Australia)
for the development of the module for Clebsch-Gordan and 3j, 6j symbols
and the module for spherical and solid harmonics was very much appreciated.


\section{Table of Operators}

\fbox{
\begin{tabular}{r l}\\
Function & Operator \\\\
%\hline
$\left( { n \atop m } \right)$ & {\tt Binomial(n,m)} \\
Bernoulli($n$) or $ B_n $ & {\tt Bernoulli(n)} \\
Euler($n$) or $ E_n $ & {\tt Euler(n)} \\
$S_n^{(m)}$ & {\tt Stirling1(n,m)} \\
${\bf S}_n^{(m)}$ & {\tt Stirling2(n,m)} \\
$\Gamma(z)$ & {\tt Gamma(z)} \\
$(a)_k$ & {\tt Pochhammer(a,k)} \\
$\psi(z)$ & {\tt Psi(z)} \\
$\psi^{(n)}(z)$ & {\tt Polygamma(n,z)} \\
$Riemann's \zeta(z)$ & {\tt Zeta(z)} \\
$J_\nu(z)$ & {\tt BesselJ(nu,z)}\\
$Y_\nu(z)$ & {\tt BesselY(nu,z)}\\
$I_\nu(z)$ & {\tt BesselI(nu,z)}\\
$K_\nu(z)$ & {\tt BesselK(nu,z)}\\
$H^{(1)}_\nu(z)$ & {\tt Hankel1(n,z)}\\
$H^{(2)}_\nu(z)$ & {\tt Hankel2(n,z)}\\
$B(z,w)$ & {\tt Beta(z,w)}\\
${\bf H}_{\nu}(z)$ & {\tt StruveH(nu,z)}\\
${\bf L}_{\nu}(z)$ & {\tt StruveL(n,z)}\\
$s_{a,b}(z)$ & {\tt Lommel1(a,b,z)}\\
$S_{a,b}(z)$ & {\tt Lommel2(a,b,z)}\\
$Ai(z)$ & {\tt Airy\_Ai(z)}\\
$Bi(z)$ & {\tt Airy\_Bi(z)}\\
$Ai'(z)$ & {\tt Airy\_Aiprime(z)}\\
$Bi'(z)$ & {\tt Airy\_Biprime(z)}\\
$M(a, b, z)$ or $_1F_1(a, b; z)$ or $\Phi(a, b; z)$ &
{\tt KummerM(a,b,z)} \\
$U(a, b, z)$ or $z^{-a}{_2F_0(a, b; z)}$ or $\Psi(a, b; z)$ &
{\tt KummerU(a,b,z)}\\
$M_{\kappa,\mu}(z)$ & {\tt WhittakerM(kappa,mu,z)}\\
$W_{\kappa,\mu}(z)$ & {\tt WhittakerW(kappa,mu,z)}\\
$B_n(x)$ & {\tt BernoulliP(n,x)}\\
$E_n(x)$ & {\tt EulerP(n,x)}\\
$C_n^{(\alpha)}(x)$ & {\tt GegenbauerP(n,alpha,x)}\\
$H_n(x)$ & {\tt HermiteP(n,x)}\\
$L_n(x)$ & {\tt LaguerreP(n,x)}\\
$L_n^{(m)}(x)$ & {\tt LaguerreP(n,m,x)}\\
$P_n(x)$ & {\tt LegendreP(n,x)}\\
$P_n^{(m)}(x)$ & {\tt LegendreP(n,m,x)}\\
\end{tabular}}

\fbox{
\begin{tabular}{r l}\\
Function & Operator \\\\
%\hline
$P_n^{(\alpha,\beta)} (x)$ & {\tt JacobiP(n,alpha,beta,x)}\\
$U_n(x)$ & {\tt ChebyshevU(n,x)}\\
$T_n(x)$ & {\tt ChebyshevT(n,x)}\\
$Y_n^{m}(x,y,z,r2)$ & {\tt SolidHarmonicY(n,m,x,y,z,r2)}\\
$Y_n^{m}(\theta,\phi)$ & {\tt SphericalHarmonicY(n,m,theta,phi)}\\
$\left( {j_1 \atop m_1} {j_2 \atop m_2}
 {j_3 \atop m_3} \right)$ & {\tt ThreeJSymbol(\{j1,m1\},\{j2,m2\},\{j3,m3\})}\\
$\left( {j_1m_1j_2m_2 | j_1j_2j_3 - m_3} \right)$ &
{\tt Clebsch\_Gordan(\{j1,m1\},\{j2,m2\},\{j3,m3\})}\\
$\left\{ {j_1 \atop l_1} {j_2 \atop l_2}
 {j_3 \atop l_3} \right\}$  & {\tt SixJSymbol(\{j1,j2,j3\},\{l1,l2,l3\})}\\
\\
$Si(z)$ & {\tt Si(z) }\\
$si(z)$ & {\tt s\_i(z) }\\
$Ci(z)$ & {\tt Ci(z) }\\
$Shi(z)$ & {\tt Shi(z) }\\
$Chi(z)$ & {\tt Chi(z) }\\
$erf(z)$ & {\tt erf(z) }\\
$erfc(z)$ & {\tt erfc(z) }\\
$Ei(z)$ & {\tt Ei(z) }\\
$dilog(z)$ & {\tt dilog(z)} \\
$C(x)$ & {\tt Fresnel\_C(x)} \\
$S(x)$ & {\tt Fresnel\_S(x)} \\
\\
$sn(u|m)$ & {\tt Jacobisn(u,m)}\\
$dn(u|m)$ & {\tt Jacobidn(u,m)}\\
$cn(u|m)$ & {\tt Jacobicn(u,m)}\\
$cd(u|m)$ & {\tt Jacobicd(u,m)}\\
$sd(u|m)$ & {\tt Jacobisd(u,m)}\\
$nd(u|m)$ & {\tt Jacobind(u,m)}\\
$dc(u|m)$ & {\tt Jacobidc(u,m)}\\
$nc(u|m)$ & {\tt Jacobinc(u,m)}\\
$sc(u|m)$ & {\tt Jacobisc(u,m)}\\
$ns(u|m)$ & {\tt Jacobins(u,m)}\\
$ds(u|m)$ & {\tt Jacobids(u,m)}\\
$cs(u|m)$ & {\tt Jacobics(u,m)}\\
$F(\phi|m)$ & {\tt EllipticF(phi,m)}\\
$K(m)$ & {\tt EllipticK(m)}\\
$E(\phi|m) or E(m)$ & {\tt EllipticE(phi,m) or EllipticE(m)}\\
$H(u|m), H_1(u|m), \Theta_1(u|m), \Theta(u|m)$ & {\tt EllipticTheta(a,u,m)}\\
$\theta_1(u|m), \theta_2(u|m), \theta_3(u|m), \theta_4(u|m)$
   & {\tt EllipticTheta(a,u,m)}\\
$Z(u|m)$ & {\tt Zeta\_function(u,m)}
\end{tabular}}

\bibliography{specfn}
\bibliographystyle{plain}

\end{document}