########################################################################
# Programming examples from Chapter 1 "Maple V Programming Guide",
# by Monagan et al.
#
# Note that Maple 6 and Maple V.5 use " " to delimit strings and ` ` to 
# delimit names; the Programming Guide (Maple V.4) uses ` ` to 
# delimit both.
########################################################################

########################################################################
# 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
########################################################################

y := x^3 - 2*x + 1;

yp := diff(y, x);

plot( [y, yp], x=-1..1 ); 

plotdiff := proc(y,x,a,b)
   local yp;
   yp := diff(y,x);
   plot( [y, yp], x=a..b );
end;

plotdiff( cos(x), x, 0, 2*Pi );

########################################################################
# 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(100);

add(n, n=1..100);

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

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

ABS(3); ABS(-2.3);

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;
   elif type(x, `*`) then
      map(ABS, x);
   else
      'ABS'(x)
   fi;
end;

ABS(a * b);

########################################################################
#  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 DOUBLE-QUOTES (") to delimit a string.
########################################################################
SUM("hello world");

########################################################################
#  Another definition of'SUM' with manual 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;

SUM("hello world");

########################################################################
#  Another definition of'SUM' with built-in type-checking.
########################################################################

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

SUM("hello world");

########################################################################
#   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::posint)
   local q;
   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::posint)
   local   a;
   a := 3 * n + 1;
   divideby2( a );
end;

checkconjecture := proc(x::posint)
   local count, n;
   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::nonnegint)
   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::nonnegint)
   option remember;
   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::nonnegint)
   local temp, fnew, fold, i;
   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::nonnegint)
   option remember;
   if n < 2 then
      RETURN(n);
   fi;
   Fibonacci(n-1) + Fibonacci(n-2);
end;

Fibonacci(10);

Fibonacci(-1);
   
########################################################################
#  Basic Data Structures
########################################################################

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

nops(X);

X[2];

add( i, i = X);

average := proc(X::list)
   local  n, total;
   n := nops(X);
   if n = 0 then ERROR("empty list") fi;
   total := add(i, i=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 ... 
########################################################################

W := a, b, c;

Y, W, 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::anything, L::{list,set} )
   local i;
   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::anything, L::{list,set} )
   local i;
   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::list, w::anything)
   local x;
   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.
#
#  Note that use of 'name' types rather than 'string' as in Programming
#  Guide.
########################################################################

BinarySearch := proc(D::list(name), w::name, s::integer, f::integer)
   local m;
   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::polynom(constant,x))
   local R, points;
   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::polynom, x::name) 
   local   i,  c,  height;
   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::polynom, x::name) 
   local   S;
   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::nonnegint, x::name)
   option remember;
   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::nonnegint, x::name)
   if n = 0 then RETURN( exp(x) ) fi;
   x^n * exp(x) - n * IntExpMonomial(n-1, x);
end;

IntExpMonomial(5, x);

collect(%, exp(x));

########################################################################
#  Note that the value computed in the text for
#  
#  IntExpPolynomial( (x^2 + 1)*(1 - 3*x) , x );
#
#  is incorrect.
########################################################################

IntExpPolynomial := proc(p::polynom, x::name)
   local   i, result;
   result := add( coeff(p, x, i)*IntExpMonomial(i, x),
                  i = 0..degree(p, x) );
   collect(result, exp(x));
end;

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

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

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::anything, x::name)
   local i;
   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);

IntExpMonomial(n, x);

diff(%,x);

simplify(%);