SICP Problem 1.28

In this post, I’ll go through some fun I had last week with problem 1.28 from SICP. The problem statement asks you to implement the Miller-Rabin primality test in Scheme. I had not read this book and picked it up because of how many times I’ve seen it recommended. So far, I have to say it has been a fun ride; I’m enjoying it a lot!

I started with that implementation, but then I sort of got into implementing my own ‘prime?’ function (a function that determines whether an integer is a prime number or not). I decided to limit it to implement a function that works for numbers up to 1 billion. This restriction is entirely arbitrary, but it’s big enough to have fun with it.

Important disclaimer: the code in this article is a fun little exploration. If you need to test primality in real life, I recommend looking at this paper before deciding on how to go with your implementation.

Let’s first recall the Miller-Rabin primality test. Given an integer odd \(n > 3\), you can pick a value \(a\) such that \(2 \leq a < n\). Since \(n\) is odd, we know that $n-1$ is even, and therefore can be decomposed into \(n – 1 = d \times 2^r\). The test checks if \(a^d \mod n = 1\) or if \(a^{d2^i} mod n = -1\) for \(i \in \{0, \ldots, r – 1\}\).

If the conditions hold, we think the number may be a prime; otherwise, we know it’s not a prime. If we test enough values of \(a\), we can be sure whether a number is prime or not. Sadly, it is not practical to check too many values of \(a\), so in practice, you test a few and assume some probability of being wrong.

From the method above, the first thing we need to do is decompose \(n-1\) into \(d\times 2^r\). We will start with that function:

; We decompose n (assuming n is even) as (d r) where n = d \times 2^r.
(define (decompose n r)
  (cond
    ((= 1 (remainder n 2)) (list n r))
    (else (decompose (/ n 2) (+ 1 r)))))

Now, we need to be able to compute \(a^d \mod n\), for this, we can write:

(define (mod-exp base exp n)
  (cond
    ((= exp 0) 1)
    ((= 1 (remainder exp 2))
     (remainder (* base (mod-exp base (- exp 1))) n))
    (else (remainder (sqr (mod-exp base (/ exp 2))) n))))

We would first compute \(a^d mod n\), test whether it’s 1 or -1, if it is, we say the number is likely to be prime. Otherwise, we proceed to the other conditions, this function for testing looks like this (I’m abusing notation a bit writing n-1):

; Tests whether squaring x r times ever hits n-1 (or -1 mod n).
(define (test-squaring x r n n-1)
  (cond
    ((= r 1) false)
    (else (let ((x-squared (remainder (sqr x) n)))
            (if (= x-squared n-1)
                true
                (test-squaring x-squared (- r 1) n n-1))))))

Now, let’s put this together into one function:

; We assume n is greater than 2 and odd.
(define (miller-rabin-test n witnesses)
  ; We decompose n (assuming n is even) as (r d) where n = r \times 2^d.
  (define (decompose n r)
    (cond
      ((= 1 (remainder n 2)) (list n r))
      (else (decompose (/ n 2) (+ 1 r)))))

  (let (
        ; We need this a few times.
        (n-1 (- n 1))
        ; We compute a list (d r) such that n-1 = d\times 2^r.
        (power-decomposition (decompose (- n 1) 0)))
    ; Compute base^exp mod n.
    (define (mod-exp base exp)
      (cond
          ((= exp 0) 1)
          ((= 1 (remainder exp 2))
           (remainder (* base (mod-exp base (- exp 1))) n))
          (else (remainder (sqr (mod-exp base (/ exp 2))) n))))
    ; Tests whether squaring x r times ever hits n-1 (or -1 mod n).
    (define (test-squaring x r)
      (cond
        ((= r 1) false)
        (else (let ((x-squared (remainder (sqr x) n)))
                (if (= x-squared n-1)
                    true
                    (test-squaring x-squared (- r 1)))))))
    ; This test whether a^d mod n is 1 or -1 mod n.
    ; If not, then we need to test the squarings of this number up to a^(n-1)
    ; (i.e., r times).
    (define (test-d a)
      (let ((a-pow-d-mod-n (mod-exp a (first power-decomposition))))
        (cond
          ((or (= a-pow-d-mod-n  1) (= a-pow-d-mod-n n-1)) true)
          (else
           (test-squaring a-pow-d-mod-n (second power-decomposition))))))
    ; Our test is whether all witnesses pass.
    (andmap test-d (filter (lambda (x) (< x n)) witnesses))))

If the conditions hold, then we still don’t know, but the chances of the number being prime are high. We can test different values of \(a\) to increase our certainty. If we examine all possible values of \(a\), we will be sure. However, this is not that interesting, since the process of testing all values of \(a\) is slower than just checking all divisors.

We can take our first try at writing our own ‘prime?’ function:

; My own implementation of 'prime?'
(define (my-prime? n)
  (let (
        ; Our potential whitnesses are 2, 3, and 5. If they test the number as prime,
        ; we assume the number is prime.
        (witnesses '(2 3 5)))
    (cond
      ((= 1 n) false)
      ((= 2 n) true)
      ((= 0 (remainder n 2)) false)
      (else (miller-rabin-test n witnesses)))))

We know this won’t always work, but how good is it? We can use the standard ‘prime?’ function to test and get exceptions:

; Checks 'my-prime?' vs the standard's library 'prime?' function up to n.
(define (check-prime n)
  (define (check-prime-it n limit)
    (let ((v n))
      (cond
        ((= n limit) true)
        (else (and
               (if (eq? (prime? v) (my-prime? v)) true
                   (begin (print v) false))
               (check-prime-it (+ n 1) limit))))))
  (check-prime-it 1 n))

When we call this with \(n = 1,000,000,000\) we get:

> (check-prime 1000000000)
25326001#f

Great, let’s remove that bad case:

; My own implementation of 'prime?'
(define (my-prime? n)
  (let (
        ; Our potential whitnesses are 2, 3, and 5. If they test the number as prime,
        ; we assume the number is prime.
        (witnesses '(2 3 5)))
    (cond
      ((= 1 n) false)
      ((= 2 n) true)
      ((= 25326001 n) false)
      ((= 0 (remainder n 2)) false)
      (else (miller-rabin-test n witnesses)))))

We run this again:

> (check-prime 1000000000)
161304001#f

We could keep doing this, but luckily, there is a nice paper (see page 10) that did this already, so we know there are only three bad examples for our selection of witnesses:

• 25326001
• 161304001
• 960946321

I still ran the last one to gain confidence in my implementation of the test. After adding the three exceptions, our function works nicely up to 1 billion. From the paper and the Wikipedia page, we know that for up to 1 billion, the recommendation seems to be using three witnesses. Can we stretch this? What if we pick only 2, but end up with a reasonable list of exceptions? This is precisely the next thing we will do. We start by defining a function that given a limit for \(n\) and a list of witnesses, tells us all the witnesses that would fail our test (we would call primes when they are not).

(define (exceptions n witnesses)
  (for/list 
    ([i (in-range 11 n 2)] 
     #:when (not (eq? (prime? i) (miller-rabin-test i witnesses)))) 
    i))

We explore numbers from 11 to n, jumping by 2 since we don’t need to explore smaller numbers or even numbers. Using this function, I built the list of exceptions for \(a \in \{2,\ldots,7\}\).

I want to find the pair of values for \(a\) that minimize the number of exceptions. This is the pair of exceptions whose intersection is smaller. The code for intersecting the lists is:

(define (intersect-merge a b)
  (cond
    ((empty? a) '())
    ((empty? b) '())
    (else (let
              ((fst-a (first a))
               (fst-b (first b)))
            (cond
              ((= fst-a fst-b) (cons fst-a (intersect-merge (rest a) (rest b))))
              ((< fst-a fst-b) (intersect-merge (rest a) b))
              (else (intersect-merge a (rest b))))))))

After defining the list of exceptions for each possible \(a \in \{2,\ldots,7\}\), we do the following:

(define exceptions-nr 
  (list known-exceptions-2 
        known-exceptions-3 
        known-exceptions-4 
        known-exceptions-5 
        known-exceptions-6 
        known-exceptions-7))

(define (zip l1 l2)
  (map list l1 l2))

(define (intersection-length a b)
  (length (intersect-merge a b)))

(define (show-combination l1 i1 l2 i2)
  (display i1) 
  (display " ") 
  (display i2) 
  (display " ") 
  (display (intersection-length l1 l2)) 
  (newline))

(define (combine l l1 i1 index)
  (map 
    (lambda (x) (show-combination l1 i1 (first x) (+ (second x) index))) 
    (zip l (range (length l)))))

(define (explore l s)
  (cond
    ((empty? l) empty)
    (else (combine (rest l) (first l) s (+ 1 s)))))

(explore exceptions-nr 2)

Where each value in exceptions-nr is calculated using the exceptions function:

(define known-exceptions-2 (exceptions 1000000000 '(2)))

This tells us that the best pair, among the ones we tested, is \(a \in \{2, 7\}\). There are 43 exceptions that we would have to write down to avoid doing one extra test. It may be that if we explore further, a better pair is available.

And this is where I stop, now into chapter 2.