Part 2: Functions in Python

A distinguishing property of programming languages is that the programmer can create their own functions. Creating a function is like teaching the computer a new trick. Typically a function will receive some data as input, will perform an algorithm involving the input data, and will output data when the algorithm terminates.

In this part, we explore Python functions. We also explore control statements, which allow a program to behave in different ways for different inputs. We also introduce the while loop, a loop whose repetition can be more carefully controlled than a for loop. As an application of these techniques, we implement the Euclidean algorithm as a Python function in a few ways, to effectively find the GCD of integers and solve linear Diophantine equations. This complements Chapter 1 of An Illustrated Theory of Numbers.

Getting started with Python functions

A function in Python is a construction which takes input data, performs some actions, and outputs data. It is best to start with a few examples and break down the code. Here is a function square. Run the code as usual by pressing shift-Enter when the code block is selected.

In [ ]:
def square(x):
    answer = x * x
    return answer

When you run the code block, you probably didn't see anything happen. But you have effectively taught your computer a new trick, increasing the vocabulary of commands it understands through the Python interpreter. Now you can use the square command as you wish.

In [ ]:
square(12)
In [ ]:
square(1.5)

Let's break down the syntax of the function declaration, line by line.

def square(x):
    answer = x * x
    return answer

The first line begins with the Python reserved word def. (So don't use def as a variable name!). The word def stands for "define" and it defines a function called square. After the function name square comes parentheses (x) containing the argument x. The arguments or parameters of a function refer to the input data. Even if your function has no arguments, you need parentheses. The argument x is used to name whatever number is input into the square function.

At the end of the function declaration line is a colon : and the following two lines are indented. As in the case of for loops, the colon and indentation are signals of scope. Everything on the indented lines is considered the scope of the function and is carried out when the function is used later.

The second line answer = x * x is the beginning of the scope of the function. It declares a variable answer and sets the value to be x * x. So if the argument x is 12, then answer will be set to 144. The variable answer, being declared within the scope of the function, will not be accessible outside the scope of the function. It is called a local variable.

The last line return answer contains the Python reserved word return, which terminates the function and outputs the value of the variable answer. So when you apply the function with the command square(1.5), the number 1.5 is the argument x, and answer is 2.25, and that number 2.25 becomes the output.

A function does not have to return a value. Some functions might just provide some information. Here is a function which displays the result of division with remainder as a sentence with addition and multiplication.

In [ ]:
def display_divmod(a,b):
    quotient = a // b # Integer division
    remainder = a % b #
    print "%d = %d (%d) + %d"%(a,quotient,b,remainder)
In [ ]:
display_divmod(23,5)

Notice that this function has no return line. The function terminates automatically at the end of its scope.

The function also uses a slick Python string manipulation called string substitution. This has changed between Python 2.x and 3.x, and as usual we are following Python 2.x syntax.

String substitution allows you to place variable numbers in a string by putting placeholders like %d within the string, and then placing the list of variables after the string (as a tuple). Here are some examples to get used to the syntax, since we will use this often in this lesson.

In [ ]:
print "My favorite number is %d"%(17)  # The % symbol substitutes 17 for %d.
In [ ]:
print "%d + %d = %d"%(13,12,13+12)

The placeholder %d is meant to hold the place of an integer. For floating point numbers, use %f instead. This will impose a bit of rounding before printing. One can modify the placeholder to change the amount of digits displayed. Here is the official reference for placeholder syntax and string substitutions. We will only use the most basic features.

In [ ]:
print "pi is approximately %f"%(3.14159265) # This uses a default precision for display.
In [ ]:
print "pi is approximately %.3f"%(3.14159265) # This sets 3 digits after the decimal point.

Exercises

  1. What are the signals of scope in Python?

  2. How is the symbol % used in Python? Can you think of three ways?

  3. Write a function called area_circle, which takes one argument radius. The function should return the area of the circle, as a floating point number. Then add one line to the function, using string substitutions, so that it additionally prints a helpful sentence of the form "The area of a circle of radius 1.0 is 3.14159." (depending on the radius and the area it computes).

  4. Can you think of a reason you might want to have a function with no arguments?

In [ ]:
#  Use this space to work on the Exercises.  
#  Remember that you can add a new cell above/below by clicking to the left of a cell,
#  (the cell will have a blue bar at the left) and then pressing "a" or "b" on the keyboard.

Control statements

It is important for a computer program to behave differently under different circumstances. The simplest control statements, if and its relative else, can be used to tell Python to carry out different actions depending on the value of a boolean variable. The following function exhibits the syntax.

In [ ]:
def is_even(n):
    if n%2 == 0:
        print "%d is even."%(n)
        return True
    else:
        print "%d is odd."%(n)
        return False
In [ ]:
is_even(17)
In [ ]:
is_even(1000)

The broad syntax of the function should be familiar. We have created a function called is_even with one argument called n. The body of the function uses the control statement if n%2 == 0:. Recall that n%2 gives the remainder after dividing n by 2. Thus n%2 is 0 or 1, depending on whether n is even or odd. Therefore the boolean n%2 == 0 is True if n is even, and False if n is odd.

The next two lines (the first print and return statements) are within the scope of the if <boolean>: statement, as indicated by the colon and the indentation. The if <boolean>: statement tells the Python interpreter to perform the statements within the scope if the boolean is True, and to ignore the statements within the scope if the boolean is False.

Putting it together, we can analyze the code.

if n%2 == 0:
        print "%d is even."%(n)
        return True

If n is even, then the Python interpreter will print a sentence of the form n is even. Then the interpreter will return (output) the value True and the function will terminate. If n is odd, the Python interpreter will ignore the two lines of scope.

Often we don't just want Python to do nothing when a condition is not satisfied. In the case above, we would rather Python tell us that the number is odd. The else: control statement tells Python what to do in case the if <boolean>: control statement receives a False boolean. We analyze the code

else:
        print "%d is odd."%(n)
        return False

The print and return commands are within the scope of the else: control statement. So when the if statement receives a false signal (the number n is odd), the program prints a sentence of the form n is odd. and then returns the value False and terminates the function.

The function is_even is a verbose, or "talkative" sort of function. Such a function is sometimes useful in an interactive setting, where the programmer wants to understand everything that's going on. But if the function had to be called a million times, the screen would fill with printed sentences! In practice, an efficient and silent function is_even might look like the following.

In [ ]:
def is_even(n):
    return (n%2 == 0)
In [ ]:
is_even(17)

A for loop and an if control statement, used together, allow us to carry out a brute force search. We can search for factors in order to check whether a number is prime. Or we can look for solutions to an equation until we find one.

One thing to note: the function below begins with a block of text between a triple-quote (three single-quotes when typing). That text is called a docstring and it is meant to document what the function does. Writing clear docstrings becomes more important as you write longer programs, collaborate with other programmers, and when you want to return months or years later to use a program again. There are different style conventions for docstrings; for example, here are Google's docstring conventions. We take a less formal approach.

In [ ]:
def is_prime(n):
    '''
    Checks whether the argument n is a prime number.
    Uses a brute force search for factors between 1 and n.
    '''
    for j in range(2,n):  # the list of numbers 2,3,...,n-1.
        if n%j == 0:  # is n divisible by j?
            print "%d is a factor of %d."%(j,n)
            return False
    return True

An important note: the return keyword terminates the function. So as soon as a factor is found, the function terminates and outputs False. If no factor is found, then the function execution survives past the loop, and the line return True is executed to terminate the function.

In [ ]:
is_prime(91)
In [ ]:
is_prime(101)

Try the is_prime function on bigger numbers -- try numbers with 4 digits, 5 digits, 6 digits. Where does it start to slow down? Do you get any errors when the numbers are large? Make sure to save your work first, just in case this crashes your computer!

In [ ]:
# Experiment with is_prime here.

There are two limiting factors, which we study in more detail in the next lesson. These are time and space (your computer's memory space). As the loop of is_prime goes on and on, it might take your computer a long time! If each step of the loop takes only a nanosecond (1 billionth of a second), the loop would take about a second when executing is_prime(1000000001). If you tried is_prime on a much larger number, like is_prime(2**101 - 1), the loop would take longer than the lifetime of the Earth.

But if you try such a big number, you'll hit a more immediate problem with space. Namely, the range(2,n) attempts to store the entire list of numbers [2,3,4,...,n-1] in the memory of your computer. Your computer has some (4 or 8 or 16, perhaps) gigabytes of memory (RAM). A gigabyte is a billion bytes, and a byte is enough memory to store a number between 0 and 255. (More detail about this later!). So a gigabyte will not even hold a billion numbers, and your computer will probably run out of memory if you attempt to try is_prime(n) when n is close to a billion.

We will address both of these problems in the next lesson, to some extent. The memory (space) problem is easier, and we will be able to adapt the is_prime(n) function to deal with n in the trillions. To go beyond this, to work with n of size $10^{100}$ or even $10^{1000}$, we will find completely different methods in a later lesson.

Exercises

  1. Create a function my_abs(x) which outputs the absolute value of the argument x. (Note that Python already has a built-in abs(x) function).

  2. Modify the is_prime function so that it prints a message Number too big and returns None if the input argument is bigger than one million. (Note that None is a Python reserved word. You can use the one-line statement return None.)

  3. Write a Python function thrarity which takes an argument n, and outputs the string threeven if n is a multiple of three, or throdd is n is one more than a multiple of three, or thrugly if n is one less than a multiple of three. Example: thrarity(31) should output throdd and thrarity(44) should output thrugly. Hint: study the if/elif syntax at the official Python tutorial

  4. Write a Python function sum_of_squares(n) which finds and prints a pair of natural numbers $x$, $y$, such that $x^2 + y^2 = n$. The function should use a brute force search.

In [ ]:
#  Use this space for your solutions to the questions.

While loops and implementation of the Eucidean algorithm

We almost have all the tools we need to implement the Euclidean algorithm. The last tool we will need is the while loop. We have seen the for loop already, for iterating over a range of numbers. The Euclidean algorithm involves repetition, but there is no way to know in advance how many steps it will take. The while loop allows us to repeat a process as long as a boolean value (sometimes called a flag) is True. The following countdown example illustrates the structure of a while loop.

In [ ]:
def countdown(n):
    current_value = n
    while current_value > 0:  # The condition (current_value > 0) is checked before every instance of the scope!
        print current_value
        current_value = current_value - 1
In [ ]:
countdown(10)

The while loop syntax begins with while <boolean>: and the following indented lines comprise the scope of the loop. If the boolean is True, then the scope of the loop is executed. If the boolean is True again afterwards, then the scope of the loop is executed again. And again and again and so on.

This can be a dangerous process! For example, what would happen if you made a little typo and the last line of the while loop read current_value = current_value + 1? The numbers would increase and increase... and the boolean current_value > 0 would always be True. Therefore the loop would never end. Bigger and bigger numbers would scroll down your computer screen.

You might panic under such a circumstance, and maybe turn your computer off to stop the loop. Here is some advice for when your computer gets stuck in such a neverending loop:

  1. Back up your work often. When you're programming, make sure everything else is saved just in case.
  2. Save your programming work (use "Save and checkpoint" under the "File" menu) often, especially before running a cell with a loop for the first time.
  3. If you do get stuck in a neverending loop, click on "Kernel... Interrupt". This will often unstick the loop and allow you to pick up where you left off.
  4. On a Mac, you might try a "Force Quit" of the Python process, using the Activity Manager.

Now, if you're feeling brave, save your work, change the while loop so that it never ends, and try to recover where you left off. But be aware that this could cause your computer to freeze or behave erratically, crashing your browser, etc. Don't panic... it won't break your computer permanently.

The neverending loop causes two problems here. One is with your computer processor, which will be essentially spinning its wheels. This is called busy waiting, and your computer will essentially be busy waiting forever. The other problem is that your loop is printing more and more lines of text into the notebook. This could easily crash your web browser, which is trying to store and display zillions of lines of numbers. So be ready for problems!

The Euclidean algorithm with a while loop

The Euclidean Algorithm is a process of repeated division with remainder. Beginning with two integers a (dividend) and b (divisor), one computes quotient q and remainder q to express a = qb + r. Then b becomes the dividend and r becomes the divisor, and one repeats. The repetition continues, and the last nonzero remainder is the greatest common divisor of a and b. It is the subject of Chapter 1 of An Illustrated Theory of Numbers.

We implement the Euclidean algorithm in a few variations. The first will be a verbose version, to show the user what happens at every step. We use a while loop to take care of the repetition.

In [ ]:
def Euclidean_algorithm(a,b):
    dividend = a
    divisor = b
    while divisor != 0:   # Recall that != means "is not equal to".
        quotient = dividend // divisor
        remainder = dividend % divisor
        print "%d = %d (%d) + %d"%(dividend, quotient, divisor, remainder)
        dividend = divisor  
        divisor = remainder
In [ ]:
Euclidean_algorithm(133, 58)
In [ ]:
Euclidean_algorithm(1312331323, 58123123)

This is excellent if we want to know every step of the Euclidean algorithm. If we just want to know the GCD of two numbers, we can be less verbose. We carefully return the last nonzero remainder after the while loop is concluded. This last nonzero remainder becomes the divisor when the remainder becomes zero, and then it would become the dividend in the next (unprinted) line. That is why we return the (absolute value) of the dividend after the loop is concluded. You might insert a line at the end of the loop, like print dividend, divisor, remainder to help you track the variables.

In [ ]:
def GCD(a,b):
    dividend = a # The first dividend is a.
    divisor = b # The first divisor is b.
    while divisor != 0:   # Recall that != means "not equal to".
        quotient = dividend // divisor
        remainder = dividend % divisor
        dividend = divisor  
        divisor = remainder
    return abs(dividend)  #  abs() is used, since we like our GCDs to be positive.

Note that the return dividend statement occurs after the scope of the while loop. So as soon as the divisor variable equals zero, the funtion GCD returns the dividend variable and terminates.

In [ ]:
GCD(111,27)
In [ ]:
GCD(111,-27)

We can refine our code in a few ways. First, note that the quotient variable is never used! It was nice in the verbose version of the Euclidean algorithm, but plays no role in finding the GCD. Our refined code reads

def GCD(a,b):
    dividend = a  
    divisor = b  
    while divisor != 0:   # Recall that != means "not equal to".
        remainder = dividend % divisor
        dividend = divisor  
        divisor = remainder
    return abs(dividend)

Now there are two slick Python tricks we can use to shorten the code. The first is called multiple assignment. It is possible to set the values of two variables in a single line of code, with a syntax like below.

In [ ]:
x,y = 2,3  # Sets x to 2 and y to 3.

This is particular useful for self-referential assignments, because as for ordinary assignment, the right side is evaluated first and then bound to the variables on the left side. For example, after the line above, try the line below. Use print statements to see what the values of the variables are afterwards!

In [ ]:
x,y = y,x #  Guess what this does!
In [ ]:
print "x =",x
print "y =",y

Now we can use multiple assignment to turn three lines of code into one line of code. For the remainder variable is only used temporarily before its value is given to the divisor variable. Using multiple assignment, the three lines

remainder = dividend % divisor
    dividend = divisor  
    divisor = remainder

can be written in one line,

dividend, divisor = divisor, dividend % divisor # Evaluations on the right occur before any assignments!

Our newly shortened GCD function looks like this.

def GCD(a,b):
    dividend = a  
    divisor = b  
    while divisor != 0:   # Recall that != means "not equal to".
        dividend, divisor = divisor, dividend % divisor
    return abs(dividend)

The next trick involves the while loop. The usual syntax has the form while <boolean>:. But if while is followed by a numerical type, e.g. while <int>:, then the scope of the while loop will execute as long as the number is nonzero! Therefore, the line

while divisor != 0:

can be replaced by the shorter line

while divisor:

This is truly a trick. It probably won't speed anything up, and it does not make your program easier to read for beginners. So use it if you prefer communicating with experienced Python programmers! Here is the whole function again.

def GCD(a,b):
    dividend = a  
    divisor = b  
    while divisor:   # Executes the scope if divisor is nonzero.
        dividend, divisor = divisor, dividend % divisor
    return abs(dividend)

The next optimization is a bit more dangerous for beginners, but it works here. In general, it can be dangerous to operate directly on the arguments to a function. But in this setting, it is safe, and makes no real difference to the Python interpreter. Instead of creating new variables called dividend and divisor, one can manipulate a and b directly within the function. If you do this, the GCD function can be shortened to the following.

In [ ]:
def GCD(a,b):
    while b:   # Recall that != means "not equal to".
        a, b = b, a % b
    return abs(a)
In [ ]:
# Try it out.  Try it on some big numbers and see how quick it runs!

This code is essentially optimal, if one wishes to execute the Euclidean algorithm to find the GCD of two integers. It almost matches the GCD code in a standard Python library. It might be slightly faster than our original code -- but there is a tradeoff here between execution speed and readability of code. In this and the following lessons, we often optimize enough for everyday purposes, but not so much that readability is lost.

Exercises

  1. Modify the is_prime function by using a while loop instead of for j in range(2,n):. Hint: the function should contain the lines j = 2 and while j < n: and j = j + 1 in various places. Why might this be an improvement from the for loop?

  2. Modify the Euclidean_algorithm function to create a function which returns the number of steps that the Euclidean algorithm requires, i.e., the number of divisions-with-remainder.

  3. Create a function which carries out division with minimal remainder. In other words, given integers a,b, the function expresses a = q(b) + r, where r is a positive or negative integer of magnitude bounded by b/2. Use such a function to create a new Euclidean algorithm function which uses minimal remainder.

  4. What GCD(a,b) function do you think strikes the best balance between efficiency and readability?

In [ ]:
# Use this space to work on the exercises.

Solving the linear Diophantine equation

In Chapter 1 of An Illustrated Theory of Numbers, we not only used the Euclidean algorithm to find the GCD of two integers, but also to solve the linear Diophantine equation $ax + by = c$. On paper, this required us to perform the Euclidean algorithm, then "work backwards" to carefully solve a series of linear equations. This process is repetetive and error-prone... perfect for a computer.

So here we develop a function solve_LDE(a,b,c) which will describe all integer solutions $x,y$ to the equation $ax + by = c$.

The idea of the algorithm is to keep track of "hops" and "skips" throughout the Euclidean algorithm. A general step in the Euclidean algorithm looks like u = q(v) + r. The remainder can then be expressed by the formula r = u - q(v). If u and v can be built from hops and skips, then r can be built from hops and skips. How many? Just tally the hops and skips to find:

`r_hops = u_hops - q (v_hops)` and `r_skips = u_skips - q (v_skips)`.

This sort of tallying is what makes the algorithm below work. The function below does not introduce any new programming concepts, but it assembles many ideas together.

In [ ]:
def hop_and_skip(a,b):
    '''
    Takes two integer arguments a,b, and prints a sentence of the form
    GCD(a,b) = x(a) + y(b).  The method is the Euclidean algorithm,
    tallying hops (units of a) and skips (units of b) along the way.
    '''
    u = a # We use u instead of dividend.
    v = b # We use v instead of divisor.
    u_hops, u_skips = 1,0 # u is built from one hop (a) and no skips, for now.
    v_hops, v_skips = 0,1 # v is built from no hops and one skip (b), for now.
    while v != 0:   # We could just write "while v:"
        q = u // v  # q stands for quotient.
        r = u % v  # r stands for remainder.  So u = q(v) + r.
        
        r_hops = u_hops - q * v_hops  # Tally hops
        r_skips = u_skips - q * v_skips  # Tally skips
        
        u,v = v,r  # The new dividend,divisor is the old divisor,remainder.
        u_hops, v_hops = v_hops, r_hops # The new u_hops, v_hops is the old v_hops, r_hops
        u_skips, v_skips = v_skips, r_skips # The new u_skips, v_skips is the old v_skips, r_skips
    
    print "%d = %d(%d) + %d(%d)"%(u,u_hops,a,u_skips,b)
In [ ]:
hop_and_skip(102,45)

Try out the hop_and_skip code on some integers of your choice. Does it behave correctly? Check the results, using Python as a calculator. Does it run quickly for large integers?

In [ ]:
#  Experimentation space here.

To conclude this lesson, we put everything together to create a long(ish) function to solve linear Diophantine equations. We want this function to be smart enough to respond when an equation has no solutions, and to describe all solutions when they exist.

The first part of the function solve_LDE is the same as the hop and skip function above. But rather than expressing $GCD(a,b)$ as $ax + by$, the function uses the GCD to determine the existence and the general form of a solution to $ax + by = c$. The formula for the general form comes from An Illustrated Theory of Numbers, Chapter 1, Corollary 1.25.

In [ ]:
def solve_LDE(a,b,c):
    '''
    Describes all of the solutions to the linear Diophantine equation
    ax + by = c.  There are either no solutions or infinitely many solutions.
    Prints a description of the solution set, and returns None if there are no solutions
    or returns a single solution if one exists.
    '''  
    u = a # We use u instead of dividend.
    v = b # We use v instead of divisor.
    u_hops, u_skips = 1,0 # u is built from one hop (a) and no skips.
    v_hops, v_skips = 0,1 # v is built from no hops and one skip (b).
    while v != 0:   # We could just write while v:
        q = u // v  # q stands for quotient.
        r = u % v  # r stands for remainder.  So u = q(v) + r.
        
        r_hops = u_hops - q * v_hops  # Tally hops
        r_skips = u_skips - q * v_skips  # Tally skips
        
        u,v = v,r  # The new dividend,divisor is the old divisor,remainder.
        u_hops, v_hops = v_hops, r_hops # The new u_hops, v_hops is the old v_hops, r_hops
        u_skips, v_skips = v_skips, r_skips # The new u_skips, v_skips is the old v_skips, r_skips
    
    g = u # The variable g now describes the GCD of a and b.
    
    if c%g == 0:  # When GCD(a,b) divides c...
        d = c/g
        x = d * u_hops
        y = d * u_skips  # Now ax + by = c is a specific solution!
        print "%d x + %d y = %d if and only if "%(a, b, c)
        print "x = %d + %d n and y = %d - %d n, for some integer n."%(x,b/g,y,-a/g)
        return x,y
    else:  # When GCD(a,b) does not divide c...
        print "There are no solutions to %d x + %d y = %d,"%(a,b,c)
        print "because GCD(%d, %d) = %d, which does not divide %d."%(a,b,g,c)
        return None
In [ ]:
solve_LDE(102,45,3)
In [ ]:
solve_LDE(72,100,17)

Exercises

  1. Solve problems 4-7 of Chapter 1 of An Illustrated Theory of Numbers using the solve_LDE function.

  2. Write an LCM function, using the previous GCD function and the GCD-LCM product formula, Theorem 1.23 of An Illustrated Theory of Numbers.

  3. Sometimes it is important to find not the integer solutions, but the positive integer solutions to a Diophantine equation. Modify the solve_LDE function to create a solve_LDE_positive(a,b,c) function. The output of the function should be all pairs of positive integers $x$, $y$, such that $ax + by = c$ (if any pairs exist), and a helpful message if no pairs exist (and a return None should be used in this case).

In [ ]:
#  Use this space to work on the exercises.