Copyright (C) 1997 Rich Neidinger, John Swallow, and Todd Will.  Free for
distribution to college and university instructors for personal,
non-commercial use.  If these notebooks are used in a course, the authors
request $20 per student.

17. Recursion II:  Dynamic Programming

Last revision:  September 25 1996

In this section we continue our study of recursion, presenting two more
examples, and discussing the relationships between recursion and dynamic
programming, a method which may increase the speed of our program at the
expense of using additional space.

Another Recursive Example:  The Euclidean Algorithm

Recall that the greatest common divisor of two integers "a" and "b" is the
largest (positive) integer which divides each evenly (with remainder 0).
Mathematica gives us a way to compute greatest common divisors (GCD's), the
"GCD[ ]" command.

GCD[72,30]

Mathematica also provides a remainder function which divides the first
parameter by the second and returns the remainder:  the "Mod[ ]" function.
(The "Mod[ ]" function is useful in determining whether one integer is
divisible by another.  For example, "Mod[a,3]==0" tests whether or not "a"
is evenly divisible by 3.)

Mod[72,30]

And, for completeness, we also have the integer quotient function,
"Quotient[ ]", which provides the quotient, which is the number of times
the first parameter divides the second, ignoring the remainder.

Quotient[72,30]

The purpose of this section is to use recursion to program our own GCD function.

A famous theorem from arithmetic implies recursion

A famous mathematical theorem states that GCD[a,b] = GCD[b,r], where r =
Mod[a,b].  What is interesting about the theorem from a recursive
standpoint is that it provides a way to "simplify" the two parameters, so
that the function can be invoked again.  How?  First, notice that the "Mod[
]" function returns a remainder, an integer between 0 and the value of the
second argument.  Now suppose that "b" is less than "a".  Then "b,r", the
new pair, is a pair of integers both less than or equal to "b", so that the
maximum of the pair has decreased, from "a" to "b".  Suppose, on the other
hand, that "a" is less than "b".  Then the new pair "b,r" is equal to
"b,a", since "a=Mod[a,b]" if "b" is greater than "a".  However, now the
pair has the first element greater than the second element, so that the
first case applies, and the next pair will have a smaller maximum value.
Finally, suppose that "a" equals "b".  Then the next pair is "b,0", and we
know that the GCD of any number and 0 is the number itself, which is in
this case "b".  Putting these facts together, we see that repeatedly
applying the rule "GCD[a,b]=GCD[b,Mod[a,b]]" reduces the maximum value of
the pair of parameters to GCD.  Since the two parameters are always
nonnegative integers, at some point we much reach a scenario in which one
parameter is zero, and the GCD at that point is easily found.  In other
words, we have not only a way to simplify parameters in order to invoke the
GCD function again, but a base case:  "GCD[b,0]=b" for any nonnegative
integer "b".  As an example, consider the following sequence of parameter
replacements:  GCD[15,6] = GCD[6,3] = GCD[3,0] = 3.

Now we are in a position to write a recursive GCD function:

myGCD[a_Integer?Positive, 0] := a
myGCD[a_Integer?Positive, b_Integer?Positive] :=
   myGCD[b, Mod[a,b]]

myGCD[15,6]

Recall that the "Trace[ ]" command can be useful with recursive functions,
and that executing "Trace[ ]" with the name of a function as the second
argument causes "Trace[ ]" to announce only each call of the function.

Trace[myGCD[15,6],myGCD]

myGCD[72,30]

Dynamic Programming; Making Recursion Run Faster

Fibonacci numbers begin as 1,1, and then continue as the sum of the
preceding two numbers:  2, 3, 5, 8, 13, 21, 34, and so on.  A simple
recursive implementation of a Fibonacci number-finder would be the
following:

Clear[fib]
fib[1]=1;
fib[2]=1;
fib[n_]:=fib[n-1]+fib[n-2]

We can find the 7th Fibonacci number, then by executing "fib[7]":

fib[7]

There are some important time issues to be considered with recursive
programs of this sort.  To highlight these, let's ask "Trace[ ]" to tell us
how many function calls to "fib[ ]" it makes when it considers "fib[7]".

fibUse = Trace[fib[7],fib[_Integer]]

Sort[Flatten[fibUse]]

We see that "fib[7]" and "fib[6]" are each computed once, that "fib[5]" is
computed twice, "fib[4]" three times, "fib[3]" five times, and "fib[2]"
eight times.   A lot of work goes into repeatedly computing "fib[5]"!
There is an obvious improvement here:  we should save the value of "fib[n]"
once we find it. We can do this by modifying our current definition only
slightly, in the next input cell.

The change is so small it might be hard to see.  In the last line of the
cell below we are not only computing "fib2[n]" but are also setting up a
special case of "fib2[n]", as follows.  By using an assignment statement
("fib2[n]=..."), we are instructing Mathematica that whenever "fib2[n]" is
executed with the same value again, it should use the value we've specified
instead of using the recursive definition of "fib2[ ]".  It will be just as
if we had specified the special value as we did with "fib2[1]=1".

Clear[fib2]
fib2[1]=1;
fib2[2]=1;
fib2[n_]:= fib2[n] = fib2[n-1]+fib2[n-2]

For example, when we execute "fib2[7]", "fib2[7]" is computed and "fib2[7]"
is set equal to the computed value, and the next time we request the value
of "fib2[7]", the value is not recomputed using recursion, but instead is
simply looked up.

fib2[7]

?fib2

Now let's return to the time issue.  We'll clear each of the two variants
of our Fibonacci-number function, redefine them, and time how long it takes
to compute the 20th Fibonacci number each way.

Clear[fib]
fib[1]=1;
fib[2]=1;
fib[n_]:=fib[n-1]+fib[n-2]

Clear[fib2]
fib2[1]=1;
fib2[2]=1;
fib2[n_]:= fib2[n] = fib2[n-1]+fib2[n-2]

Timing[fib[20]]

Timing[fib2[20]]

Fairly dramatic.  Even more dramatic are the next calls:

Timing[fib[20]]

Timing[fib2[20]]

Dynamic Programming for the Euclidean Algorithm?

The section above provided us a first look at using dynamic programming to
solve a programming problem, which is a method whereby the program (a)
saves certain results for recall as they are computed, but (b) does not
simply compute all of the possible results ahead of time.  Dynamic
programming is a via media between two extremes, one in which, say, the
first 1000 Fibonacci numbers are computed ahead of time and then are looked
up when requested, the other in which no values are saved as they are
computed but are always recomputed when requested.  In the first extreme,
we find a large use of time at the beginning of the program, in order to
compute all of the possible values first, and a large use of space (in the
computer's memory) to store all of the possible results.  At the other
extreme, we finda small space burden but a probable chance that time will
be wasted by recomputing results already computed before. Dynamic
programming uses time judiciously, using memory as needed to store
requested results for later, but not (pre-)computing any results which we
may not ever need.

Is dynamic programming an appropriate method for implementing the Euclidean
Algorithm?  In other words, would implementing such a method balance time
and space constraints?  We first realize that under special circumstances
the choice may be clear.  For instance, if the computer in question has
very little memory or if we are entirely satisfied with the amount of time
the computer takes to compute the GCD of any two numbers we might be
interested in, we would not want to use dynamic programming.  Can we answer
the question more generally, however?  The difference between our recursive
code for the Euclidean Algorithm and our recursive code for Fibonacci
numbers lies in the fact that the definition of "fib[n]" forces the
computation (or recall) of every "fib[x]" for integers x < n, while
"myGCD[a,b]" forces the computation (or recall) of only certain values of
the parameters of "myGCD[ ]".  Put another way, computing "fib[15]" will
always benefit from knowing the values which "fib[8]" was forced to
compute, while computing "myGCD[10,15]" does not use any of the values
found during the computation of "myGCD[8,16]".  We might conclude, roughly,
that if we will be computing many, many GCDs of integers whose values fall
in a relatively small range (say 1 to 1000), using dynamic programming
would prove a good idea, since we would need to store at most 1000 * 1000 =
1,000,000 integers, taking up 1 megabyte or less, and the frequency of the
computations would insure that at some point the algorithm would begin
recalling, instead of computing, values of "myGCD[ ]".  If the range of the
parameters is very large, or if we will be computing relatively few
instances of "myGCD[ ]", then dynamic programming may have very little
benefit, either forcing us to set aside a large amount of memory, as in the
first case, or storing values which may never be recalled again, as in the
second case.