Primality Testing

One classic application of Monte Carlo algorithms is primality testing. Primality testing means to test a number to determine whether it’s a prime number. Finding prime numbers is not just a trivial math concept; it’s important in many contexts, particularly in the field of cryptography.

In case you need a quick refresher, a prime number is one that cannot be divided by any other whole number (save for the same number itself or the number 1) without leaving a remainder. The number 17, for example, is prime because every number we try dividing it by will leave a remainder. The fact that 17 can be divided by 1 and 17 doesn’t rule 17 out from being prime, since every number can be divided by itself or the number 1.

Another way of defining a prime number is that a prime number is one in which no two whole numbers (that are both greater than 1) multiplied together produce that number.

The opposite of a prime number is a composite number. A composite number can be divided by another number (that isn’t 1) without leaving a remainder. The number 15 is a composite number because it can be divided by 3 or 5 without leaving a remainder. Said another way: because we can find two whole numbers that multiply together to produce 15, namely the numbers 3 and 5, the number 15 is composite.

Perhaps the most straightforward way to test whether a number is prime is by using brute force. That is, we’ll try dividing the number by every other smaller whole number. So, if our number is 11, we’ll first try dividing by 2, and then by 3, and then by 4, and so on. As soon as we find a division that doesn’t produce a remainder, we would conclude that 11 is composite. However, if we try all the possible divisions and can’t get a remainder (which indeed is the case with the number 11), we’ll know that the number 11 is prime.

In Volume 1, Chapter 3, we wrote a function that does this. Here it is:

	def is_prime(number):
	for i in range(2, number):
	if number % i == 0:
	return False

	return True

Here, our function accepts the number we’re testing for primality. We run a loop that divides the number by every integer i from 2 up until the number itself. If we find one such division that leaves no remainder, we return False to indicate that the number is not prime. However, if we get through the entire loop without finding a “remainderless” division, we return True since the number must be prime.

Improving Brute-Force Primality Testing

We can optimize our is_prime function in a few ways. As things stand now, we take N steps to run this primality test, if we consider number to be N. That could take a long time if we’re testing whether the number 20,988,936,657,440,586,486,151,264,256,610,222,593,863,921 is prime! (It is.)

One simple trick we can do to speed things up a bit is to first check if number is even. If it is, number is certainly composite since it can be divided by 2. (The exception is if number is itself 2; then it’s prime.)

With even numbers out of the way, our loop only has to divide number by odd integers. Here’s the code:

	def is_prime(number):
	if number == 2:
	return True

	if number % 2 == 0:
	return False

	i = 3
	while i < number:
	if number % i == 0:
	return False

	i += 2

	return True

With this easy trick of skipping even numbers, we’ve effectively shaved off half our search time.

Now, if we think about things a bit more (it’s always good to think about things), we’ll find that we can speed up our primality test even further. To help paint the picture, I’m going to refer to some division jargon—the type of stuff we learned about in fourth grade but subsequently forgot:

In the division operation 36 / 9 = 4, 36 is the dividend, 9 is the divisor, and 4 is the quotient.

An important property of division is that if 36 / 9 = 4, then it’s also true that 36 / 4 = 9. This is because 4 * 9 and 9 * 4 both equal 36. In other words, a division equation holds true even after we swap the divisor with the quotient.

Now, let’s say that I want to test whether the number 37 is prime. Let me begin:

	Dividend Divisor Quotient Remainder
	37 / 2 = 18 1
	37 / 3 = 12 1
	37 / 4 = 9 1
	37 / 5 = 7 2
	37 / 6 = 6 1

Here, I tried dividing 37 by the numbers 2 through 6. The quotients range from 18 down to 6. Now, watch what happens when I divide 37 by 7 and 8:

	Dividend Divisor Quotient Remainder
	37 / 7 = 5 2
	37 / 8 = 4 5

Note that at this point, the quotients we’re getting now are the same as the divisors from our first bunch of division computations. What this means is that I don’t have to bother dividing 37 by 7 or 8 to see whether 37 is prime. After trying to divide 37 by 6, from here on in, I’m only going to get quotients that lie in the range of divisors that I already tried previously. Accordingly, I know that if I couldn’t find any remainderless quotients earlier, there’s no way I’m going to find any remainderless quotients now. And so I can conclude that 37 is prime.

This inflection point occurs when the divisor reaches the square root of the dividend. As a reminder, because 6 squared is 36, we say that 6 is the square root of 36.

When I divide 36 by 2 and 3 and so on, the quotients are still new numbers we haven’t encountered before as divisors. This changes after I divide 36 by its square root, which is 6. Henceforth, all the quotients we’ll get by dividing by 7 and on will all lie in the range of divisors we’ve already tried before. The same applies to the 37 example since the 6 is the approximate square root of 37.

With this all in mind, we can now optimize the is_prime function to only run its loop up until the square root of number:

	import math


	def is_prime(number):
	if number == 2:
	return True

	if number % 2 == 0:
	return False

	i = 3
	while i <= math.sqrt(number):
	if number % i == 0:
	return False

	i += 2

	return True

This algorithm belongs to a category of Big O notation we haven’t encountered before in this book: O(sqrt(N)). An algorithm is described as O(sqrt(N)) when the algorithm takes sqrt(N) steps when there are N elements of data. In our case, the algorithm takes only half that number of steps since we’re only dealing with odd numbers. Although this is truly sqrt(N)/2 steps, we drop the constant, leaving us with O(sqrt(N)).

These are some pretty clever optimizations, so we should definitely pat ourselves on the back. However, even this optimized algorithm is still no match for massive numbers. The square root of a massive number may be a massive number itself. As such, our primality test will still be impossibly slow.

And that’s where Monte Carlo algorithms come in.