1.2.6 Example: Testing for Primality - SICP Comparison Edition

Permalink copied!

1.2.6

This section describes two methods for checking the primality of an integer $n$, one with order of growth $\Theta(\sqrt{n})$, and a probabilistic algorithm with order of growth $\Theta(\log n)$. The exercises at the end of this section suggest programming projects based on these algorithms.

Searching for divisors

Since ancient times, mathematicians have been fascinated by problems concerning prime numbers, and many people have worked on the problem of determining ways to test if numbers are prime. One way to test if a number is prime is to find the number's divisors. The following program finds the smallest integral divisor (greater than 1) of a given number $n$. It does this in a straightforward way, by testing $n$ for divisibility by successive integers starting with 2.

function smallest_divisor(n) {
    return find_divisor(n, 2);
}
function find_divisor(n, test_divisor) {
    return square(test_divisor) > n
           ? n
           : divides(test_divisor, n)
           ? test_divisor
           : find_divisor(n, test_divisor + 1);
}
function divides(a, b) {
    return b % a === 0;
}

We can test whether a number is prime as follows: $n$ is prime if and only if $n$ is its own smallest divisor.

function is_prime(n) {
    return n === smallest_divisor(n);
}

The end test for find_divisor is based on the fact that if $n$ is not prime it must have a divisor less than or equal to $\sqrt{n}$.

[1]

This means that the algorithm need only test divisors between 1 and $\sqrt{n}$. Consequently, the number of steps required to identify $n$ as prime will have order of growth $\Theta(\sqrt{n})$.

The Fermat test

The $\Theta(\log n)$ primality test is based on a result from number theory known as Fermat's Little Theorem.

[2]

Fermat's Little Theorem: If $n$ is a prime number and $a$ is any positive integer less than $n$, then $a$ raised to the $n$th power is congruent to $a$ modulo $n$.
フェルマーの小定理： $n$ が素数で、$a$ が $n$ より小さい任意の正の整数であるとき、$a$ の $n$ 乗は $a$ と法 $n$ で合同です。

(Two numbers are said to be congruent modulo $n$ if they both have the same remainder when divided by $n$. The remainder of a number $a$ when divided by $n$ is also referred to as the remainder of $a$ modulo $n$, or simply as $a$ modulo $n$.)

If $n$ is not prime, then, in general, most of the numbers $a < n$ will not satisfy the above relation. This leads to the following algorithm for testing primality: Given a number $n$, pick a random number $a < n$ and compute the remainder of $a^n$ modulo $n$. If the result is not equal to $a$, then $n$ is certainly not prime. If it is $a$, then chances are good that $n$ is prime. Now pick another random number $a$ and test it with the same method. If it also satisfies the equation, then we can be even more confident that $n$ is prime. By trying more and more values of $a$, we can increase our confidence in the result. This algorithm is known as the Fermat test.

To implement the Fermat test, we need a function that computes the exponential of a number modulo another number:

function expmod(base, exp, m) {
    return exp === 0
           ? 1
           : is_even(exp)
           ? square(expmod(base, exp / 2, m)) % m
           : (base * expmod(base, exp - 1, m)) % m;
}

This is very similar to the fast_expt function of section 1.2.4. It uses successive squaring, so that the number of steps grows logarithmically with the exponent.

[3]

The Fermat test is performed by choosing at random a number $a$ between 1 and $n-1$ inclusive and checking whether the remainder modulo $n$ of the $n$th power of $a$ is equal to $a$. The random number $a$ is chosen using the primitive function math_random, which returns a nonnegative number less than 1. Hence, to obtain a random number between 1 and $n-1$, we multiply the return value of math_random by $n-1$, round down the result with the primitive function math_floor, and add 1:

function fermat_test(n) {
    function try_it(a) {
        return expmod(a, n, n) === a;
    }
    return try_it(1 + math_floor(math_random() * (n - 1)));
}

The following function runs the test a given number of times, as specified by a parameter. Its value is true if the test succeeds every time, and false otherwise.

function fast_is_prime(n, times) {
    return times === 0
           ? true
           : fermat_test(n)
           ? fast_is_prime(n, times - 1)
           : false;
}

Probabilistic methods

The Fermat test differs in character from most familiar algorithms, in which one computes an answer that is guaranteed to be correct. Here, the answer obtained is only probably correct. More precisely, if $n$ ever fails the Fermat test, we can be certain that $n$ is not prime. But the fact that $n$ passes the test, while an extremely strong indication, is still not a guarantee that $n$ is prime. What we would like to say is that for any number $n$, if we perform the test enough times and find that $n$ always passes the test, then the probability of error in our primality test can be made as small as we like.

Unfortunately, this assertion is not quite correct. There do exist numbers that fool the Fermat test: numbers $n$ that are not prime and yet have the property that $a^n$ is congruent to $a$ modulo $n$ for all integers $a < n$. Such numbers are extremely rare, so the Fermat test is quite reliable in practice.

[4]

There are variations of the Fermat test that cannot be fooled. In these tests, as with the Fermat method, one tests the primality of an integer $n$ by choosing a random integer $a < n$ and checking some condition that depends upon $n$ and $a$. (See exercise 1.28 for an example of such a test.) On the other hand, in contrast to the Fermat test, one can prove that, for any $n$, the condition does not hold for most of the integers $a < n$ unless $n$ is prime. Thus, if $n$ passes the test for some random choice of $a$, the chances are better than even that $n$ is prime. If $n$ passes the test for two random choices of $a$, the chances are better than 3 out of 4 that $n$ is prime. By running the test with more and more randomly chosen values of $a$ we can make the probability of error as small as we like.

The existence of tests for which one can prove that the chance of error becomes arbitrarily small has sparked interest in algorithms of this type, which have come to be known as probabilistic algorithms. There is a great deal of research activity in this area, and probabilistic algorithms have been fruitfully applied to many fields.

[5]

Exercise 1.21

Use the smallest_divisor function to find the smallest divisor of each of the following numbers: 199, 1999, 19999.

smallest_divisor(199);
// smallest_divisor(1999);
// smallest_divisor(19999);

Exercise 1.22

Assume a primitive function get_time of no arguments that returns the number of milliseconds that have passed since 00:00:00 UTC on Thursday, 1 January, 1970.

[6]

The following timed_prime_test function, when called with an integer $n$, prints $n$ and checks to see if $n$ is prime. If $n$ is prime, the function prints three asterisks

[7]

followed by the amount of time used in performing the test.

function timed_prime_test(n) {
    display(n);
    return start_prime_test(n, get_time());
}

function start_prime_test(n, start_time) {
    return is_prime(n)
           ? report_prime(get_time() - start_time)
           : false;
}

function report_prime(elapsed_time) {
    display(" *** ");
    display(elapsed_time);
    return true;
}

Using this function, write a function search_for_primes that checks the primality of consecutive odd integers in a specified range. Use your function to find the three smallest primes larger than 1000; larger than 10,000; larger than 100,000; larger than 1,000,000. Note the time needed to test each prime. Since the testing algorithm has order of growth of $\Theta(\sqrt{n})$, you should expect that testing for primes around 10,000 should take about $\sqrt{10}$ times as long as testing for primes around 1000. Do your timing data bear this out? How well do the data for 100,000 and 1,000,000 support the $\sqrt{n}$ prediction? Is your result compatible with the notion that programs on your machine run in time proportional to the number of steps required for the computation?

function search_for_primes(start, times) {
    return times === 0
           ? true
           : start > 2 && start % 2 === 0
           ? search_for_primes(start + 1, times)
           // if we get true, it's a prime
           : timed_prime_test(start)
           ? search_for_primes(start + 2, times - 1)
           : search_for_primes(start + 2, times);
}

The timing data quite clearly supports the $\sqrt{n}$ for prediction that is sufficiently big, such as 100,000 and 1,000,000.

Exercise 1.23

The smallest_divisor function shown at the start of this section does lots of needless testing: After it checks to see if the number is divisible by 2 there is no point in checking to see if it is divisible by any larger even numbers. This suggests that the values used for test_divisor should not be 2, 3, 4, 5, 6, … but rather 2, 3, 5, 7, 9, …. To implement this change, declare a function next that returns 3 if its input is equal to 2 and otherwise returns its input plus 2. Modify the smallest_divisor function to use next(test_divisor) instead of test_divisor + 1. With timed_prime_test incorporating this modified version of smallest_divisor, run the test for each of the 12 primes found in exercise 1.22. Since this modification halves the number of test steps, you should expect it to run about twice as fast. Is this expectation confirmed? If not, what is the observed ratio of the speeds of the two algorithms, and how do you explain the fact that it is different from 2?

function next(input) {
    return input === 2 
           ? 3 
           : input + 2;
}

function find_divisor(n, test_divisor) {
    return square(test_divisor) > n
           ? n
           : divides(test_divisor, n)
           ? test_divisor
           : find_divisor(n, next(test_divisor));
}

The ratio of the speeds of the two algorithms is not exactly 2, but this might be due to hardware / network issues. It is about 1.5 times faster compared to previous solution.

Exercise 1.24

Modify the timed_prime_test function of exercise 1.22 to use fast_is_prime (the Fermat method), and test each of the 12 primes you found in that exercise. Since the Fermat test has $\Theta(\log n)$ growth, how would you expect the time to test primes near 1,000,000 to compare with the time needed to test primes near 1000? Do your data bear this out? Can you explain any discrepancy you find?

function timed_prime_test(n) {
    display(n);
    return start_prime_test(n, get_time());
}
function start_prime_test(n, start_time) {
    return fast_is_prime(n, math_floor(math_log(n)))
           ? report_prime(get_time() - start_time)
           : true;
}
function report_prime(elapsed_time) {
    display(" *** ");
    display(elapsed_time);
}

The time to test primes near 1,000,000 using fast_is_prime is about 4 ms, 4 times the time needed to test primes near 1,000. This is faster compared to 8 ms that we achieved if we use the is_prime. However, despite being 4 times slower, this fact cannot lead us to believe that it has a greater growth than $\Theta(\log n)$, as it should be tested with greater numbers to gain a more accurate understanding of the growth of the function.

Exercise 1.25

Alyssa P. Hacker complains that we went to a lot of extra work in writing expmod. After all, she says, since we already know how to compute exponentials, we could have simply written

function expmod(base, exp, m) {
    return fast_expt(base, exp) % m;
}

Is she correct? Would this function serve as well for our fast prime tester? Explain.

Alyssa's suggestion is correct at first sight: her expmod function computes $\textit{base}^{\textit{exp}}$ and then finds its remainder modulo $m$, as required in the Fermat test.

However, for large bases, Alyssa's method will quickly bump into limitations because JavaScript uses 64 bits to represent numbers, following the double-precision floating point standard. When the numbers become so large that they cannot be represented precisely any longer in this standard, the results become unreliable. Even worse, the method might exceed the largest number that can be represented in this standard, and the computation leads to an error.

For small bases, however, Alyssa's method may be even faster than the original expmod function, because it will carry out only one single remainder operation.

Exercise 1.26

Louis Reasoner is having great difficulty doing exercise 1.24. His fast_is_prime test seems to run more slowly than his is_prime test. Louis calls his friend Eva Lu Ator over to help. When they examine Louis's code, they find that he has rewritten the expmod function to use an explicit multiplication, rather than calling square:

function expmod(base, exp, m) {
    return exp === 0
           ? 1
           : is_even(exp)
           ? (  expmod(base, exp / 2, m)
              * expmod(base, exp / 2, m)) % m
           : (base * expmod(base, exp - 1, m)) % m;
}

I don't see what difference that could make, says Louis. I do. says Eva. By writing the function like that, you have transformed the $\Theta(\log n)$ process into a $\Theta(n)$ process. Explain.

Eva is correct: by evaluating the expression:

(expmod(base, exp / 2, m) * expmod(base, exp / 2, m)) % m

the expression expmod(base, exp / 2, m) is evaluated twice at each step in the computation when the exponent is even, eliminating the benefit of the fast exponentiation algorithm—which halves the exponent when the exponent is even—therefore eliminating the feature of the algorithm that makes it faster.

Exercise 1.27

Demonstrate that the Carmichael numbers listed in footnote 4 really do fool the Fermat test. That is, write a function that takes an integer $n$ and tests whether $a^n$ is congruent to $a$ modulo $n$ for every $a < n$, and try your function on the given Carmichael numbers.

function carmichael(n) {
    function expmod(base, exp, m) {
        return exp === 0
               ? 1
               : is_even(exp)
               ? square(expmod(base, exp / 2, m)) % m
               : (base * expmod(base, exp - 1, m)) % m;
    }
    function fermat_test(n, a) {
        return expmod(a, n, n) === a;
    }
    function iter(n, i) {
        return i === n
               ? true
               : fermat_test(n, i)
               ? iter(n, i + 1)
               : false;
    }
    return iter(n, 2);
}

Exercise 1.28

One variant of the Fermat test that cannot be fooled is called the Miller–Rabin test (Miller 1976; Rabin 1980). This starts from an alternate form of Fermat's Little Theorem, which states that if $n$ is a prime number and $a$ is any positive integer less than $n$, then $a$ raised to the $(n-1)$st power is congruent to 1 modulo $n$. To test the primality of a number $n$ by the Miller–Rabin test, we pick a random number $a < n$ and raise $a$ to the $(n-1)$st power modulo $n$ using the expmod function. However, whenever we perform the squaring step in expmod, we check to see if we have discovered a nontrivial square root of 1 modulo $n$, that is, a number not equal to 1 or $n-1$ whose square is equal to 1 modulo $n$. It is possible to prove that if such a nontrivial square root of 1 exists, then $n$ is not prime. It is also possible to prove that if $n$ is an odd number that is not prime, then, for at least half the numbers $a < n$, computing $a^{n-1}$ in this way will reveal a nontrivial square root of 1 modulo $n$. (This is why the Miller–Rabin test cannot be fooled.) Modify the expmod function to signal if it discovers a nontrivial square root of 1, and use this to implement the Miller–Rabin test with a function analogous to fermat_test. Check your function by testing various known primes and non-primes. Hint: One convenient way to make expmod signal is to have it return 0.

function random(n) {
    return math_floor(math_random() * n);
}
function miller_rabin_test(n) {
    function expmod(base, exp, m) {
        return exp === 0
               ? 1
               : is_even(exp)
               ? square(trivial_test(expmod(base,
                                            exp / 2,
                                            m), 
                                     m))
                 % m
               : (base * expmod(base, exp - 1, m)) 
                 % m;
    }
    function trivial_test(r, m) {
        return r === 1 || r === m - 1
               ? r
               : square(r) % m === 1
               ? 0
               : r;
    }
    function try_it(a) {
        return expmod(a, n - 1, n) === 1;
    }
    return try_it(1 + random(n - 1));
}
function do_miller_rabin_test(n, times) {
    return times === 0
           ? true
           : miller_rabin_test(n)
           ? do_miller_rabin_test(n, times - 1)
           : false;
}

[1]

If $d$ is a divisor of $n$, then so is $n/d$. But $d$ and $n/d$ cannot both be greater than $\sqrt{n}$.

[2]

Pierre de Fermat (1601–1665) is considered to be the founder of modern number theory. He obtained many important number-theoretic results, but he usually announced just the results, without providing his proofs. Fermat's Little Theorem was stated in a letter he wrote in 1640. The first published proof was given by Euler in 1736 (and an earlier, identical proof was discovered in the unpublished manuscripts of Leibniz). The most famous of Fermat's results—known as Fermat's Last Theorem—was jotted down in 1637 in his copy of the book Arithmetic (by the third-century Greek mathematician Diophantus) with the remark I have discovered a truly remarkable proof, but this margin is too small to contain it. Finding a proof of Fermat's Last Theorem became one of the most famous challenges in number theory. A complete solution was finally given in 1995 by Andrew Wiles of Princeton University.

[3]

The reduction steps in the cases where the exponent $e$ is greater than 1 are based on the fact that, for any integers $x$, $y$, and $m$, we can find the remainder of $x$ times $y$ modulo $m$ by computing separately the remainders of $x$ modulo $m$ and $y$ modulo $m$, multiplying these, and then taking the remainder of the result modulo $m$. For instance, in the case where $e$ is even, we compute the remainder of $b^{e/2}$ modulo $m$, square this, and take the remainder modulo $m$. This technique is useful because it means we can perform our computation without ever having to deal with numbers much larger than $m$. (Compare exercise 1.25.)

[4]

Numbers that fool the Fermat test are called Carmichael numbers, and little is known about them other than that they are extremely rare. There are 255 Carmichael numbers below 100,000,000. The smallest few are 561, 1105, 1729, 2465, 2821, and 6601. In testing primality of very large numbers chosen at random, the chance of stumbling upon a value that fools the Fermat test is less than the chance that cosmic radiation will cause the computer to make an error in carrying out a correct algorithm. Considering an algorithm to be inadequate for the first reason but not for the second illustrates the difference between mathematics and engineering.

[5]

One of the most striking applications of probabilistic prime testing has been to the field of cryptography. Although it is computationally infeasible to factor an arbitrary 300-digit number as of this writing (2021), the primality of such a number can be checked in a few seconds with the Fermat test. This fact forms the basis of a technique for constructing unbreakable codes suggested by Rivest, Shamir, and Adleman (1977). The resulting RSA algorithm has become a widely used technique for enhancing the security of electronic communications. Because of this and related developments, the study of prime numbers, once considered the epitome of a topic in pure mathematics to be studied only for its own sake, now turns out to have important practical applications to cryptography, electronic funds transfer, and information retrieval.

[6]

This date is called the UNIX epoch and is part of the specification of functions that deal with time in the UNIX$^{\textrm{TM}}$ operating system.

[7]

The primitive function display returns its argument, but also prints it. Here " *** " is a string, a sequence of characters that we pass as argument to the display function. Section 2.3.1 introduces strings more thoroughly.

< Previous

Next >

1.2.6

Example: Testing for Primality