SICP Exercise 1.28: The Miller-Rabin test

From SICP section 1.2.6 Example: Testing for Primality

(Note: The explanation of the Miller-Rabin test given in SICP is slightly different from every other source I've read. I don't know if this was done purposefully by the authors for the sake of simplicity, or if the algorithm has simply evolved since SICP's publication. For the sake of consistency, I'm solving the exercise using the explanation given in the text. If anyone is interested, Programming Praxis has a short explanation of the algorithm and an implementation in Scheme. Naturally, you can also read about it on Wikipedia.)

Exercise 1.28 introduces a variation on the Fermat test that isn't fooled by Carmichael numbers called the Miller-Rabin test. We start with an alternate form of Fermat's Little Theorem. In exercise 1.24 we saw that if n is prime and a is any positve integer < n,

aⁿ ≡ a (mod n)

If we divide both sides of the congruency by a, we get

a^n-1 ≡ 1 (mod n)

When we check this congruency using the expmod procedure, at the squaring step we need to also check to see if we've discovered a "non-trivial square root of 1 modulo n," i.e., a number not equal to 1 or (n - 1) whose square is equal to 1 modulo n. It's possible to prove that if a non-trivial square root of 1 exists, then n is not prime. It's also possible to prove that if n is composite, then "for at least half the numbers a < n, computing a^n-1 in this way will reveal a nontrivial square root of 1 modulo n." (I'm not sure if this is just a very conservative estimate or if it was the best information available at the time that SICP was written. Every other source I've read claims that at least 3/4 of the bases you choose will reveal a nontrivial square root of 1 modulo n.)

We're asked to modify the expmod procedure to check for non-trivial square roots of 1 modulo n, and to use it to implement the Miller-Rabin test with a procedure analogous to the fermat-test that we wrote before. We can check the new procedure by testing various known primes and non-primes.

The first thing we need to do is modify the square procedure so that it checks for a non-trivial square root of 1 modulo n. We'll write a new procedure called square-check that checks for two conditions:

1. a number not equal to 1 or (n - 1)
2. whose square is equal to 1 modulo n

If both of these conditions are met, we'll signal that we've found a non-trivial square root by returning 0. If these conditions are not met, we'll perform the remainder and squaring steps that expmod did before. Here's the complete code listing:

(require (lib "27.ss" "srfi"))

(define (square-check x m)
   (if (and (not (or (= x 1) (= x (- m 1))))
            (= (remainder (* x x) m) 1))
       0
       (remainder (* x x) m)))

(define (even? n)
   (= (remainder n 2) 0))

(define (expmod base exp m)
   (cond ((= exp 0) 1)
         ((even? exp)
          (square-check (expmod 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-integer (- n 1)))))

As with the fermat-test procedure from exercise 1.27, prime numbers should always pass. Composite numbers should have a high likelihood of failing, even the Carmichael numbers (so if you repeatedly test these values enough times, eventually you will see some false positives). You can test several primes and composites in a Scheme interpreter. Here's the output I got when I tested the first six Carmichael numbers.

> (miller-rabin-test 561)
#f
> (miller-rabin-test 1105)
#f
> (miller-rabin-test 1729)
#f
> (miller-rabin-test 2465)
#f
> (miller-rabin-test 2821)
#f
> (miller-rabin-test 6601)
#f

Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

"Goes to" Considered Harmful

A few months ago the following question was asked (half jokingly) on Stack Overflow,

What is the name of this operator: “-->”?

The question includes the following C++ code:

#include <stdio.h>
int main()
{
    int x = 10;
    while( x --> 0 ) // x goes to 0
    {
        printf("%d ", x);
    }
}

The output of this program is to simply count down from 9 to 0, printing each value. As many others pointed out, "-->" isn't a single operator at all, but two operators smashed together. The condition of the while loop above can be rewritten with proper spacing as:

while( x-- > 0 ) // x-- is greater than  0
{
    ...
}

Joshua Bloch commented on the goes-to operator on Twitter yesterday.



He links to the results of a Google Code search for instances of "-->" in the wild. The results are specific to the C programming language, but "-->" turns up in similar searches for C++ and Java as well (and probably any other language that allows it).

What's the harm?

So what's so harmful about this bit of cleverness? The main problem is that it's too clever. Any time you use standard operators in a non-standard way (even if the language specification doesn't strictly forbid it) you're negatively impacting the readability of your code. In this specific case, the "-->" operator isn't documented anywhere. It's just not reasonable to expect other programmers to know what it means.

On top of that, this situation isn't best handled by a while loop to begin with. If you know ahead of time how many iterations your loop needs to make, use a for loop. Kernighan and Richie covered this a very long time ago in The C Programming Language. It's still true today.

The for is preferable when there is a simple initialization and increment, since it keeps the loop control statements close together and visible at the top of the loop.

This, of course, is also true for a simple decrementing loop. The meaning of the following code should be obvious to any junior programmer.

int i;
for( i = 9; i >= 0; i-- )
{
    printf("%d ", i);
}

Writing code that's "too clever" for other programmers to understand immediately just obfuscates your code needlessly. You'll seem much cleverer to your colleagues if you write your code in a clear and readable style to begin with. Besides that, using the common idioms of your language is probably the quickest way to reduce your WTFs/minute score at your next code review.

SICP Exercise 1.27: Carmichael numbers

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.27 asks us to demonstrate that the first six Carmicheal numbers (561, 1105, 1729, 2465, 2821, and 6601) really do fool the Fermat test. We're to write a procedure that takes an integer n and tests whether aⁿ is congruent to a modulo n for every value less than n, then test the procedure on each of the Carmichael numbers listed.

The fermat-test procedure that we first saw in exercise 1.24 is a good starting point.

(define (fermat-test n)
   (define (try-it a)
     (= (expmod a n n) a))
   (try-it (+ 1 (random-integer (- n 1)))))

Let' start by modifying this procedure so that we pass it the value to test, instead of generating a random value.

(define (fermat-test n a)
   (= (expmod a n n) a))

Now we just need to define a new procedure that will call this one for every value less than the one we want to test. If fermat-test passes for all values tested we should return true, but if it fails for any value we should return false. Here is the complete code listing for the exercise.

(define (square x)
   (* x x))

(define (even? n)
   (= (remainder n 2) 0))

(define (expmod base exp m)
   (cond ((= exp 0) 1)
         ((even? exp)
          (remainder (square (expmod base (/ exp 2) m))
                      m))
        (else
          (remainder (* base (expmod base (- exp 1) m))
                     m))))

(define (fermat-test n a)
   (= (expmod a n n) a))

(define (fermat-full n)
   (define (iter a)
     (cond ((= a 1) #t)
           ((not (fermat-test n a)) #f)
           (else (iter (- a 1)))))
   (iter (- n 1)))

Prime numbers will legitimately pass this test, and Carmichael numbers will fool it. All other composites should fail. You can test out a few primes and composites in a Scheme interpreter to make sure it works as described. Here's the output for the first six Carmichael numbers.

> (fermat-full 561)
#t
> (fermat-full 1105)
#t
> (fermat-full 1729)
#t
> (fermat-full 2465)
#t
> (fermat-full 2821)
#t
> (fermat-full 6601)
#t

Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

SICP Exercise 1.26: Explicit multiplication vs. squaring

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.26 asks us to consider yet another variation on the fast-prime? procedure from exercise 1.24. This time the expmod procedure has been modified to use an explicit multiplication instead of calling square:

(define (expmod base exp m)
   (cond ((= exp 0) 1)
         ((even? exp)
          (remainder (* (expmod base (/ exp 2) m)
                        (expmod base (/ exp 2) m))
                     m))
         (else
          (remainder (* base (expmod base (- exp 1) m))
                     m))))

This change causes fast-prime? to run more slowly than the original prime? test by transforming the O(log n) process into a O(n) process. We're asked to explain this transformation.

A cursory inspection of the code makes it seem like the explicit multiplication will cause twice as many calls to expmod because each input argument to * will be evaluated separately, instead of only once when expmod is passed as the single argument to square. This isn't enough to account for the reported slow down.

Let's take a closer look at the process generated by each version of expmod. (Since the only difference between the two versions of expmod is in the branch where exp is even, I'm using powers of 2 for that argument to better illustrate the difference in the resulting processes. I should point out that this is a worst case scenario.)

Here is what the process generated by the original expmod procedure using square might look like:

(expmod 2 8 2)
   (expmod 2 4 2)
       (expmod 2 2 2)
           (expmod 2 1 2)
               (expmod 2 0 2)
               1
           0
       0
   0
0

This is a pretty straightforward linear recursive process. Now let's look at the process for expmod when explicit multiplication is used. (This is only part of the process diagram, but it's enough to see what's going on.)


Right away, it's easy to see what went wrong here. By replacing the call to square with explicit multiplication, we've transformed expmod from a linear recursive process (like the factorial example) into a tree recursion (like the original Fibonacci example). This causes the number of recursive calls to expmod to grow exponentially instead of simply doubling. What was once a O(log n) process is now O(log 2ⁿ), which simplifies to O(n).


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

SICP Exercise 1.25: A closer look at expmod

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.25 asks us to take a closer look at the expmod procedure that we used in exercise 1.24 to compute the exponential of a number modulo another number:

(define (expmod base exp m)
   (cond ((= exp 0) 1)
         ((even? exp)
          (remainder (square (expmod base (/ exp 2) m))
                     m))
         (else
          (remainder (* base (expmod base (- exp 1) m))
                     m))))

We're asked to consider a more straightforward implementation of the same function. This version makes use of the fast-expt procedure from section 1.2.4, which I've also included here.

(define (expmod base exp m)
   (remainder (fast-expt base exp) m))

(define (fast-expt b n)
   (cond ((= n 0) 1)
         ((even? n) (square (fast-expt b (/ n 2))))
         (else (* b (fast-expt b (- n 1))))))

Will this version perform just as well in the fast prime tester from exercise 1.24? The easiest way to answer that question is to drop it in and test it out. Starting out with relatively small values, you can see that the new procedure doesn't perform very well at all.

> (timed-prime-test 1999)

1999 *** 138.0
> (timed-prime-test 10007)

10007 *** 681.0
> (timed-prime-test 100003)

100003 *** 29059.0

These values are several orders of magnitude smaller than those used in the previous exercise, but the new procedure is taking a much longer time to complete. Why is that?

The answer is buried in a footnote to the expmod code (#46).

The reduction steps in the cases where the exponent e is greater than 1 are based on the fact that, for any integers x, y, and m, we can find the remainder of x times y modulo m by computing separately the remainders of x modulo m and y modulo m, multiplying these, and then taking the remainder of the result modulo m. For instance, in the case where e is even, we compute the remainder of be/2 modulo m, square this, and take the remainder modulo m. This technique is useful because it means we can perform our computation without ever having to deal with numbers much larger than m.

The important point is that the original expmod procedure uses successive squaring to perform its computations without ever having to deal with numbers larger than m. Simple arithmetic with very large numbers is much slower than arithmetic with 32-bit integers. That's because the time complexity for arithmetic operations is based on the number of bits in the operands. The intermediate results during computation in the fast-expt procedure get very big, very quickly (the final result is growing exponentially, after all). Since we're only interested in the remainder, modular arithmetic provides a much sleeker solution, as explained in the footnote.


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

SICP Exercise 1.24: The Fermat Test

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.24 asks us to once again modify the timed-prime-test procedure from exercise 1.22, this time using fast-prime? (which uses the Fermat test), then test the improvement using the timing statistics we gathered before.

The Fermat test is a probabilistic method for testing primality, meaning that if a given number passes the Fermat test, the best we can say is that the number is very probably prime. It's based on Fermat's Little Theorem which states:

If n is a prime number and a is any positive integer less than n, then a raised to the nth power is congruent to a modulo n.

So to test n for primality, we select a number, a, which is less than n, and raise it to the nth power. We then divide aⁿ by n to get the remainder. If n is prime, then the remainder will be equal to a, the base we selected.

Fermat's Little Theorem says that this property will always hold true if n is prime, but it doesn't say anything about when n isn't prime. If aⁿ mod n is not equal to a, then we know for certain that n is not prime. However, there are values of n and a that will pass the Fermat test even when n is not prime. If we test enough values of a for a given n, we can increase our confidence that n is prime. Unfortunately, there are extremely rare values of n called Carmichael numbers that pass the Fermat test for any value of a. There are variations on the Fermat test that cannot be fooled. We'll look at one of those variations, the Miller-Rabin test, in a later exercise.

The complete code listing below shows what modifications were necessary to use the Fermat test in the timed-prime-test procedure. The fast-prime? procedure takes two arguments, n, and times. The first argument is the number to test, the second argument tells the procedure how many random values to test with. The more random values we use, the higher our confidence that n is prime, so we should test a lot of values. I've (rather arbitrarily) chosen to test 100 random values.

One other change that you may notice is the inclusion of a library module in the first line of code. When testing the original code listing from the book, I found that Scheme's primitive random procedure has a limit of 4,294,967,087. This won't do for the values we're testing, so I replaced it with the random-integer procedure, which doesn't have this limitation. I found out about this procedure and library through the Scheme Cookbook.

(require (lib "27.ss" "srfi"))

(define (timed-prime-test n)
   (newline)
   (display n)
   (start-prime-test n (current-inexact-milliseconds)))

(define (start-prime-test n start-time)
   (cond ((fast-prime? n 100)
          (report-prime (- (current-inexact-milliseconds) start-time)))))

(define (report-prime elapsed-time)
   (display " *** ")
   (display elapsed-time))

(define (square x)
   (* x x))

(define (even? n)
   (= (remainder n 2) 0))

(define (expmod base exp m)
   (cond ((= exp 0) 1)
         ((even? exp)
          (remainder (square (expmod base (/ exp 2) m))
                     m))
         (else
          (remainder (* base (expmod base (- exp 1) m))
                     m))))

(define (fermat-test n)
   (define (try-it a)
     (= (expmod a n n) a))
   (try-it (+ 1 (random-integer (- n 1)))))

(define (fast-prime? n times)
   (cond ((= times 0) true)
         ((fermat-test n) (fast-prime? n (- times 1)))
         (else false)))

Test Results

We need to run timed-prime-test again with the new modifications using the same set of inputs that we found in the previous exercises, then compare the results to see how much we improved the run time. The following table compares the original data we gathered in exercise 1.22, the improved algorithm we used in exercise 1.23 (improved time column), and the new results we got using the Fermat test (all three time columns are in milliseconds). The ratio column is the ratio of the improved values from exercise 1.23 to the fermat time column.

prime	original time	improved time	fermat time	ratio
10000000019	141	103	178	0.58
10000000033	172	101	187	0.54
10000000061	154	112	173	0.65
100000000003	516	310	181	1.71
100000000019	493	295	186	1.59
100000000057	527	305	189	1.61
1000000000039	1627	861	189	4.56
1000000000061	1559	884	195	4.53
1000000000063	1549	873	197	4.43
10000000000037	5014	2671	208	12.84
10000000000051	4932	2668	211	12.64
10000000000099	4855	2697	208	12.97
100000000000031	15745	8531	233	36.61
100000000000067	16022	8264	231	35.77
100000000000097	15861	8530	225	37.91
1000000000000037	48950	26346	244	108.0
1000000000000091	48836	26210	255	102.8
1000000000000159	49008	26256	257	102.2

Analysis

The timing numbers from the Fermat test start out looking pretty poor compared to what we've seen previously. This has mostly to do with the completely arbitrary number of random values I chose to test each prime with. As we increase the size of the numbers we're testing, you can see that the time using the Fermat test barely increases at all. We can verify that the time increase is logarithmic by observing that as the numbers under test increase by a factor of 10, the ratio column increases by a factor of roughly three. This logarithmic characteristic is due to the fact that the expmod procedure used in fermat-test has a logarithmic time complexity.


Related:

For links to all of the SICP lecture notes and exercises that I've done so far, see The SICP Challenge.

Solution to A Curiously Clever Chess Puzzle

In yesterday's post I asked two questions about the following chess position.



First, how is this a legal position, and second, how can white guarantee a checkmate in four moves or fewer?

This puzzle is much more commonly known as Lord Dunsany's Chess Problem, although he had several others. I first ran across this particular problem by Lord Dunsany in Martin Gardner's My Best Mathematical and Logic Puzzles.

I'll let Mr. Gardner explain the first part of the puzzle:

The key to Lord Dunsany's chess problem is the fact that the black queen is not on a black square as she must be at the start of a game. This means that the black king and queen have moved, and this could have happened only if some black pawns have moved. Pawns cannot move backward, so we are forced to conclude that the black pawns reached their present positions from the other side of the board! With this in mind, it is easy to discover that the white knight on the right has an easy mate in four moves.

If that seems like a dirty, dirty trick to you, that's because it is. Chess diagrams are traditionally drawn from the perspective of the white player (white pieces at the bottom of the diagram in the starting position), so this puzzle really forces you to check your assumptions and "think outside the box."

As Gardner mentioned, the mate in four is easy once you've mentally turned the board around. First, white moves his knight in front of his king, threatening mate in two more moves. Black can delay white's plan by one move by moving his own knight out into the bishop's file.



When white advances his knight on the bishop's file, black can protect the mating square with his own knight.



Next, white can capture the knight with his queen. After that, black is defenseless to stop the white knight from delivering checkmate on the fourth move. The black king is hemmed in by his own pieces.

All chess diagrams are courtesy of the Online Chess Diagram Editor.