SICP Exercise 1.28: The Miller-Rabin Test

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 procedure. 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 procedure to signal if it discovers a nontrivial square root of 1, and use this to implement the Miller-Rabin test with a procedure analogous to fermat-test. Check your procedure by testing various known primes and non-primes. Hint: One convenient way to make expmod signal is to have it return 0.

Let's start with figuring out how to test whether or not a number is a non-trivial square root of 1 modulo n. We want to produce a predicate that returns true if, and only if, the number passed to it is not 1 or n - 1 and, when you square it and modulo the result by n you get 1. This is pretty straightforward:

(define (non-trivial-sqrt? x n)
  (and (not (= x 1))
       (not (= x (- n 1)))
       (= (remainder (square x) n) 1)))

Now we're given the hint that "one convenient way to make expmod signal [if it's found a non-trival square root of 1 modulo n] is to have [expmod] return 0." How can we do that? First note that we want to apply the non-trivial-sqrt? test to the result of the recursive call, so we need to put our test in between the recursive call and the subsequent squaring of the result of the recursive call and application of modulo m to that. What happens when you square 0? You get 0. What happens when you calculate 0 modulo m? You get 0 for any m. So if the non-trivial-sqrt? test fails then we can return 0 instead of the value of the recursive call, and this value will then propagate back up the chain of calls and become the final result of the expmod call.

To do this we introduce a procedure that performs the test and returns 0 if it's found a non-trivial square root of 1 modulo n, and just returns the number tested otherwise:

(define (zero-if-non-trivial-sqrt x n)
  (if (non-trivial-sqrt? x n) 0 x))

Then we can modify expmod, simply wrapping the recursive call in a call to this. I renamed this version of expmod to mrt-expmod, as it no longer performs a true calculation of base^exp mod m:

(define (mrt-expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp) (remainder
                       (square
                         (zero-if-non-trivial-sqrt
                           (mrt-expmod base (/ exp 2) m) m)) m))
        (else (remainder (* base (mrt-expmod base (- exp 1) m)) m))))

Finally, we can wrap all this up in a test:

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

And then try it out:

> (miller-rabin-test 561)
#f
> (miller-rabin-test 1105)
#f
> (miller-rabin-test 1729)
#f
> (miller-rabin-test 2465)
#f
> (miller-rabin-test 2821)
#f
> (miller-rabin-test 6601)
#f

Of course, this is just testing a random value between 1 and n. To be truly sure that the Miller-Rabin test is not defeated by the Charmichael numbers we tested in exercise 1.27, let's check all of the numbers in the range 1 to n using:

(define (miller-rabin-test-all n)
  (define (try-it a)
    (cond ((= a n) #t)
          ((= (mrt-expmod a (- n 1) n) 1) (try-it (+ a 1)))
          (else #f)))
  (try-it 1))

And here are the results of our testing:

> (miller-rabin-test-all 561)
#f
> (miller-rabin-test-all 1105)
#f
> (miller-rabin-test-all 1729)
#f
> (miller-rabin-test-all 2465)
#f
> (miller-rabin-test-all 2821)
#f
> (miller-rabin-test-all 6601)
#f

And just for the sake of sanity:

> (miller-rabin-test-all 1)
#t
> (miller-rabin-test-all 2)
#t
> (miller-rabin-test-all 3)
#t
> (miller-rabin-test-all 4)
#f
> (miller-rabin-test-all 5)
#t
> (miller-rabin-test-all 6)
#f
> (miller-rabin-test-all 7)
#t
> (miller-rabin-test-all 8)
#f
> (miller-rabin-test-all 9)
#f