17. Recursion II:  Dynamic Programming

Last revision:  February 8 1999

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.

Mathematica operations/structures introduced: 
     recursive definitions
     GCD[]
     Mod[]
     Quotient[]
     Trace[]

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] equals GCD[b,r], where r is 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", from "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] equals 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] equals 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:

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

recurGCD[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[recurGCD[15,6],recurGCD]

recurGCD[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[0]=1;
fib[1]=1;
fib[n_]:=fib[n-1]+fib[n-2]

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

fib[6]

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[6]".

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

Sort[Flatten[fibUse]]

We see that "fib[6]" and "fib[5]" are each computed once, that "fib[4]" is computed twice, "fib[3]" three times, "fib[2]" five times, and "fib[1]" eight times.   A lot of work goes into repeatedly computing "fib[4]"!  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[0]=1".

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

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

fib2[6]

?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[0]=1;
fib[1]=1;
fib[n_]:=fib[n-1]+fib[n-2]

Clear[fib2]
fib2[0]=1;
fib2[1]=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]]

Wow.  Take a look at what Mathematica is storing for fib2.

??fib2

What if we try another option: a tail recursive Fibonacci function, fib3?

Clear[fib3]
fib3[
  n_Integer?((#>=0)&), 
  n0_Integer:1, 
  n1_Integer:1
]:= 
  Which[
    n == 0, n0,
    n == 1, n1,
    True,   fib3[n-1, n1, n0+n1]
  ]

Map[fib3, Range[0, 10]]

Timing[fib3[20]]

Hmmm, that's pretty fast too.  Take a look at what Mathematica is storing for fib3.

??fib3

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 compromise 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 find a 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 "recurGCD[a,b]" forces the computation (or recall) of only certain values of the parameters of "recurGCD[]".  Put another way, computing "fib[15]" will always benefit from knowing the values which "fib[8]" was forced to compute, while computing "recurGCD[10,15]" does not use any of the values found during the computation of "recurGCD[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 "recurGCD[]".  If the range of the parameters is very large, or if we will be computing relatively few instances of "recurGCD[]", 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.