########################################################################
# Programming examples from Chapter 1 "Maple V Programming Guide",
# by Monagan et al.
########################################################################

########################################################################
#  INPUT ERROR marks statements which will cause 'ERROR' to be 
#  called and thus halt 'read <inputfile>' in GUI.
########################################################################

########################################################################
# Getting Started
# 
# By default, the value returned by a procedure is the value of 
# the last expression computed by the procedure.
########################################################################

a := 103993 / 33102;

half := proc(x)
   evalf(x/2);
end;

half(2/3);

half(a);

half(1) + half(2);

f := proc()
   a := 103993 / 33102;
   evalf(a/2);
end;

f();

########################################################################
# Illustrating use of ':' to supress echo of procedure definition
########################################################################

g := proc()
   a := 103993 / 33102;
   evalf(a/2);
end;

g;

eval(g);

########################################################################
# Locals and Globals
########################################################################

b := 2;

g := proc()
   local b;
   b := 103993 / 33102;
   evalf(b/2);
end;

h := proc()
   global b;
   b := 103993 / 33102;
   evalf(b/2);
end;

g();

b;

h();

b;

########################################################################
# Inputs, Parameters, Arguments
########################################################################

k := proc(p,q)
   p/q;
end;

k(103993,33102);

k( 23, 0.56);

(2 + 3*I)^2;

k(2 + 3*I,");

k(1.362, 5*I);

abnorm := proc(a,b)
   sqrt(a^2 + b^2);
end;

abnorm(2, 3);

znorm := proc(z)
   sqrt( Re(z)^2 + Im(z)^2 );
end;

znorm( 2 + 3*I );


abznorm := proc(z)
   local r, i;
   r := Re(z);
   i := Im(z);
   abnorm(r, i);
end;

abznorm( 2 + 3*I );

abznorm( x + y*I );

########################################################################
# Basic Programming Constructs
########################################################################

########################################################################
# The Assignment Statement
#
# THIS EXAMPLE OMITTED since 'plot' in Maple V.3 does not support 
# multiple-plot feature ...
########################################################################

########################################################################
# The for Loop
########################################################################

########################################################################
# A simple example of a 'for' loop.
########################################################################

total := 0;
for i from 1 to 5 do
   total := total + 1;
od;

########################################################################
# A procedure for computing the sum of the first n natural numbers 
# using a 'for' loop.
########################################################################

SUM := proc(n)
   local i, total;
   total := 0;
   for i from 1 to n do
      total := total + i;
   od;
#-----------------------------------------------------------------------
#  Evaluate 'total' so that the procedure will return its value.
#-----------------------------------------------------------------------
   total;
end;

SUM(1000);

########################################################################
# NOTE: Maple V.3 does not support 'add' command
########################################################################

########################################################################
# The Conditional Statement
########################################################################

ABS := proc(x)
   if x < 0 then
      -x;
   else
       x;
   fi;
end;

ABS(3); ABS(-2.3);

########################################################################
#  INPUT ERROR
########################################################################
# ABS(q);

########################################################################
# This version checks that its argument is numeric; if it isn't the 
# procedure returns UNEVALUATED.  Note that the quotes used are 
# (regular) single quotes ('), not backquotes (`) which are used 
# to delimit strings.
########################################################################

ABS := proc(x)
   if( type(x,numeric) ) then
      if x < 0 then -x else x fi;
   else
      'ABS'(x);
   fi;
end;

ABS(q);

########################################################################
# Use of nested if-statements.
########################################################################

HAT := proc(x)
   if type(x, numeric) then
      if x <= 0 then
         0;
      else
         if x <= 1 then
            x;
         else
            if x <= 2 then
               2 - x;
            else
               0;
            fi;
         fi;
      fi;
   else
      'HAT'(x);
   fi;
end;

HAT(-1);
HAT(1/2);
HAT(3/2);
HAT(3);
HAT(q);

########################################################################
# Another version of HAT which uses 'elif' clauses.
########################################################################

HAT := proc(x)
   if type(x, numeric) then
      if x <= 0 then 0;
      elif x <= 1 then x;
      elif x <= 2 then 2 - x;
      else 0;
      fi;
   else
      'HAT'(x);
   fi;
end;

HAT(-1);
HAT(1/2);
HAT(3/2);
HAT(3);
HAT(q);

########################################################################
# Symbolic Transformations.  Note the use of the 'unassign' command
# to unassign the values of global variables.
########################################################################

unassign('a');
unassign('b');

ABS(a * b);

########################################################################
# Version of ABS which 'maps over products'.  Note that for products
# this procedure calls itself recursively.
########################################################################

ABS := proc(x)
   if type(x, numeric) then
      if x < 0 then -x else x fi;
#-----------------------------------------------------------------------
#  Note:  In Maple V.3, multiplication operator (*) must be 
#  BACK-QUOTED (`*`) NOT FORWARD-QUOTED ('*').
#-----------------------------------------------------------------------
   elif type(x, `*`) then
      map(ABS, x);
   else
      'ABS'(x)
   fi;
end;

########################################################################
#  Type Checking:  This is an area of major difference between
#  Maple V.3 and V.4.  Maple V.3 DOES NOT SUPPORT THE '::<type>'
#  construct in function definitions for automatic type checking 
#  of arguments.  You must use the 'type' function to do explicit 
#  type-checking.
########################################################################

########################################################################
#  Define 'SUM' as previously: no type-checking.
########################################################################

SUM := proc(n)
   local i, total;
   total := 0;
   for i from 1 to n do
      total := total + i;
   od;
   total;
end;

########################################################################
#  Note use of BACKQOUTES (`) to delimit a string.
########################################################################
########################################################################
#  INPUT ERROR
########################################################################
# SUM(`hello world`);

########################################################################
#  Another definition of'SUM' with type-checking.
#
#  Note the use of the maple function, 'ERROR' to print an error 
#  message and exit the procedure immediately if a bad argument
#  is detected
########################################################################

SUM := proc(n)
   local i, total;
   if not type(n, integer) then
      ERROR(`input must be an integer`);
   fi;
   total := 0;
   for i from 1 to n do total := total + i od;
   total;
end;

########################################################################
#  INPUT ERROR
########################################################################
# SUM(`hello world`);

########################################################################
#  IN THE FOLLOWING EXAMPLES, KEEP IN MIND THAT OUR VERSION OF 
#  'maple' DOES NOT SUPPORT AUTOMATIC TYPE-CHECKING OF ARGUMENTS AND
#  THUS PROCEDURES HAVE BEEN MODIFIED TO INCLUDE EXPLICIT TYPE-CHECKING
########################################################################

########################################################################
#   The while Loop
########################################################################

########################################################################
#   The 'iquo' and 'irem' commands calculate the quotient and 
#   remainder respectively using integer division.
########################################################################
iquo(8,3);

irem(8,3);

divideby2 := proc(n)
   local q;
   if not type(n, 'posint' ) then 
      ERROR(`input must be a positive integer`);
   fi;
   q := n;
   while irem(q, 2) = 0 do
      q := iquo(q, 2);
   od;
   q;
end;

divideby2(32);

divideby2(48);

########################################################################
# Modularization
########################################################################

40;

divideby2( " );

3 * " + 1;

divideby2( " );

3 * "  + 1;

divideby2( " );

iteration := proc(n)
   local   a;
   if not type(n, 'posint' ) then
      ERROR(`input must be a positive integer`);
   fi;
   a := 3 * n + 1;
   divideby2( a );
end;

checkconjecture := proc(x)
   local count, n;
   if not type(x, 'posint' ) then
      ERROR(`input must be a positive integer`);
   fi;
   count := 0;
   n := divideby2(x);
   while n > 1 do
      n := iteration(n);
      count := count + 1;
   od;
   count;
end;

checkconjecture( 40 );

checkconjecture( 4387 );

########################################################################
# Recursive Procedures 
########################################################################

Fibonacci := proc(n)
   if not type(n, 'nonnegint' ) then
      ERROR(`input must be a non-negative integer`);
   fi;
   if n < 2 then
      n;
   else 
      Fibonacci(n-1) + Fibonacci(n-2);
   fi;
end;

seq( Fibonacci(i), i = 0..15 );

Fibonacci(10);

time( Fibonacci(20) );


########################################################################
# Version of 'Fibonacci' which remembers previously computed values
########################################################################

Fibonacci := proc(n)
   option remember;
   if not type(n, 'nonnegint' ) then
      ERROR(`input must be a non-negative integer`);
   fi;
   if n < 2 then
      n;
   else 
      Fibonacci(n-1) + Fibonacci(n-2);
   fi;
end;

Fibonacci(10);

time( Fibonacci(20) );

Fibonacci(1000);

time( Fibonacci(2000) );

########################################################################
# Iterative (looping) version of 'Fibonacci'
########################################################################

Fibonacci := proc(n)
   local temp, fnew, fold, i;
   if not type(n, 'nonnegint' ) then
      ERROR(`input must be a non-negative integer`);
   fi;
   if n < 2 then
      n;
   else
      fold := 0;
      fnew := 1;
      for i from 2 to n do
         temp := fnew + fold;
         fold := fnew;
         fnew := temp;
      od;
      fnew;
   fi;
end;

Fibonacci(10);

time( Fibonacci(20) );

Fibonacci(1000);

time( Fibonacci(2000) );

########################################################################
# The RETURN command
########################################################################

Fibonacci := proc(n)
   option remember;
   if not type(n, 'nonnegint' ) then
      ERROR(`input must be a non-negative integer`);
   fi;
   if n < 2 then
      RETURN(n);
   fi;
   Fibonacci(n-1) + Fibonacci(n-2);
end;

Fibonacci(10);

########################################################################
#  INPUT ERROR
########################################################################
# Fibonacci(-1);
   
########################################################################
#  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);