Feeds:
Posts
Comments

## SICP Exercises Section 1.3.4

Exercise 1.40

Cubic creates and returns a procedure that takes an argument x and maps it to the value of x3 + ax2 + bx + c.

(define (cubic a b c)
(lambda (x)
(+ (* x x x)
(* a x x)
(* b x)
c)))


Even though this is a very simple example the smooth transition between algebraic notation and code is impressive.

Exercise 1.41

(define (double f)
(lambda (x)
(f (f x))))

(define (inc x)
(+ x 1))

(((double (double double)) inc) 5)
=> 21


The evaluation of the expression above results in 21 = 5 + 16. This might seem strange at first because doubling something successively 3 times is equivalent to multiplying it by 23 = 8. So we could expect to get 13 = 5 + 8. But this is incorrect. The outermost double procedure creates 2 copies of the 2 inner double’s resulting in a “cascade” of 4 consecutive double’s, for a total of 24 = 16 applications of the procedure f. This can be clearly seen by applying the substitution model of evaluation.

Exercise 1.42

The capability of a procedure to take another procedure(s) as argument and to return a procedure as a value makes the implementation of compose simple and elegant:

(define (compose f g)
(lambda (x)
(f (g x))))


Exercise 1.43

The returned procedure can in turn be passed as an argument to a procedure which returns a procedure. This “closure property” (in the mathematical sense) adds additional power to the language as can be seen in this and the next exercise.

(define (repeated f n)
(if (= n 1)
f
(compose f
(repeated f (- n 1)))))


Exercise 1.44

(define (average-3 x y z)
(/ (+ x y z) 3.0))

(define dx 0.001)

(define (smooth f)
(lambda (x)
(average-3 (f (- x dx))
(f x)
(f (+ x dx)))))

(define (n-fold-smooth n)
(repeated smooth n))


Exercise 1.45

As an experimental tool we can define a procedure n-th-root-test:

;; Tries to compute the n'th root of x using t average-damps.
(define (n-th-root-test n x t)
(fixed-point ((repeated average-damp t)
(lambda (y) (/ x (expt y (- n 1)))))
1.0))


It takes three arguments – n, x, and t. N determines the order of the root we want to compute, x is the number whose root we are computing, and t is the number of damping operations to be applied.

Running tests for n’s from 2 to 16 on x = 42 and taking t to be as small as possible we get:

;; Tests
(n-th-root-test 2 42 1)
(n-th-root-test 3 42 1)
(n-th-root-test 4 42 2) ; does not terminate with t = 1

(n-th-root-test 5 42 2)
(n-th-root-test 6 42 2)
(n-th-root-test 7 42 2)
(n-th-root-test 8 42 3) ; does not terminate with t = 2

(n-th-root-test 9 42 3)
(n-th-root-test 10 42 3)
(n-th-root-test 11 42 3)
(n-th-root-test 12 42 3)
(n-th-root-test 13 42 3)
(n-th-root-test 14 42 3)
(n-th-root-test 15 42 3)
(n-th-root-test 16 42 4) ; does not terminate with t = 3


A clear pattern emerges. Each time n reaches the next power of 2 one more average damping operation is needed to achieve convergence. In other words to compute the n-th root of x we need to apply ⌊ log2n⌋ damping operations. Embedding this knowledge into a new n-th-root procedure we get:

(define (n-th-root n x)
(fixed-point ((repeated average-damp
(floor (/ (log n) (log 2.0))))
(lambda (y) (/ x (expt y (- n 1)))))
1.0))


Exercise 1.46

Iterative-improve is slightly different from the procedures we saw so far. It can not return a lambda expression in the most straight forward manner. The reason is that the procedure returned is syntactically recursive (even though it generates an iterative process). So it must have a name:

(define (iterative-improve good-enough? improve)
(define (p guess)
(if (good-enough? guess)
guess
(p (improve guess))))
p)

(define (my-sqrt x)
(define (average x y)
(/ (+ x y) 2))
(define (good-enough? guess)
(< (abs (- (* guess guess) x)) 0.00001))
(define (improve guess)
(average guess (/ x guess)))
((iterative-improve good-enough? improve) x))

(define (fixed-point f first-guess)
(define (good-enough? guess)
(< (abs (- guess (f guess))) 0.00001))
(define (improve guess)
(f guess))
((iterative-improve good-enough? improve) first-guess))


Read Full Post »

## SICP Exercises Section 1.3.3

Exercise 1.35

Indeed φ is a fixed point of the transformation

$x \mapsto 1 + \dfrac{1}{x}$

This is a direct consequence of the fact that φ is a root of the equation:

$x^2 - x - 1 = 0$

Using the fixed-point procedure we get a good approximation of φ:

(fixed-point (lambda (x) (+ 1.0 (/ 1.0 x))) 1.0)
=> 1.6180327868852458


Exercise 1.36

(define (average x y)
(/ (+ x y) 2.0))

(define (tolerance) 0.00001)
(define (fixed-point f first-guess)
(define (close-enough? v1 v2)
(< (abs (- v1 v2)) (tolerance)))
(define (try guess)
(display "guess: ")
(display guess)
(newline)
(let ((next (f guess)))
(if (close-enough? guess next)
next
(try next))))
(try first-guess))

(fixed-point (lambda (x) (/ (log 1000.0) (log x))) 2.0)
=> 4.555532270803653 ; 34 guesses

(fixed-point (lambda (x) (average x (/ (log 1000.0) (log x))) 2.0))
=> 4.555537551999825 ; 9 guesses


The computation that uses average damping converges much faster. It takes only 9 guesses as opposed to the 34 needed for the non-damped one.

Exercise 1.37

Computing continued fractions turns out to be very cool. The recursive version is straight forward as usual. However the iterative implementation requires computing the continued fraction backwards! See the code below for more details:

(define (cont-frac n d k)
;; Recursive implementation
(define (rec i)
(if (> i k)
0
(/ (n i)
(+ (d i) (rec (+ i 1))))))

;; Iterative implementation
(define (iter i result)
(if (= i 0)
result
(iter (- i 1)
(/ (n i)
(+ (d i) result)))))

(iter k 0))


Experimenting a little we see that for the 12-term finite continued fraction with Ni = 1 and Di = 1 we get an approximation of 1 / φ accurate to 4 decimal places:

(cont-frac (lambda (i) 1.0)
(lambda (i) 1.0)
12)
=> 0.6180257510729613


Exercise 1.38

(define (den i)
(if (< (remainder i 3) 2)
1
(* 2 (/ (+ i 1) 3))))

(define e (+ 2 (cont-frac (lambda (i) 1.0)
den
100)))

e
=> 2.7182818284590455


This 100-term continued fraction approximation for e is accurate to 14 decimal places.

Exercise 1.39

(define (tan-cf x k)
(define (num i)
(if (= i 1)
x
(- (* x x))))
(define (den i)
(- (* 2 i) 1))
(cont-frac num den k))

(tan-cf 3.0 100)
=> -0.14254654307427775


Read Full Post »

## SICP Exercises Section 1.3.2

Exercise 1.34

Applying the substitution model of evaluation we get:

(f f) =>
(f 2) =>
(2 2)


The attempt to evaluate the last expression will result in an error since 2 is not a procedure the interpreter knows how to apply.

Read Full Post »

## SICP Exercises Section 1.3.1

Exercise 1.29

This one is just a simple mapping between the given formula and Scheme code using the higher order sum template:

(define (simpson f a b n)
(let ((h (/ (- b a) n)))
(define (simpson-term i)
(let ((y (f (+ a (* h i)))))
(cond ((= i 0) y)
((= i n) y)
((even? i) (* 2 y))
(else (* 4 y)))))
(define (simpson-next i)
(+ i 1))
(* (/ h 3.0) (sum simpson-term 0 simpson-next n))))

(simpson cube 0 1 100)
=> 0.25

(simpson cube 0 1 1000)
=> 0.25

(integral cube 0 1 0.0001)
=> 0.24999999874993412


The results speek for themselves with simpson being exact for both n = 100 and n = 1000. As can be expected the initial integral procedure is not as accurate.

Exercise 1.30

(define (sum term a next b)
(define (iter a result)
(if (> a b)
result
(iter (next a) (+ (term a) result))))
(iter a 0))


Exercise 1.31

The higher-order product procedure is very similar to sum:

(define (product term a next b)
(if (> a b)
1
(* (term a)
(product term (next a) next b))))


The only difference is that the base case returns 1 (the multiplicative identity, as opposed to the additive identity 0), and the consecutive terms are multiplied instead of added, This gives us a hint that there is an even more general abstraction waiting to be revealed.

Returning to the problem – computing factorial simply becomes:

(define (factorial n)
(product identity 1 inc n))


There are several ways to use product to compute an approximation of π based on the formula:

$\dfrac{\pi}{4} = \dfrac{2.4.4.6.6.8.\dots}{3.3.5.5.7.7.\dots}$

The one I chose groups the 4’s of the numerator with the 3’s of the denominator, the 6’s in the numerator with the 5’s in the denominator and so on. The general term of the product then becomes:

$\left(\dfrac{2i}{2i - 1}\right)^2$

for i starting at 2 and going on to some fixed number n. In addition we have to compensate for one factor of 2 in the numerator and for one factor of 2n + 1 in the denominator:

(define (compute-pi n)
(define (square x) (* x x))
(define (term-pi x)
(square (/ (* 2 x) (- (* 2 x) 1))))
(define (next-pi x)
(+ x 1))
(* (* (/ 2.0 (+ (* 2.0 n) 1))
(product term-pi 2 next-pi n))
4.0))

(compute-pi 1000)
=> 3.1408077460303945


Exercise 1.32

In this exercise we capture and express as a procedure a more general abstraction that spans both sum and product – the accumulate abstraction:

(define (accumulate combiner null-value term a next b)
(if (> a b)
null-value
(combiner (term a)
(accumulate combiner
null-value
term
(next a)
next
b))))


It begins to become apparent why it is so important for a language to be able to capture these higher-order abstractions. With accumulate it becomes trivial to express both sum:

(define (sum term a next b)
(accumulate + 0 term a next b))


and product:

(define (product term a next b)
(accumulate * 1 term a next b))


Here is the iterative version of the accumulate procedure:

(define (accumulate-iter combiner null-value term a next b)
(define (iter a result)
(if (> a b)
result
(iter (next a) (combiner (term a) result))))
(iter a null-value))


Exercise 1.33

The accumulate abstraction can be taken even further by adding a filter to determine which elements should selectively be accumulated:

(define (filtered-accumulate combiner null-value term a next b filter)
(if (> a b)
null-value
(let ((fa (filtered-accumulate combiner
null-value
term
(next a)
next
b
filter)))
(if (filter a)
(combiner (term a) fa)
fa))))


Now the sum of the squares of the prime numbers in the interval [a, b] can be expressed as:

(define (sum-squares-primes a b)
(filtered-accumulate + 0 square a inc b prime?))


And here is a procedure that computes the product of all positive integers less than n, that are relatively prime to n:

(define (rel-prime-product n)
(define (rel-prime a)
(= (gcd a n) 1))
(filtered-accumulate * 1 identity 1 inc (- n 1) rel-prime))


The ability to generalize using higher-order procedures is a powerful conceptual tool which allows us to write simple and elegant code.

Read Full Post »

## SICP Exercises Section 1.2.6

Exercise 1.21

(define (square n)
(* n n))

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

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


Exercise 1.22

Since I am using PLT Racket which does not support runtime, the tool I am relying upon in this and the next few exercises is time.

To compensate for Moore’s law and for the fact that time is measuring the CPU time not in microseconds but in milliseconds, I am running the primality test for each of the primes 10000 times. After all, it is always sound to base conclusions on data from more than just one or a few measurements.

The code and the results obtained follow:

;; Checks if n is prime in O(sqrt(n)).
(define (prime? n)
(= n (smallest-divisor n)))

;; Tests whether n is prime and if so prints it.
(define (prime-test n)
(if (prime? n)
(report-prime n)
#f))

;; Pretty prints a prime p.
(define (report-prime p)
(display "*** ")
(display p)
(newline)
#t)

;; Checks all odd numbers in the range [a, b]
;; for primallity and returns the number of
;; primes it has found.
(define (search-for-primes a b)
(cond ((> a b) 0)
((even? a) (search-for-primes (+ a 1) b))
(else (+ (if (prime-test a) 1 0)
(search-for-primes (+ a 2) b)))))

(search-for-primes 1000 1019) ; primes: 1009, 1013, 1019

(search-for-primes 10000 10037) ; primes: 10007, 10009, 10037

(search-for-primes 100000 100043) ; primes: 100003, 100019, 100037

(search-for-primes 1000000 1000037) ; primes: 1000003, 1000033, 1000037

;; Applies "proc" with argument "arg" "times" number of times.
(define (repeat proc arg times)
(if (= times 0)
'done
(begin
(proc arg)
(repeat proc arg (- times 1)))))

; 10^3
(time (repeat prime? 1009 10000))
=> cpu time: 110 real time: 109 gc time: 0

(time (repeat prime? 1013 10000))
=> cpu time: 94 real time: 94 gc time: 0

(time (repeat prime? 1019 10000))
=> cpu time: 109 real time: 109 gc time: 0

; 10^4
(time (repeat prime? 10007 10000))
=> cpu time: 344 real time: 344 gc time: 0

(time (repeat prime? 10009 10000))
=> cpu time: 328 real time: 328 gc time: 0

(time (repeat prime? 10037 10000))
=> cpu time: 329 real time: 328 gc time: 0

; 10^5
(time (repeat prime? 100003 10000))
=> cpu time: 1046 real time: 1047 gc time: 0

(time (repeat prime? 100019 10000))
=> cpu time: 1047 real time: 1047 gc time: 0

(time (repeat prime? 100043 10000))
=> cpu time: 1047 real time: 1047 gc time: 0

; 10^6
(time (repeat prime? 1000003 10000))
=> cpu time: 3297 real time: 3297 gc time: 0

(time (repeat prime? 1000033 10000))
=> cpu time: 3312 real time: 3312 gc time: 0

(time (repeat prime? 1000037 10000))
=> cpu time: 3297 real time: 3297 gc time: 0


This can be summarized in the following small table:

order 103     (109 + 94 + 109) / 3         ≈   104 msec
order 104     (344 + 328 + 328) / 3       ≈   333 msec
order 105     (1047 + 1047 + 1047) / 3 ≈ 1047 msec
order 106     (3297 + 3312 + 3297) / 3 ≈ 3302 msec

Comparing the adjacent rows we get:

333 / 104     ≈ 3.20
1047 / 333   ≈ 3.14
3302 / 1047 ≈ 3.15

The conclusion is that the “experimental data” strongly support the square root prediction (√10 ≈ 3.16).

Exercise 1.23

The modified find-divisor procedure:

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

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


Results obtained using the modified procedure:

; 10^3
(time (repeat prime? 1009 10000))
=> cpu time: 79 real time: 78 gc time: 0

(time (repeat prime? 1013 10000))
=> cpu time: 78 real time: 78 gc time: 0

(time (repeat prime? 1019 10000))
=> cpu time: 78 real time: 78 gc time: 0

; 10^4
(time (repeat prime? 10007 10000))
=> cpu time: 344 real time: 344 gc time: 94

(time (repeat prime? 10009 10000))
=> cpu time: 235 real time: 234 gc time: 0

(time (repeat prime? 10037 10000))
=> cpu time: 234 real time: 234 gc time: 0

; 10^5
(time (repeat prime? 100003 10000))
=> cpu time: 734 real time: 734 gc time: 0

(time (repeat prime? 100019 10000))
=> cpu time: 734 real time: 734 gc time: 0

(time (repeat prime? 100043 10000))
=> cpu time: 750 real time: 750 gc time: 0

; 10^6
(time (repeat prime? 1000003 10000))
=> cpu time: 2344 real time: 2344 gc time: 0

(time (repeat prime? 1000033 10000))
=> cpu time: 2344 real time: 2344 gc time: 0

(time (repeat prime? 1000037 10000))
=> cpu time: 2343 real time: 2344 gc time: 0


In summary:

order 103     (78 + 78 + 78) / 3             ≈     78 msec
order 104     (250 + 234 + 234) / 3       ≈   239 msec
order 105     (734 + 734 + 750) / 3       ≈   739 msec
order 106     (2344 + 2344 + 2344) / 3 ≈ 2344 msec

Comparing with the times from the previous exercise we get:

78 / 104       ≈ 0.75
239 / 333     ≈ 0.72
739 / 1047   ≈ 0.71
2344 / 3302 ≈ 0.71

The modification did not result in halving the runtime. This is not surprising. Indeed the old version of find-divisor performs twice as many divisibility tests as the new version. But at each step the new find-divisor is making a call to the next procedure. This results in evaluating an additional if expression, which turns out to be a substantial overhead.

Exercise 1.24

;; Runs the Fermat test 10 times.
(define (fast-prime-10? n)
(fast-prime? n 10))

; 10^3
(time (repeat fast-prime-10? 1009 10000))
=> cpu time: 1641 real time: 1640 gc time: 0

(time (repeat fast-prime-10? 1013 10000))
=> cpu time: 1703 real time: 1703 gc time: 0

(time (repeat fast-prime-10? 1019 10000))
=> cpu time: 1797 real time: 1797 gc time: 0

; 10^4
(time (repeat fast-prime-10? 10007 10000))
=> cpu time: 2219 real time: 2218 gc time: 0

(time (repeat fast-prime-10? 10009 10000))
=> cpu time: 2141 real time: 2141 gc time: 0

(time (repeat fast-prime-10? 10037 10000))
=> cpu time: 2203 real time: 2203 gc time: 0

; 10^5
(time (repeat fast-prime-10? 100003 10000))
=> cpu time: 5375 real time: 5375 gc time: 155

(time (repeat fast-prime-10? 100019 10000))
=> cpu time: 5282 real time: 5281 gc time: 187

(time (repeat fast-prime-10? 100043 10000))
=> cpu time: 5282 real time: 5281 gc time: 188

; 10^6
(time (repeat fast-prime-10? 1000003 10000))
=> cpu time: 9812 real time: 9812 gc time: 2563

(time (repeat fast-prime-10? 1000033 10000))
=> cpu time: 7375 real time: 7375 gc time: 172

(time (repeat fast-prime-10? 1000037 10000))
=> cpu time: 7672 real time: 7672 gc time: 234


In summary:

order 103     (1640 + 1703 + 1797) / 3 ≈ 1713 msec
order 104     (2218 + 2141 + 2203) / 3 ≈ 2187 msec
order 105     (5375 + 5281 + 5281) / 3 ≈ 5312 msec
order 106     (7249 + 7375 + 7672) / 3 ≈ 7432 msec

From the fact that the Fermat test has a time complexity of Θ(log n) we can predict that the time needed to test for primality a number close to 106 will be approximately twice the time to test a number close to 103, since log2(106) / log2(103) ≈ 2.

However, our empirical data show that this ratio is actually close to 4. The reason is a bit subtle.

Let’s consider the expmod procedure:

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


It involves a multiplication in which both the multiplier and the multiplicand can be as large as m – 1. The intermediate result after performing this multiplication can be close to m2. When dealing with m’s of order 103 no problems arise. But when m is of the order 106 the intermediate result can easily exceed the range of the standard 32-bit integer. In such cases the system has to do something about it. One option might be to switch to a 64-bit integer representation, another option is to switch to the representation used for arbitrary precision integers. No matter which option is chosen by George (or whoever implemented the system) some slowdown will inevitably result.

The change forced on the internal representation of integers when computing with numbers of order 106 can explain the observed anomaly.

Exercise 1.25

The problem with Ms. Hacker’s code lies in the fact, that the intermediate results in computing fast-expt can get very large very fast.

In Lisp systems integer overflow is not an issue since arbitrary precision integers are supported. But as the numbers grow the arithmetic operations on them can become very expensive.

On the other hand, the original expmod procedure never allows the intermediate results to exceed m2.

Exercise 1.26

Mr. Reasoner has fallen into a common “trap” in implementing procedures such as expmod. The mistake is that his code recursively solves the same subproblem twice. The correct way is to solve it only once, to store the result, and, when the time comes, to multiply it by itself. That is exactly what the version that uses square does. It stores the result by binding it to the formal parameter x.

The recurrence relation corresponding to Mr. Reasoner’s implementation is:

$T(exp) = 2T\left(\dfrac{exp}{2}\right)+ O(1)$

By case 1 of Master theorem:

$T(exp) = \Theta(exp)$

Exercise 1.27

The results bellow clearly demonstrate that the Carmichael numbers really do fool the Fermat test – proof by complete enumeration.

(define (pass-fermat n)
(define (pass-fermat-iter a)
(cond ((= a n) #t)
((= (expmod a n n) a) (pass-fermat-iter (+ a 1)))
(else #f)))
(pass-fermat-iter 1))

(pass-fermat 561)
=> t

(pass-fermat 1105)
=> t

(pass-fermat 1729)
=> t

(pass-fermat 2465)
=> t

(pass-fermat 2821)
=> t

(pass-fermat 6601)
=> t


Exercise 1.28

For diversity the implementation of Miller-Rabin is in Common Lisp. As the results show the Carmichael numbers don’t fool it as they did the Fermat test. (Well actually if you run the program locally there is a chance of less than 1 in 1000 that it will be tricked to believe that a Carmichael number, or any other composite number, is actually prime).

(defun square (x)
(* x x))

(defun expmod (base exp m)
(cond ((= exp 0) 1)
((evenp exp)
(let ((a (expmod base (/ exp 2) m)))
(if (and (not (= a 1))
(not (= a (- m 1)))
(= (rem (square a) m) 1))
0
(rem (square a) m))))
(t (rem (* base (expmod base (- exp 1) m)) m))))

(defun miller-rabin-test (n)
(let ((a (+ (random (- n 1)) 1)))
(if (= 1 (expmod a (- n 1) n))
t
nil)))

(defun fast-prime (n times)
(cond ((= times 0) t)
((miller-rabin-test n) (fast-prime n (- times 1)))
(t nil)))

(fast-prime 561 10)
=> NIL

(fast-prime 1105 10)
=> NIL

(fast-prime 1729 10)
=> NIL

(fast-prime 2465 10)
=> NIL

(fast-prime 2821 10)
=> NIL

(fast-prime 6601 10)
=> NIL

(fast-prime 103 10)
=> T


Read Full Post »

Older Posts »