HW: Getting started with Sage

Here is an example of the kind of programming we might be doing in Sage. The problem: Create a visualization of integer divisibility. This is a widely open-ended problem and it has many interpretations and solutions. Here is how I will choose to interpret it today: Consider all pairs of integers $(n,m)$ as points on the $xy$-plane, and highlight all the points $(m,n)$ such that $n$ is divisible by $m$.

Let’s go through the steps for doing this one at a time.

Deciding when $n$ is divisible by $m$

There are many ways to do this in Sage. One way is to use the “remainder after division” operator % provided by Python. A more Sage-y way to do it is to test whether $\frac nm$ lies in the ring of integers, which is denoted ZZ in sage.

8/2 in ZZ

Output: True

8/3 in ZZ

Output: False

Looping to create an array of pairs

We want to consider all values of $n$ between 1 and 100. Sage has a quick way of producing these numbers:

range(1,101)

Output: An array containing [1,2,3,…,100]

Notice that we had to use 101 and not 100 as our second argument, because the range function always stops just before the second argument!

Now we want to loop over all $m$ and $n$ in this range and enter the pair $(m,n)$ into an array whenever $n$ is divisible by $m$. First we create an empty array for the data. Then we use the for variable in array: construct, which performs the following instructions over and over again, with the given variable looping through all values in the given array. We also use the if condition: construct, which performs the following instructions only if the given condition holds.

divisor_list = []
for n in range(1,101):
    for m in range(1,n+1):
        if Integer(n)/Integer(m) in ZZ:
            divisor_list.append((n,m))

[Notice that instead of writing n/m we wrote Integer(n)/Integer(m). This is just a technicality that you might encounter occaisionally. The integers that are created inside ranges are actually Python integers, which get rounded down when they divide, and that is not what we want! Writing Integer forces them to be Sage integers, which divide symbolically.]

Plotting the output

We can plot an array of points using the point command. Not surprisingly, there is also a handy plot command, but always remember: plot is good for functions, point is good for data points.

point(divisor_list)

Output: an interesting visualization. Try staring at it until you see patterns.

An added flourish

Sometimes it is helpful to visualize different types of data with different colors. Suppose in this example that we wanted to see all of the prime divisors in red instead of blue. We can put just the primes into a different array and plot them separately.

prime_divisor_list = []
nonprime_divisor_list = []
for n in range(1,101):
    for m in range(1,n+1):
        if Integer(n)/Integer(m) in ZZ:
            if is_prime(m):
                prime_divisor_list.append((n,m))
            else:
                nonprime_divisor_list.append((n,m))
point(nonprime_divisor_list) + point(prime_divisor_list, color="red")

Notice that to superimpose two plots we use the + operation between them.

Putting everything into a function

The code above works perfectly well for the first 100 integers, but if we want a larger picture we have to modify the code. To generalize this, we can wrap the code inside a function with 100 replaced by a parameter. We begin the code using the def function name(parameter list): construction.

def plot_divisors(k):
    prime_divisor_list = []
    nonprime_divisor_list = []
    for n in range(1,k+1):
        for m in range(1,n+1):
            if Integer(n)/Integer(m) in ZZ:
                if is_prime(m):
                    prime_divisor_list.append((n,m))
                else:
                    nonprime_divisor_list.append((n,m))
    return point(nonprime_divisor_list) + point(prime_divisor_list, color="red")

Notice that last line is preceded by the return instruction, which precedes the final value that the function produces. In this case, the final value is the plot itself. Now we can obtain a plot of any size using a function call.

plot_divisors(150)

Output: a similar plot but going all the way up to 150 instead of just 100.

Your turn

Your task will be the following: Create a visualization of how common it is for a number to be prime.. Once again, this is wildly open-ended, so if you’d like you may use these more specific instructions: Create a graph of the points $(n,m)$ where $m$ is the number of primes that are less than or equal to $n$.

Getting started

Log in to the sage notebook at sagenb.org using your google email account. When you type a command, you have to use [shift]+[enter] to run them.

It will help you to start with the following warm-up: Write a function that takes one parameter, $n$, and returns an array containing all of the prime numbers that are less than or equal to $n$.

You can find the length of an array using the len function.

Here is what my output picture looks like:

prime_count