Category: SICP

Related the classic text, ‘Structure and Interpretation of Computer Programs’.

Solution to SICP Exercise 1.26

Structure and Interpretation of Computer Programs

Solution to Exercise 1.26:

Because Louis’s version makes two calls to expmod for the even case instead of one, his process comes to resemble a binary tree of depth log n. The original creates a process that resembled a chain of length log n. When we count the number of nodes in Louis’s binary tree we get n – 1, bringing it back to Θ(n).

Solution to SICP Exercise 1.25

Structure and Interpretation of Computer Programs

Solution to Exercise 1.25:

Both the original and Alyssa’s version of expmod give the same result, but Alyssa’s version has much poorer performance. Here is a table comparing times spent (in milliseconds) calculating fast-prime? with the two different versions:

Prime number Original Alyssa’s
1009 608 2937
1013 578 3094
1019 687 2877
10007 874 93900
10009 875 92301
10037 749 94035

Times for Alyssa’s version are much higher than those for the original. Why is that?

As I couldn’t figure out how to get DrScheme’s profiler to measure time spent in built-in procedures such as remainder and *, I have to guess where the inefficiency is.

The original expmod calculates its result incrementally, calling remainder many times on relatively small values. Alyssa’s expmod, on the other hand, builds up a huge intermediate value, which is then passed to remainder once. My guess: because the numbers have grown so large in Alyssa’s version, the built-in operations, remainder, *, and / are consuming much more time than they do in the original.

Moral of the story: arbitrary precision arithmetic can be expensive.

Solution to SICP Exercise 1.24

Structure and Interpretation of Computer Programs

Solution to Exercise 1.24:

We need to modify start-prime-test to use fast-prime?. Here are my modifications:

(define (start-prime-test i n start-time)
(if (fast-prime? n 10)
(if (= i 0)
(report-prime
(time-difference
(current-time time-process)
start-time))
(start-prime-test (- i 1) n start-time))))

The inevitable problem when switching to fast-prime? is deciding what value to use for times. I’ve somewhat arbitrarily decided to use a constant value of 10 for now, which gives the following output:

Welcome to DrScheme, version 209.
Language: Textual (MzScheme, includes R5RS).

1009 *** 0s 2190000ns
1013 *** 0s 2340000ns
1019 *** 0s 7650000ns
(# # #)

10007 *** 0s 2810000ns
10009 *** 0s 2820000ns
10037 *** 0s 2810000ns
(# # #)

100003 *** 0s 5940000ns
100019 *** 0s 5470000ns
100043 *** 0s 5780000ns
(# # #)

1000003 *** 0s 6410000ns
1000033 *** 0s 3600000ns
1000037 *** 0s 3750000ns
(# # #)

As log1000=3 and log1000000=6, we should expect the processing times for the primes near 1000000 to be twice that of the primes near 1000.

The average time for the primes near 1000000 is 4586667ns. For those near 1000, 4060000.

These are very close; definitely not the factor of two we expected. What’s going on?

Luckily DrScheme has a built in profiler. Perhaps it can shed some light on what’s going on.

Here is the breakdown of time spent in each procedure for the values near 1000 (in milliseconds).

FunctionNumber 1009 1013 1019 Average
start-prime-test 672 766 797 745
fast-prime? 641 751 749 713.67
try-it 563 719 608 630
expmod 516 687 593 598.67
square 78 187 109 124.67
fermat-test 0 0 32 10.67
report-prime 0 0 0 0
timed-prime-test 0 0 0 0

And here are the numbers for the values near 1000000.

FunctionNumber 1000003 1000033 1000037 Average
start-prime-test 1546 1578 1719 1614.33
fast-prime? 1530 1562 1704 1598.67
try-it 1499 1436 1580 1505
expmod 1483 1420 1533 1478.67
square 284 392 235 303.67
fermat-test 15 32 15 20.67
report-prime 0 0 0 0
timed-prime-test 0 0 0 0

Comparing the results:

Function Ratio Deviation from expectations
start-prime-test 2.17 8%
fast-prime? 2.24 12%
try-it 2.39 19%
expmod 2.47 23%
square 2.44 22%
fermat-test 1.94 -3%

As you can see, the measurements aquired by the profiler match our expectations quite closely.

Lesson of the day: DrScheme’s profiler give much more accurate measurements than its implementation of SRFI 19.

Solution to SICP Exercise 1.23

Structure and Interpretation of Computer Programs

Solution to Exercise 1.23:

Here is the definition of a next procedure.

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

To use it, we need to modify find-divisor like so:

(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)))))

When I rerun the code from exercise 1.22, I get the following output:


1009 *** 0s 310000ns
1013 *** 0s 320000ns
1019 *** 0s 150000ns
(# # #)

10007 *** 0s 780000ns
10009 *** 0s 780000ns
10037 *** 0s 940000ns
(# # #)

100003 *** 0s 2500000ns
100019 *** 0s 7500000ns
100043 *** 0s 2500000ns
(# # #)

1000003 *** 0s 2190000ns
1000033 *** 0s 7810000ns
1000037 *** 0s 2190000ns
(# # #)

Compared to the results of exercise 1.22:

Number Ex 1.22 Time Ex 1.23 Time Ratio
1009 460000 310000 1.48
1013 470000 320000 1.47
1019 320000 150000 2.13
10007 1250000 780000 1.6
10009 1410000 780000 1.81
10037 1090000 940000 1.16
100003 5940000 2500000 2.38
100019 3910000 7500000 0.52
100043 4060000 2500000 1.62
1000003 7500000 2190000 3.42
1000033 1870000 7810000 0.24
1000037 2340000 2190000 1.07

If I throw away the top two and bottom two ratios, the average ratio is 1.54. This is close to but less than the expected speed up of 2. If the numbers are at all reliable — and I have my doubts — the extra 0.46 can be attributed to a combination of the procedure call to next (which could take longer than the call to + because it is user-defined versus built-in) and the additional if, with its test for equality to 2.

Solution to SICP Exercise 1.22

Structure and Interpretation of Computer Programs

Solution to Exercise 1.22:

DrScheme does not have a built-in runtime procedure that I could find, so I modified the code for timed-prime-test to use SRFI 19, like so

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

(define (timed-prime-test n)
(newline)
(display n)
(start-prime-test n (current-time time-process)))

(define (start-prime-test n start-time)
(if (prime? n)
(report-prime
(time-difference
(current-time time-process)
start-time))))

(define (report-prime elapsed-time)
(display " *** ")
(display (time-second elapsed-time))
(display "s ")
(display (time-nanosecond elapsed-time))
(display "ns"))

I implemented seach-for-primes as follows:

(define (search-for-next-prime starting-at)
(if (prime? starting-at)
starting-at
(search-for-next-prime (+ starting-at 2))))

(define (search-for-primes find-n starting-at)
(if (= find-n 0)
null
(let ((next-prime (search-for-next-prime starting-at)))
(cons next-prime
(search-for-primes (- find-n 1) (+ next-prime 2))))))

I could then run some tests:

(define (time-prime-tests primes)
(map timed-prime-test primes))

(time-prime-tests (search-for-primes 3 1001))
(time-prime-tests (search-for-primes 3 10001))
(time-prime-tests (search-for-primes 3 100001))
(time-prime-tests (search-for-primes 3 1000001))

This gave the following output:


1009 *** 0s 0ns
1013 *** 0s 0ns
1019 *** 0s 0ns
(# # #)

10007 *** 0s 0ns
10009 *** 0s 0ns
10037 *** 0s 0ns
(# # #)

100003 *** 0s 0ns
100019 *** 0s 0ns
100043 *** 0s 0ns
(# # #)

1000003 *** 0s 0ns
1000033 *** 0s 0ns
1000037 *** 0s 0ns
(# # #)

As you can see, all of the tests for prime took no time at all. This is obviously wrong.

I have a feeling that all these processes complete in less time than a single tick of the process clock (which, on Windows, appears to have a resolution in milliseconds).

Let’s magnify the results by calling prime? in timed-prime-test 1000 times instead of just once. Here’s the modified definition:

(define (timed-prime-test n)
(newline)
(display n)
(start-prime-test 1000 n (current-time time-process)))
(define (start-prime-test i n start-time)
(if (prime? n)
(if (= i 0)
(report-prime
(time-difference
(current-time time-process)
start-time))
(start-prime-test (- i 1) n start-time))))

Now I get the following output:


1009 *** 0s 460000ns
1013 *** 0s 470000ns
1019 *** 0s 320000ns
(# # #)

10007 *** 0s 1250000ns
10009 *** 0s 1410000ns
10037 *** 0s 1090000ns
(# # #)

100003 *** 0s 5940000ns
100019 *** 0s 3910000ns
100043 *** 0s 4060000ns
(# # #)

1000003 *** 0s 7500000ns
1000033 *** 0s 1870000ns
1000037 *** 0s 2340000ns
(# # #)

Alright! We have some numbers. Unfortunately, when I run the tests again, I get totally different numbers. Grr! I’m going to consider these numbers representative of many runs, but they probably aren’t. Let’s start with some analysis.

Let’s start by averaging the times for each group. We can multiply this by √10 to obtain an expected time for the following group, which we can compare against the measured time.

Group Average Time Expected time (based on lower group) Difference
~1000 416667 — —
~10000 1250000 1317616 -5%
~100000 4636667 3952847 17%
~1000000 3903333 14662427 -73%

With a difference of over 70% on the final group, this data clearly does not support the √n hypothesis, but it doesn’t disprove it, either. For a definitive answer, we need a larger sample group. We need data for more values of n and for several runs. I’m not going to do that here because this post is already too long, and I’m sure you have better things to do.

Solution to SICP Exercise 1.20

Structure and Interpretation of Computer Programs

Solution to Exercise 1.20:

;;; Normal Order

; remainder count: 0
(gcd 206 40)

; rc: 0
(if (= 40 0)
206
(gcd 40 (remainder 206 40)))

; rc: 0
(gcd 40 (remainder 206 40))

; rc: 0
(if (= (remainder 206 40) 0)
40
(gcd (remainder 206 40)
(remainder 40 (remainder 206 40))))

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

; rc: 1
(gcd (remainder 206 40)
(remainder 40 (remainder 206 40)))

; rc: 1
(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)))))

; rc: 2
(if (= (remainder 40 6) 0)
(remainder 206 40)
(gcd (remainder 40 (remainder 206 40))
(remainder (remainder 206 40)
(remainder 40
(remainder 206 40)))))

; rc: 3
(if (= 4 0)
(remainder 206 40)
(gcd (remainder 40 (remainder 206 40))
(remainder (remainder 206 40)
(remainder 40
(remainder 206 40)))))

; rc: 3
(gcd (remainder 40 (remainder 206 40))
(remainder (remainder 206 40)
(remainder 40
(remainder 206 40))))

; rc: 3
(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))))))

; rc: 4
(if (= (remainder 6
(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))))))

; rc: 5
(if (= (remainder 6
(remainder 40
6))
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))))))

; rc: 6
(if (= (remainder 6 4) 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))))))

; rc: 7
(if (= 2 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))))))

; rc: 7
(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)))))

; rc: 7
(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)))))))

; rc: 14 (there were 7 calls to remainder in the = form above)
(if (= 0 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)))))))

; rc: 14
(remainder (remainder 206 40)
(remainder 40
(remainder 206 40)))

; rc: 15
(remainder 6
(remainder 40
(remainder 206 40)))

; rc: 16
(remainder 6 (remainder 40 6))

; rc: 17
(remainder 6 4)

; rc: 18
2

For normal order, remainder is performed 18 times.

;;; Applicative Order

; rc: 0
(gcd 206 40)

; rc: 0
(if (= 40 0)
206
(gcd 40 (remainder 206 40)))

; rc: 0
(gcd 40 (remainder 206 40))

; rc: 1
(gcd 40 6)

; rc: 1
(if (= 6 0)
40
(gcd 6 (remainder 40 6)))

; rc: 1
(gcd 6 (remainder 40 6))

; rc: 2
(gcd 6 4)

; rc: 2
(if (= 4 0)
6
(gcd 4 (remainder 6 4)))

; rc: 2
(gcd 4 (remainder 6 4))

; rc: 3
(gcd 4 2)

; rc: 3
(if (= 2 0)
4
(gcd 2 (remainder 4 2)))

; rc: 3
(gcd 2 (remainder 4 2))

; rc: 4
(gcd 2 0)

; rc: 4
(if (= 0 0)
2
(gcd 0 (remainder 2 0)))

; rc: 4
2

For applicative order, remainder is only performed four times.