Introduction to Algorithms
Week 2 – Notes

primality testing

As we know from mathematics, a prime number is an integer greater than 1 whose only factors are 1 and itself. For example, 2, 7, 47 and 101 are all prime. A integer greater than 1 that is not prime is said to be composite.

Prime numbers were known to the ancient Greeks, and mathematicians have been interested in prime numbers for millennia. With the advent of digital computers, prime numbers have taken on a new importance since large primes are commonly used in cryptographic algorithms.

We would now like to write a program that tests whether a given number is prime. To do this, we'll use a simple algorithm called trial division, which means dividing by each possible factor in turn. Here is a naive implementation of trial division:

n = int(input('Enter n: '))

prime = True

for i in range(2, n):  # loop from 2 .. (n - 1)
    if n % i == 0:
        prime = False
        break

if prime:
    print('n is prime')
else
    print('not prime')

This works fine, but is inefficient because it may need to test all integers from 2 up to (n – 1). When n is large, this can take a long time.

Actually for a given n we need test only the values up to √n, i.e. the square root of n. To see this, consider the following fact:

Fact. If n is composite, then it has a prime factor that is less than or equal to √n.

Proof. If n is composite, then n = ab for some integers 1 < a, b < n. Suppose that a > √n and b > √n. Then ab > √n ⋅ √n = n, which is a contradiction since ab = n. So either a ≤ √n or b ≤ √n.

It follows that if we have tested all the values from 2 through √n and none of them divide n, then n cannot be composite, so it is prime.

Here's a more efficent version of our program that stops at √n:

from math import sqrt

n = int(input('Enter n: '))

prime = True
for i in range(2, int(sqrt(n)) + 1):  # loop from 2 .. sqrt(n)
    if n % i == 0:
        prime = False
        break

if prime:
    print('n is prime')
else:
    print('not prime')

This program is far more efficient. For example, suppose n is a prime number that is approximately 1,000,000. Then we will now check only √n ≈ 1,000 rather than n ≈ 1,000,000 possible factors, which is 1000 times as efficient!

Further optimizations are possible. For example, if n is odd, then it cannot be divisible by any even number. So if we first check that n is not 2 (which is prime) or any other even number (which is trivially composite), we can then loop only over the possible factors 3, 5, 7, and so on up to √n. This will check only about half as many possible factors as the code above. You may wish to implement this optimization as an exercise.

prime factorization

The Fundamental Theorem of Arithmetic states that

Every positive integer has a unique prime factorization.

For example:

You can see a proof of this theorem (which is not completely trivial) in an introductory number theory course.

We can use trial division to factor an integer. This is similar to our primality testing algorithm, which also used trial division. Given an integer N, we first try dividing N by 2, then by 3 and so on. If n is actually divisible by a number such as 3, then we must repeatedly attempt to divide by that same number. since the same prime factor may appear multiple times in the factorization.

Here is a first program for prime factorization:

n = int(input('Enter n: '))

f = 2
while n > 1:
    if n % f == 0:
        print('factor: ', f)
        n = n // f
    else:
        f += 1

Here's the program in action:

$ py factor.py
Enter n: 60
factor:  2
factor:  2
factor:  3
factor:  5

Study the program to understand how it works. The program tries dividing by all possible values. Once n is 1, we have divided out all factors, so the factorization is complete.

Note that the program will attempt to divide by non-prime factors, such as when i = 4. But n will never be divisible by such values. That's because any non-prime factor i is itself the product of two (or more) smaller primes, and we have already divided those primes out of n, so n cannot possibly be divisible by i.

The program works, but is inefficient because it potentially tests all values from 1 to n. Let's modify it so that we only test for factors up to sqrt(n), just like in our second primality testing program above:

from math import sqrt

n = int(input('Enter n: '))

i = 2
s = ''
while i <= int(sqrt(n)):
    if n % i == 0:  # n is divisible by i
        s += str(i) + ' * '
        n = n // i
    else:
        i += 1

s += str(n)     # append the last factor
print(s)

This version of the program also prints out all the factors on a single line:

Enter n: 644
2 * 2 * 7 * 23

example: Fermat numbers

As an interesting example involving prime numbers, consider the Fermat numbers Fn, defined as

Fn = 2(2^n) + 1

We have

All of the numbers F0 … F4 are prime. Seeing this, the famous French mathematician Pierre de Fermat conjectured in 1650 that all Fn are prime.

He was wrong, since Leonhard Euler showed in 1732 that F5 is composite:

In fact, no more prime numbers Fn have been discovered past F4! Many mathematicians believe that in fact no Fn with n > 4 is prime, though nobody has yet been able to prove this.

Let's use the program we wrote above to factor F5 and F6:

$ py factor.py
Enter n: 4_294_967_297
641 * 6700417

$ py factor.py
Enter n: 18_446_744_073_709_551_617    
274177 * 67280421310721

This is impressive. However, factoring

F7 = 340,282,366,920,938,463,463,374,607,431,768,211,457

using this program would take a very long time (e.g. a year or more), since its smallest factor is a 17-digit number!

greatest common divisor

For integers a and b, the greatest common divisor of a and b, also written as gcd(a, b), is the largest integer that evenly divides both a and b. For example:

The greatest common divisor is useful in various situations, such as simplifying fractions. For example, suppose that we want to simplify 35/15. gcd(35, 15) = 5, so we can divide both the numerator and denominator by 5, yielding 7/3.

Naively, we can find the gcd of two values by trial division:

a = int(input('Enter a: '))
b = int(input('Enter b: '))

for i in range(min(a, b), 0, -1):
  if a % i == 0 and b % i == 0:
    break

print('The gcd is', i)

But this may take time O(N), where N = min(a, b). We'd like to be more efficient, especially since in some situations (e.g. cryptography) we may want to take the gcd of numbers that are very large.

Another way to find the greatest common divisor of two numbers is to factor the numbers into primes. For every prime p, if the prime factorization of a includes pi, and the prime factorization of b includes pj, then the prime factorization of gcd(a, b) will include pmin(i, j). (In other words, the prime factorization of gcd(a, b) is the intersection of the prime factorizations of a and b.)

For example consider finding gcd(252, 90). 252 = 9 * 28 = 22 * 32 * 71. And 90 = 9 * 10 = 21 * 32 * 51. So the common primes are 21 * 32 = 18, and we have gcd(252, 90) = 18.

Euclid's algorithm

Euclid’s algorithm is a much more efficient way to find the gcd of two numbers than prime factorization. It is based on the fact that for all positive integers a and b, gcd(a, b) = gcd(b, a mod b). (We will prove this below.) So, for example,

gcd(252, 90)
= gcd(90, 72)
= gcd(72, 18)
= gcd(18, 0)
= 18

We can implement Euclid's algorithm in Python like this:

a = int(input('Enter a: '))
b = int(input('Enter b: '))

while b > 0:
  a, b = b, a % b

print('The gcd is', a)

Notice that the line "a, b = b, a % b" is a simultaneous assignment. You can read it like this: the new value of a is the previous value of b, and the new value of b is the previous value of (a % b).

In all the examples of gcd(a, b) above we had a > b. But note that this function works even when a < b ! For example, if we call gcd(5, 15), then on the first iteration of the while loop we have a % b = 5 mod 15 = 5. So then we assign a = 15 and b = 5, and further iterations continue as if we had called gcd(15, 5).

Above, we said that Euclid's algorithm works because of this mathematical fact: for all integers a and b,

gcd(a, b) = gcd(b, a mod b)

Let's now prove that. First, let's introduce the notion of a linear combination of two integers a and b, which is any integer of the form ia + jb for some integers i and j. For example, 26 is a linear combination of 3 and 5, because 2 · 3 + 4 · 5 = 26. To put it differently, a linear combination of a and b is any number that we can form by adding up some number of a's plus some number of b's.

The notation d | n means that d divides n, i.e. n is a multiple of d. Notice that if d | a and d | b, then d also divides any linear combination of a and b. For example, any linear combination of 6 and 10 must be divisible by 2, because 2 | 6 and 2 | 10.

We'd like to show that the sets {a, b} and {b, a mod b} have the same common divisors. Suppose that d | b and d | a mod b. We have already seen this equation that relates a, b, a div b and a mod b:

a = b · (a div b) + (a mod b)

We see that a is a linear combination of b and (a mod b). So we know that d | a. And so if d is a divisor of both b and (a mod b), then it's a divisor of both a and b.

Now suppose that d | a and d | b. We can rearrange the equation above to to

a mod b = a – b · (a div b)

We see that (a mod b) is a linear combination of a and b. So we know that d | (a mod b). And so if d is a divisor of both a and b, then it's a divisor of both b and (a mod b).

So the sets {a, b} and {b, a mod b} have the same common divisors. Therefore the sets have the same greatest common divisors, and so gcd(a, b) = gcd(b, a mod b).

Our implementation of Euclid's algorithm is quite fast. To see it in action, let's add a print statement so that we can see the values of a and b as the algorithm runs:

a = int(input('Enter a: '))
b = int(input('Enter b: '))

while b > 0:
    print('a =', a, '; b =', b)
    a, b = b, a % b

print(a)

Now let's try running it on two large integers:

$ py gcd.py 
Enter a: 384753412340000
Enter b: 937845943000
a = 384753412340000 ; b = 937845943000
a = 937845943000 ; b = 236575710000
a = 236575710000 ; b = 228118813000
a = 228118813000 ; b = 8456897000
a = 8456897000 ; b = 8239491000
a = 8239491000 ; b = 217406000
a = 217406000 ; b = 195469000
a = 195469000 ; b = 21937000
a = 21937000 ; b = 19973000
a = 19973000 ; b = 1964000
a = 1964000 ; b = 333000
a = 333000 ; b = 299000
a = 299000 ; b = 34000
a = 34000 ; b = 27000
a = 27000 ; b = 7000
a = 7000 ; b = 6000
a = 6000 ; b = 1000
1000

Notice how quickly a and b decrease as the algorithm runs. In the next lecture, we'll see how to express the running time of Euclid's algorithm in a mathematically precise way.

As a final note on this topic, suppose that we want to compute the greatest common divisor of three numbers, i.e. gcd(a, b, c). This is easy: we can just compute

gcd(a, b, c) = gcd(a, gcd(b, c))

least common multiple

We have seen how to compute the greatest common divisor of two integers. A related concept is the least common multiple of two integers p and q, which is the smallest integer that is divisible by both p and q. For example,

lcm(60, 90) = 180

As with the gcd, we can find the least common multiple of two numbers by computing their prime factorizations. For every prime p, if the prime factorization of a includes pi, and the prime factorization of b includes pj, then the prime factorization of gcd(a, b) will include pmax(i, j). (In other words, the prime factorization of gcd(a, b) is the union of the prime factorizations of a and b.)

How can we compute a least common multiple efficiently? Here is a useful fact: for all integers and b,

gcd(a, b) · lcm(a, b) = a · b

Why is this true? First observe that for any integers i and j, trivially

min(i, j) + max(i, j) = i + j

Now, given integers a and b, consider any prime p. Suppose that the prime factorization of a contains pi, and the prime factorization of b contains pj. Then we have already seen that the prime factorization of gcd(a, b) contains pmin(i, j), and the prime factorization of lcm(a, b) contains pmax(i, j). Therefore the prime factorization of the product gcd(a, b) · lcm(a, b) will contain pmin(i, j) + max(i, j) = pi + j, which is also present in the prime factorization of the product a · b. Since this is true for all primes in the factorizations, we have gcd(a, b) · lcm(a, b) = a · b.

Now we can rearrange the formula above to give us lcm(a, b):

lcm(a, b) = a · b / gcd(a, b)

We can compute the gcd efficiently using Euclid's algorithm, so this formula gives us an efficient way to compute the lcm.

Similarly to the gcd, if we want the least common multiple of several numbers we can just compute

lcm(a, b, c) = lcm(a, lcm(b, c))

limits

You should all have seen the concept of limits in your high school mathematics courses. In computer science we are often concerned with limits as a variable goes to infinity. For example:

orders of growth

We will say that a function f(n) grows faster than a function g(n) if the limit of f(n) / g(n) as n approaches ∞ is ∞. In other words, the ratio f(n) / g(n) becomes arbitrarily large as n increases. For example, n2 grows faster than 10n, since

lim n→∞ (n2 / 10n) = lim n→∞ (n / 10) = ∞

Here is a list of some common functions of n in increasing order. Each function in the list grows faster than all of its predecessors: