Francisco Claude-Faust

SICP Problem 1.28

· 4 min read

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:

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.

← all notes