Bill the Lizard

SICP Exercise 1.28: The Miller-Rabin test

2010-03-04T21:18:00.003-05:00

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:

a number not equal to 1 or (n - 1)
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

2010-02-27T17:08:00.012-05:00

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

2010-02-23T19:23:00.005-05:00

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

2010-02-23T15:49:00.007-05:00

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

2010-02-21T12:13:00.003-05:00

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

2010-02-20T16:28:00.010-05:00

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

2010-02-14T22:44:00.015-05:00

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.

A Curiously Clever Chess Puzzle

2010-02-14T21:32:00.007-05:00

It's been some time since I've posted any logic puzzles, and I don't remember ever posting a chess puzzle, so today I'll kill two birds with one stone. In fact, the following position represents two puzzles in one. (You don't have to be a grandmaster at chess to solve it, so even if you only know the basic rules of the game you should give it a try.)

The first part of the puzzle requires more logical thinking than chess skill. Explain how this highly unlikely position could possibly be reached in a game. You don't have to show all the moves, but you should be able to convince yourself that it is a legal position.

Second, it's white's turn to move. Show how white can checkmate black in at most four moves.

I'll give the answer and the original source of this puzzle in an upcoming post. (Update: the full solution is up.)

The diagram image is courtesy of the Online Chess Diagram Editor.

SICP Exercise 1.23: Improved Prime Test

2010-02-08T17:02:00.005-05:00

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.23 asks us to add a modification to the smallest-divisor procedure, then test the improvement using the timing statistics we gathered in exercise 1.22.

We start with the observation that smallest-divisor inefficiently finds the divisor of n by checking every candidate value from 2 to √n. The number of candidates can easily be cut nearly in half by simply not checking even numbers greater than 2.

Our first task is to define a procedure next that returns 3 if its input is equal to 2 and otherwise returns its input plus 2. The code is as straightforward as the description:

(define (next x)
  (if (= x 2) 3 (+ x 2)))

Our next task is to modify the find-divisor procedure to use (next test-divisor) instead of (+ test-divisor 1). This is a straight substitution, and no other changes to the code from exercise 1.22 are necessary.

(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
        ((divides? test-divisor n) test-divisor)
        (else (find-divisor n (next test-divisor)))))

Finally, we need to run timed-prime-test with these modifications using the same set of inputs that we found in the previous exercise, then compare the results to see if we really did cut the run time in half. The following table compares the original data to new timing data gathered using the improved algorithm. (Values in the new time column are averaged from three runs of the procedure.)

prime	old time (ms)	new time (ms)	ratio
10000000019	141	103	1.37
10000000033	172	101	1.70
10000000061	154	112	1.38
100000000003	516	310	1.67
100000000019	493	295	1.67
100000000057	527	305	1.73
1000000000039	1627	861	1.89
1000000000061	1559	884	1.76
1000000000063	1549	873	1.77
10000000000037	5014	2671	1.88
10000000000051	4932	2668	1.85
10000000000099	4855	2697	1.80
100000000000031	15745	8531	1.85
100000000000067	16022	8264	1.94
100000000000097	15861	8530	1.86
1000000000000037	48950	26346	1.86
1000000000000091	48836	26210	1.86
1000000000000159	49008	26256	1.87

Analysis

We're seeing a clear improvement in the new procedure, but it's not quite as fast as we expected. The first thing that needs to be explained in this data is the fact that the first three values shows very little performance gain, the next three a little more, then fairly consistent results for the remaining data. I think this can be explained by other processes running on the computer. Measuring shorter runs of the procedure (those in the 100-500 millisecond range) is going to be much more sensitive to measurement error due to being interrupted by background processes. These errors will be a less significant proportion of the total run time for longer runs of the procedure.

We're also seeing that the procedure is only running approximately 1.85 times as fast, instead of the expected factor of 2. This may be explained by the fact that we replaced a primitive operation, (+ test-divisor 1), by a user-defined operation, (next test-divisor). Each time that user-defined operation is called, an extra if must be evaluated (to check if the input is 2). Other than this small discrepancy, I think the improvement is quite good for such a small change to the code.

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.22: Timed Prime Test

2010-02-05T17:36:00.018-05:00

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.22 introduces a new primitive called runtime...

Most Lisp implementations include a primitive called runtime that returns an integer that specifies the amount of time the system has been running (measured, for example, in microseconds).

Unfortunately, this turns out not to be the case for Scheme. After a bit of research, I found what I think to be a suitable replacement: current-inexact-milliseconds. It's not as precise as I'd like, but there are a couple of workarounds to this problem that I'll explain shortly.

The exercise goes on to to introduce the timed-prime-test procedure, which prints its input n and tests to see if n is prime. If n is prime, the procedure prints three asterisks followed by the amount of time used in performing the test. Here is the modified version that I'll be using, including all of the required support procedures from the book:

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

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

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

(define (smallest-divisor n)
   (find-divisor n 2))

(define (find-divisor n test-divisor)
   (cond ((> (square test-divisor) n) n)
         ((divides? test-divisor n) test-divisor)
         (else (find-divisor n (+ test-divisor 1)))))

(define (divides? a b)
   (= (remainder b a) 0))

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

(define (prime? n)
   (= n (smallest-divisor n)))

Our task (finally) is to use timed-prime-test to write a procedure named search-for-primes that checks the primality of consecutive odd integers in a specified range. We're instructed to use this procedure to find the three smallest primes larger than 1000, 10000, 100000, and 1000000. Since the prime? procedure has an order of growth of O(√n), we should expect a run time increase of a factor of √10 at each step. We'll use the output from timed-prime-test to check this assumption.

The following procedure checks the primality of consecutive odd integers in the range from start to end.

(define (search-for-primes start end)
   (if (even? start)
       (search-for-primes (+ start 1) end)
       (cond ((< start end) (timed-prime-test start)
                            (search-for-primes (+ start 2) end)))))

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

The procedure initially checks to see if the start input argument is even. If it is, the procedure simply starts over with the next (odd) integer. Next the procedure checks to see if start is less than end. If so, it performs a timed-prime-test on start, then recursively calls itself with the next odd integer as a starting value. This process will continue until start exceeds end.

With a little bit of trial and error we can find start and end values that surround the first three primes larger than our target values. Here's a sample run starting at 1000.

> (search-for-primes 1000 1020)

1001
1003
1005
1007
1009 *** 1.0
1011
1013 *** 1.0
1015
1017
1019 *** 0.0

We can tell from this output that our timer's resolution isn't fine enough to give us meaningful results on modern desktop hardware (even if it is cheap). I see two things we can do to work around this problem. First, we could ignore the suggested inputs from the exercise and just keep increasing the input values until we get meaningful measurements. Second, we could run prime? in a loop and measure how long it takes to test a value 1000 (or maybe more) times. The code is ready as-is for the first option, so that's what I'm going to try. (If anyone implements the second option, please leave a comment so we can compare results.)

The first consistent results I got were by using a start value of 10¹⁰. It took an average of 154 milliseconds to test each of the first three primes after that value. The following table shows details of the time required to test the primality of numbers spread over several orders of magnitude.

start	primes	time(ms)	avg time(ms)	ratio
10¹⁰	10000000019	141	155.67	-
	10000000033	172
	10000000061	154
10¹¹	100000000003	516	512	3.289
	100000000019	493
	100000000057	527
10¹²	1000000000039	1627	1578.33	3.082
	1000000000061	1559
	1000000000063	1549
10¹³	10000000000037	5014	4933.67	3.126
	10000000000051	4932
	10000000000099	4855
10¹⁴	100000000000031	15745	15876	3.218
	100000000000067	16022
	100000000000097	15861
10¹⁵	1000000000000037	48950	48931.33	3.082
	1000000000000091	48836
	1000000000000159	49008

The last column shows the ratio of the average time for each row to the average time for the preceding row. These ratios stay very close to √10, which is approximately 3.162. These results do seem to verify that programs on my machine run in time proportional to the number of steps required for the computation.

Related:

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

Unsolved: Perfect Numbers and Mersenne Primes

2010-02-01T17:11:00.001-05:00

In an earlier post I talked about perfect numbers, and how it was still unknown whether there are any odd perfect numbers. There's another unsolved problem surrounding the perfect numbers, and this one intersects with a special set of numbers that you may have heard of, the Mersenne primes.

Mersenne numbers and Mersenne primes

Mersenne numbers take their name from the French monk and mathematician Father Marin Mersenne, even though they were studied as long ago as the times of Euclid. You construct Mersenne numbers by subtracting 1 from a power of 2.

M_n = 2ⁿ - 1

The following table details the construction of the first sixteen Mersenne numbers.

n	2ⁿ	2ⁿ - 1
1	2	1
2	4	3
3	8	7
4	16	15
5	32	31
6	64	63
7	128	127
8	256	255
9	512	511
10	1024	1023
11	2048	2047
12	4096	4095
13	8192	8191
14	16384	16383
15	32768	32767
16	65536	65535

If you look closely at the highlighted fields you may notice what mathematicians have known for centuries. A Mersenne number (2ⁿ - 1) cannot be prime unless the exponent n is prime. This is often expressed with a slightly different notation.

M_p = 2^p - 1

What wasn't known in the time of the ancient Greeks was whether or not the exponent being prime meant that M_p would always be prime. It wasn't until the 16th century, when it was discovered that 2¹¹ - 1 = 2047 is not prime, that this question was finally settled (2047 = 23 * 89). In order for M_p to be prime, it is necessary but not sufficient for the exponent p to also be prime.

Constructing Perfect Numbers

Perfect numbers have also been studied since the times of the early Greek mathematicians. It was Euclid who first discovered that the first four perfect numbers are given by the formula

P = 2^p-1 * (2^p - 1)

6 = 2¹ * (2² - 1)
28 = 2² * (2³ - 1)
496 = 2⁴ * (2⁵ - 1)
8128 = 2⁶ * (2⁷ - 1)

It wasn't until 1849 that it was shown by Euler (in a paper published after his death) that all even perfect numbers are of the form 2^p-1 * (2^p - 1), where (2^p - 1) is a Mersenne prime.

Since every even perfect number corresponds to a Mersenne prime (and vice versa), searching for perfect numbers is the same as searching for Mersenne primes. This leads us to two unsolved problems.

It isn't known if there are infinitely many even perfect numbers.

It isn't known if there are infinitely many Mersenne primes.

(In fact, given the 1-to-1 correspondence between even perfect numbers and Mersenne primes, it really wouldn't be incorrect to say that either one of these problems is just a restatement of the other.)

The largest prime number yet discovered as of this writing is the Mersenne prime 2^43,112,609 - 1, which has a starggering 12,978,189 digits. The corresponding perfect number, 2^43,112,608 * (2^43,112,609 - 1), has 25,956,377 digits.

Additional reading

In this article I've barely scratched the surface of prime number research. There is a very large network of people and computers dedicated to the search for Mersenne primes. For more information (and to possibly join in the search) have a look at the Great Internet Mersenne Prime Search (GIMPS).

Interesting Miscellany

2010-01-30T12:34:00.004-05:00

Here are some of the links I've found interesting enough to tweet or retweet in that past several weeks.

Programming
A beautiful obfuscated Pi calculator
Programming language trends of 2009
Tweetable game of life in MATLAB
Factorization of RSA768
Calculating Sines

Science
Leonardo da Vinci's resume
Richard Dawkins: Growing up in the universe
Richard Hamming: You and your research

Math
A quick little puzzle (in probability)
Euler Book Prize (A list to watch in the future. Since the site isn't updated for 2010 yet, the winner was Euler's Gem by David S. Richeson.)
A cryptic math advertisement

Miscellaneous
Succed Blog
RockYou.com's 20 most common passwords
50 Books Every Geek Should Read

SICP Exercise 1.21: Smallest Divisor

2010-01-30T09:44:00.004-05:00

From SICP section 1.2.6 Example: Testing for Primality

Exercise 1.21 asks us to use the smallest-divisor procedure to find the smallest divisor of each of the following numbers: 199, 1999, 19999.

The smallest-divisor procedure was given earlier in the section.

(define (smallest-divisor n)
   (find-divisor n 2))

(define (find-divisor n test-divisor)
   (cond ((> (square test-divisor) n) n)
         ((divides? test-divisor n) test-divisor)
         (else (find-divisor n (+ test-divisor 1)))))

(define (divides? a b)
   (= (remainder b a) 0))

This procedure makes use of the square procedure that was defined earlier. You'll need to include it if you want to run the code in an interpreter.

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

The smallest-divisor procedure finds the smallest integral divisor (greater than 1) of a given number n in a straightforward way, by testing to see if n is divisible by each integer from 2 to √n. Since the number of steps is between 1 and √n, the order of growth of the algorithm is O(√n).

We can find the solutions to the exercise by simply running the code in an interpreter.

> (smallest-divisor 199)
199
> (smallest-divisor 1999)
1999
> (smallest-divisor 19999)
7

As you can see, the first two numbers are prime because they are their own smallest divisor, but the third is divisible by 7.

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.20: GCD

2010-01-26T13:06:00.006-05:00

From SICP section 1.2.5 Greatest Common Divisors

Exercise 1.20 asks us to consider the following iterative procedure that uses Euclid's Algorithm to comput the GCD:

(define (gcd a b)
   (if (= b 0)
       a
       (gcd b (remainder a b))))

We're asked to illustrate the process generated by this procedure when using normal-order and applicative-order evaluation rules to evaluate (gcd 206 40). How many remainder operations are actually performed using each evaluation model?

We learned the rules for both applicative-order and normal-order evaluation back in SICP section 1.1.5 The Substitution Model for Procedure Application. We learned that under normal-order evaluation rules, the interpreter fully expands a procedure before reducing it. Operand expressions are substituted for parameters until an expression involving only primitive operators is reached. Under applicative-order rules, arguments are evaluated before operators are applied. This can often avoid multiple unnecessary evaluations of the same expression, which is one of the reasons why Lisp uses applicative-order evaluation.

In section 1.1.6 Conditional Expressions and Predicates (see exercise 1.5) we're told to assume that the evaluation rule for the special form if is the same whether we use normal or applicative order evaluation. "The predicate expression is evaluated first, and the result determines whether to evaluate the consequent or the alternative expression." We'll continue using this assumption until we're told otherwise.

Now let's put all of this together to illustrate both processes for (gcd 206 40).

Normal-order evaluation - fully expand and then reduce.

(gcd 206 40)

(if (= 40 0)
    206
    (gcd 40 (remainder 206 40)))

(gcd 40 (remainder 206 40))

(if (= (remainder 206 40) 0)
    40
    (gcd (remainder 206 40) (remainder 40 (remainder 206 40))))

(if (= 6 0) ;; 1 remainder evaluated
    40
    (gcd (remainder 206 40) (remainder 40 (remainder 206 40))))

(gcd (remainder 206 40) (remainder 40 (remainder 206 40)))

(if (= (remainder 40 (remainder 206 40)) 0)
    (remainder 206 40)
    (gcd (remainder 40 (remainder 206 40))
         (remainder (remainder 206 40)
                    (remainder 40 (remainder 206 40)))))

At this point both of the nested remainder operations that were substituted in for b in the (= b 0) procedure need to be evaluated. I'll combine these and show them as one step.

(if (= 4 0) ;; 2 remainders evaluated
    (remainder 206 40)
    (gcd (remainder 40 (remainder 206 40))
         (remainder (remainder 206 40)
                    (remainder 40 (remainder 206 40)))))

(gcd (remainder 40 (remainder 206 40))
     (remainder (remainder 206 40)
                (remainder 40 (remainder 206 40))))

(if (= (remainder (remainder 206 40)
                  (remainder 40 (remainder 206 40)))
       0)
       (remainder 40 (remainder 206 40))
       (gcd (remainder (remainder 206 40)
                       (remainder 40 (remainder 206 40)))
            (remainder (remainder 40 (remainder 206 40))
                       (remainder (remainder 206 40)
                                  (remainder 40 (remainder 206 40))))))

Here again, multiple nested remainder operations that were substituted in for b in the (= b 0) procedure need to be evaluated. I'll combine all four and show them as one step.

(if (= 2 0) ;; 4 remainders evaluated
    (remainder 40 (remainder 206 40))
    (gcd (remainder (remainder 206 40)
                    (remainder 40 (remainder 206 40)))
         (remainder (remainder 40 (remainder 206 40))
                    (remainder (remainder 206 40)
                               (remainder 40 (remainder 206 40))))))

(gcd (remainder (remainder 206 40)
                (remainder 40 (remainder 206 40)))
     (remainder (remainder 40 (remainder 206 40))
                (remainder (remainder 206 40)
                           (remainder 40 (remainder 206 40)))))

(if (= (remainder (remainder 40 (remainder 206 40))
                  (remainder (remainder 206 40)
                             (remainder 40 (remainder 206 40)))) 0)
    (remainder (remainder 206 40)
               (remainder 40 (remainder 206 40)))
    (gcd (remainder (remainder 40 (remainder 206 40))
                    (remainder (remainder 206 40)
                               (remainder 40 (remainder 206 40))))
         (remainder (remainder (remainder 206 40)
                               (remainder 40 (remainder 206 40)))
                    (remainder (remainder 40 (remainder 206 40))
                               (remainder (remainder 206 40)
                                          (remainder 40
                                                     (remainder 206 40)))))))

There are seven nested remainder operations to evaluate at this point.

(if (= 0 0) ;; 7 remainders evaluated
    (remainder (remainder 206 40)
               (remainder 40 (remainder 206 40)))
    (gcd (remainder (remainder 40 (remainder 206 40))
                    (remainder (remainder 206 40)
                               (remainder 40 (remainder 206 40))))
         (remainder (remainder (remainder 206 40)
                               (remainder 40 (remainder 206 40)))
                    (remainder (remainder 40 (remainder 206 40))
                               (remainder (remainder 206 40)
                                          (remainder 40
                                                     (remainder 206
                                                                40)))))))

(remainder (remainder 206 40)
           (remainder 40 (remainder 206 40)))

Finally, we can evaluate these four remainder operations to come up with the answer (GCD) and the total.

2 ;; 18 total remainders evaluated

Applicative-order evaluation - evaluate arguments and then apply operands. This one is a little bit easier to follow. With each substitution, we just evaluate the operands that we need before applying the operators at each step.

(gcd 206 40)

(if (= 40 0)
    206
    (gcd 40 (remainder 206 40))))

(gcd 40 (remainder 206 40))

(gcd 40 6) ;; 1st remainder evaluated

(if (= 6 0)
    40
    (gcd 6 (remainder 40 6)))

(gcd 6 (remainder 40 6))

(gcd 6 4) ;; 2nd remainder evaluated

(if (= 4 0)
    6
    (gcd 4 (remainder 6 4)))

(gcd 4 (remainder 6 4))

(gcd 4 2) ;; 3rd remainder evaluated

(if (= 2 0)
    4
    (gcd 2 (remainder 4 2)))

(gcd 2 (remainder 4 2))

(gcd 2 0) ;; 4th remainder evaluated

(if (= 0 0)
    2
    (gcd 0 (remainder 2 0)))

2

As you can see, using applicative-order evaluation, remainder was only evaluated four times, compared to 18 times using normal-order evaluation.

Related:

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

Unsolved: Odd Perfect Numbers

2010-01-23T11:10:00.005-05:00

A perfect number is a number that is the sum of its proper divisors (i.e., including 1 but excluding itself).

For example, the smallest perfect number is 6.

1 + 2 + 3 = 6

The next two perfect numbers are

1 + 2 + 4 + 7 + 14 = 28
1 + 2 + 4 + 8 + 16 + 31 + 62 + 124 + 248 = 496

See sequence A000396 in OEIS for a list of the first ten perfect numbers.

On a related note, a number that is smaller than the sum of its proper divisors is an abundant number, while a number greater than the sum of its proper divisors is called deficient. The smallest abundant number is 12, whose proper divisors sum to 16.

All of the perfect numbers that have been discovered so far are even.

The problem is to find an odd perfect number, or prove that no odd perfect numbers exist.

If there is an odd perfect number, we already know a lot about it. It must be at least 300 digits long, with at least 75 prime factors, at least 9 of them distinct. Its largest prime factor must be greater than 100,000,000. Its second largest prime factor is greater than 10,000, and its third largest is greater than 100.

SICP Exercise 1.19: Computing Fibonacci numbers

2010-01-22T15:43:00.004-05:00

From SICP section 1.2.4 Exponentiation

Exercise 1.19 asks us to complete a procedure for computing Fibonacci numbers in a logarithmic number of steps. The following code is given:

(define (fib n)
   (fib-iter 1 0 0 1 n))

(define (fib-iter a b p q count)
   (cond ((= count 0) b)
         ((even? count)
          (fib-iter a
                    b
                    <??>      ; compute p'
                    <??>      ; compute q'
                    (/ count 2)))
         (else (fib-iter (+ (* b q) (* a q) (* a p))
                         (+ (* b p) (* a q))
                         p
                         q
                         (- count 1)))))

We're reminded of the transformation of the state variables a and b in the original fib-iter procedure from section 1.2.2: a ← a + b and b ← a. If these state changes are labeled transformation T, then applying T repeatedly for n iterations starting with a = 1 and b = 0 produces the pair a = Fib(n + 1) and b = Fib(n). So the Fibonacci numbers are produced by the nth power of the transformation T, or Tⁿ, starting with the pair (1, 0).

We are then asked to consider the family of transformations T_pq which transforms the pair (a, b) according to the following rules:

a ← bq + aq + ap
b ← bp + aq

We can verify by quick substitution that the original transformation T is just a special case of T_pq, where p = 0 and q = 1.

a ← b(1) + a(1) + a(0)
a ← b + a

b ← b(0) + a(1)
b ← a

We are asked to show that if we apply T_pq twice, the effect is the same as using a single transformation T_p'q' of the same form, and compute p' and q' in terms of p and q. This will give us an explicit way to square these transformations, which we can use to compute Tⁿ using successive squaring, just like the fast-expt procedure from exercise 1.16.

We can apply T_pq twice by defining new variables and using substitution. Let's define a₁ and b₁ as the results of applying transformation T_pq once

a₁ = bq + aq + ap
b₁ = bp + aq

The next step is to define a₂ and b₂ and apply the tranformation a second time, this time using a₁ and b₁ in place of a and b.

a₂ = b₁q + a₁q + a₁p
b₂ = b₁p + a₁q

Now that we have a system of equations defined, we can use substitution on our way to simplifying.

a₂ = (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
b₂ = (bp + aq)p + (bq + aq + ap)q

The second equation is shorter, so it should be easier to manipulate. Remember, we're trying to find p' and q', so we need to rewrite the equation to fit the form

b₂ = bp' + aq'

where p' and q' can be computed in terms of q and p.

b₂ = (bp + aq)p + (bq + aq + ap)q
= (bpp + apq) + (bqq + aqq + apq)
= bpp + apq + bqq + aqq + apq
= (bpp + bqq) + (2apq + aqq)
= b(pp + qq) + a(2qp + qq)

From this we can see that p' and q' can be computed using the following equations:

p' = p² + q²
q' = 2pq + q²

Manipulating the equation for a₂ in the same way, we can verify those results. This time we're trying to fit the form

a₂ = bq' + aq' + ap'

The groupings required for this manipulation are made even easier by the fact that we now already know the formulas for p' and q'.

a₂ = (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
= (bpq + aqq) + (bqq + aqq + apq) + (bpq + apq + app)
= bpq + aqq + bqq + aqq + apq + bpq + apq + app
= (bpq + bpq + bqq) + (apq + apq + aqq) + (app + aqq)
= b(pq + pq + qq) + a(pq + pq + qq) + a(pp + qq)
= b(2pq + qq) + a(2pq + qq) + a(pp + qq)

Now that we've verified the formulas for p' and q' in terms of p and q, we can use them to complete the procedure we were given.

(define (fib n)
   (fib-iter 1 0 0 1 n))

(define (fib-iter a b p q count)
   (cond ((= count 0) b)
         ((even? count)
          (fib-iter a
                    b
                    (+ (* p p) (* q q))     ; compute p'
                    (+ (* 2 p q) (* q q))   ; compute q'
                    (/ count 2)))
         (else (fib-iter (+ (* b q) (* a q) (* a p))
                         (+ (* b p) (* a q))
                         p
                         q
                         (- count 1)))))

You can test several known values in an interpreter to verify that this actually works.

> (fib 0)
0
> (fib 1)
1
> (fib 2)
1
> (fib 3)
2
> (fib 5)
5
> (fib 10)
55
> (fib 20)
6765

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.18: Iterative multiplication

2010-01-14T14:26:00.003-05:00

From SICP section 1.2.4 Exponentiation

Exercise 1.18 asks us to use the results from the two previous exercises (1.16 and 1.17) to devise a procedure that generates an iterative process for multiplying two integers in terms of adding, doubling and halving, and which uses a logarithmic number of steps.

The key idea from exercise 1.16 was to keep an additional state variable a, and manipulate the state transition so that the product (a * bⁿ) remained unchanged. We can modify our solution to exercise 1.17 to use the same idea. Since we're now dealing with multiplication instead of exponentiation, the invariant quantity will be (a * b + p) where a and b are the two operands, and p will be our state variable. Initially p will be 0, but by the end of the iterative process it will hold the product (a * b).

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

(define (double x)
   (+ x x))

(define (halve x)
   (/ x 2))

(define (mult-iter a b p)
   (cond ((= 0 b) p)
         ((even? b) (mult-iter (double a) (halve b) p))
         (else (mult-iter a (- b 1) (+ a p)))))

(define (mult a b)
   (mult-iter a b 0))

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.17: Multiplication by repeated addition

2010-01-13T14:36:00.005-05:00

From SICP section 1.2.4 Exponentiation

Exercise 1.17 asks us to design a multiplication procedure analogous to fast-expt that uses a logarithmic number of steps. We're told that in addition to adding, we can use procedures for doubling or halving an integer argument.

(define (double x)
   (+ x x))

(define (halve x)
   (/ x 2))

We'll also use the same even? procedure we've seen before:

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

The following multiplication procedure that takes a linear number of steps is given:

(define (* a b)
   (if (= b 0)
       0
       (+ a (* a (- b 1)))))

This procedure defines multiplication in terms of the algorithm we all learned in grade school: add a to itself b times.

The exercise also makes reference to the following fast-expt procedure, which we can use as a guide:

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

This procedure makes use of the following observations:

bⁿ = (b^n/2)² if n is even
bⁿ = b * b^n-1 if n is odd

We can define similar rules for integer multiplication.

a * b = 2 * (a * b/2) if be is even
a * b = a + a * (b - 1) if b is odd

The first rule tells us when to use our procedures double and halve. The second rule just lets us take an odd b and make it even in preparation for the next iteration. Translating these rule into a Scheme procedure give us:

(define (fast-mult a b)
   (cond ((= b 0) 0)
         ((= b 1) a)
         ((even? b) (double (fast-mult a (halve b))))
         (else (+ a (fast-mult a (- b 1))))))

The procedure takes a logarithmic number of steps. If you double the size of the input parameter b, the procedure takes only one more step to compute the product.

Related:

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

Guest article on Programming Praxis

2010-01-12T08:26:00.003-05:00

I was excited to be asked to be the guest author for an article on one of my favorite programming blogs, Programming Praxis. If you're not familiar with the site, Programming Praxis provides weekly practice exercieses to "sharpen your saw" on. It was inspired in part by Project Euler, but the problems on Programming Praxis focus more on developing programming skills and usually require less advanced math than those on Project Euler.

To follow up earlier exercises on Calculating Pi and Logarithms, my article is on Calculating Sines using two different methods, the Taylor series and the Triple-angle formula. The sample solutions are written in Scheme, but readers are encouraged to submit solutions in any language they choose (see the HOWTO on posting source code to post your own solution). Since I'm still a relative newbie to Scheme (I'm currently still in chapter one of SICP), I'll be interested to see what alternate solutions you can come up with.

The article is also currently featured on Planet Scheme.

SICP Exercise 1.16: Fast Exponentiation

2010-01-11T16:07:00.007-05:00

From SICP section 1.2.4 Exponentiation

Exercise 1.16 asks us to design an exponentiation procedure that evolves an iterative process using successive squaring and uses a logarithmic number of steps. A hint tells us to use the observation that (b^n/2)² = (b²)^n/2 and to keep, along with the exponent n and base b, an additional state variable a, and to define the state transformation in such a way that the product a * bⁿ is unchanged from state to state. The value of a should start at 1, and should contain the final answer by the end of the process.

Earlier in section 1.2.4, we were given the procedure:

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

which tests whether an integer is even. This will come in handy again for this exercise, since we'll need to adjust a by a different amount when n is even than when n is odd.

We'll also use the following simple procedure for squaring a number:

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

With those building blocks in place, we can define our iterative procedure.

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

We call the procedure expt-iter with the arguments base b, exponent n, and the state variable a. The first thing we do is check to see if the exponent is 0. If so, we've reached the end of the procedure and just return the value of a. Next we check to see if n is even. If it is, we use the observation that (b^n/2)² = (b²)^n/2 to transform state information by squaring b and dividing n by 2. Finally, if n is odd we transform state information by reducing n by 1 and multiplying a by the base b.

The only thing that's left to do is define a procedure that calls expt-iter with the correct initial values.

(define (fast-expt b n)
   (expt-iter b n 1))

Here are a few sample runs of the procedure:

> (fast-expt 2 2)
4
> (fast-expt 2 3)
8
> (fast-expt 2 10)
1024
> (fast-expt 3 3)
27
> (fast-expt 5 5)
3125

You can try running the procedure with large values of n to verify that it really is fast.

Related:

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

Math visualization: Multiplying Two Binomials

2010-01-10T08:39:00.006-05:00

In my earlier post Math visualization: (x + 1)², I showed how to draw a picture that proved the equality

(x + 1)² = x² + 2x + 1

An alert reader commented that (x + 1)² is a specific instance of the binomial squared

(x + a)² = x² + 2ax + a²

Here's a visualization of that equality:

It turns out that the illustration can be generalized a little bit further to cover multiplying two binomials that form a rectangle (not just a perfect square).

(x + a)(x + b) = x² + ax + bx + ab

If you want to see how this idea can be extended into the third dimension, check out Montessori World's page on the Binomial Cube.

Unsolved: Twin Primes Conjecture

2010-01-04T21:54:00.011-05:00

Prime numbers, as most school children can tell you, are those numbers that are evenly divisible by only themselves and 1. The first few primes are

2, 3, 5, 7, 11, 13, 17...

We've known that there are infinitely many primes since c. 300 BC when it was proven by Euclid of Alexandria. Euclid used a very simple and elegant proof by contradiction.

Euclid's proof of the infinitude of primes

First, assume that there are a finite number of primes, p1, p2, p3, ..., pn.

Now let

Q = (p1 * p2 * ... * pn) + 1

That is, Q is equal to all of the primes multiplied together plus one.

By the Fundamental Theorem of Arithmetic, Q is either prime or it can be written as the product of two or more primes. However, none of the primes in our list evenly divides Q. If any prime in our list did evenly divide Q, then that same prime would also evenly divide 1, since

Q - (p1 * p2 * ... * pn) = 1

This contradicts the assumption that we had listed all the primes. So no matter how many primes we start with in our list, there must always be more primes. (Note that this proof does not claim that Q itself is prime, just that there must be some prime not in the initial list.)

Twin primes

Twin primes are those pairs of numers that have a difference of two, and that are both prime.

3, 5
5, 7
11, 13
17, 19
...

The Twin prime conjecture, which dates back to the 18th century, simply states that there are infinitely many twin primes.

Given the simplicity of Euclid's proof of the infinitude of primes, it's tempting to hope for an equally simple proof to the Twin prime conjecture. Needless to say, such a proof has not been found.

Other facts about twin primes:

Other than (3, 5), all twin primes have the form (6n - 1, 6n + 1).

In 1919 Viggo Brun showed that the sum of the reciprocals of the twin primes converges to a definite number, now known as Brun's constant (approximately 1.902160578).

In 1994, while in the process of estimating Brun's constant by calculating the twin primes up to 10¹⁴, Thomas Nicely discovered the infamous Pentium bug.

The largest known twin primes (as of January 2010) are a pair of 100,355 digit primes with the values 65516468355 * 2³³³³³³ ± 1.

Math visualization: (x + 1)2

2009-12-26T23:27:00.007-05:00

You may remember this from high-school algebra (or perhaps earlier for some). Expand (x + 1)² using the FOIL method.

(x + 1)² = (x + 1)(x + 1)
= x² + x + x + 1
= x² + 2x + 1

A neat way to visualize this equality, and hopefully help remember the factorization of the resulting polynomial, is to look at how the pieces fit together.

When they're arranged in this way, it's easy to see that the pieces squeeze together to form a larger square that has a length of (x + 1) on a side, proving the equality.

Related posts:

Six Visual Proofs
Multiplying Two Binomials

SICP Exercise 1.15: Calculating sines

2009-12-24T23:00:00.007-05:00

From SICP section 1.2.3 Orders of Growth

Exercise 1.15 shows an interesting procedure for computing the sine of an angle (in radians) by combined use of the approximation sin x ~ x (if x is sufficiently small) and the trigonometric identity

sin r = 3 * sin(r/3) - 4 * sin³(r/3)

which is used to reduce the argument of sin. The code provided implements the computation:

(define (cube x) (* x x x))

(define (p x) (- (* 3 x) (* 4 (cube x))))

(define (sine angle)
   (if (not (> (abs angle) 0.1))
       angle
       (p (sine (/ angle 3.0)))))

The exercise goes on to ask the following two questions:

a. How many times is the procedure p applied when (sine 12.15) is evaluated?

b. What is the order of growth in space and number of steps (as a function of a) used by the process generated by the sine procedure when (sine a) is evaluated?

Looking at how the process executes will help us answer both questions.

(sine 12.15)
(p (sine 4.05))
(p (p (sine 1.35)))
(p (p (p (sine 0.45))))
(p (p (p (p (sine 0.15)))))
(p (p (p (p (p (sine 0.05))))))
(p (p (p (p (p 0.05)))))
(p (p (p (p 0.1495))))
(p (p (p 0.4351345505)))
(p (p 0.9758465331678772))
(p -0.7895631144708228)
-0.39980345741334

As long as the angle stays greater than 0.1, the process keeps recursing. There are five calls to the procedure p before that happens.

Complexity

The sine procedure is recursive, but it makes only a single call to itself, so it isn't tree recursive and its complexity isn't nearly as bad as some of the procedures we've seen in previous exercises. The complexity for this procedure in both space and number of steps is actually better than linear.

Space - The amount of space required by this procedure would increase linearly as the input size is tripled. Only the calls to procedure p are deferred, and we can triple the size of the input before an additional call would be made.

Number of steps - It's also easy to see that we can triple the value of the starting input angle and only one more step would be required by this procedure.

The complexity of both the space and the number of steps required by this procedure can be described as O(log 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.14: Counting change

2009-12-22T23:50:00.016-05:00

From SICP section 1.2.3 Orders of Growth

Exercise 1.14 asks us to draw the tree illustrating the process generated by the count-change procedure presented in section 1.2.2 in making change for 11 cents.

(define (count-change amount)
   (cc amount 5))

(define (cc amount kinds-of-coins)
   (cond ((= amount 0) 1)
         ((or (< amount 0) (= kinds-of-coins 0)) 0)
         (else (+ (cc amount
                      (- kinds-of-coins 1))
                  (cc (- amount
                         (first-denomination kinds-of-coins))
                      kinds-of-coins)))))

(define (first-denomination kinds-of-coins)
   (cond ((= kinds-of-coins 1) 1)
         ((= kinds-of-coins 2) 5)
         ((= kinds-of-coins 3) 10)
         ((= kinds-of-coins 4) 25)
         ((= kinds-of-coins 5) 50)))

Running the code with the specified input shows us how many ways there are to make change for 11 cents with five coins.

> (count-change 11)
4
>

A quick scan of the code shows that the cc procedure contains two recursive calls to itself, resulting in the tree-shaped process structure mentioned in the question. Leaf nodes of the tree are reached when the base cases are met, that is when the amount is less than or equal to zero, or when the kinds of coins remaining reaches zero.

When a node splits, its left-hand child is passed the same amount and one fewer kinds of coins. The right-hand child is passed an amount that is reduced by the highest-denomination coin available to its parent node, and the same number of kinds of coins. (Click the images enlarge.)

I've color-coded the nodes of the tree. White nodes are those leaf nodes that evaluate to 0, blue nodes are those leaf nodes that evaluate to 1 and contribute to the final answer. The yellow nodes are those that branch recursively.

Orders of growth

The exercise also asks us to find the orders of growth of the space and number of steps used by the count-change process as the amount to be changed increases. (Note that we are not really concerned with the orders of growth as the number of kinds of coins grows, only the amount to be changed.)

Space - As we saw in the Fibonacci procedure in section 1.2.2, the space required by the count-change procedure grows linearly with the input. The space required is proportional to the maximum depth of the tree. At any point in the computation we only need to keep track of the nodes above the current node.

Since the height of the tree is proportional to the amount to be changed, the order of growth of the space required by the count-change procedure is O(n).

Number of steps - The diagram illustrates that this procedure, much like the recursive procedure for computing Fibonacci numbers that we saw before, in not very efficient. Much of the computation being done is repeated.

If we look closely at a call to cc where kinds-of-coins is equal to one, we can see each call generates two more calls until we reach a terminal node.

If we define the function T(a, k) to be the number of calls to cc generated by calling (cc n k), where n is the amount and k is the number of kinds of coins, then

T(n, 1) = 2n + 1

or in Big-O notation,

T(n, 1) = O(n)

Now let's take a look at a partial graph of a call to (cc 50 2).

There are two interesting things to point out here. First, there are n/5 calls to cc made where k = 2 kinds of coins (look down the right hand side of the diagram above). This is because each call is made with 5 cents less in the amount argument until a terminal node is reached. The second interesting thing is that each of these n/5 calls to (cc n 2) generates an entire (cc n 1) subtree.

That means that

T(n, 2) = (n/5) * 2n + 1

and because of the multiplication of terms

T(n, 2) = O(n²)

Similarly, if we look at the graph of (cc 50 3) we can see n/10 calls to cc where k = 3, each of which generates its own (cc n 2) subtree.

From this we find that

T(n, 3) = O(n³)

Finally, it's easy enough to show that we'll get the same pattern when k = 4 and k = 5. The n/50 nodes generated when k = 5 each generate n/25 nodes with k = 4, each of which generates a node where k = 3.

And so the final answer we arrive at for the time complexity of the count-change procedure with five kinds of coins is

T(n, 5) = O(n⁵)

Related:

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