SICP Exercise 1.24: The Fermat Test

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.24 asks us to once again modify the timed-prime-test procedure from exercise 1.22, this time using fast-prime? (which uses the Fermat test), then test the improvement using the timing statistics we gathered before.

The Fermat test is a probabilistic method for testing primality, meaning that if a given number passes the Fermat test, the best we can say is that the number is very probably prime. It's based on Fermat's Little Theorem which states:

If n is a prime number and a is any positive integer less than n, then a raised to the nth power is congruent to a modulo n.

So to test n for primality, we select a number, a, which is less than n, and raise it to the nth power. We then divide aⁿ by n to get the remainder. If n is prime, then the remainder will be equal to a, the base we selected.

Fermat's Little Theorem says that this property will always hold true if n is prime, but it doesn't say anything about when n isn't prime. If aⁿ mod n is not equal to a, then we know for certain that n is not prime. However, there are values of n and a that will pass the Fermat test even when n is not prime. If we test enough values of a for a given n, we can increase our confidence that n is prime. Unfortunately, there are extremely rare values of n called Carmichael numbers that pass the Fermat test for any value of a. There are variations on the Fermat test that cannot be fooled. We'll look at one of those variations, the Miller-Rabin test, in a later exercise.

The complete code listing below shows what modifications were necessary to use the Fermat test in the timed-prime-test procedure. The fast-prime? procedure takes two arguments, n, and times. The first argument is the number to test, the second argument tells the procedure how many random values to test with. The more random values we use, the higher our confidence that n is prime, so we should test a lot of values. I've (rather arbitrarily) chosen to test 100 random values.

One other change that you may notice is the inclusion of a library module in the first line of code. When testing the original code listing from the book, I found that Scheme's primitive random procedure has a limit of 4,294,967,087. This won't do for the values we're testing, so I replaced it with the random-integer procedure, which doesn't have this limitation. I found out about this procedure and library through the Scheme Cookbook.

(require (lib "27.ss" "srfi"))

(define (timed-prime-test n)
   (newline)
   (display n)
   (start-prime-test n (current-inexact-milliseconds)))

(define (start-prime-test n start-time)
   (cond ((fast-prime? n 100)
          (report-prime (- (current-inexact-milliseconds) start-time)))))

(define (report-prime elapsed-time)
   (display " *** ")
   (display elapsed-time))

(define (square x)
   (* x x))

(define (even? n)
   (= (remainder n 2) 0))

(define (expmod base exp m)
   (cond ((= exp 0) 1)
         ((even? exp)
          (remainder (square (expmod base (/ exp 2) m))
                     m))
         (else
          (remainder (* base (expmod base (- exp 1) m))
                     m))))

(define (fermat-test n)
   (define (try-it a)
     (= (expmod a n n) a))
   (try-it (+ 1 (random-integer (- n 1)))))

(define (fast-prime? n times)
   (cond ((= times 0) true)
         ((fermat-test n) (fast-prime? n (- times 1)))
         (else false)))

Test Results

We need to run timed-prime-test again with the new modifications using the same set of inputs that we found in the previous exercises, then compare the results to see how much we improved the run time. The following table compares the original data we gathered in exercise 1.22, the improved algorithm we used in exercise 1.23 (improved time column), and the new results we got using the Fermat test (all three time columns are in milliseconds). The ratio column is the ratio of the improved values from exercise 1.23 to the fermat time column.

prime	original time	improved time	fermat time	ratio
10000000019	141	103	178	0.58
10000000033	172	101	187	0.54
10000000061	154	112	173	0.65
100000000003	516	310	181	1.71
100000000019	493	295	186	1.59
100000000057	527	305	189	1.61
1000000000039	1627	861	189	4.56
1000000000061	1559	884	195	4.53
1000000000063	1549	873	197	4.43
10000000000037	5014	2671	208	12.84
10000000000051	4932	2668	211	12.64
10000000000099	4855	2697	208	12.97
100000000000031	15745	8531	233	36.61
100000000000067	16022	8264	231	35.77
100000000000097	15861	8530	225	37.91
1000000000000037	48950	26346	244	108.0
1000000000000091	48836	26210	255	102.8
1000000000000159	49008	26256	257	102.2

Analysis

The timing numbers from the Fermat test start out looking pretty poor compared to what we've seen previously. This has mostly to do with the completely arbitrary number of random values I chose to test each prime with. As we increase the size of the numbers we're testing, you can see that the time using the Fermat test barely increases at all. We can verify that the time increase is logarithmic by observing that as the numbers under test increase by a factor of 10, the ratio column increases by a factor of roughly three. This logarithmic characteristic is due to the fact that the expmod procedure used in fermat-test has a logarithmic time complexity.


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.