Thursday, March 4, 2010

SICP Exercise 1.28: The Miller-Rabin test

From SICP section 1.2.6 Example: Testing for Primality

(Note: The explanation of the Miller-Rabin test given in SICP is slightly different from every other source I've read. I don't know if this was done purposefully by the authors for the sake of simplicity, or if the algorithm has simply evolved since SICP's publication. For the sake of consistency, I'm solving the exercise using the explanation given in the text. If anyone is interested, Programming Praxis has a short explanation of the algorithm and an implementation in Scheme. Naturally, you can also read about it on Wikipedia.)

Exercise 1.28 introduces a variation on the Fermat test that isn't fooled by Carmichael numbers called the Miller-Rabin test. We start with an alternate form of Fermat's Little Theorem. In exercise 1.24 we saw that if n is prime and a is any positve integer < n,

an ≡ a (mod n)

If we divide both sides of the congruency by a, we get

an-1 ≡ 1 (mod n)

When we check this congruency using the expmod procedure, at the squaring step we need to also check to see if we've discovered a "non-trivial square root of 1 modulo n," i.e., a number not equal to 1 or (n - 1) whose square is equal to 1 modulo n. It's possible to prove that if a non-trivial square root of 1 exists, then n is not prime. It's also possible to prove that if n is composite, then "for at least half the numbers a < n, computing an-1 in this way will reveal a nontrivial square root of 1 modulo n." (I'm not sure if this is just a very conservative estimate or if it was the best information available at the time that SICP was written. Every other source I've read claims that at least 3/4 of the bases you choose will reveal a nontrivial square root of 1 modulo n.)

We're asked to modify the expmod procedure to check for non-trivial square roots of 1 modulo n, and to use it to implement the Miller-Rabin test with a procedure analogous to the fermat-test that we wrote before. We can check the new procedure by testing various known primes and non-primes.

The first thing we need to do is modify the square procedure so that it checks for a non-trivial square root of 1 modulo n. We'll write a new procedure called square-check that checks for two conditions:
  1. a number not equal to 1 or (n - 1)
  2. whose square is equal to 1 modulo n
If both of these conditions are met, we'll signal that we've found a non-trivial square root by returning 0. If these conditions are not met, we'll perform the remainder and squaring steps that expmod did before. Here's the complete code listing:
(require (lib "" "srfi"))

(define (square-check x m)
(if (and (not (or (= x 1) (= x (- m 1))))
(= (remainder (* x x) m) 1))
(remainder (* x x) m)))

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

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

(define (miller-rabin-test n)
(define (try-it a)
(= (expmod a (- n 1) n) 1))
(try-it (+ 2 (random-integer (- n 2)))))

As with the fermat-test procedure from exercise 1.27, prime numbers should always pass. Composite numbers should have a high likelihood of failing, even the Carmichael numbers (so if you repeatedly test these values enough times, eventually you will see some false positives). You can test several primes and composites in a Scheme interpreter. Here's the output I got when I tested the first six Carmichael numbers.
> (miller-rabin-test 561)
> (miller-rabin-test 1105)
> (miller-rabin-test 1729)
> (miller-rabin-test 2465)
> (miller-rabin-test 2821)
> (miller-rabin-test 6601)

Since the Miller-Rabin test is a refinement of the Fermat test that we saw earlier, it's natural to assume that you can substitute our new procedure for the old one. This is a correct assumption. To upgrade the fast-prime? procedure we saw in exercise 1.24, just change the call to fermat-test into a call to miller-rabin-test as follows:
(define (fast-prime? n times)
(cond ((= times 0) true)
((miller-rabin-test n) (fast-prime? n (- times 1)))
(else false)))


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


kami said...

hi bill, i am a Stack Overflow user, i tried to end one of my questions but no luck, i saw you are the last editer on that question, could you tell me how to do it? the question link is

Bill the Lizard said...

I went ahead and deleted it for you. I'm a moderator on Stack Overflow, so you can click the "flag" link below the tags, then select "Requires moderator attention" from the list that appears whenever a question or answer on that site needs special attention.

kami said...

Thank you very much! also for the SO internal contact method.

Yonatan Erez said...

Hi Bill,

It seems like your miller-rabin-test implementation only tests a single random integer. On SICP, when they implemented the fermat-test, they also made the wrapping (fast-prime?) procedure that runs the fermat-test numerous times on random integers.

I think you should add a similar procedure to the post.

Yonatan Erez said...

BTW, I'm taking a Computer Science course that's based on SICP this semester, and your notes are really helpful for me. Keep up the good work and thanks. :)

Bill the Lizard said...

Thanks Yonatan, I took your suggestion and added the new version of fast-prime? to this post.

Thanks also for telling me about your class. It's good to hear that someone is getting some use out of these notes. I've been pretty busy with a new job, but I'm almost done with my notes for the next lecture and section of the text. I should have those posted soon, and some new exercise solutions to follow.

Ramakrishnan VU3RDD said...

Thanks for the nice discussions on sicp. I am also working through the book and it is very informative to look at the way others have done the exercises.

If the rand-integer generates 0, then base becomes 1. In that case, wouldn't the algorithm give false positive?

Bill the Lizard said...

Ramakrishnan VU3RDD,
Yes, it absolutely would return a false positive in that case. To see it fail, just replace the call to the random number generator in last line of the code listing with the value 1:

(try-it 1))

I've fixed the code listing so it now generates random numbers in the range 2 to (n-1).

Nice catch, and thanks a lot for reading!

Anonymous said...

Hi Bill, in order to be sure that n is prime you have to test more than half the numbers 0 < a < n. Therefore you have to be able to control the "as", as in your current version even if you chose a high value for the times parameter of fast-primes? the random value generated for a in each round could always be the same.

Dan Corneanu said...
This comment has been removed by the author.
Dan Corneanu said...

Hello, very nice discussion on SICP.

I think there might be an error at the beginning of your statement. You say "If we divide both sides of the congruency by a, we get

a^(n-1) ≡ 1 (mod n)"

I do not believe this is generally true. I think that the equivalence between a^n ≡ a (mod n) and a^(n-1) ≡ 1 (mod n) holds true only if n is prime and a does not divide n. From your statement, one could conclude that he/she could actually divide both sides of a congruency (any congruency) by the same number and it will still hold true. I do not think this is the case.

Am I terribly wrong?

Anonymous said...


mt said...

The way Miller-Rabin is described in SICP seems slightly misleading to me. They say " 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." But there's no nontrivial square root of 1 modulo 9 for example. You have to run a Fermat test *and* check for nontrivial square roots at the same time; if the Fermat test fails OR if you find a nontrivial square root then n is composite.

Vishy Ranganath said...

Here is my implementation:

I run it in a loop to test for all numbers less than the number being tested so that we get to see that the Miller-Rabin test really cannot be fooled while the Fermat test can be.

One clarifying question I have: when we test for non-trivial square roots of unity modulo n, we only should be looking in the range between 1 and n-1 right? Otherwise it won't work? For example: 7 is a prime number. 36 is a non-trivial square-root of unity modulo 7. (Because 36*36 = 1296 which when divided by 7 yields a remainder of 1.)