## Monday, February 8, 2010

### SICP Exercise 1.23: Improved Prime Test

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
100000000191411031.37
100000000331721011.70
100000000611541121.38
1000000000035163101.67
1000000000194932951.67
1000000000575273051.73
100000000003916278611.89
100000000006115598841.76
100000000006315498731.77
10000000000037501426711.88
10000000000051493226681.85
10000000000099485526971.80
1000000000000311574585311.85
1000000000000671602282641.94
1000000000000971586185301.86
100000000000003748950263461.86
100000000000009148836262101.86
100000000000015949008262561.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.

#### 1 comment:

Mark Amery said...

Testing this using the Guile interpreter, I got pretty similar results to you - the time taken to test a prime was reduced by a factor of almost - but nonetheless noticeably less than - 2.

I wasn't as quick as you to reach the conclusion that this was due to (next test-divisor) being more expensive to calculate than (+ test-divisor 1). The first possibility I considered - I wonder if you thought about and dismissed it as too implausible to be worth testing? - was that some other part of (find-divisor) ran significantly faster when testing even numbers than odd ones, and so skipping the even numbers had eliminated less than half of the true workload.

Testing revealed that your explanation here is the correct one, though, rather than my first guess. (find-divisor) is pretty simple, consisting solely of calls to (square), (divides?) and (next), and we can demonstrate (by timing lots of calls on a loop) that (square) and (divides?) both handle odd and even numbers equally fast (at least on my machine).

We can also demonstrate compellingly that the extra cost imposed by having to call (next) is at fault by redefining (next) to return (+ x 1) when x!=3 instead of (+ x 2), and then comparing the results of our benchmark against those from Ex-1.23 before (next) was introduced. Since we're then iterating over exactly the same numbers, any time difference can be definitively blamed upon us having changed the procedure that we use to increment test-divisor.

As you would expect, it turns out (on my machine) that the time taken to test each prime is greater using this modified version of (next) than it was in Ex-1.22, and the absolute time difference is pretty much exactly double (since (next) is called twice as many times) the difference between the actual and naively-expected results in this exercise.