From SICP section
1.2.6 Example: Testing for PrimalityExercise 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
n 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
n 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 ResultsWe 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 |
AnalysisThe 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.