Why does this Miller-Rabin Procedure in Scheme works when the code seems to be wrong?

Question

I am working through SICP. In exercise 1.28 about the Miller-Rabin test. I had this code, that I know is wrong because it does not follow the instrcuccions of the exercise.

(define (fast-prime? n times)           
  (define (even? x)             
    (= (remainder x 2) 0))          
  (define (miller-rabin-test n)         
    (try-it (+ 1 (random (- n 1)))))        
  (define (try-it a)                
    (= (expmod a (- n 1) n) 1))         
  (define (expmod base exp m)           
    (cond ((= exp 0) 1) 
      ((even? exp)              
       (if (and (not (= exp (- m 1))) (= (remainder (square exp) m) 1)) 
         0
         (remainder (square (expmod base (/ exp 2) m)) m)))
      (else         
        (remainder (* base (expmod base (- exp 1) m)) m))))
  (cond ((= times 0) true)          
    ((miller-rabin-test n) (fast-prime? n (- times 1))) 
    (else false)))              

In it I test if the square of the exponent is congruent to 1 mod n. Which according to what I have read, and other correct implementations I have seen is wrong. I should test the entire number as in:

...
        (square 
         (trivial-test (expmod base (/ exp 2) m) m))
...

The thing is that I have tested this, with many prime numbers and large Carmicheal numbers, and it seems to give the correct answer, though a bit slower. I don't understand why this seems to work.

Answer 1

Your version of the function "works" only because you are lucky. Try this experiment: evaluate (fast-prime? 561 3) a hundred times. Depending on the random witnesses that your function chooses, sometimes it will return true and sometimes it will return false. When I did that I got 12 true and 88 false, but you may get different results, depending on your random number generator.

> (let loop ((k 0) (t 0) (f 0))
    (if (= k 100) (values t f)
      (if (fast-prime? 561 3)
          (loop (+ k 1) (+ t 1) f)
          (loop (+ k 1) t (+ f 1)))))
12
88

I don't have SICP in front of me -- my copy is at home -- but I can tell you the right way to perform a Miller-Rabin primality test.

Your expmod function is incorrect; there is no reason to square the exponent. Here is a proper function to perform modular exponentiation:

(define (expm b e m) ; modular exponentiation
  (let loop ((b b) (e e) (x 1))
    (if (zero? e) x
      (loop (modulo (* b b) m) (quotient e 2)
            (if (odd? e) (modulo (* b x) m) x)))))

Then Gary Miller's strong pseudoprime test, which is a strong version of your try-it test for which there is a witness a that proves the compositeness of every composite n, looks like this:

(define (strong-pseudoprime? n a) ; strong pseudoprime base a
  (let loop ((r 0) (s (- n 1)))
    (if (even? s) (loop (+ r 1) (/ s 2))
      (if (= (expm a s n) 1) #t
        (let loop ((r r) (s s))
          (cond ((zero? r) #f)
                ((= (expm a s n) (- n 1)) #t)
                (else (loop (- r 1) (* s 2)))))))))

Assuming the Extended Riemann Hypothesis, testing every a from 2 to n-1 will prove (an actual, deterministic proof, not just a probabilistic estimate of primality) the primality of a prime n, or identify at least one a that is a witness to the compositeness of a composite n. Michael Rabin proved that if n is composite, at least three-quarters of the a from 2 to n-1 are witnesses to that compositeness, so testing k random bases demonstrates, but does not prove, the primality of a prime n to a probability of 4^−k. Thus, this implementation of the Miller-Rabin primality test:

(define (prime? n k)
  (let loop ((k k))
    (cond ((zero? k) #t)
          ((not (strong-pseudoprime? n (random (+ 2 (- n 3))))) #f)
          (else (loop (- k 1))))))

That always works properly:

> (let loop ((k 0) (t 0) (f 0))
    (if (= k 100) (values t f)
      (if (prime? 561 3)
          (loop (+ k 1) (+ t 1) f)
          (loop (+ k 1) t (+ f 1)))))
0
100

I know your purpose is to study SICP rather than to program primality tests, but if you're interested in programming with prime numbers, I modestly recommend this essay at my blog, which discusses the Miller-Rabin test, among other topics. You should also know there are better (faster, less likely to report erroneous result) primality tests available than randomized Miller-Rabin.

Answer 2

It seems to me, you still got correct answer, because in each iteration of expmod you check conditions for previous iteration. You could try to debug exp value using display function inside expmod. Really, your code is not very different from this one.