Algorithms intro: implementing a math proof

What's a nice way of showing a visual example to better understand how a math proof can be implemented as an algorithm in code? Here we are showing an example of implementing an algorithm, comparing iterative and recursive algorithms for solving a method; and validating the results with profiling.

In problem solving, a great methodology to use is: Understand, Plan, Execute and Review.
We will use a bit of geometry to understand the problem, and then plan how to solve it; then we will implement code to execute that plan and review how things worked out.

Understand

The Euclidean algorithm calculates the greatest common divisor (GCD) of two natural numbers a and b. (or x and y).

We can use some geometry to visualize this, check out this video of it in action. From that, it looks like a matter of:

1. Fill with squares along the shortest side of the rectangle
2. For the remaining rectangle, fill with squares again, until smallest square is found.

Plan

Lets write down the steps in more detail:

• fill rectangle with squares
• find shortest sid
• divide long side by square side, get remainder
• if remainder is 0, smallest square is found
• for remainder, fill rectangle with squares

Execute

To model the visualization, lets make a Rectangle object, and use that with a function. Lets start with what we know, the Rectangle:

class Rectangle(object):
    def __init__(self, x, y):
        self.x = x
        self.y = y
    def __str__(self):
        return "%s, %s" % (self.x, self.y)

Now lets plan how to fill the rectangle just like the visualization happens, the next rectangle is filled, and so on until the end. This style of implementation is recursive.

• fill_rectangle_with_squares(2,4)
• fill_rectangle_with_squares(2,2)
• 2 is smallest square

We keep calling the fill_rectangle_with_squares function, until the x and y sides are equal. This is the stopping condition.

def fill_rectangle_with_squares(rectangle):
    print("Fill " + str(rectangle))
    # check the sides, if they aren't equal make a new rectangle with the remaining
    if rectangle.x == rectangle.y:
        print('smallest square found %s ' % str(rectangle))
        return rectangle.x
    elif rectangle.x > rectangle.y:
        return fill_rectangle_with_squares(Rectangle(rectangle.x - rectangle.y, rectangle.y))
    elif rectangle.y > rectangle.x:
        return fill_rectangle_with_squares(Rectangle(rectangle.x, rectangle.y - rectangle.x))

Now we call the function with an assert to test the result

To test to see if this function returns what we want, we use assert to check. We want to assert that a 2 by 4 rectangle is anything other than 2. Lets test some others as well

assert fill_rectangle_with_squares(Rectangle(2,2)) == 2
assert fill_rectangle_with_squares(Rectangle(140,68)) == 4

assert fill_rectangle_with_squares(Rectangle(24,12)) == 12

Review

Lets keep reviewing our solution with other sizes, like: Rectangle(143,62) and Rectangle(270,192). These work fine, but what happens when larger values are involved?

Rectangle(2400000223334219423423424234240,1945303434230234234242322)

When this would fail would vary from machine to machine, but on my home laptop I get a:

RuntimeError: maximum recursion depth exceeded


This means the call stack limit is reached, and this limit is bound by the design of the python language itself. Some languages like LISP don't have this limitation, it's just ends up using what resources it has. The limit is imposed because the amount of callable memory is a finite space on the machine. If you are using a large amount of recursion depth, and not getting the results you need then there is way to only bound by time, and the running of the CPU.

How can we make this better? Lets implement the algorithm using iteration.

def fill_rectangle_with_squares_iterative(x,y):
    rectangle = Rectangle(x,y)
    print('fill %s with squares ' % str(rectangle))
    found = False
    while not found:
        found = rectangle.x == rectangle.y
        if rectangle.x > rectangle.y:
            rectangle.x = rectangle.x - rectangle.y
        elif rectangle.y > rectangle.x:
            rectangle.y = rectangle.y - rectangle.x
    print('smallest square found %s ' % str(rectangle))

So far, so good, and the tests work with fill_rectangle_with_squares_iterative(24000000,19453003)

Can the iterative version take the massive difficult rectangle that blew the call stack with the recursive version?

fill_rectangle_with_squares_iterative(2400000223334219423423424234240,1945303434230234234242322)

It does, but takes a couple seconds.

Optimize

I know it took a couple seconds because the console lets me know how long any code execution takes. What's really happening here? To find out more I need to profile the code which will will in the second part of this post.

Before we do though, lets take a crack at optimizing the code even farther.

When you look at the iterative function, see that we are just modifying the Rectangle objects x and y values. We can just keep track of the x and y and not bother making the Rectangle object at all. We could maybe make one object and pass it instead of making a new object for each call, but lets go the step further. It's only 2 variables that makes it concise.

We used the Rectangle as a way to visualize the problem, but nice terse notation is possible and ends up being very clear, so lets boil down the description to the original:

Given two (natural) numbers not prime to one another, to find their greatest common measure.
(The Elements: Book VII: Proposition 2)

def greatest_common_divisor(x,y):
    while not x == y:
        if x > y:
            x = x - y
        elif y > x:
            y = y - x
    return x
    
print('greatest common divisor found %s' % (str(greatest_common_divisor(24,10))))

That's much more straightforward, and can probably be more concise, but I'd like to stop here to keep it readable and easier to understand.

Wrap it up

Lets review what we did in the post. we understood the problems and used some geometry to show a solution. Since it was Euclids proof, it seems like the right thing to do. We then planned a solution using recursion as it matched our real world model.
This worked to a point, but we reviewed our options and implemented an iterative method that could handle much larger rectangles. Once we had that solved we optimized for the computational abilities of the language and the machine we are using.

Our last implementation works for all the tests above, and certainly seems to run faster. How much faster?
We will profile all these functions and find out in the next post:  http://jseller.blogspot.ca/2016/11/algorithms-intro-profiling-python-code.html