SICP exercise 1.28 - miller-rabin primality test From SICP


  Exercise 1.28: 
  
  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 \mod n\$. 
  
  To test the
  primality of a number \$n\$ by the Miller-Rabin test, we pick a random
  number \$a &lt; n\$ and raise \$a\$ to the \$n-1\$-st power \$\mod 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
\mod n\$,” that is, a number not equal to \$1\$ or \$n-1\$ whose square is
  equal to \$1 \mod 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 &lt; n\$, computing \$a^{n-1}\$ in this way will reveal a
  nontrivial square root of \$1 \mod 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\$.


According to the book this is a probabalistic algorithm. My tests gave mostly correct results.

Here's my code.

;; modified expmod procedure
(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (if (and (not (= base 1)) 
                  (not (= base (- m 1))) 
                  (not (= exp (- m 1))) 
                  (= (remainder base m) 1))
             0
             (remainder 
              (expmod (square base) (/ exp 2) m)
              m)))
        (else
         (remainder 
          (* base (expmod base (- exp 1) m))
          m))))

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

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


How can I improve this code and make it faster?
To ensure I understood the SICP exercise correctly, I wrote my own version of the Miller-Rabin algorithm. You can look at my code and compare against yours.

With this in mind, here's my correctness-related feedback:


Your test expression, (and (not (= base 1)) (not (= base (- m 1))) (not (= exp (- m 1))) (= (remainder base m) 1)), seems to muddle up the relationship between the previous and current values. Really, what you want to find out is whether the previous value is between 2 and (- m 2) (mod m) and the current value is 1 (mod m).
Your (remainder (expmod (square base) (/ exp 2) m) m) is not optimal. All values returned by expmod should already be mod m so doing a further remainder call is pointless. You should, however, clamp the result of square to mod m.


Other style notes:


Not sure why you didn't just use a one-liner for the miller-rabin-test function:

(define (miller-rabin-test n)
  (= (expmod (+ 1 (random (- n 1))) (- n 1) n) 1))


or, if you're trying to break up the expression a bit,

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


No need for an internal function!
You can rewrite the main loop as a single or/and expression:

(define (fast-prime? n times)
  (or (zero? times)
      (and (miller-rabin-test n)
           (fast-prime? n (- times 1)))))