\documentstyle[fullpage,maple,shadow,11pt]{article}
%% <SHAREFILE=program/program/program.tex >
%% <DESCRIBE>
%%                "Programming in Maple: The Basics"
%%                An introduction to programming in Maple with exercises.
%%                This gives more insight and examples than the Language
%%                Reference Manual does and covers some other topics.
%%                AUTHOR: Michael Monagan, monagan@inf.ethz.ch
%% </DESCRIBE>

\date{}
\title{Programming in Maple: The Basics.}
\author{Michael Monagan\\
Institut f\"ur Wissenschaftliches Rechnen\\
ETH-Zentrum, CH-8092 Z\"urich, Switzerland\\
monagan@inf.ethz.ch}

\begin{document}

{
\Large
\maketitle
}


\begin{abstract}
This is a tutorial on programming in Maple.
The aim is to show how you can write simple programs in Maple for
doing numerical calculations, linear algebra, and programs for
simplifying or transforming symbolic expressions or mathematical formulae.
It is assumed that the reader is familiar with using Maple
interactively as a calculator.
\end{abstract}

\tableofcontents
\section{Introduction}

A few words to those who are familiar with other programming languages.

\begin{itemize}
\item
Maple is a procedural programming language.  It also includes a number
of functional programming constructs.  If you have written programs in
Basic, Pascal, Algol, C, Lisp, or Fortran, you should be able to write
numerical programs in Maple very quickly.
\item
Maple is not strongly typed like C and Pascal.  No declarations
are required.  Maple is more like Basic and Lisp in this respect.
However types exist.  Type checking is done at run time and must
be programmed explicitly.
\item
Maple is interactive and the programming language is interpreted.
Maple is not suitable for running numerically intensive programs
because of the interpreter overhead.  Though it is suitable for
high-precision numerical calculations and as a tool for generating
numerical codes.
\end{itemize}

This document is based on Maple version V Release 3.
Maple development
continues.  New versions come out every one or two years which
contain not only changes to the mathematical capabilities of Maple,
but also changes to the programming language and user interface.
For this reason, I have tried to steer you away from constructs in the
language which are being, or are likely to be removed, and have mentioned
some things that will be in the language in future versions of Maple.

The {\em Maple V Language Reference Manual} is the main reference for
programming in Maple.  It is published by Springer-Verlag.  The ISBN
number is 0-387-87621-3 .
Other useful sources of information include
the Release Notes for Maple V Release 3, available from Waterloo
Maple Software, 
the {\em First Leaves: Tutorial
Introduction to Maple}, also published by Springer-Verlag,
and the extensive on-line documentation which is accessible
with the \verb+?+ command in Maple.
Part of the reason for writing this document is that some of the important
things I have stated here are not mentioned elsewhere, or they are buried
in the documentation.

\newpage
\subsection{Evaluation}

Before I begin, I want to point out the most important difference
between Maple and traditional programming languages.
If an identifier has not been assigned a value, then it stands for itself.
It is a {\em symbol}.  Symbols are used to represent unknowns in equations,
variables in polynomials, summation indices, etc.
Consider the Maple assignment statement

\begin{maple}
> p := x^2+4*x+4;

                                     2
                               p := x  + 4 x + 4
\end{maple}

Here the identifier $p$ has been assigned the formula $x^2 + 4 x + 4$.
The identifier $x$ has not been assigned a value, it is just a symbol,
an unknown.  The identifier $p$ has been assigned a value.  It is now
like a programming variable, and its value can be used in subsequent
calculations just like a normal programming variable.
What is the value of $p$?

\begin{maple}
> p;
                                   2
                                  x  + 4 x + 4
\end{maple}

\noindent
It is the formula $x^2 + 4 x + 4$.  What is the value of $x$?

\begin{maple}
> x;
                                       x
\end{maple}

\noindent
It is the symbol $x$.
Because a variable can be assigned a value which contains symbols,
the issue of {\em evaluation} arises.  Consider

\begin{maple}
> x := 3;
                                     x := 3

> p;
                                       25
\end{maple}

Here we assigned $x$ the value 3 and asked Maple to print
out the value of $p$.  What should Maple print out?  Should it print out
the polynomial $x^2+4 x + 4$ or should it evaluate the polynomial,
i.e. compute $3^2 + 4 \times 3 + 4$ and return the result 25?
We see that Maple does the latter.
This issue of evaluation will come up from time to time as it may
affect the efficiency and semantics of a program.
This difference between Maple and traditional programming languages,
where identifiers can be used for both programming variables and mathematical
unknowns is nice.  But be careful not to mix the two.
Many problems that users encounter have to do with using identifiers
for both symbols and programming variables.
For example, what happens if we now try to compute $\int p {\rm d} x$?

\begin{maple}
> int(p,x);
Error, (in int) wrong number (or type) of arguments
\end{maple}

An error has occurred in the integration function \verb+int+.
Here we are thinking of $x$ as a symbol, the integration variable.
But $x$ was previously assigned the integer $3$.
Maple evaluates the arguments to the \verb+int+ function
and tries to integrate 25 with respect to 3.
It doesn't make sense to integrate with respect to 3!
How does one convert $x$ back into a symbol?
In Maple one {\em unassigns} the variable $x$ by doing

\begin{maple}
> x := 'x';
                                     x := x

> int(p,x);
                                   3      2
                              1/3 x  + 2 x  + 4 x
\end{maple}

\subsection{Expressions: Sums, Products, Powers, Functions}

In Maple, mathematical formulae, e.g. things
like $\sin (x+\pi/2)$, and $x^3 y^2 - 2/3$ are called {\em expressions}.  
They are made up of symbols, numbers, arithmetic
operators and functions.  Symbols are things like \verb+sin, x, y, Pi+ etc.
Numbers include $12, \frac{2}{3}, 2.1$ etc.
The arithmetic operators are
\verb:+: (addition),
\verb:-: (subtraction),
\verb:*: (multiplication),
\verb:/: (division),
and \verb:^: (exponentiation).  And examples of functions include
\verb+sin(x)+, \verb+f(x,y)+, \verb+min(x1,x2,x3,x4)+.
For example, the formula $p ~=~ x^2 y+3 x^3 z+2$, which is a polynomial,
is input in Maple as

\begin{maple}
> p := x^2*y + 3*x^3*z + 2;
                                   2        3
                             p := x  y + 3 x  z + 2
\end{maple}

\noindent
and the formula $\sin (x+\pi/2) e^{-x}$ is input as

\begin{maple}
> sin(x+Pi/2)*exp(-x);
                                cos(x) exp(- x)
\end{maple}

\noindent
Notice that Maple simplified $sin(x+\pi/2)$ to $cos(x)$ for us.
Formulae in Maple are represented as {\em expression trees}
or {\em DAGs} (Directed Acyclic Graphs) in computer jargon.
When we program Maple functions to manipulate formulae, we are
basically manipulating expression trees.  The three basic
routines for examining these expression trees are
\verb+type, op+ and \verb+nops+ .
The \verb+type+ function

\begin{center}
\verb+type+( $f$, $t$ )
\end{center}

\noindent
returns the value \verb+true+ if the expression $f$ is of type $t$.
The basic types are {\tt string}, {\tt integer}, {\tt fraction},
{\tt float}, \verb:`+`:, \verb:`*`:, \verb:`^`:, and {\tt function}.
The \verb+whattype+ function is also useful for printing out the
type of an expression.
For example, our polynomial $p$ is a sum of 3 terms.  Thus

\begin{maple}
> type( p, integer );
                                     false

> whattype(p);
                                       +

> type( p, `+` );
                                      true
\end{maple}

\noindent
Note the use of the back quote character \verb+`+ here.
Back quotes are used for
strings in Maple which contain funny characters like \verb+/,.+ etc.
Back quotes are not the same as forward quotes (the apostrophe) \verb+'+
or double quotes \verb+"+.

\noindent
Two other numerical types are \verb+rational+, and \verb+numeric+.
The type \verb+rational+ refers to the rational numbers, i.e. integers
and fractions.  The type \verb+float+ refers to floating point
numbers, i.e. numbers with a decimal point in them.  The type
\verb+numeric+ refers to any of these kinds of numbers, i.e.
numbers of type \verb+rational+ or \verb+float+.
Maple users will have noticed that Maple distinguishes between
exact rational numbers and approximate numbers.
The presence of a decimal point is significant!  Consider

\begin{maple}
> 2/3;
                                      2/3

# That is a rational number, the following is a floating point number
> 2/3.0;
                                  .6666666667
\end{maple}

\vspace{2mm}
\sbox{
\parbox{6in}{
ATTENTION Engineers: many of you are used to using and writing Fortran or C
software.  In Fortran and C you always use decimal numbers.  I.e. you
have {\em learned} to always write 0.5 for 1/2.  In Maple, if you are typing
in a formula, you probably mean 1/2.
If you are typing in some data, you probably mean 0.5 .
1/2 is not the same as 0.5 in Maple.
}}
\vspace{2mm}

\noindent
Our example $p$ is a sum of terms.
How many terms does it have?  The \verb+nops+ function returns the {\bf n}umber
of {\bf op}erand{\bf s} of an expression.
For a sum, this means the number of terms in the sum.
For a product, it means the number of terms in the product.
Hence

\begin{maple}
> nops(p);
                                       3
\end{maple}

\noindent
The \verb+op+ function is used to extract one of the {\bf op}erands of an
expression.  It has the syntax

\begin{center}
\verb+op+( $i$, $f$ )
\end{center}

\noindent
meaning extract the $i$'th operand of the expression $f$ where $i$
must be in the range 1 to the \verb+nops+ of $f$.  In our example,
this means extract the $i$'th term of the sum $f$.

\begin{maple}
> op(1,p);
                                       2
                                      x  y

> op(2,p);
                                        3
                                     3 x  z

> op(3,p);
                                       2
\end{maple}

\noindent
The \verb+op+ function can also be used in the following way

\begin{center}
\verb+op+( $i..j$, $f$ ).
\end{center}

\noindent
This returns a sequence of operands of $f$ from $i$ to $j$.  For example

\begin{maple}
> op(1..3,p);
                                 2       3
                                x  y, 3 x  z, 2
\end{maple}

\noindent
A useful abbreviation is ${rm op}(f)$ which is equivalent to
${\rm op}(1..{\rm nops}(f),f)$ which means create a sequence
of all the operands of $f$.  What about the second term of the sum $p$?
It is a product of 3 factors

\begin{maple}
> type(op(2,p),`*`);        # It is a product

                                      true

> nops(op(2,p));            # It has 3 factors

                                       3

> op(1,op(2,p));            # Here is the first factor

                                       3

> op(op(2,p));              # Here is a sequence of all 3 factors

                                        3
                                    3, x , z
\end{maple}


\subsection*{Indexed names and Functions}

Maple has two kinds of variables or {\em names} as Maple calls them.
There are strings like \verb+x+, \verb+sin+, \verb+Pi+, which are of
type \verb+string+ and indexed names or subscripted variables
like $A_1, A_{i,j}, A_{i_j}$ which are of type \verb+indexed+.
These examples are input \verb+A[1]+, \verb+A[i,j]+, \verb+A[i][j]+ in Maple.
Most functions in Maple accept both kinds of variables.
The Maple type \verb+name+ means either a string or a subscript.
Example

\begin{maple}
> type(a,string);
                                      true

> type(a,name);
                                      true

> whattype(A[1]);
                                     indexed

> type(A[1],indexed);
                                      true

> type(A[1],name);
                                      true
\end{maple}

\noindent
If $f$ is an indexed name, then the \verb+nops+ function returns the
number of indices, and the \verb+op+ function returns the $i$'th index.
Also ${\rm op}(0,f)$ returns the name of the index.  Example

\begin{maple}
> nops(A[i,j]);
                                       2

> op(1,A[i,j]);
                                       i

> op(0,A[i,j]);
                                       A

> nops(A[i][j]);
                                       1

> op(1,A[i][j]);
                                       j

> op(0,A[i][j]);
                                      A[i]
\end{maple}


\noindent
Functions work very similarly to indexed names.
The syntax for a function call is

\begin{center}
$f( x_1, x_2, \ldots, )$
\end{center}

\noindent
where $f$ is the name of the function, and $x_1, x_2, \ldots$ are the arguments.
The \verb+nops+ function returns the number of arguments,
and the \verb+op+ function returns the $i$'th argument.
Also op$(0,f)$ returns the name of the function.
Example

\begin{maple}
> nops(f(x,y,z));
                                       3

> op(1..3,f(x,y,z));
                                    x, y, z

> op(0,f(x,y,z));
                                       f
\end{maple}

\noindent
We can now create and pick apart any Maple formula.
Exercise 6 asks you to write a Maple program that picks apart
a Maple formula.
We conclude here with some examples showing how this can be
done interactively.

\begin{maple}
> f := sin(x[1])^2*(1-cos(Pi*x[2]));

                                  2
                         sin(x[1])  (1 - cos(Pi x[2]))

> whattype(f);
                                       *

> op(1,f);
                                            2
                                   sin(x[1])

> op(1,op(1,f));
                                   sin(x[1])

> op(1,op(1,op(1,f)));
                                      x[1]

\end{maple}

\subsection{Statements: Assignment, Conditional, Loops}
The Maple syntax for the assignment, if, for and while statements
is taken from Algol 60.  The assignment statement looks like

\begin{quote}
{\em name} := {\em expr}
\end{quote}


\noindent
where {\em expr} is any expression and {\em name} is a variable name.
We want to mention here an evaluation problem that you will sooner
or later run into.  It arises because variables can be assigned expressions
which contain symbols as well as numbers.
Consider the assignment to $p$ earlier.

\begin{maple}
> p := x^2+4*x+4;

                                     2
                               p := x  + 4 x + 4
\end{maple}

\noindent
Here we assigned the name $p$ to a formula that explicitly contains
the symbol $x$.
What happens now if we use the same name $p$ on both sides of the
assignment statement?

\begin{maple}
> p := p+x;
                                     2
                               p := x  + 5 x + 4
\end{maple}

\noindent
Nothing spectacular.  But let us see how the evaluation rules operate here.
The name $p$ is assigned the result of evaluating and also simplifying the
right hand side $p+x$.  In evaluating $p+x$,
the value of $p$ is the polynomial $x^2 + 4 x + 4$ and the
value of $x$ is just $x$.
The result of evaluation is therefore the expression $(x^2 + 4 x + 4) + x$.
Next this is simplified to the final result $x^2 + 5 x + 4$.
Now comes the big question.  What happens if $p$ was not already assigned
a value?  Let's see what happens.  Using $q$ instead of $p$

\begin{maple}
> q := q+x;

                     Warning: Recursive definition of name

                                   q := q + x
\end{maple}

Well, you may say, obviously the user forgot to assign $q$ a value first.
Indeed that is probably the case, but let us continue and see what
happens to Maple.  Maple certainly allows $q$ to not have a value just
like $x$ in the previous example.

The above statement resulted in a warning from Maple.
The name $q$ is now assigned a formula which involves $q$.
What happens now if we try to evaluate $q$?   We would have to evaluate $q+x$.
But to evaluate $q+x$ we have to evaluate $q$ again.  Therein lies
an infinite evaluation loop.  On my system, if I try to evaluate
$q$, Maple dies with a rather spectacular crash

\begin{maple}
> q := q+x^2;
Segmentation fault
\end{maple}

\vspace{2mm}
\sbox{
\parbox{6in}{
ATTENTION: a recursive definition of a name like $q$ in the above
example will when evaluated typically result in crash as Maple runs
out of Stack space.  On some systems Maple will die.  On other systems,
Maple is able to stop just before it dies and report that there
is a problem.  If Maple ever crashes on you, there is a good chance
that this is the reason.  Maple will also die of course if there
is a true infinite loop in an algorithm.
}}
\vspace{2mm}

\noindent
To recover from such a warning, simply unassign the variable $q$, i.e.
do {\tt q := 'q';}.  Note that Maple does not detect all recursive
assignments because this would be too expensive to do in general.

In a conventional programming language, this problem of a recursive
definition cannot arise because all variables must have values.
If they haven't been assigned a value, this is really an error.
Depending on the language, variables either get a default value,
or they get whatever junk happens to be present in memory at the time,
or in a few cases, the language attempts to detect this and issue an error.

The conditional statement in Maple has the following syntax.

\begin{quote}
{\bf if} {\em expr} {\bf then} {\em statseq}

\hspace{1cm}
	[ {\bf elif} {\em expr} {\bf then} {\em statseq} ]*

\hspace{1cm}
	[ {\bf else} {\em statseq} ]

{\bf fi}
\end{quote}

\noindent
where {\em statseq} is a sequence of statements separated by semi-colons,
[] denotes an optional part, and * denotes a part which can
be repeated zero or more times.  A typical {\bf if} statement would be

\begin{quote}
\begin{verbatim}
if x < 0 then -1 elif x = 0 then 0 else 1 fi
\end{verbatim}
\end{quote}

\noindent
The {\bf for} statement has two variations.  The first one is

\begin{quote}
[ {\bf for} {\em name} ] [ {\bf from} {\em expr} ] [ {\bf by} {\em expr} ] [ {\bf to} {\em expr} ] [ {\bf while} {\em expr} ]

{\bf do} {\em statseq} {\bf od}
\end{quote}

\noindent
Thus each of the {\bf for}, {\bf from}, {\bf by}, {\bf to}, {\bf while} clauses can be omitted.
If omitted, the default values are a dummy variable, 1, 1, $\infty$ and true
respectively.  A typical example of a {\bf for} loop is

\begin{quote}
\begin{verbatim}
for i to 10 do print(i^2) od;
\end{verbatim}
\end{quote}

\noindent
If one omits the {\bf for}, {\bf from}, {\bf by} and {\bf to}
clauses we have a so called {\bf while} loop.

\begin{quote}
\begin{verbatim}
i := 10^10+1;
while not isprime(i) do i := i + 2 od;
print(i);
\end{verbatim}
\end{quote}

\noindent
Combining the {\bf for} and {\bf while} loops is sometimes nice, e.g. this example
of searching for the first prime $> 10^{10}$ could be done as

\begin{quote}
\begin{verbatim}
for i from 10^10+1 by 2 while not isprime(i) do od;
print(i);
\end{verbatim}
\end{quote}

\noindent
The above example illustrates that the value of the loop index can be
accessed after the loop has terminated.

The second variation on the {\bf for} loop is the so called {\bf for-in} loop.
It is really just a short hand for this loop which occurs very often

\begin{quote}
\begin{verbatim}
for i to nops(s) do f(op(i,s)) od;
\end{verbatim}
\end{quote}

\noindent
where $s$ may be a sum, but in general, any Maple expression or
data structure such as a list or set -- see next section.
One can code this as a {\bf for-in} loop as follows

\begin{quote}
\begin{verbatim}
for i in s do f(i) od;
\end{verbatim}
\end{quote}

\noindent
The syntax for the {\bf for-in} loop is

\begin{quote}
[ {\bf for} {\em name} ] [ {\bf in} {\em expr} ] [ {\bf while} {\em expr} ]

{\bf do} {\em statseq} {\bf od}
\end{quote}

\section{Data Structures}

More complicated programs involve manipulating and storing data.
How we represent out data affects the algorithms that we write, and how fast
our programs will run.  Maple has a good set of data structures.
The ones we will look at here are sequences, lists (or vectors), sets,
tables (or hash tables), and arrays.  Maple does not have records or linked lists.
We'll say how these can be implemented in Maple at the end of this section.

\subsection{Sequences}

A sequence is a sequence of expressions separated by commas.  For example

\begin{maple}
> s := 1,4,9,16,25;
                              s := 1, 4, 9, 16, 25

> t := sin,cos,tan;
                               t := sin, cos, tan
\end{maple}

\noindent
A sequence of sequences simplifies into one sequence, that is,
sequences are {\em associative}.  For example

\begin{maple}
> s := 1,(4,9,16),25;
                              s := 1, 4, 9, 16, 25

> s,s;
                        1, 4, 9, 16, 25, 1, 4, 9, 16, 25
\end{maple}

\noindent
The special symbol NULL is used for the empty sequence.
Sequences are used for many purposes.  The next section shows how
lists and sets are constructed from sequences.  Here we note that
function calls are really constructed from sequences.
For example, the {\tt min} and {\tt max} functions in Maple take an
arbitrary number of values as arguments, i.e. a sequence of arguments

\begin{maple}
> max(s);

                                       25

> min(s,0,s);

                                       0
\end{maple}

\vspace{2mm}
\sbox{
\parbox{6in}{
ATTENTION: the {\tt op} and {\tt nops} functions cannot be applied to sequences.
This is because the sequence itself becomes the arguments to a function.
Thus the calls {\tt op(s);} and {\tt nops(s);} are equivalent to
{\tt op(1,4,9,16,25);} and {\tt nops(1,4,9,16,25)} respectively,
which result in an error.
Put a sequence into a list first if you want to use {\tt op} or {\tt nops}.
}}
\vspace{2mm}

The \verb+seq+ function is extremely useful for creating sequences.
It comes in two flavours, corresponding to the two kinds of {\bf for} loops.
The first one is

\begin{center}
\verb+seq+( $f(i)$, $i=m..n$ ).
\end{center}

\noindent
The way {\tt seq} works is just as if you had programmed the following loop.

\begin{quote}
\begin{verbatim}
s := NULL;
for i from m by 1 to n do s := s, f(i) od;
\end{verbatim}
\end{quote}

\noindent
For example

\begin{maple}
> seq( i^2, i=1..5 );
                                1, 4, 9, 16, 25

> s := NULL; for i from 1 to 5 do s := s, i^2 od;

                                     s :=

                                     s := 1

                                   s := 1, 4

                                  s := 1, 4, 9

                                s := 1, 4, 9, 16

                              s := 1, 4, 9, 16, 25
\end{maple}

\noindent
Notice that \verb+seq+ is more efficient than the {\bf for} loop because
it does not create the intermediate sequences.
The other flavour of \verb+seq+ is

\begin{center}
\verb+seq+( $f(i)$, $i=a$ ).
\end{center}

\noindent
This is equivalent to

\begin{quote}
\begin{verbatim}
seq( f(op(i,a)), i=1..nops(a) )
\end{verbatim}
\end{quote}

\noindent
Here are a couple of interesting examples.  The \verb+coeff+ function
computes the coefficient of the term of degree $i$ of a polynomial in $x$.
The \verb+D+ function in Maple is the derivative operator.

\begin{maple}
> a := 3*x^3+y*x-11;
                                      3
                              a := 3 x  + y x - 11

> seq( coeff(a,x,i), i=0..degree(a,x) );

                                  -11, y, 0, 3

> seq( D(f), f=[sin,cos,tan,exp,ln] );

                                         2
                      cos, - sin, 1 + tan , exp, a -> 1/a
\end{maple}

\subsection{Lists and Sets}

Lists, sets and functions are constructed from sequences.
A list is a data structure for collecting objects together.
Square brackets are used to create lists, for example

\begin{maple}
> l := [x,1,1-z,x];
                             l := [x, 1, 1 - z, x]

> whattype(l);

                                      list
\end{maple}

\noindent
The empty list is denoted by [].
Sets can also be used to collect objects together.  The difference
between lists and sets is that duplicates are removed from sets.
Squiggly brackets are used for sets, for example

\begin{maple}
> s := {x,1,1-z,x};
                               s := {1, 1 - z, x}

> whattype(s);
                                      set
\end{maple}

\noindent
The empty set is denoted by \{\}.
The \verb+nops+ function returns the number of elements of a list or set
and the \verb+op+ function extracts the $i$'th element.
You can also use the subscript notation to access the $i$'th element of
a sequence, list or set.  For example,

\begin{maple}
> op(1,s);
                                       1

> s[1];
                                       1

> op(1..3,s);
                                  1, 1 - z, x

> s[1..3];

                                  1, 1 - z, x
\end{maple}

\noindent
Here is a loop that prints true if a list or set contains the element $x$,
and false otherwise.

\begin{quote}
\begin{verbatim}
for i to nops(s) while s[i] <> x do od;
if i > nops(s) then print(false) else print(true) fi;
\end{verbatim}
\end{quote}

\noindent
The built in \verb+member+ function does this for you.
\verb+member+($x$, $s$)
returns true if the element $x$ is in the list or set $s$.
You can append a new element to a list l by doing

\begin{quote}
\begin{verbatim}
l := [op(l),x];
\end{verbatim}
\end{quote}

\noindent
You can remove the $i$'th element of a list by doing

\begin{quote}
\begin{verbatim}
l := [ l[1..i-1], l[i+1..nops(l)] ];
\end{verbatim}
\end{quote}

\noindent
Or, better, use the \verb+subsop+ function

\begin{quote}
\begin{verbatim}
l := subsop(i=NULL,l);
\end{verbatim}
\end{quote}

\noindent
The \verb+subsop+ function note can be used for any type of expression.
The set operators are \verb+union+, \verb+intersect+ and
\verb+minus+, e.g.

\begin{maple}
> t := {u,x,z};
                                 t := {x, z, u}

> s union t;
                                  {x, z, y, u}
\end{maple}

The Maple user will probably have noticed that Maple orders
the elements of sets in a strange order which seems to be unpredictable.
The algorithms that Maple uses to remove duplicates, and to do set union,
intersection and difference all work by first sorting the elements of the
input sets.  Maple sorts by {\em machine address}, i.e. in the order
that the elements occur in computer memory.  Since this depends on where
they are put in memory, this is why the order is strange and unpredictable.
The reason Maple sets are implemented in this way is to make the
set operations very fast.

\subsection{Tables}

Tables or hash tables are extremely useful for writing efficient programs.
A table is a one-to-many relation between two discrete sets of data.
For example, here is a table of colour translations for English
to French and German.

\begin{maple}
> COLOUR[red] := rouge,rot;

                           COLOURS[red] := rouge, rot

> COLOUR[blue] := bleu,blau;

                          COLOURS[blue] := bleu, blau

> COLOUR[yellow] := jaune,gelb;

                         COLOURS[yellow] := jaune, gelb
\end{maple}

\noindent
The domain of the \verb+COLOUR+ table is the name of the colour in English.
The range of the table is a sequence of the names of the colours in
French and German.  In general, the domain and range values
can be sequences of zero or more values.
The domain values in the table are called the {\em keys} or {\em indices}.
The Maple \verb+indices+ function returns a sequence of them.
The range values in the table are called the {\em values} or {\em entries}.
The Maple \verb+entries+ function returns a sequence of them.
For example

\begin{maple}
> indices(COLOUR);

                            [red], [yellow], [blue]

> entries(COLOUR);

                   [rouge, rot], [jaune, gelb], [bleu, blau]
\end{maple}

\noindent
Note, the order that the {\tt indices} and {\tt entries} functions return
the table indices and entries is not necessarily the same order as the
order in which they were entered into the table.  This is because Maple makes
use of {\em hashing} to make table searching very fast and as a consequence,
the order in which the entries were made is lost.  However, there is a one to
one correspondence between the output of the {\tt indices}
and the {\tt entries} functions.

What can one do with a table?
Given a table {\em key} or {\em index},
we can look up the corresponding entry very quickly.
That is, the operation

\begin{maple}
> COLOUR[red];

                                   rouge, rot
\end{maple}

\noindent
returns the French and German names for the colour red quickly.  How quickly?
Approximately in constant time, no matter how large the table is.
Even if our table contains 1000 different colours,
the time to search the table will be very fast.
That is the basic use of a table.
Information is entered into a table by assignment, and table look up
is done by table subscript.  What else can one do with a table?
We can test if an entry is in the table using the {\tt assigned} function
and we can remove entries from a table by unassignment.
For example

\begin{maple}
> assigned(COLOUR[blue]);
                                      true

> COLOUR[blue] := 'COLOUR[blue]';

                          COLOUR[blue] := COLOUR[blue]

> assigned(COLOUR[blue]);
                                     false

> print(COLOUR);
                           table([
                               red = (rouge, rot)
                               yellow = (jaune, gelb)
                           ])

\end{maple}

\subsection{Arrays}

A one-dimensional array is created using the \verb+array+ command

\begin{center}
array( $m$..$n$ );
\end{center}

\noindent
This creates an array with indices $m, m+1, \ldots, n$.
Often $m$ is 1.  Entries can be inserted into the array by assignment
as for tables, e.g.

\begin{quote}
\begin{verbatim}
v := array(1..n);
v[1] := a;
for i from 2 to n do v[i] := a*v[i-1] mod n od;
\end{verbatim}
\end{quote}
 
One-dimensional arrays are like tables with one index which is restricted to
be in a fixed integer range.  Arrays are more efficient to access than tables
and range checking on the indices is done.
Here is an example of a one-dimensional array used to sort
a sequence of numbers.  For instance, suppose we are given an array
$v$ = \verb+array(1..n)+ as above.
The code presented here
sorts $v$ into ascending order using the {\em bubblesort} algorithm:

\begin{quote}
\begin{verbatim}
for i to n-1 do
    for j from i+1 to n do
        if v[i] > v[j] then temp := v[i]; v[i] := v[j]; v[j] := temp fi
    od
od;
\end{verbatim}
\end{quote}

Another application of one-dimensional arrays is as an intermediate
data structure.  For example, suppose we represent a polynomial
$a(x) = \sum _{i=0} ^ m$ as a list $[a_0, \ldots, a_m]$ of coefficients.
Suppose we are given also a polynomial $b$ of degree $n$.  We can
use an array $c$ to compute the product $a \times b$ as follows

\begin{quote}
\begin{verbatim}
m := nops(a)-1; # degree of a
n := nops(b)-1; # degree of b
c := array(0..m+n); # allocate storage for the product
for i from 0 to m+n do c[i] := 0 od;
for i from 0 to m do
    for j from 0 to n do
        c[i+j] := c[i+j] + a[i+1]*b[j+1]
    od
od:
[seq(c[i],i=0..n+m)]; # put the product in a list
\end{verbatim}
\end{quote}

Two-dimensional arrays, and higher dimensional arrays work similarly.
A two-dimensional array is created by the command

\begin{center}
array( $c..d$, $m$..$n$ );
\end{center}

\noindent
As an example, suppose we are given the vector $v$ of
symmetric polynomials in $3$ variables $x_1, x_2, x_3$.

\begin{maple}
v := array(1..4):
v[1] := 1:
v[2] := x[1] + x[2] + x[3]:
v[3] := x[1]*x[2] + x[1]*x[3] + x[2]*x[3]:
v[4] := x[1]*x[2]*x[3]:
\end{maple}

\noindent
Let us construct a two-dimensional array $M$, where
$M_{i,j}$ is the derivative of $v_i$ wrt $x_j$.

\begin{maple}
> M := array(1..4,1..3):
> for i to 4 do for j to 3 do M[i,j] := diff(v[i],x[j]) od od:
> M;
                                       M

> eval(M);
                   [      0            0            0      ]
                   [                                       ]
                   [      1            1            1      ]
                   [                                       ]
                   [ x[2] + x[3]  x[1] + x[3]  x[1] + x[2] ]
                   [                                       ]
                   [  x[2] x[3]    x[1] x[3]    x[1] x[2]  ]
\end{maple}

\vspace{2mm}
\sbox{
\parbox{6in}{
ATTENTION: Notice that the value of \verb+M+ is just the name of the
array \verb+M+.  Evaluation rules for arrays and tables are special.
The reason is so that one can have arrays with unassigned entries.
We will not give details about this here.  For the moment,
whenever you want to print an array or table or return an array or
table from a procedure, use the \verb+eval+ function. }
}
\vspace{2mm}

For further information on arrays see \verb+?array+ .
Note also that one-dimensional arrays indexed from 1 are used to
represent vectors and two-dimensional arrays indexed from 1 are
used to represent matrices.  See the section on Matrices and Vectors.

\subsection{Records}

Maple doesn't explicitly have a record data structure like Pascal's
{\tt record} or C's {\tt struct}.
By a record data structure we mean a data structure for keeping together
a heterogeneous collection of objects, i.e. a set of objects not
necessarily of the same type.

An example of a place where you would like a record data structure would be
in choosing a data structure to represent a quarternion.
A quarternion is a number of the form $a + b i + c j + d k$ where $a, b, c, d$
are real numbers.  To represent a quarternion, we need to store only the
four quantities $a, b, c$ and $d$.
Another example would be a data structure to represent 
the factorization of a polynomial in ${\bf Q}[x]$.
The factorization of $a(x)$ looks like

\[
a(x) ~ = ~ u \times f_1 ^ {e_1} \times \ldots \times f_n ^ {e_n}
\]

\noindent
where each of the factors $f_i \in {\bf Q}[x]$ is monic and irreducible.
We need to store the factors $f_i$ and the exponents $e_i$ and the
unit $u \in {\bf Q}$.

There are several possibilities for representing records in Maple.
The simplest, and most obvious, is to use a list.  I.e. we would
represent the quarternion $a + b i + c j + d k$ as the list $[a,b,c,d]$,
and the factorization of $a(x)$ as the list $[u,f]$ where $f$ is a
list of lists of the form $[f_i,e_i]$.
We use subscripts to refer to a component of the data structure,
e.g. the unit part of a factorization would be given by {\tt a[1]}.
You can use the {\em macro} facility to define an identifier
to be equal to a constant if you prefer to use a symbol to reference
a component as shown in the following example.

\begin{maple}
> a := [-1/2,[[x+1,2],[x-1,1]]];

                     a := [-1/2, [[x + 1, 2], [x - 1, 1]]]

> macro(unit=1,factors=2,base=1,exponent=2);
> a[unit];
                                      -1/2

> a[factors][1][base];
                                     x + 1

> a[unit]*a[factors][1][base]^a[factors][1][exponent]*a[2][2][1]^a[2][2][2];

                                          2
                             - 1/2 (x + 1)  (x - 1)
\end{maple}

A second possibility for representing records in Maple is to use a
function call.  I.e. we could represent $a + b i + c j + d k$
as {\tt QUARTERNION(i,j,k,l)}.
An advantage of this representation is that you we can tell Maple how
to do various operations on functions.
We will go into details later.
Here we shall only mention that you can define
how to {\em pretty print} a function.
The example below will show what we mean

\begin{maple}
> QUARTERNION(2,3,0,1);

                             QUARTERNION(2, 3, 0, 1)

> `print/QUARTERNION` := proc(a,b,c,d)  a + b*'i' + c*'j' + d*'k' end:
> QUARTERNION(2,3,0,1);

                                  2 + 3 i + k

\end{maple}

\noindent
Here we have defined a printing procedure or subroutine.  This routine
is called once for each different {\tt QUARTERNION} function in a result from Maple
prior to displaying the result.  Note the use of quotes in the
procedure, because we want to see the identifiers $i, j$, and $k$ in the
output, and not the value of the variables $i,j,k$ which we might be using
for something else.

A third possibility for representing a record is to think of it as
a multivariate polynomial in the field names, and to store the values
in the coefficients.  This is quite useful when the fields are
numerical and you wish to be able to do arithmetic on the fields.
For example, we could represent a quarternion as a polynomial in the
variables $i,j,k$ as in the output representation above.  I.e.

\begin{maple}
> z1 := 2 + 3*i + k;

                               z1 := 2 + 3 i + k

> z2 := 2 - 3*i + 2*j + 2*k;

                           z2 := 2 - 3 i + 2 j + 2 k

> coeff(z1,i); # the coefficient in i

                                       3

> z1 + z2;

                                 4 + 3 k + 2 j
\end{maple}

\noindent
Although this looks nice, we don't recommend using the names $i, j$, or $k$
because you will use them for {\bf for} loop variables!

\subsection{Linked Lists}

Maple lists are not linked lists.  Maple lists are really arrays of
pointers to their entries.  Maple lists differ from Maple arrays in
that they are read only, i.e. you can't assign to a component of a Maple list.
Linked lists are recursive data structures.
The difference can be seen by studying a simple example.
Consider representing a polynomial $a(x) = \sum _ {i=0} ^ n a_i x^i$
in Maple.  One possibility would be to store the coefficients
$a_0, \ldots a_n$ in a Maple list.
For instance, we could represent the polynomial $p = x^4+3 x^2+2 x+11$ as follows

\begin{quote}
\begin{verbatim}
[11, 2, 3, 0, 1]
\end{verbatim}
\end{quote}

\noindent
The degree of the polynomial is the number of entries in the list minus one.
Alternatively we could use a linked list

\begin{quote}
\begin{verbatim}
[ 1, [ 0, [ 3, [ 2, [ 11, NIL ]]]]]
\end{verbatim}
\end{quote}

\noindent
We see that the linked list is a recursive data structure.
It is either list of two values, traditionally called the {\em CAR}
and {\em CDR} -- terminology from the Lisp programming language,
or it is the special value NIL signifying the empty linked list.
The first field contains a data value, in our case,
a coefficient.  The second field is a pointer to another linked list
which contains the rest of the data values.
Note that in order to compute the degree of a polynomial represented
by a linked list, we have to compute the depth of the linked list.
If $p$ is a linked list, you could do this using the loop

\begin{quote}
\begin{verbatim}
for n from 0 while p <> NIL do p := p[2] od;
\end{verbatim}
\end{quote}

Now, why would we represent a polynomial using a linked list instead
of a Maple list?  The principle reason is that we can put a new
entry onto the front of a linked list in constant time instead of linear time.
For example, suppose we wanted to add the term $5 x^5$ to our polynomial $p$.
In the Maple list representation we must create a
new list with 6 entries by doing

\begin{maple}
[op(p),5];
\end{maple}

\noindent
This requires at least 6 words of storage for the new Maple list.
Don't be fooled.  The {\tt op} call is short for
{\tt op(1..nops(p),p)} which creates a sequence of all the entries
in the Maple list $p$.  This takes constant time and uses no storage
but now, we have sequence {\tt (11, 2, 3, 0, 1), 5} which
results in the new larger sequence {\tt 11, 2, 3, 0, 1, 5} being created.
This is where 6 words of storage get allocated.
In general, adding $a_{n+1} x^{n+1}$ to a polynomial of degree $n$
takes $O(n)$ time and storage.  But what about the linked list?
To add the term $5 x^5$ we do

\begin{maple}
[5,p];
\end{maple}

\noindent
which only needs to create a new Maple list of length two,
hence constant storage.
The time will also be linear time if $p$ is a global variable,
but constant time if $p$ is a local variable or parameter inside a procedure.  
This evaluation detail will be explained in the section on procedures.
But for now assume that the running time for this operation
is also $O(1)$ for linked lists.

\vspace{2mm}
\sbox{
\parbox{6in}{
ATTENTION Lisp Programmers: You can't change the components of
a Maple list.  There can be no equivalent of {\tt REPLACA} and
{\tt REPLACDR}.  Maple lists, sets, sequences, and functions are
read only data structures.  Only arrays and tables can be changed.
}}

\section{Maple Procedures}

\subsection{Parameters, Local Variables, RETURN, ERROR}

A Maple procedure has the following syntax

\begin{quote}
{\bf proc} ( {\em nameseq } )

\hspace{1cm}
[ {\bf local } { \em nameseq } ; ]

\hspace{1cm}
[ {\bf global} { \em nameseq } ; ]

\hspace{1cm}
[ {\bf options} { \em nameseq } ; ]

\hspace{1cm}
{\em statseq }

{\bf end}
\end{quote}

\noindent
where {\em nameseq} is a sequence of symbols separated by commas,
and {\em statseq} is a sequence of statements separated by semicolons.
Here is a simple procedure which, given $x,y$, computes $x^2+y^2$.

\begin{quote}
\begin{verbatim}
proc(x,y) x^2 + y^2 end
\end{verbatim}
\end{quote}

This procedure has two parameters $x$ and $y$.
It has no local variables, no global variables, no options,
and only one statement.
The value returned by the procedure is $x^2 + y^2$.
In general the value returned by a procedure is the last value computed
unless there is an explicit return statement.  An example of a procedure with
an explicit return statement is the following \verb+MEMBER+ procedure.
\verb+MEMBER(x,L)+ returns \verb+true+ if $x$ is in the list $L$,
false otherwise.

\begin{quote}
\begin{verbatim}
MEMBER := proc(x,a) local v;
    for v in L do if v = x then RETURN(true) fi od;
    false
end;
\end{verbatim}
\end{quote}

\noindent
The \verb+MEMBER+ procedure has a local variable $v$ so that it does
not interfere with the global user variable $v$.
Variables that appear in a procedure which are neither parameters nor
locals and are not explicitly declared global are in this new version of
Maple given a default declaration local if either the following is true
\begin{itemize}
\item
A variable or suscripted variable on the left hand side of an assignment
\item
A variable or suscripted variable in a for or seq loop index
\end{itemize}
Part of the reason for this is to catch likely errors.  For example,

\begin{maple}
> proc(x,n) s := 1; for i to n do s := s+x^i od; s end;

  proc(x,n) local s,i; s := 1; for i to n do  s := s+x^i od; s end
\end{maple}
Notice that Maple has declared s and i to be local.  If you want variables
to really be global, you should declare them explicitly using the
{\tt global} statement.

The \verb+ERROR+ function can be used to generate an error message
from within a procedure.
For example, the \verb+MEMBER+ routine should check that the argument
really is a list and issue an appropriate error message otherwise.

\begin{quote}
\begin{verbatim}
MEMBER := proc(x,L) local v;
    if not type(L,list) then ERROR(`2nd argument must be a list`) fi;
    for v in L do if v = x then RETURN(true) fi od;
    false
end;
\end{verbatim}
\end{quote}


\subsection{Tracing Procedure Execution: {\tt printlevel}}

Here is a Maple procedure which computes the greatest common
divisor of two non-negative integers using the Euclidean algorithm.
For example, the greatest common divisor of the integers 21 and 15 is 3 because
3 is the largest integer that divides both 21 and 15.  By the way, the
Euclidean algorithm is one of the oldest known algorithms in Mathematics.
It dates back to around 300 BC.

\begin{quote}
\begin{verbatim}
GCD := proc(a,b) local c,d,r;
    c := a;
    d := b;
    while d <> 0 do r := irem(c,d); c := d; d := r od;
    c
end;
\end{verbatim}
\end{quote}

The \verb+irem+ function computes the {\bf i}nteger {\bf rem}ainder
of two integers.  How does the GCD routine really work?
The simplest tool for looking at the execution of a procedure is
the {\tt printlevel} facility.  The \verb+printlevel+ variable is a global
variable that is initially assigned 1.
If you set it to a higher value, a trace of all assignments, and procedure
entries and exits is printed.
Let's see what happens when we compute {\tt GCD(21,15);}

\begin{maple}
> printlevel := 100;
> GCD(21,15);
{--> enter GCD, args = 21, 15
                                    c := 21

                                    d := 15

                                     r := 6

                                    c := 15

                                     d := 6

                                     r := 3

                                     c := 6

                                     d := 3

                                     r := 0

                                     c := 3

                                     d := 0

                                       3

<-- exit GCD (now at top level) = 3}
                                       3

\end{maple}

We see that the input arguments to the GCD procedure are displayed
together with the value returned.  The execution of each assignment statement
is also displayed.

An interesting point about our Maple GCD procedure is that
this routine works for integers of any size
because Maple uses arbitrary precision integer arithmetic.
For example, here is the greatest common divisor between 100! and $2^{100}$

\begin{maple}
> GCD(100!,2^100);
                         158456325028528675187087900672
\end{maple}

\noindent
Our \verb+GCD+ procedure could also have been been written recursively
and should also include type checking on the input parameters
in this way

\begin{quote}
\begin{verbatim}
GCD := proc(a:integer,b:integer)
    if b = 0 then a else GCD(b,irem(a,b)) fi
end;
\end{verbatim}
\end{quote}

The recursive version is simpler and easier to understand.
Let us see a trace of the recursive version to see how the {\rm printlevel}
facility can show the flow of computation.

\begin{maple}
> GCD(15,21);
{--> enter GCD, args = 15, 21
{--> enter GCD, args = 21, 15
{--> enter GCD, args = 15, 6
{--> enter GCD, args = 6, 3
{--> enter GCD, args = 3, 0
                                       3

<-- exit GCD (now in GCD) = 3}
                                       3

<-- exit GCD (now in GCD) = 3}
                                       3

<-- exit GCD (now in GCD) = 3}
                                       3

<-- exit GCD (now in GCD) = 3}

<-- exit GCD (now at top level) = 3}
                                       3

\end{maple}

\subsection{Arrow Operators}
For procedures which compute only a formula, there is an alternative
syntax called the arrow syntax.  This mimics the syntax for functions
often used in algebra.  For functions of one parameter the syntax is

\begin{center}
{\em symbol} \verb+->+ [{\bf local} {\em nameseq};] {\em expr} .
\end{center}

\noindent
For 0 or more parameters, parameters are put in parentheses i.e.

\begin{center}
( {\em nameseq} ) \verb+->+ [{\bf local} {\em nameseq};] {\em expr}
\end{center}

\noindent
The example which computes $x^2+y^2$ can be written more succinctly as

\begin{quote}
\begin{verbatim}
(x,y) -> x^2+y^2;
\end{verbatim}
\end{quote}

In Maple V Release 2, the syntax has been extended to allow the body
of the procedure to be an {\bf if} statement.  So you can define piecewise
functions.  For example
\begin{quote}
\begin{verbatim}
x -> if x<0 then 0 elif x<1 then x elif x<2 then 2-x else 0 fi;
\end{verbatim}
\end{quote}

\subsection{Scope Rules: Parameters, Locals, Globals}

Maple supports nested procedures.  For example, you can write

\begin{quote}
\begin{verbatim}
f1 := proc(x) local g; g := x -> x+1; x*g(x) end;
\end{verbatim}
\end{quote}

Procedure $f_1$ has a local variable $g$ which is a procedure.
$f_1$ computes $x*(x+1)$.  However, nested parameters and local variables
do not use nested scoping rules.  E.g. the above procedure is not
equivalent to this one

\begin{quote}
\begin{verbatim}
f2 := proc(x) local g; g := () -> x+1; x*g() end;
\end{verbatim}
\end{quote}

\noindent
because the reference to $x$ in the $g$ procedure does not refer to
the parameter $x$ in $f_2$.  It refers to the global variable
$x$.  Consider these examples 

\begin{maple}
> f1(a);

                                   a (a + 1)

> f2(a);

                                   a (x + 1)

> x := 7;

                                     x := 7

> f2(a);

                                      8 a
\end{maple}

One similarly cannot refer to local variables in outer scopes.

\vspace{2mm}
\sbox{
\parbox{6in}{
ATTENTION Computer Scientists:  Although Maple supports nested procedures,
which can be returned as function values, it does not support nested lexical
scopes, so you cannot return closures directly.
For example, {\tt f -> x -> f(x+1)} does NOT define a shift operator.
This is being looked at for a future version of Maple.
}}

\subsection{Evaluation Rules: Actual and Formal Parameters, Locals, Globals}

Consider the function call {\tt f(x1,x2,...,xn)}.  The execution
of this function call proceeds as follows.  The function name {\em f}
is evaluated.  Next the arguments {\em x1, x2, ..., xn} are evaluated
from left to right.  Then if {\em f} evaluated to a procedure, the
procedure is executed on the evaluated arguments.  There are only 6
exceptions to this, including {\tt eval}, {\tt assigned}, and {\tt seq}.
Now, what about the evaluation of variables inside procedures?
Consider the \verb+procedure+

\begin{quote}
\begin{verbatim}
f := proc() local p; global x; p := x^2+4*x+4; x := 5; p end
\end{verbatim}
\end{quote}

\noindent
Here $p$ is a local variable and there are no parameters.
And x is a global variable.
What are the evaluation rules for parameters, locals, and global variables?
In the introduction we considered the problem

\begin{maple}
> p := x^2+4*x+4;
                                     2
                               p := x  + 4 x + 4

> x := 3;
                                     x := 3

> p;
                                       25
\end{maple}

Here $p$ and $x$ are global variables.  Global variables are evaluated fully,
i.e. recursively, hence the result is $25$.
What if $p$ is a local variable as in our procedure \verb+f+ above?
I.e. what is the result of executing

\begin{quote}
\begin{verbatim}
f();
\end{verbatim}
\end{quote}

And what if $p$ is a parameter as in the following procedure

\begin{quote}
\begin{verbatim}
> x := 'x'; # Make x a symbol first

                                     x := x

> g := proc(p) global x; x := 5; p end:
> g(x^2+4*x+4);

                                   2
                                  x  + 4 x + 4
\end{verbatim}
\end{quote}

For reasons of efficiency and desirability, the Maple designers have
decided that local variables and parameters evaluate one level,
i.e. the value of $p$ in the above two examples is the
polynomial $x^2 + 4 x + 4$ not the value $49$.
Full evaluation only occurs for global variables.
The \verb+eval+ function can be used to get full evaluation for local
variables and parameters, and one level evaluation of global variables
should you ever need it.  For example

\begin{quote}
\begin{verbatim}
> x := 'x'; # Make x a symbol first

                                     x := x

> g := proc(p) global x; x := 5; eval(p) end:
> g(x^2+4*x+4);

                                       49

\end{verbatim}
\end{quote}


\subsection{Recurrence Equations and Option Remember}

The Fibonacci numbers $F_n$ are defined by the linear recurrence
$F_0 = 0$, $F_1 = 1$, and $F_n = F_{n-1} + F_{n-2}$.
This can be coded directly by

\begin{quote}
\begin{verbatim}
F := proc(n) 
    if n = 0 then 0 elif n = 1 then 1 else F(n-1)+F(n-2) fi
end;
\end{verbatim}
\end{quote}

\noindent
Here are the first few Fibonacci numbers

\begin{maple}
> seq( F(i), i=0..10 );

                      0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55
\end{maple}

\noindent
However this is a not an efficient way to compute the Fibonacci numbers.
In fact, you will never be able to compute
\verb+F(100)+ using this procedure even on the fastest computer.
If you count the number of calls to the F procedure, you will see that
it is called repeatedly on the same arguments.
It is clear that one should remember the previous two values when computing
the next value.  This could be done in a loop in the following way

\begin{quote}
\begin{verbatim}
F := proc(n) local fnm1,fnm2,f;
    if n = 0 then RETURN(0) fi;
    fnm2 := 0;
    fnm1 := 1;
    for i to n do f := fnm1 + fnm2; fnm2 := fnm1; fnm1 := f od;
    fnm1
end;
\end{verbatim}
\end{quote}

\noindent
Another way to code this is to use {\em option remember}.
This option is used to store values as they are computed
so that they can be used when they are needed.
Consider

\begin{quote}
\begin{verbatim}
F := proc(n) option remember;
    if n = 0 then 0 elif n = 1 then 1 else F(n-1)+F(n-2) fi
end;
\end{verbatim}
\end{quote}

\noindent
This program computes $F(100)$ quite quickly.
Each Maple procedure has an associated {\em remember} table.
The table index is the arguments and the table entry is the function value.
When $F$ is called with $n$, Maple first looks up $F$'s {\em remember}
table to see if $F(n)$ has already been computed.
If it has, it returns the result from $F$'s remember table.
Otherwise, it executes the code for the procedure $F$, and
automatically stores the pair $n, F(n)$ in $F$'s remember table.

We also illustrate the possibility of
explicitly saving values in a remember table by using
the so called {\em functional assignment}.
This is more flexible than the remember option because it allows
one to save only selected values in the remember table.

\begin{quote}
\begin{maple}
F := proc(n) F(n) := F(n-1)+F(n-2) end;
F(0) := 0;
F(1) := 1;
\end{maple}
\end{quote}

\subsection{Types and Map}
The \verb+type+ function can be used to code a routine that does
different things depending on the type of the input.
For example, the \verb+DIFF+ routine below differentiates expressions
which are polynomials in the given variable $x$.

\begin{quote}
\begin{verbatim}
DIFF := proc(a:algebraic,x:name) local u,v;
    if type(a,numeric) then 0
    elif type(a,name) then if a = x then 1 else 0 fi
    elif type(a,`+`) then map( DIFF, a, x )
    elif type(a,`*`) then u := op(1,a); v := a/u; DIFF(u,x)*v + DIFF(v,x)*u
    elif type(a,anything^integer) then
        u := op(1,a); v := op(2,a); v*DIFF(u,x)*u^(v-1)
    else ERROR(`don't know how to differentiate`,a)
    fi
end;
\end{verbatim}
\end{quote}

\noindent
Types are used in the {\tt DIFF} procedure for two different purposes.
The first usage is for type checking.  The inputs must be an algebraic
expression or a formula, and a name for the differentiation variable.
The second usage is to examine the type of the input -- is it a sum or product --
and decide what differentiation rule is to be used.
The \verb+DIFF+ example shows the use of the \verb+map+ function,
a very useful function, which we will now explain.  It has the following syntax

\begin{center}
map( $f$, $a$, $x_1$, ..., $x_n$ )
\end{center}


\noindent
The meaning is to apply the function $f$ to the operands of the expression $a$
passing the additional arguments $x_1$, ..., $x_n$ to $f$.
In our {\tt DIFF} procedure, there is one additional argument {\tt x}.
Often, there are no additional arguments.
Formally, this is equivalent to computing the sequence

\begin{quote}
\begin{verbatim}
seq( f( op(1,a), x1, ..., xn ), i=1..nops(a) );
\end{verbatim}
\end{quote}

\noindent
and then creating from this sequence a value of the same type as the type of $f$.
Here are some examples

\begin{maple}
> p := x^3+2*x+1;
                                     3
                               p := x  + 2 x + 1

> map( F, p );
                                3
                             F(x ) + F(2 x) + F(1)

> map( x -> x^2, p );
                                  6      2
                                 x  + 4 x  + 1

> map( degree, p, x );
                                       4
\end{maple}

The {\tt DIFF} function also shows the use of a {\em structured} type.
The types {\tt anything}, {\tt name}, \verb:`+`: and \verb+`*`+ are
simple types.  The type \verb+anything^integer+ is a {\em structured} type.
It means that the value must be a power, and the base can be anything,
i.e. any type, but the exponent must be an integer.  It is equivalent to writing

\begin{quote}
\begin{verbatim}
if type(a,`^`) and type(op(2,a),integer) then
\end{verbatim}
\end{quote}

\noindent
Structured types allow you to replace long type tests with concise tests.
Let us illustrate another common case.
%The type {\tt name=algebraic..algebraic} is the type of the second argument
%for definite integration, and definite summation, as in
%
%\begin{maple}
%> int( exp(-2*t)*ln(t), t=0..infinity );
%
%                            - 1/2 ln(2) - 1/2 gamma
%
%> sum( 1/(i^2-1), i=2..infinity );
%
%                                      3/4
%\end{maple}
%
%\noindent
%Of course, the integration and summation commands also allow indefinite
%integration and summation, i.e.
%
%
%\begin{maple}
%> int( exp(-2*x)*ln(x), x );
%
%                     - 1/2 exp(- 2 x) ln(x) + 1/2 Ei(- 2 x)
%
%> sum( 1/(n^2-1), n );
%
%                                           1
%                              - 1/n - -----------
%                                      2 n (n - 1)
%\end{maple}
%
%\noindent
%So the type of the second argument can be either just a
%variable name or, the above type.  We do this as follows
%\verb+{ name, name=algebraic..algebraic }+.
%The use of set braces is for denoting alternatives.
%
Many routines take a set or list of names or equations as arguments.
For example, the {\tt solve} command allows one to solve a set
of equations for a set of unknowns.  For example

\begin{maple}
> solve( {x+y=2, x-y=3}, {x,y} );

                              {y = -1/2, x = 5/2}
\end{maple}
A set of zero or more equations
can be tested with the type: \verb+set(equation)+, and
as set of zero or more unknowns with the type \verb+set(name)+.
But the {\tt solve} command also allows a set of algebraic
formulae which are implicitly equated to zero, i.e. the example
above could have been input this way

\begin{maple}
> solve( {x+y-2, x-y-3}, {x,y} );

                              {y = -1/2, x = 5/2}
\end{maple}

\noindent
Hence the type should of the first argument should be either
a set of equations, or a set of algebraic expressions,
i.e. the type \verb+{set(algebraic),set(equation)}+.
Notice that this is not the same as \verb+set({algebraic,equation})+.
Why?

Further information about the various different types can be obtained
from the on-line help system under ?type .

\subsection{Variable Number of Arguments: {\tt args} and {\tt nargs}}

It is possible in Maple for a function to take a variable number of parameters.
An example of such a function is the \verb+max+ function.
Here is an initial attempt to code up this function

\begin{quote}
\begin{verbatim}
MAX := proc(x1) local maximum,i;
    maximum := x1;
    for i from 2 to nargs do
        if args[i] > maximum then maximum := args[i] fi
    od;
    maximum
end;
\end{verbatim}
\end{quote}

\noindent
The special variable {\tt nargs} is the number of arguments, and the
variable {\tt args} is a sequence of the arguments, hence \verb+args[i]+
is the $i$'th argument.  

\subsection{Returning Unevaluated}

The \verb+MAX+ procedure that we have just written works for numerical
arguments only.  If you try the Maple function \verb+max+ you will see
that it also works for symbolic arguments.  Consider

\begin{quote}
\begin{verbatim}
> MAX(1,2,x);
Error, (in MAX) cannot evaluate boolean
> max(1,2,x);

                                   max(2, x)
\end{verbatim}
\end{quote}

\noindent
Maple cannot execute the procedure \verb+MAX+ because it cannot compute
whether \verb+args[i] < maximum+ for a non-numeric value.
We want \verb+MAX+ of some numbers to compute the answer, but otherwise,
we want \verb+MAX(x,y)+ to stay as \verb+MAX(x,y)+ so we can compute with
\verb+MAX(x,y)+ symbolically just
like \verb+sin(x)+ stays as \verb+sin(x)+.

Also, our \verb+MAX+ procedure that we wrote only works for
numbers of type \verb+numeric+.
It would be nice if \verb+MAX+ knew that $\pi > \sqrt 2$ for instance.

To help us write such a \verb+MAX+ function, we will make use of
the \verb+signum+ function which provides a more powerful comparison
of two real values.
The \verb+signum+ function returns $-1$ if it can
show that the $x < 0$, +1 if $x \ge 0$, otherwise it returns
{\em unevaluated}, i.e. it returns \verb+signum(x)+.
For example


\begin{maple}
> signum(sqrt(2)-1);
                                       1

> signum(sqrt(2)-Pi);
                                       -1

> signum(a-b);
                                 signum(a - b)
\end{maple}


\noindent
Let us employ the signum function to make our \verb+MAX+ function smarter
and also let our \verb+MAX+ function handle symbolic arguments.

\begin{quote}
\begin{verbatim}
MAX := proc() local a,i,j,n,s;
    n := nargs;
    # First, put the arguments in an array
    a := array(1..n); 
    for i to n do a[i] := args[i] od;
    # Compare a[i] with a[j] for 1 <= i < j <= n
    i := 1;
    while i < n do
        j := i+1;
        while j <= n do
            s := signum(a[i]-a[j]);
            if s = 1 then # i.e. a[i] >= a[j]               
                a[j] := a[n]; n := n-1;
            elif s = -1 then # i.e. a[i] < a[j]
                a[i] := a[j]; a[j] := a[n]; j := n; n := n-1; i := i-1;
            else # cannot determine the sign
                j := j+1
            fi
        od;
        i := i+1
    od;
    if n = 1 then RETURN(a[1]) fi;
    'MAX'( seq(a[i], i=1..n) );
end;
\end{verbatim}
\end{quote}

\noindent
What is most interesting about the above code is the last line.
The back quotes \verb+'+ are used to prevent the \verb+MAX+ function
call from executing as otherwise it would go into an infinite loop.
Instead, the unevaluated function call {\tt MAX}( ... ) is returned,
indicating that the maximum could not be computed.  However,
some simplifications may have taken place.  For example

\begin{maple}
> MAX( x, 1, sqrt(2), x+1 );
                                            1/2
                                MAX(x + 1, 2   )
\end{maple}

\subsection{Simplifications and Transformation Rules}

Often one wants to introduce simplifications which can be described
algebraically or by transformation rules.  For instance, given a function
$f$, we may know that $f$ is both commutative and associative.
max is in fact such a function, i.e. it is true that
$\max(a,b) = \max(b,a)$ and
$\max(a,\max(b,c)) = \max(\max(a,b),c) = \max(a,b,c)$.
How can we implement these properties in Maple?
What we want is a {\em canonical} way for writing
expressions involving the max function.
We can implement commutativity by sorting the arguments.
For associativity we can unnest any nested max calls.
I.e. we would transform both $\max(\max(a,b),c)$ and $\max(a,\max(b,c)$ into
$\max(a,b,c)$.  Actually this also implements $\max(\max(a)) = \max(a)$, i.e.
max is idempotent.  Here is a \verb+MAX+ function to do this

\begin{quote}
\begin{verbatim}
MAX := proc() local a;
    a := [args];
    a := map( flatten, a, MAX ); # unnest nested MAX calls
    'MAX'( op(sort(a)) );
end;
flatten := proc(x,f)
    if type(x,function) and op(0,x) = f then op(x) else x fi
end;
\end{verbatim}
\end{quote}

\noindent
For example
\begin{maple}
> MAX(a,MAX(c,b),a);

                                MAX(a, a, b, c)
\end{maple}

We see that we should also recognize the property that
$\max(a,a) = a$.  To do this, instead of putting the arguments in a
list, we will put them in a set so that duplicates are removed.
Also, since sets are sorted automatically, we can remove the call to sort.
Hence we have

\begin{quote}
\begin{verbatim}
MAX := proc() local a;
    a := {args};
    a := map( flatten, a, MAX );
    'MAX'( op(a) );
end;
\end{verbatim}
\end{quote}

\begin{maple}
> MAX(a,MAX(c,b),a);

                                  MAX(a, b, c)
\end{maple}

The reader may be a little puzzled as to just what our {\tt MAX} procedure
is doing.  We have seen earlier that if we assign a positive integer
to the {\tt printlevel} variable, we get a trace of all statements executed.
However, often the output from this simple tracing facility is too much.
In this case, we would also get the output from the {\tt flatten} procedure.
The {\tt trace} function can be used instead to selectively trace procedures.
Let's use it to trace the {\tt MAX} procedure

\begin{maple}
> trace(MAX);

                                      MAX

> MAX(a,MAX(b,a),c);
{--> enter MAX, args = b, a

                                  a := {a, b}

                                  a := {a, b}

                                   MAX(a, b)

<-- exit MAX (now at top level) = MAX(a,b)}
{--> enter MAX, args = a, MAX(a,b), c

                             a := {a, c, MAX(a, b)}

                                 a := {a, c, b}

                                  MAX(a, c, b)

<-- exit MAX (now at top level) = MAX(a,c,b)}

                                  MAX(a, c, b)
\end{maple}


\subsection{Optional Arguments and Default Values}

Many Maple routines accept optional arguments.  This is often used to
allow the user to use default values instead of having to specify
all parameters.  Examples are the functions {\tt plot}, {\tt factor},
{\tt collect}, and {\tt series}.  Let us consider the {\tt degree} function.
The {\tt degree} function computes the degree of a univariate polynomial in
one variable, and for multivariate polynomials, the total degree.
For example

\begin{maple}
> p := x^3+2*x+1;

                                     3
                               p := x  + 2 x + 1

> degree(p);

                                       3

> q := 3*x^2*y+2*y^2-x*z+7;

                                  2        2
                          q := 3 x  y + 2 y  - x z + 7

> degree(q);

                                       3
\end{maple}

\noindent
Sometimes you will want to compute the degree in a specific variable, say $x$.
This can be done by specifying an optional second argument to the degree
function, namely, the variable.  For example

\begin{maple}
> degree(p,x);

                                       3

> degree(q,x);

                                       2
\end{maple}

\noindent
How would we code the degree function?  Let us assume that the input
is a formula and if an optional second argument is given, it
is a name or set of names for the variables.
We would write

\begin{quote}
\begin{verbatim}
DEGREE := proc(a:algebraic,x:{name,set(name)}) local s,t;
    if nargs = 1 then # determine the variable(s) for the user
        s := indets(a); # the set of all the variables of a
        if not type(s,set(name)) then ERROR(`input not a polynomial`) fi;
        DEGREE(a,s)
    elif type(a,constant) then 0
    elif type(a,name) then
        if type(x,name) then if a = x then 1 else 0 fi
        else if member(a,x) then 1 else 0 fi
        fi
    elif type(a,`+`) then max( seq( DEGREE(t,x), t=a ) )
    elif type(a,`*`) then
        s := 0;
        for t in a do s := s + DEGREE(t,x) od;
        s
    elif type(a,algebraic^integer) then DEGREE(op(1,a),x) * op(2,a)
    else ERROR(`cannot compute degree`)
    fi
end;
\end{verbatim}
\end{quote}

The {\tt indets} function used here returns a set of all the 
{\em indeterminates} (or variables) that appear in the input.
We leave it to the reader to study each rule that is being
used here, and the order in which the cases are done.

\subsection{Returning Results Through Parameters}

Many functions in Maple return more than one value.  Of course, it
is always possible to return more than one value in a sequence or list.
However, it is also possible to return values through parameters like
in other programming languages, and often, this is more convenient.
For example, consider the {\tt divide} function in Maple
which does polynomial long division.  The call {\tt divide(a,b)} returns
true if and only if the polynomial $b$ divides the polynomial $a$ with
no remainder, e.g.

\begin{maple}
> divide(x^3-1,x-1);

                                      true
\end{maple}

\noindent
But usually, if $b$ divides $a$, one wants to do something with the
quotient $q$.  In Maple, this can be done by giving the {\tt divide}
function a third parameter which is a name which will be assigned
the quotient if $b$ divides $a$, e.g.

\begin{maple}
> if divide(x^3-1,x-1,'q') then print(q) fi;

                                    2
                                   x  + x + 1
\end{maple}

\noindent
Notice the use of quotes here to pass to the divide function the name $q$
and not the value of $q$.  Let us consider another example and study how
to write a program which assigns a value to an optional parameter.  Consider
our {\tt MEMBER} function which tested to see if a value $x$ appears in
a list $L$.  Let us modify our function such that {\tt MEMBER(x,L,'p')}
still returns whether $x$ appears in the list $L$, and in addition, assigns
the name $p$ the position of the first appearance of $x$ in $L$.

\begin{quote}
\begin{verbatim}
MEMBER := proc(x,L:list,p:name) local i;
    for i to nops(L) do
        if x = L[i] then
            if nargs = 3 then p := i fi;
            RETURN(true)
        fi
    od;
    false
end;
\end{verbatim}
\end{quote}

\noindent
Here is an example

\begin{maple}
> MEMBER(4,[1,3,5],'position');

                                     false

> position;

                                    position

> MEMBER(3,[1,3,5],'position');

                                      true

> position;

                                       2
\end{maple}

\noindent
We see that the effect of the assignment to the formal parameter $p$
inside the {\tt MEMBER} procedure is that the actual parameter {\tt position}
is assigned.

\section{Programming in Maple}

\subsection{Matrix and Vector Computation in Maple}

A vector in Maple is represented by a one-dimensional
array indexed from 1, and a matrix is represented by a two-dimensional
array, with row and column indices starting from 1.
Here is one way to create a 5 by 5 Hilbert matrix.
Recall that a Hilbert matrix is a symmetric matrix whose (i,j)'th
entry is $1/(i+j-1)$.

\begin{maple}
> H := array(1..5,1..5):
> for i to 5 do for j to 5 do H[i,j] := 1/(i+j-1) od od;
> H;

                                       H

> eval(H);
                          [  1   1/2  1/3  1/4  1/5 ]
                          [                         ]
                          [ 1/2  1/3  1/4  1/5  1/6 ]
                          [                         ]
                          [ 1/3  1/4  1/5  1/6  1/7 ]
                          [                         ]
                          [ 1/4  1/5  1/6  1/7  1/8 ]
                          [                         ]
                          [ 1/5  1/6  1/7  1/8  1/9 ]
\end{maple}

\vspace{2mm}
\sbox{
\parbox{6in}{
ATTENTION: Notice that the value of \verb+H+ is just the name of the
matrix \verb+H+.  Evaluation rules for arrays, hence matrices and vectors
are special.  The reason is technical.  For the moment, whenever you
want to print a matrix or vector use the \verb+eval+ function. }
}
\vspace{2mm}

The \verb+linalg+ package contains many functions for computing with
vectors and matrices in Maple.  This matrix could also have been created
in the following way using the \verb+matrix+ command in the linear
algebra package.

\begin{quote}
\begin{verbatim}
linalg[matrix](5,5,(i,j) -> 1/(i+j-1));
\end{verbatim}
\end{quote}

Maple can compute the determinant and inverse of the matrix H and many
other matrix operations.  See \verb+?linalg+ for a list of the operations.

Here we show a program for doing a so called row reduction on a matrix
using Gaussian elimination.  \verb+GaussianElimination(A,'r')+
computes the reduced matrix, an upper triangular matrix.  It also
has an optional 2nd argument which, if given, is assigned the
rank of the matrix.

\begin{quote}
\begin{verbatim}
GaussianElimination := proc(A:matrix(rational),rank:name)
local m,n,i,j,B,r,c,t;

    m := linalg[rowdim](A); # the number of rows of the matrix
    n := linalg[coldim](A); # the number of columns of the matrix

    B := array(1..m,1..n);
    for i to m do for j to n do B[i,j] := A[i,j] od od;

    r := 1; # r and c are row and column indices
    for c to m while r <= n do

        for i from r to n while B[i,c] = 0 do od;  # search for a pivot
        if i <= n then
            if i <> r then # interchange row i with row r
                for j from c to m do
                    t := B[i,j]; B[i,j] := B[r,j]; B[r,j] := t
                od
            fi;
            for i from r+1 to n do
                if B[i,c] <> 0 then
                    t := B[i,c]/B[r,c];
                    for j from c+1 to m do B[i,j] := B[i,j]-t*B[r,j] od;
                    B[i,c] := 0
                fi
            od;
            r := r + 1                # go to next row
        fi
    od;                       # go to next column

    if nargs>1 then rank := r-1 fi;
    eval(B)

end:
\end{verbatim}
\end{quote}

\noindent
The type {\tt matrix(rational)} specifies a matrix (a 2-dimensional
array in Maple) whose entries are all rational numbers.

\subsection{Numerical Computation in Maple}

Floating point numbers in Maple can be input either as decimal
numbers or directly using the \verb+Float+ function

\begin{center}
\verb+Float(m,e)+  =  $m*10^e$
\end{center}

\noindent
where the mantissa $m$ is an integer of any size, but
the exponent $e$ is restricted to be a machine size integer,
typically 31 bits.  For example, the number 3.1 and be input as
either {\tt 3.1} or {\tt Float(31,-1)}.
The \verb+op+ function can be used to
extract the mantissa and exponent of a floating point number.
Note: the floating point constant 0.0 is treated specially and
is simplified to 0, i.e. the integer 0 automatically.

Here is a uniform random number generator that generates uniform
random numbers with exactly 6 decimal digits precision in the
range $[0,1)$.

\begin{maple}
> UniformInteger := rand(0..10^6-1):
> UniformFloat := proc() Float(UniformInteger(),-6) end:
> seq( UniformFloat(), i=1..6 );

              .669081, .693270, .073697, .143563, .718976, .830538
\end{maple}

\noindent
Note, the built-in Maple function \verb+rand+ does
not return a random number in the given range.
It returns instead a Maple procedure (a random number generator)
which when called returns random integers in the given range.

Floating point arithmetic is done in decimal with rounding.
The precision used is controlled by the global variable \verb+Digits+ which
has 10 as its default value.  It can be set to any value however.
The \verb+evalf+ function is used to evaluate an exact
symbolic constant to a floating point approximation.
For example

\begin{maple}
> Digits := 25:
> sin(1.0);
                           .8414709848078965066525023

> sin(1);
                                     sin(1)

> evalf(");
                           .8414709848078965066525023
\end{maple}

Maple knows about the elementary functions and many special functions
such as $J_v (x)$, $\Gamma (x)$, and $\zeta (x)$.  In Maple these are
\verb+BesselJ(v,x)+ \verb+GAMMA(x)+ and \verb+Zeta(x)+ respectively.
They are all computed to high precision by summing various series.

The model of floating point arithmetic used is that the
relative error of the result is less than $10^{1-{\rm Digits}}$.
This is a stronger model than that used by hardware implementations
of floating point arithmetic and it requires that intermediate calculations
to be done at higher precision to obtain such accurate results.
Here is an example of summing a Taylor series.
Suppose we want to compute the error function

\[
{\rm erf}(x) ~=~ \frac{2}{\sqrt{\pi}} \int _0 ^x e^{-t^2} {\rm d}t
\]

\noindent
for small $x$.  If you are unfamiliar with the error function, take a moment
to plot it in Maple.
We can make use of the Taylor series for ${\rm erf}(x)$ about $x=0$

\[
\frac{\sqrt {\pi}}{2} {\rm erf}(x) ~=~
\sum _{n=0} ^ \infty \frac{ (-1)^n x^{2 n + 1} } {n! (2 n + 1) } ~=~ 
x-{\frac {x^{3}}{3}}+{\frac {x^{5}}{10}}-{\frac {x^{7}}{42}}+{\frac {x
^{9}}{216}}+ \ldots
\]

\noindent
for computing ${\rm erf}(x)$ for small $x$, $-1 < x < 1$
(with no error checking) as follows

\begin{quote}
\begin{verbatim}
ErfSmall := proc(a) local x,x2,result,term,sumold,sumnew;
    x := evalf(a); # evaluate the input at Digits precision
    Digits := Digits + 2; # add some guard digits
    sumold := 0;
    term := x;
    sumnew := x;
    x2 := x^2;
    for n from 1 while sumold <> sumnew do
        sumold := sumnew;
        term := - term * x2 / n;
        sumnew := sumold + term / (2*n + 1);
    od;
    result := evalf( 2/sqrt(Pi)*sumnew );
    Digits := Digits-2;
    evalf(result) # round the result to Digits precision
end;
\end{verbatim}
\end{quote}

\noindent
For large $x$, this routine will be inaccurate and inefficient.
See the exercises for how to compute ${\rm erf}(x)$ for large $x$ .

For a second example, let us consider writing a Newton iteration
to compute the square root of an integer to 50 Digits.
Recall that the Newton iteration for solving an equation $f(x) = 0$ is

\[
x_{k+1} ~=~ x_k - \frac{f(x_k)}{f'(x_k)}
\]

\noindent
In our case, to compute $\sqrt a$ we want to solve $x^2 - a = 0$.
So the iteration is

\[
x_{k+1} ~=~ ( x_k + \frac{a}{x_k} ) / 2
\]

\noindent
The following routine will do the job.

\begin{quote}
\begin{verbatim}
SqrtNewton := proc(a:integer) local xk, xkm1;
    if a < 0 then ERROR(`square root of a negative integer`) fi;
    Digits := 55; # add some guard digits
    xkm1 := 0;
    xk := evalf(a/2); # initial floating point approximation
    while abs(xk-xkm1) > abs(xk)*10^(-50) do 
        xkm1 := xk;
        print(xk);
        xk := (xk + a/xk)/2;
    od;
    Digits := 50;
    evalf(xk); # round the result to 50 Digits
end;
\end{verbatim}
\end{quote}

\noindent
We have put in the print statement to show the result of each iteration.
Here is what happens

\begin{maple}
> SqrtNewton(2);

                                       1.

            1.500000000000000000000000000000000000000000000000000000

            1.416666666666666666666666666666666666666666666666666667

            1.414215686274509803921568627450980392156862745098039216

            1.414213562374689910626295578890134910116559622115744045

            1.414213562373095048801689623502530243614981925776197429

            1.414213562373095048801688724209698078569671875377234002

            1.414213562373095048801688724209698078569671875376948073

            1.414213562373095048801688724209698078569671875376948073
\end{maple}

It is known that a Newton iteration will converge quadratically
provided $f'(x_k)$ is not close to zero and $x_k$ is sufficiently
close to a root of $f$, as illustrated in the above example.
For very high precision though, i.e. Digits $>$ 1000, one does not want
to do every iteration at full precision as the early steps are not
accurate to full precision.  Why do all that work when you
are only getting a few digits correct?
Modify the Newton iteration used in the SqrtNewton procedure
appropriately so that it doubles the number of Digits at each step.
How much faster does this make the iteration run?

Because Maple's floating point model is implemented in software,
it is much slower than hardware floating point arithmetic.
Maple also has a function called \verb+evalhf+ for
{\tt eval}uating in {\tt h}ardware {\tt f}loating point arithmetic.
This function uses the builtin C library routines for floating
point arithmetic.  Consequently it is much faster than Maples software
floats, but it is still slower than hardware floats because
it is not compiled.  See \verb+?evalhf+ for details.

\subsection{Computing with Polynomials in Maple}

Computing with polynomials and also rational functions is Maple's forte.
Here is a program that computes the Euclidean norm of a polynomial.
I.e. given the polynomial $a(x) = \sum _{i=0} ^ n a_i x^i$, it computes
the $ \sqrt { \sum _{i=0} ^n a_i^2 } $.

\begin{quote}
\begin{verbatim}
EuclideanNorm := proc(a)
    sqrt( convert( map( x -> x^2, [coeffs( expand(a) )] ), `+` ) )
end;
\end{verbatim}
\end{quote}

\noindent
Reading this one liner inside out, the input polynomial is first expanded.
The \verb+coeffs+ function returns a sequence of the coefficients
which we have put in a list.
Each element of the list is squared yielding a new list.
Then the list of squares is converted to a sum.

Why is the polynomial expanded?
The \verb+coeff+ and \verb+coeffs+ functions insist on having the input
polynomial {\em expanded} because you cannot compute the
coefficient(s) otherwise.  The following example should clarify

\begin{maple}
> p := x^3 - (x-3)*(x^2+x) + 1;

                               3             2
                         p := x  - (x - 3) (x  + x) + 1

> coeffs(p);
Error, invalid arguments to coeffs
> expand(p);
                                    2
                                 2 x  + 3 x + 1

> coeffs(expand(p));
                                    1, 3, 2

> EuclideanNorm(p);
                                       1/2
                                     14
\end{maple}

\vspace{2mm}
\sbox{
\parbox{6in}{
ATTENTION: If you want to compute the coefficients of a polynomial
in some variable, you should always expand the polynomial in that
variable first.  You can use the \verb+expand+ function but this
expands the polynomial in all variables.  Alternatively, you can use
the \verb+collect+ function.  \verb+collect(p,x)+ expands a polynomial
$p$ in $x$ only.  See \verb+?collect+ for other details and for
working with multivariate polynomials. }
}
\vspace{2mm}

The \verb+EuclideanNorm+ procedure works for multivariate
polynomials in the sense
that it computes the square root of the sum of the squares of the numerical
coefficients.  E.g. given the polynomial $p = u x^2 + y^2 + v$,
the \verb+EuclideanNorm+ procedure returns $\sqrt 3$.
However, you may want to view this polynomial as a polynomial
in $x,y$ whose coefficients are symbolic
coefficients in $u,v$.
We really want to be able to tell the \verb+EuclideanNorm+
routine what the polynomial variables are.  We can do that by specifying an
additional optional parameter, the variables, as follows

\begin{quote}
\begin{verbatim}
EuclideanNorm := proc(a,v:{name,set(name),list(name)})
    if nargs = 1 then
        sqrt( convert( map( x -> x^2, [coeffs(expand(a))] ), `+` ) )
    else
        sqrt( convert( map( x -> x^2, [coeffs(expand(a),v)] ), `+` ) )
    fi
end;
\end{verbatim}
\end{quote}

\noindent
The type \verb+{name,set(name),list(name)}+ means that the
2nd parameter $v$ may be a single variable, or a set of variables,
or a list of variables.
Notice that the \verb+coeffs+ function itself accepts this argument
as a 2nd optional argument.
Finally, our routine doesn't ensure that the input is a polynomial.
Let's add that to

\begin{quote}
\begin{verbatim}
EuclideanNorm := proc(a,v:{name,list(name),set(name)})
    if nargs = 1 then
        if not type(a,polynom) then
            ERROR(`1st argument is not a polynomial`,a) fi;
        sqrt( convert( map( x -> x^2, [coeffs(expand(a))] ), `+` ) )
    else
        if not type(a,polynom(anything,v)) then
            ERROR(`1st argument is not a polynomial in`,v) fi;
        sqrt( convert( map( x -> x^2, [coeffs(expand(a),v)] ), `+` ) )
    fi
end;
\end{verbatim}
\end{quote}

\noindent
The type \verb+polynom+ has the following general syntax

\begin{center}
polynom( $R$, $X$ )
\end{center}

\noindent
which means a polynomial whose coefficients are of type $R$ in the
variables $X$.  An example is \verb+polynom(rational,x)+ which specifies
a univariate polynomial in $x$ with rational coefficients,
i.e. a polynomial in $Q[x]$.  If $R$ and $X$ are not specified, as in the
first case above, the expression must be a polynomial in all
its variables.

Basic functions for computing with polynomials are {\tt degree}, {\tt coeff},
{\tt expand}, {\tt divide}, and {\tt collect}.
There are many other functions for computing with polynomials
in Maple, including facilities for polynomial division, greatest common
divisors, resultants, etc.  Maple can also factor polynomials over
different number fields including finite fields.  See the on-line help
for {\tt ?polynom} for a list of facilities.

%We include here a simple program for computing the greatest common divisor
%of two univariate polynomials $a,b$ using the Euclidean algorithm.
%Note, Maple does not use the Euclidean algorithm because this turns
%out to be very inefficient for polynomials of high degree.
%However, this example summarizes many of the things we have mentioned
%so far.
%
%\begin{quote}
%\begin{verbatim}
%GCD := proc(a,b,x) local c,d,r;
%    if not type(x,name) then ERROR(`3rd argument (variable) must be a name`) fi;
%    if not type([a,b],'polynom(anything,x)') then
%        ERROR(`1st and 2nd arguments must be polynomials in `,x)
%    fi;
%    c := collect(a,x,recursive,normal);
%    d := collect(a,x,recursive,normal);
%    while d <> 0 do r := rem(c,d,x); c := d; d := r od;
%    if c <> 0 then # Make c monic
%        c := collect(c/lcoeff(c,x),x,recursive,normal);
%    fi;
%    c
%end;
%\end{verbatim}
%\end{quote}
%
%\noindent
%The use of {\tt collect} here is used to simplify the coefficients
%of a polynomial, in particular, to recognize the zero polynomial.
%The call {\tt collect(a,x,recursive,normal)} means group together
%the coefficients of like powers of $x$, and apply the normal function
%to them to simplify them.

To illustrate facilities for computing with polynomials over finite fields,
we conclude with a program to compute the first {\em primitive trinomial} of
degree $n$ over $GF(2)$ if one exists .
That is, we want to find an irreducible polynomial $a$ of the form $x^n + x^m + 1$ such
that $x$ is a primitive element in $GF(2)[x]/(a)$.  The {\tt iquo} function
computes the {\bf i}nteger {\bf quo}tient of two integers.

\begin{quote}
\begin{verbatim}
trinomial := proc(n:integer) local i,t;
    for i to iquo(n+1,2) do
        t := x^n+x^i+1;
        if Primitive(t) mod 2 then RETURN(t) fi;
    od;
    FAIL
end;
\end{verbatim}
\end{quote}


\subsection{Reading and Saving Procedures: {\tt read} and {\tt save}}

You can write one or two line Maple programs interactively.
But for larger programs you will want to save them in a file.
Typically, one would use an editor to write the Maple program and
save it into a file, then read the program into Maple before using it.
A program can be read into Maple using the read command.
For example, if we have written a Maple procedure \verb+MAX+ in the
file {\tt MAX}, in Maple we read this file into Maple by doing

\begin{quote}
\begin{verbatim}
read MAX;
\end{verbatim}
\end{quote}

You can save your Maple procedures or any formulae that you have computed
in your Maple session in a file from inside Maple using the \verb+save+
statement, which has the form

\begin{center}
\verb+save f1, f2, ..., filename;+
\end{center}

\noindent
This saves the values of the variables \verb+f1, f2, ...+
in text format in the file \verb+filename+.
You can also save Maple data in internal format or the so called ``.m''
format.  This format is more compact and can be read faster by Maple.
One simply appends ``.m'' to the filename as follows

\begin{center}
\verb+save f1, f2, ..., `filename.m`;+
\end{center}

\noindent
This saves the values of \verb+f1, f2, ...+ in the file \verb+filename.m+
You can then read them back in to Maple using the {\tt read} command:

\begin{quote}
\begin{verbatim}
read `filename.m`;
\end{verbatim}
\end{quote}

\subsection{Debugging Maple Programs}

The simplest debugging tool is the {\em printlevel} facility.
{\tt printlevel} is a global variable that is initially assigned 1.
If you set it to a higher value, a trace of all assignments, procedure
entry and exits are printed.
The higher the value of {\tt printlevel}, the more levels of
procedure execution that will be traced.
Often, however, the output from the printlevel facility will be too much.
The {\tt trace} function allows you to trace the execution of
the specified functions only.
Examples of these two tools have already been given in this document.
We mention here one further tool.  If you set the {\tt printlevel}
variable to 3 or higher, then if a run-time error occurs, Maple will print
a stack trace of the calling sequence at the time the error occurs.
Specifically, it will print the arguments to all procedures currently
being executed, and the values of the local variables
and the statement being executed in the procedure where the error occured.
This is simplest illustrated by an example.

\begin{maple}
> f := proc(x) local y; y := 1; g(x,y); end:
> g := proc(u,v) local s,t; s := 0; t := v/s; s+t end:
> printlevel := 4:
> f(3);
Error, (in g) division by zero
        executing statement: t := v/s
        locals defined as: s = 0, t = t
        g called with arguments: 3, 1
        f called with arguments: 3
\end{maple}

\subsection{Interfacing with other Maple Facilities}

We have shown how one can tell Maple properties about a function $f$ by
coding them as a procedure.  Suppose instead you wish to teach Maple to
differentiate a formula involving $f$.
You may want to teach Maple how to evaluate $f$ numerically so that $f$
can be plotted, or how to simplify expressions involving $f$ etc.
What do you need to do? 
Many Maple routines have interfaces that allow you to teach these
routines about your $f$ function.
These include {\tt diff}, {\tt evalf}, {\tt expand}, {\tt combine},
{\tt simplify}, {\tt series}, etc.
To teach Maple how to differentiate a function
$W$ one writes a routine called \verb+`diff/W`+.
If the \verb+diff+ routine is called with an expression $f(x)$ which contains
$W(g)$ then the \verb+diff+ routine will invoke \verb+`diff/W`(g,x)+ to
compute the derivative of $W(g)$ with respect to $x$.
Suppose we know that $W'(x) = W(x)/(1+W(x))$.
Then we can write the following procedure which
explicitly codes the chain rule.

\begin{quote}
\begin{verbatim}
`diff/W` := proc(g,x) diff(g,x) * W(g)/(1+W(g)) end;
\end{verbatim}
\end{quote}

\noindent
Hence we have

\begin{maple}
> diff(W(x),x);
                                      W(x)
                                    --------
                                    1 + W(x)

> diff(W(x^2),x);
                                          2
                                     x W(x )
                                  2 ---------
                                           2
                                    1 + W(x )
\end{maple}

As a second example, suppose we want to manipulate symbolically
the Chebyshev polynomials of the first kind $T_n(x)$.  We can represent these in Maple
as \verb+T(n,x)+.  Suppose also that for a particular value of $n$
we want to expand \verb+T(n,x)+ out as a polynomial in $x$.
We can tell Maple's \verb+expand+ function how to do this by writing
the routine `expand/T`.
When \verb+expand+ sees \verb+T(n,x)+,
it will invoke \verb+`expand/T`(n,x)+.
Recall that $T_n(x)$ satisfies the linear recurrence $T_n(0) = 1$,
$T_1(x) = x$, and $T_n(x) = 2 x T_{n-1}(x) - T_{n-2}(x)$.
Hence we can write

\begin{quote}
\begin{verbatim}
`expand/T` := proc(n,x) option remember;
    if n = 0 then 1 
    elif n = 1 then x
    elif not type(n,integer) then T(n,x) # can't do anything
    else expand(2*x*T(n-1,x) - T(n-2,x))
    fi
end;
\end{verbatim}
\end{quote}

\noindent
This routine is recursive, but because we used the remember option,
this routine will compute $T(100,x)$ quite quickly.
Here are some examples

\begin{maple}
> T(4,x);
                                   T(4, x)

> expand(T(4,x));
                                   4      2
                                8 x  - 8 x  + 1

> expand(T(100,x));

             2                                                           100
   1 - 5000 x  ... output deleted ...  + 633825300114114700748351602688 x
\end{maple}

One can also tell Maple how to evaluate your own function $f$ numerically by
defining a Maple procedure \verb+`evalf/f`+ such that
\verb+`evalf/f`(x)+ computes $f(x)$ numerically.
For example, suppose we wanted to use the Newton iteration shown
earlier for computing the square root of a numerical value.
Our function might look like

\begin{quote}
\begin{verbatim}
`evalf/Sqrt` := proc(a) local x,xk,xkm1;
    x := evalf(a); # evaluate the argument in floating point
    if not type(a,numeric) then RETURN( Sqrt(x) ) fi;
    if x<0 then ERROR(`square root of a negative number`) fi;
    Digits := Digits + 3; # add some guard digits
    xkm1 := 0;
    xk := evalf(x/2); # initial floating point approximation
    while abs(xk-xkm1) > abs(xk)*10^(-Digits) do
        xkm1 := xk;
        xk := (xk + x/xk)/2;
    od;
    Digits := Digits - 3;
    evalf(xk); # round the result to Digits precision
end;
\end{verbatim}
\end{quote}

\begin{maple}
> x := Sqrt(3);
                                  x := Sqrt(3)

> evalf(x);
                                  1.732050808

> Digits := 50;
                                  Digits := 50

> evalf(Sqrt(3));

              1.7320508075688772935274463415058723669428052538104

\end{maple}

\subsection{Calling External Programs}

Often one wants to get Maple to call an external program, which
might have been written in C or Fortran.
There are really two reasons we might want to do this.
Obviously, if Maple cannot so something, and another program can,
it may make sense to use the other program, instead of reimplementing
an algorithm in Maple.  The other reason is efficiency.
Although Maple is efficient at symbolic computations, it is not generally
efficient at machine precision numerical computations.  E.g. you
may want to use a Fortran library routine to compute numerical eigenvectors.

Communication of data in Maple must be done via files.  
The Maple program must write the data needed by the external program
out to an input file.
Maple can start running the external program
with the {\tt system} command.
The external program must read the data from the input file
and must write its results into an output file.
After the program has finished executing, Maple reads the results from
the output file back into Maple.
A sketch of the Maple code to do this would be

\begin{quote}
\begin{verbatim}
interface(quiet=true);
writeto(input);
... # write any data into the file input
writeto(terminal);
interface(quiet=false);
system(...); # execute the external program
read output;
... # continue processing in Maple
\end{verbatim}
\end{quote}

Let us look at this more closely.  The first command,
{\tt interface(quiet=true);} turns off all diagnostics from Maple, i.e.
bytes used messages, warnings, etc. so that they will not be
written into the {\tt input} file.
This is reset by {\tt interface(quiet=false);}
after the data has been written to the {\tt input} file.

The next command {\tt writeto(input);} opens the file named {\tt input}
for writing.  All Maple output from this point will go into this file.
Output will typically be created by the {\tt lprint}
or the {\tt printf} commands.
The command {\tt writeto} overwrites the file.
If you want to just append some data to what is
already in the file, use the command {\tt appendto} instead of {\tt writeto}.

After the data has been written to the file, the file is closed
implicitly by resetting output back to the terminal with the command
{\tt writeto(terminal);}.

The {\tt system} command is used to execute the external program.
The exact call that the system command makes will depend on the system.
Under Unix, the call might look something like

\begin{quote}
\begin{verbatim}
system(`foo < in > out`);
\end{verbatim}
\end{quote}

The external program {\tt foo} is executed on the input file {\tt in}
and its results are written to the  file {\tt out}.
Note that the system command returns a value indicating whether the external
program terminated normally.  Under Unix, it returns 0 for normal
termination.

Finally, we read the results in the {\tt out} file back into
Maple using the {\tt read} command if the data is in Maple's format,
or the {\tt readline}, {\tt sscanf}, or {\tt readdata} commands
described below for reading arbitrary data.

In this mode, Maple calls the external program like a subroutine.
After reading the results back into Maple, Maple continues execution.
This is actually a very simple means for calling an external program,
and easy to debug.  You can check the files to make sure that
they have the data in the format expected.
However, the disadvantage is that communicating
data via a file of text is not the most efficient approach.
In many applications this
inefficiency will not matter if the work done by the external
program is significant.

\subsection{File Input/Output of Numerical Data}

We wish now to discuss in more detail how to output data from Maple
into a file suitable for reading by another program, and how the data
should be formatted by the external program so that Maple can read
it in.

In Maple V Release 2 the file I/O facilities now include the
the routines {\tt printf}, {\tt sscanf} and {\tt readline}.
The routine {\tt printf} is used for general formatted printing.
It is modelled after the printf routine from C and is the same except
that it also accepts the \%a option for algebraic output.
Here is a typical example of using {\tt printf} to print floating point data.

\begin{maple}
> printf(`A float in fixed point: %f and scientific notation: %e\n`, x, x );
A float in fixed point: 0.577216 and scientific notation: 5.772157e-01

\end{maple}
The {\tt printf} function takes a string as it's first argument which
specifies text to be printed and how the values, which are
given as additional arguments are to be printed using \% control sequences.   
\%e control sequences.  The escape sequence \verb+\n+ is used to print a new line.
Other useful control sequences are \%d and \%s for printing integers and
strings, and also the control sequence \%a for printing a Maple algebraic
expression.  See {\tt ?printf} for other options.

The {\tt readline} routine is used to read a single line of text from
a text file as a string.  The {\tt sscanf} routine can be used to scan
a string and extract any numbers.  It is the inverse of {\tt printf}.
So to read in some data, one can read each line of a file using {\tt readline}
and for each line use {\tt sscanf} to extract data.
The utility routine {\tt readdata} has been written to read files of
text which contain columns of numerical data separated by blanks or tabs.
The {\tt readdata} routine returns a list of the data, or a list of lists
of the data if more than one column is specified.
Suppose the file {\tt foo} looks like this
\begin{quote}
\begin{verbatim}
2.1    3.2   1.2
2.4    3.1   3.4
3.9    1.1   5.5
\end{verbatim}
\end{quote}
\noindent
Here is the result of reading in the data.
\begin{maple}
> readlib(readdata): # load the readdata function from the library
> readdata(foo,2); # read the first two columns of the data

                  [[2.100000000000000, 3.200000000000000],

                     [2.400000000000000, 3.100000000000000],

                     [3.900000000000000, 1.100000000000000]]

\end{maple}

\subsection{Fortran and C output}

The {\tt fortran} and {\tt C} commands
generate output of Maple formulae in a format suitable for
a Fortran or C compiler.  This is useful when you have developed a
formula in Maple, or perhaps a vector or matrix of formulae,
and you wish to evaluate these formulae in a Fortran or C subroutine.
How do you translate the Maple formulae into Fortran or C?
The Maple functions {\tt fortran} and {\tt C} are written for this purpose.
Here is an example.  Suppose we have created the following polynomial
in Maple which, by the way, is an approximation to the complementary
error function ${\rm erfc}(x) = 1-{\rm erf}(x)$ on the range [2,4] accurate
to 5 decimal digits.

\begin{maple}
 f :=  - 3.902704411 x + 1.890740683 - 1.714727839 x  + 3.465590348 x

                        7                4                 5                  6
      - .0003861021174 x  + .5101467996 x  - .09119265524 x  + .009063185478 x
\end{maple}

It is not important to know what the complementary error function is for
the purpose of this example though.
It is in fact related to the Normal distribution in statistics.
Neither is it important here to know how we created the approximation $f$.
We needed a rough approximation to this function in the given
range because our Fortran and C libraries did not have this
function built in.  For the interested reader, we used the command
\verb:chebyshev(erfc(x),x=2..4,10^(-5): to create a Chebyshev series
approximation and used our {\tt `expand/T`} routine that we wrote earlier
to convert it to a polynomial.
To evaluate the approximation $f$ above efficiently,
we want to write the polynomial in Horner form.
Then we want to generate Fortran code.  We could do

\begin{maple}
> h := convert(f,horner);

h := 1.890740683 + ( - 3.902704411 + (3.465590348 + ( - 1.714727839

 + (.5101467996 + ( - .09119265524 + (.009063185478 - .0003861021174 x) x) x) x

) x) x) x

> fortran(h);
      t0 = 0.1890741E1+(-0.3902704E1+(0.346559E1+(-0.1714728E1+(0.510146
     #8E0+(-0.9119266E-1+(0.9063185E-2-0.3861021E-3*x)*x)*x)*x)*x)*x)*x
\end{maple}

Maple has generated two lines of Fortran code complete with continuation
character since the formula is longer than one line.  The floating point
numbers have automatically been translated to single precision Fortran {\tt E}
notation and have been truncated to 7 decimal digits.
And Maple has put the result in the variable {\tt t0} for us.
Let us now output C code into a file {\tt temp.c} assigned to
the variable $r$.  Note, the {\tt C} function must first
be loaded into Maple.

\begin{maple}
> readlib(C):
> C([r = h],filename=`temp.c`);

\end{maple}


\noindent
If we look at the file `temp.c` we find that it contains

\begin{verbatim}
      r = 0.1890740683E1+(-0.3902704411E1+(0.3465590348E1+(-0.1714727839E1+(
0.5101467996+(-0.9119265524E-1+(0.9063185478E-2-0.3861021174E-3*x)*x)*x)*x)*x)*
x)*x;
\end{verbatim}

The reader can study the help pages for {\tt ?fortran} and {\tt ?C} for
additional information capabilities and options to these commands.

\newpage

\section{Exercises}

\begin{enumerate}

\item
Write a Maple procedure called {\tt remove1} where
{\tt remove1(x,L)} removes the first occurrence of the value $x$
from the list $L$.  If $x$ is not in the list $L$, return FAIL.

\item
Write a Maple procedure called \verb+variance+ that computes
the variance of a list of numbers.  I.e. \verb+variance(x)+ computes

\[
\frac{1}{n} ~ \sum _{i=1}^n ~ (x_i - \mu )^2
\]

\noindent
where $n$ is the number of elements of the list $x$ and $\mu$ is the
average of the numbers in list $x$.
Your routine should print an error if the list is empty.

\item

Write a Maple procedure that computes the Frobenius norm of an $m$ by $n$ matrix $A$.
The Frobenius norm is $\sqrt { \sum _{i=1} ^ m \sum _{j=1} ^n |A_{i,j}|^2 }$.

\item

Write a Maple procedure which sorts an array of numbers
faster than the {\em bubblesort} algorithm in the section on arrays,
e.g. by using shell sort, quick sort, merge sort, or heap sort or
some other sorting algorithm.
Modify your Maple procedure to accept an optional 2nd argument $f$,
a boolean function to be used for comparing two
elements of the array.

\item

Write a Maple procedure which computes the Fibonacci polynomials $F_n(x)$.
These polynomials satisfy the linear recurrence
$F_0(x) = 1$, $F_1(x) = x$, and $F_n(x) = x F_{n-1} (x) + F_{n-2} (x)$.
Compute and factor the first 10 Fibonacci polynomials.
Can your program compute $F_{50}(x)$?

\item
Write a Maple procedure called \verb+structure+ that when given a Maple
expression returns  the structure of the expression as a tree, namely

\begin{itemize}
\item
for atomic objects (integers and strings) just return the object
\item
otherwise, return a list with the type of the object as the first element,
and the remaining elements the structure of the operands
of the object.
\end{itemize}

For example: \verb|structure(x^3*sin(x)/cos(x));| should return
\begin{verbatim}
    [*, [^, x, 3], [function, sin, x], [^, [function, cos, x], -1]]
\end{verbatim}

Your routine should handle the Maple types
        integer, fraction, float, `\verb|+|`, `\verb|*|`, `\verb|^|`,
        string, indexed,
        function, range, equation, set, list, series and table.
Test your function on the following input

\begin{verbatim}
   Int( exp(-t)*sqrt(t)*ln(t), t=0..infinity ) =
   int( exp(-t)*sqrt(t)*ln(t), t=0..infinity );
\end{verbatim}

\item 
Extend the \verb+DIFF+ procedure to differentiate general powers,
the functions $\ln, \exp, \sin, \cos$
and to return unevaluated instead of giving an error when
it does not know how to differentiate the given input.
Extend \verb+DIFF+ to handle repeated derivatives, i.e. how to
\verb+DIFF+ the function \verb+DIFF+.
For example, \verb+DIFF(DIFF(f(x,y),x),y)+ should output the same
result as \verb+DIFF(DIFF(f(x,y),y),x)+.

\item

Write a routine {\tt comb} that when given a set of values and a non-negative integer
$n$ returns a list of all combinations of the values of size $n$.
For example

\begin{maple}
> comb( {a,b,c,d}, 3 );

                    { {a,b,c}, {a,b,d}, {a,c,d}, {b,c,d} } 
\end{maple}


\noindent
Modify your routine to work for lists which may contain duplicates, e.g.

\begin{maple}
> comb( [a,b,b,c], 2 );

                         [[a,b], [a,c], [b,b], [b,c]]
\end{maple}

\item
The Maple function \verb+degree+ computes the total degree of a polynomial
in a given variable(s).  For example, the polynomial

\begin{maple}
> p := x^3*y^3 + x^4 + y^5:
> degree(p,x); # degree in x

                                       4

> degree(p,{x,y}); # total degree in x and y

                                       6
\end{maple}

We showed how one could code the degree function in Maple with
our {\tt DEGREE} function.  However, the \verb+degree+ function
and our {\tt DEGREE} function only work for integer exponents.
Extend the {\tt DEGREE} so that it computes the degree of a general
expression, which may have rational or even symbolic exponents, e.g.

\begin{maple}
> f := x^(3/2) + x + x^(1/2);
                                    3/2        1/2
                              f := x    + x + x

> DEGREE(f,x);
                                      3/2

> h := (x^n + x^(n-1) + x^2) * y^m;

                                 n    (n - 1)    2   m
                          h := (x  + x        + x ) y

> DEGREE(h,x);
                                   max(n, 2)

> DEGREE(h,{x,y});
                                 max(n, 2) + m

\end{maple}

\item
Write Maple procedures to add and multiply univariate polynomials 
$a(x) = \sum _{i=0} ^m a_i x^i$ and $b(x) = \sum _{i=0}^n b_i x^i$
which are represented as

\begin{enumerate}
\item
a Maple list of coefficients i.e. $[a_n, a_{n-1}, \ldots, a_1, a_0]$

\item
a linked list of coefficients i.e.
$[a_n, [a_{n-1}, \ldots, [a_1, [a_0, NIL]] \ldots ]]$.
\end{enumerate}

\item
The data structures in the above exercise are inefficient for sparse
polynomials because zeros are represented explicitly.  E.g. consider
the polynomial $x^{10}+x+1$.  Represented as a Maple list it would be

\begin{center}
{\tt [10,0,0,0,0,0,0,0,1,1]}
\end{center}

\noindent
and as a linked list it is

\begin{center}
{\tt [10,[0,[0,[0,[0,[0,[0,[0,[0,[1,[1,NIL]]]]]]]]]]]}.
\end{center}

\noindent
A more efficient representation for sparse polynomials is to store
just the non-zero terms.  Let us store the $k$'th non-zero
term $a_k x^{e_k}$ as $[a_k,e_k]$ where $e_k > e_{k+1} \ge 0$.
Now write procedures to add and multiply sparse univariate polynomials
which are represented as

\begin{enumerate}
\item
a Maple list of non-zero terms

\item
a linked list of non-zero terms
\end{enumerate}

\item
The Gaussian integers ${\bf Z}[i]$ are the set of complex numbers with
integer coefficients,
i.e. numbers of the form $a + b i$ where $a,b \in {\bf Z}$ and $i = \sqrt {-1}$.
What is most interesting about the Gaussian integers is that
they form a Euclidean domain, i.e. they can be factored into primes,
and you can use the Euclidean algorithm to compute
the GCD of two Gaussian integers.
Given a Gaussian integer $x = a + b i$ let
$||x|| = a^2 + b^2$ be the norm of $x$.
Write a Maple procedure \verb+REM+ that computes the remainder $r$ of two
non-zero Gausian integers $x,y$ such that $x = q y + r$ where $||r|| < ||y||$.
Use your remainder routine to compute the GCD of two Gaussian integers.

\item

Write a Maple n-ary procedure \verb+GCD+ that does integer and
symbolic simplification of GCD calculations in ${\bf Z}$.
E.g. on input of \verb+GCD(b,GCD(b,a),-a)+,
your procedure should return \verb+GCD(a,b)+.
For integer inputs, your routine should compute the GCD directly.
For symbolic inputs, your \verb+GCD+ procedure should know the following

\begin{enumerate}
\item
GCD(a,b) = GCD(b,a) (GCD's are commutative)
\item
GCD(GCD(a,b),c) = GCD(a,GCD(b,c)) = GCD(a,b,c) (GCD's are associative)
\item
GCD(0,a) = GCD(a) (the identity is 0)
\item
GCD(a) = GCD($-$a) = abs(a)
\end{enumerate}

\item

Write a Maple procedure \verb+monomial(v,n)+ which
computes a list of the monomials of total degree $n$ in the
list of variables $[v_1, v_2, \ldots , v_m]$.
For example, \verb+monomial([u,v],3)+ would
return the list $[ u^3, u^2 v, u v^2, v^3 ]$.
Hint: consider the Taylor series expansion of the product

\[
\prod _{i=1} ^m \frac{1}{1-v_i t}
\]

in $t$ to order $t^n$ -- see Maple's \verb+taylor+ command.

\item

Given a sequence of points $(x_0, y_0), (x_1, y_1), \ldots , (x_n, y_n)$
where the $x_i$'s are distinct, the Maple library function \verb+interp+
computes a polynomial of degree $\le n$ which interpolates the points
using the Newton interpolation algorithm.  Write a Maple procedure which
computes a natural cubic spline for the points.  A natural cubic spline is a
piecewise cubic polynomial where each interval $x_i \le x \le x_{i+1}$ is
defined by the cubic polynomial

\[
f_i(x) ~=~ a_i x^3 + b_i x^2 + c_i x + d_i  ~~~~ 1 \le i \le n
\]

where the $4 n$ unknown coefficients are uniquely determined by the
following $4 n$ conditions

\[
\begin{array}{r c l r c l l l }
f_i(x_{i-1}) & = & y_{i-1}, & f_i(x_i) & = & y_i, & & i = 1 \ldots n \\
f'_i(x_i) & = & f'_{i+1}(x_{i-1}), & f''_i(x_i) & = & f''_{i+1}(x_{i-1}), & & i = 1 \ldots n-1 \\
f''_1(x_0) & = & 0, & f''_n(x_n) & = & 0
\end{array}
\]

These conditions mean that the resulting piecewise polynomial is $C^2$
continuous.
Write a Maple procedure which on input of the points as a list of
$(x_i, y_i)$'s in the form \verb+[x0, y0, x1, y1, ..., xn, yn]+,
and a variable, outputs a list of the segment
polynomials \verb+[f1, f2, ..., fn]+.
This requires that you create the segment polynomials with unknown
coefficients, use Maple to compute the derivatives and solve the
resulting equations.
For example

\begin{maple}
> spline([0,1,2,3],[0,1,1,2],x);

          3               3      2                      3      2
  [- 1/3 x  + 4/3 x, 2/3 x  - 3 x  + 13/3 x - 1, - 1/3 x  + 3 x  - 23/3 x + 7]
\end{maple}

\item

Design a data structure for representing a piecewise function.
The representation for the piecewise cubic spline above does not interface with
other facilities in Maple.  Let us represent a piecewise function
as an unevaluated function instead as follows.  Let

\[
{\rm IF}( c_1, f_1, ..., c_{n-1}, f_{n-1}, f_n )
\]

\noindent
mean

\begin{center}
\begin{tabular}{l l}
if $c_1$	& then $f_1$ \\
if $c_2$	& then $f_2$ \\
$\ldots$ \\
        & else $f_n$
\end{tabular}
\end{center}

\begin{itemize}
\item
Write a Maple procedure called \verb+IF+ that simplifies
such a piecewise function for conditions which are purely numerical.
\item
Write a procedure called `evalf/IF` which evaluates an IF expression
in floating point, hence allowing IF expressions to be plotted.
\item
Write a procedure called `diff/IF` which differentiates an IF expression.
\end{itemize}

E.g. your procedures should result in the following behaviour.

\begin{maple}
> IF( x<0, sin(x), cos(x) );

                           IF(x < 0, sin(x), cos(x))
> diff(",x);
                          IF(x < 0, cos(x), - sin(x))

> IF( Pi/3<0, sin(Pi/3), cos(Pi/3) );
                                             1/2
                         IF(1/3 Pi < 0, 1/2 3   , 1/2)
> evalf(");
                                  .5000000000
\end{maple}

Modify your \verb+IF+ procedure to simplify nested IF expressions,
and to handle constant conditions like the example above.
E.g. your IF procedure should result in the following behaviour.

\begin{maple}
> IF( x<0, 0, IF( x<1, x^2, IF(x>2, 0, 1-x^2) ) );

                                        2                 2
                   IF(x < 0, 0, x < 1, x , 2 < x, 0, 1 - x )

> IF( Pi/3<0, sin(Pi/3), cos(Pi/3) );

                                      1/2
\end{maple}

\item
The following iteration is known as the Halley iteration

\[
x _ {k+1} = x_k - f(x_k)/f'(x_k) / \left( 1 - \frac{f(x_k) f''(x_k)}{2 f'(x_k)^2} \right) ~ .
\]

Program this iteration in Maple and check that it converges cubically.
Prove that the convergence is cubic.

\item

The command \verb+dsolve+ solves ordinary differential equations
analytically, although not all ODE's can be solved in closed form.
But it is always possible to obtain a series solution which may be useful.
We are given the ODE

\[
y'(x) = f(x,y(x))
\]

and the initial condition $y(0) = y_0$ and we want to find the first
few terms in the series

\[
y(x) = \sum _{k=0} ^ \infty y_k x^k ~ .
\]

Write a Maple procedure which on input of $f(x,y(x))$ and the initial
condition $y_0$ constructs a linear system of equations to solve.
I.e. let

\[
y(x) = y_0 + \sum _{k=1} ^n y_i x^i
\]

substitute this finite sum into the ODE, equate coefficients and solve
for the unknowns $y_i$.  Note, you will want to use the \verb+taylor+ command
to truncate the result to order $n$.  Test your Maple procedure on
the following ODE

\[
y'(x) = 1 + x^2 - 2 x y(x) + y(x)^2, ~~ y(0) = 1
\]

You should get the following series

\[
y(x)=1+2\,x+x^{2}+x^{3}+x^{4}+x^{5}+O(x^{6})
\]

Compute the solution analytically using Maple's \verb+dsolve+ command
and check that your series solution agrees with Maple's analytic solution.

Develop a Newton iteration to solve for the series $y(x)$ which
converges quadratically.  I.e. the $k$'th approximation

\[
x _ {k+1} = x_k - f(x_k)/f'(x_k)
\]

is accurate to $O(x^{2^k})$.  Thus the iteration starts with $x_0 = y_0$
and computes $x_1 = y_0 + y_1 x$ then $x_2 = y_0 + y_1 x + y_2 x^2 + y_3 x^3$
and so on.  Write a Maple procedure to compute the series.

\item
Modify the GaussianElimination procedure so that it does a complete
row reduction of a rectangular matrix of rational numbers
reducing the matrix to Reduced Row Echelon Form.
Now modify your procedure to work with a matrix of formulae.
Use the \verb+simplify+ function to simplify the intermediate results.

\item
The use of the Taylor series for computing ${\rm erf}(x)$
\[
\frac{\sqrt {\pi}}{2} {\rm erf}(x) ~=~
\sum _{n=0} ^ \infty \frac{ (-1)^n x^{2 n + 1} } {n! (2 n + 1) }
\]
becomes increasingly inefficient for large $x$ because it converges slowly.
It also becomes increasingly inaccurate because of cancellation of positive
and negative terms.  Write a Maple procedure which sums the terms
of the asymptotic series for ${\rm erf}(x)$
\[
\frac{\sqrt {\pi}}{2} {\rm erf}(x) ~=~
1 - \frac{ 2 \sqrt {\pi} }{ x } e ^ {- 2 x }
\sum _{n=1} ^ \infty \frac{ (-1)^n {1 \times 3 \times \ldots \times {2 n - 1}}}
{2^n x^{2 n} }
\]
for large $x$ until it converges.  How large must $x$ be before
this series converges at \verb+Digits+ decimal digits of precision?
Combine the use of both the Taylor series and the asymptotic series
to write a procedure that computes ${\rm erf}(x)$ for all $x$ accurate
to \verb+Digits+ decimal digits of precision.

\end{enumerate}
\end{document}
