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

This entry was posted in Programming. Bookmark the permalink.

### 3 Responses to SICP Problem 1.28

1. Dengke Ng says:

SICP is one of the greatest books about “the art of computing” I have read.
But I also haven’t finished all of it yet (for many years). lol

BTW, is webgraphs.recoded.cl down or anything? I am currently working on web graph compression (maybe an old topic for you) and seeking for implementation of Re-Pair to reproduce the results as described in your paper Fast and Compact Web Graph Representations.

• fclaude says:

Hey! The webgraphs page is down at the moment, but I’m trying to recover it too. Will probably have a follow up on the status of that this weekend :).

• fclaude says:

The site is up again! https://webgraphs.recoded.cl

This site uses Akismet to reduce spam. Learn how your comment data is processed.