# 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 < 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 < 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:

1. 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).

2. 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:

1. 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!

2. 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)))))

• I don't understand correctness-related feedback 2. You said I should "clamp the result of square to mod m". Can you explain what you mean by clamp? And you said "All values returned by expmod should already be mod m so doing a further remainder call is pointless.", but isn't calling the remainder exactly why the return value is mod m? – morbidCode Jun 3 '16 at 4:07
• @morbidCode Clamping x to mod m means (modulo x m), so that the result is always between 0 and m-1. The return value of expmod should, by contract, always be clamped to mod m, so the remainder call in your (remainder (expmod (square base) (/ exp 2) m) m) is redundant. – Chris Jester-Young Jun 3 '16 at 4:11
• @morbidCode "You should, however, clamp the result of square to mod m" means use (expmod (remainder (square base) m) (/ exp 2) m) instead. Or, as in my version, (remainder (square (expmod base (/ exp 2) m)) m). Notice how, in both cases, the shape of the expression is (remainder (square ...)), not (remainder (expmod ...)). – Chris Jester-Young Jun 3 '16 at 4:12
• I used (square base) in my solution so I can check base mod 1 inside the procedure. How does (remainder (square (expmod base (/ exp 2) m)) m) gave the correct result? Note I haven't seen your code yet. I'll look at it when I've improved mine. – morbidCode Jun 3 '16 at 4:33
• @morbidCode (expmod base (/ exp 2) m) is the recursive case. You then square the result to get what you need, then you mod m so you can compare the result to 1. – Chris Jester-Young Jun 3 '16 at 4:39