2
\$\begingroup\$

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?

\$\endgroup\$

1 Answer 1

3
\$\begingroup\$

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)))))
    
\$\endgroup\$
9
  • \$\begingroup\$ 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? \$\endgroup\$ Jun 3, 2016 at 4:07
  • 1
    \$\begingroup\$ @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. \$\endgroup\$ Jun 3, 2016 at 4:11
  • \$\begingroup\$ @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 ...)). \$\endgroup\$ Jun 3, 2016 at 4:12
  • \$\begingroup\$ 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. \$\endgroup\$ Jun 3, 2016 at 4:33
  • \$\begingroup\$ @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. \$\endgroup\$ Jun 3, 2016 at 4:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.