########################################################################
#  Basic Data Structures
########################################################################

X := [1.3, 5.3, 11.2, 2.1, 2.1];

nops(X);

X[2];

########################################################################
#  As stated previously, Maple V.3 does not have the 'add' command.
#  To remedy this shortcoming, we define an 'add' procedure which 
#  adds the numbers in a list.
########################################################################

add := proc(l)
   local  i,  total;
   if not type(l, 'list') then
      ERROR(`input must be a list`);
   fi;
   total := 0;
   for i from 1 to nops(l) do 
      total := total + l[i];
   od;
   total;
end:

add(X);


########################################################################
#  Now we can continue with the 'average' example.
########################################################################

average := proc(X)
   local  n, total;
   if not type(X, 'list') then
       ERROR(`input must be a list`);
   fi;
   n := nops(X);
   if n = 0 then ERROR(`empty list`) fi;
   total := add(X);
   total / n;
end:

average(X);

average( [ a, b, c ] );

########################################################################
#  Extracting a sequence from a list.
########################################################################

Y := X[];

########################################################################
#  Selecting elements and ranges of elements from a sequence
########################################################################
Y[3];

Y[2..4];

########################################################################
#  Note that a sequence of sequences is automatically "flattened"
#  into a single sequence ... Also note that in Maple V.3 'W'
#  is protected (stands for Lambert's W function), so we need to 
#  use a different variable name:
########################################################################

WW := a, b, c;

Y, WW, Y;

########################################################################
#  ... whereas a list of lists remains a list of lists.
########################################################################

[ X, [a, b, c], X ];

########################################################################
#  Defining a set by enclosing a sequence in braces.  
########################################################################

Z := { Y };

########################################################################
#  Recall that a Maple set, as in mathematics, is an unordered collection 
#  of distinct objects.  Thus, whereas the sequence Y has 5 elements,
#  the set Z has only four.
########################################################################

nops(Z);

########################################################################
#  Examples of the use of 'seq' to construct sequences, lists and 
#  sets.
########################################################################

seq( i^2 , i = 1..5 );

unassign(`f`);
seq( f(i) , i = X );

[ seq( { seq(i^j , j =1..3) }, i=-2..2 ) ];

########################################################################
#  Creating a sequence using a loop: NULL is the empty sequence.
########################################################################

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

########################################################################
#  A MEMBER procedure: once again, recall that we must do explicit
#  type checking of the procedure arguments.  Note that the first 
#  argument can be anything, so there is no need to check its type.
########################################################################

MEMBER := proc( a, L )
   local i;
   if not (type(L,'list') or type(L,'set')) then 
      ERROR(`second argument must be a list or set`);
   fi;
   for i from 1 to nops(L) do
      if a = L[i] then RETURN(true) fi;
   od;
   false;
end:

MEMBER( 3, [1,2,3,4,5,6] );

########################################################################
#  MEMBER rewritten to use the special form of the for loop which 
#  automatically iterates over elements of set or list.
########################################################################

MEMBER := proc( a, L )
   local i;
   if not (type(L,'list') or type(L,'set')) then
      ERROR(`second argument must be a list or set`);
   fi;
   for i in L do
      if a = L[i] then RETURN(true) fi;
   od;
   false;
end:

MEMBER( x, {1, 2, 3, 4} );
MEMBER( 1, {1, 2, 3, 4} );


########################################################################
#  Binary Search
########################################################################

########################################################################
#  First, a brute-force search
########################################################################

Search := proc(Dictionary, w)
   local x;
   if not type(Dictionary,'list') then
      ERROR(`first argument must be a list`);
   fi;
   for x in Dictionary do
      if x = w then RETURN(true) fi;
   od;
   false;
end:

Dictionary := [ induna, ion, logarithm, meld ];

Search(Dictionary,ion);
Search(Dictionary,anion);

########################################################################
#  The Binary-Search algorithm assumes that entires in the Dictionary 
#  are in ascending lexicographical order.
# 
#  Observe how we can check that the dictionary is of type 
#  'list of strings'.
########################################################################

BinarySearch := proc(D, w, s, f)
   local m;
   if not type(D,'list(string)') then 
      ERROR(`first argument must be a list of strings`) 
   fi;
   if not type(w,'string') then ERROR(`second argument must be a string`) fi;
   if not type(s,'integer') then ERROR(`third argument must be an integer`) fi;
   if not type(f,'integer') then ERROR(`fourth argument must be an integer`) fi;
   if s > f then RETURN(false) fi;     # entry was not found
   m := iquo(s+f+1, 2);                # midpoint of D
   if w = D[m] then
      true;
   elif lexorder(w, D[m]) then     
      BinarySearch(D, w, s, m-1);
   else
      BinarySearch(D, w, m+1, f);
   fi;
end:

BinarySearch( Dictionary, hedgehogs, 1, nops(Dictionary) );

BinarySearch( Dictionary, logarithm, 1, nops(Dictionary) );

BinarySearch( Dictionary, meoldious, 1, nops(Dictionary) );

########################################################################
#  Plotting the Roots of a Polynomial
########################################################################

plot( [ [0,0], [1,2], [-1,2]], style=point,color=black);

y := x^3 - 1;

########################################################################
#  Use 'fsolve' to find roots numerically.  Type '?fsolve' for more 
#  information.
########################################################################

R := [ fsolve(y = 0, x, complex) ];

########################################################################
#  Converting list of complex numbers into a list of points in the 
#  plane. 
########################################################################

points := map( z -> [Re(z), Im(z)], R );

plot( points, style=point );

########################################################################
#  Input to rootplot should be polynomial in x with constant 
#  coefficients; again, we can use 'type' to see whether 'p'
#  is of this type.
########################################################################

rootplot := proc(p)
   local R, points;
   if not type(p,'polynom(constant,x)') then
      ERROR(`input must be polynomial in x with constant coefficients`);
   fi;
   R := [ fsolve(p, x, complex) ];
   points := map( z -> [Re(z), Im(z)], R );
   plot( points, style = point );
end:

rootplot( x^6 + 3*x^5 + 5*x + 10 );

########################################################################
#  Generate a random polynomial
########################################################################

y := randpoly(x, degree=100);

rootplot(y);

########################################################################
#  Computing with Formulae.
########################################################################

########################################################################
#  The Height of a Polynomial
########################################################################

HGHT := proc(p, x) 
   local   i,  c,  height;
   if not type(p,'polynom') then
      ERROR(`first argument must be a polynomial`);
   fi;
   if not type(x,'name') then
      ERROR(`second argument must be a name`);
   fi;
   height := 0;
   for i from 0 to degree(p, x) do
      c := coeff(p, x, i);
      height := max(height, abs(c));
   od;
   height;
end:

p := 32*x^6 - 48*x^4 + 18*x^2 - 1;

HGHT(p,x);

coeffs(p, x);

S := map( abs, {"} );

max( S[] );

HGHT := proc(p, x) 
   local   S;
   if not type(p,'polynom') then
      ERROR(`first argument must be a polynomial`);
   fi;
   if not type(x,'name') then
      ERROR(`second argument must be a name`);
   fi;
   S := { coeffs(p, x) };
   S := map( abs, S );
   max( S[] ):
end:

p := randpoly(x, degree=100);

HGHT(p, x);

map(f, p);

map(abs, p);

coeffs( " );

max( " );

p := randpoly(x , degree=50) * randpoly(x, degree=99);

max( coeffs( map(abs, expand(p)) ) );

########################################################################
#  The Chebyshev Polynomials, T_n(x).
########################################################################

T := proc(n, x)
   option remember;
   if not type(n,'nonnegint') then
      ERROR(`first argument must be a non-negative number`);
   fi;
   if not type(x,'name') then
      ERROR(`second argument must be a name`);
   fi;
   if   n = 0 then
      RETURN(1);
   elif n = 1 then
      RETURN(x);
   fi;
   2 * x * T(n-1, x) - T(n-2, x);
end:

T(4,x);

expand(");

########################################################################
#  Integration by parts.
########################################################################

int( x^2 * exp(x), x );

int( x^3 * arcsin(x), x);

int( u(x) * v(x), x ) = 
   u(x) * int( v(x), x ) - int( diff(u(x),x) * int(v(x), x), x );

diff(",x);

evalb(");

IntExpMonomial := proc(n, x)
   if not type(n,'nonnegint') then
      ERROR(`first argument must be a non-negative number`);
   fi;
   if not type(x,'name') then
      ERROR(`second argument must be a name`);
   fi;
   if n = 0 then RETURN( exp(x) ) fi;
   x^n * exp(x) - n * IntExpMonomial(n-1, x);
end;

IntExpMonomial(5, x);

collect(", exp(x));

########################################################################
#  Recall that we don't have the 'add' command, so we do the 
#  sum over the monomials using a loop.
#
#  Also, in Maple V.3, 'coeff' will not extract the coefficients 
#  of an unexpanded polynomial.
# 
#  Finally, note that the value computed in the text for
#  
#  IntExpPolynomial( (x^2 + 1)*(1 - 3*x) , x );
#
#  is incorrect.
########################################################################

IntExpPolynomial := proc(p, x)
   local   i, result,  pexpanded;
   if not type(p,'polynom') then
      ERROR(`first argument must be a polynomial`);
   fi;
   if not type(x,'name') then
      ERROR(`second argument must be a name`);
   fi;
   result := 0;
   pexpanded := expand(p);         # Expand input poly so that coeff works
   for i from 0 to degree(pexpanded,x) do
      result := result + coeff(pexpanded, x, i) * IntExpMonomial(i, x);
   od;
   collect(result, exp(x));
end:

IntExpPolynomial( (x^2 + 1)* (1 - 3*x), x );

########################################################################
#  COMPUTING WITH SYMBOLIC PARAMETERS
########################################################################

########################################################################
#  INPUT ERROR
########################################################################
# IntExpPolynomial( a*x^n, x );

IntExpPolynomial(x, x);
IntExpPolynomial(x^2, x);
IntExpPolynomial(x^3, x);

assume(n, integer);
additionally(n >= 0);

type(n,integer);

is(n,integer), is(n >=0);

IntExpMonomial := proc(n, x)
   local i;
   if not type(x,name) then 
      ERROR(`second argument must be a name`);
   fi;
   if is(n,integer) and is(n >= 0) then
      n! * exp(x) * sum( ( (-1)^(n-i)*x^i )/i!, i = 0..n );
   else
      ERROR(`Expected a non-negative integer but received `,n);
   fi;
end:

IntExpMonomial(4, x);