SICP 1.37, 1.38, and 1.39: Continued Fractions

From SICP section 1.3.3 Procedures as General Methods

Exercise 1.37 explains that an infinite continued fraction is an expression of the form

The infinite continued fraction where all Ni and Di terms are equal to 1 produces 1/ϕ, where ϕ is the golden ratio described in SICP section 1.2.2.

We can approximate the value of an infinite continued fraction by truncating to a given number of terms.

We are asked to define a procedure cont-frac that takes the arguments n, d, and k. n and d are procedures of one argument (the term index i) that compute the N_i and D_i terms of the continued fraction. k is the number of terms to expand.

We can check our procedure by approximating 1/ϕ using

(cont-frac (lambda (i) 1.0)
       (lambda (i) 1.0)
       k)

for successive values of k. We're asked to show how large k must be in order to approximate 1/ϕ to 4 decimal places. As usual, we're also asked to write two versions of the procedure, one that generates a recursive process and one iterative.

Since we're given the ending value of k and it won't do to count down from k to 1 as we compute the accumulated value, we'll define a helper function frac that counts up from 1 to k. I'll show the recursive version first since that comes more naturally.

(define (cont-frac n d k)
  (define (frac i)
     (if (< i k)
         (/ (n i) (+ (d i) (frac (+ i 1))))
         (/ (n i) (d i))))
  (frac 1))

The cont-frac procedure simply defines the helper function frac then calls it with its starting value of 1. The frac procedure checks to see if i is less than k. If so, it divides the result of applying the n procedure to i by the sum of the result of applying the d procedure to i and a recursive call to frac with the next value of i. If i is equal to k, the recursion ends and one last division operation is performed.

We can check the results by using the supplied code with an arbitrary value of k.

> (cont-frac (lambda (i) 1.0)
            (lambda (i) 1.0)
            5)
0.625

This gives us a value of 1/ϕ that's in the same ball park as the expected value of about 0.61803, but we were asked to find a value of k where the result is accurate to 4 decimal places. We can just test a few more values to find the answer.

> (cont-frac (lambda (i) 1.0)
         (lambda (i) 1.0)
         8)
0.6176470588235294
> (cont-frac (lambda (i) 1.0)
         (lambda (i) 1.0)
         9)
0.6181818181818182
> (cont-frac (lambda (i) 1.0)
         (lambda (i) 1.0)
         10)
0.6179775280898876

We can define an iterative version of cont-frac by modifying the helper function. Since the iterative procedure will be carrying the intermediate result at each step as a parameter, this time we'll count down from k to 0 instead of counting up from 1 to k as we did before.

(define (cont-frac-iter n d k)
   (define (frac-iter i result)
       (if (= i 0)
           result
           (frac-iter (- i 1) (/ (n i) (+ (d i) result)))))
   (frac-iter (- k 1) (/ (n k) (d k))))

You can run cont-frac-iter with the same inputs that we used above to verify that you get the same results.



Exercise 1.38 asks us to write a program that uses our cont-frac procedure to approximate e, the base of the natural logarithm. We're to use the continued fraction expansion for e - 2 published by Euler in 1737. In this fraction all of the N_i terms are 1, and the D_i terms are

1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8,...

Since cont-frac was already written in the last exercise, this problem is reduced to writing a good function for computing the D_i terms. The series is extremely regular, except for having only one "1" at the beginning. The indices of the non-1 values in the series are 2, 5, 8, 11, 14, 17,... These values are all one less than a multiple of three, or 3i - 1. So we know to begin with that the procedure (d i) should return a 1 when (i + 1) is not divisible by 3.

Since the series is so regular, it's easy to find a formula for the remaining values as well. You could do it by plotting the indices and their values in a spreadsheet, but that's not really necessary. First, take a look at the two sequences side-by-side.

i = 2, 5, 8, 11, 14, 17, ...
d(i) = 2, 4, 6, 8, 10, 12, ...

What I notice right away is that when I add 3 to the index, the result only increases by 2. This means that I should be able to get the result by dividing the index by 3 and multiplying by 2 (after applying an offset, which is 1 in this case). In mathematical terms:

d(i) = 2(i + 1) / 3

Putting it all together in code, it looks like the following:

(define (d i)
   (if (not (= 0 (remainder (+ i 1) 3)))
       1
       (* 2 (/ (+ i 1) 3))))

Now to compute e, we can use 1 for each N_i term and our new d procedure for each D_i term in a call to cont-frac. Remember that Euler's continued fraction computed e - 2, so we need to take that into account.

(define e
   (+ 2 (cont-frac (lambda (i) 1.0) d 10)))

> e
2.7182817182817183

Exercise 1.39 asks us to define a procedure (tan-cf x k) that computes an approximation to the tangent function based on the following continued fraction (published in 1770 by Johann Heinrich Lambert).

The rule for generating the N_i term is fairly simple. If i = 1, the term is equal to x. Otherwise it's -x². We have to negate all but the first term because our cont-frac procedure adds each term and we need to subtract them in this case.

The rule for the D_i term is even simpler, since they're just the odd numbers. D_i = 2(i) - 1.

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

(define (tan-cf x k)
   (define (n k)
       (if (= k 1)
           x
           (- (square x))))
   (define (d k)
       (- (* 2 k) 1))
   (cont-frac n d k))

We can use Scheme's built-in tan function to check our work with a few common angles (remember the angle is in radians). We'll stick with k = 10 terms of the continued fraction since we had such good results with that before.

> (tan (/ pi 6))
0.5773502691896257
> (tan-cf (/ pi 6) 10)
0.5773502691896257
> (tan (/ pi 4))
0.9999999999999999
> (tan-cf (/ pi 4) 10)
1.0
> (tan (/ pi 3))
1.7320508075688767
> (tan-cf (/ pi 3) 10)
1.732050807568877

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.36: Fixed points and Average damping

From SICP section 1.3.3 Procedures as General Methods

Exercise 1.36 asks us to modify fixed-point so that it prints the sequence of approximations it generates. We can do that with the newline and display primitives we saw in exercise 1.22.

Next we're asked to find a solution to x^x = 1000 by finding a fixed point of x → log(1000) / log(x). We can use Scheme's primitive log procedure to compute natural logarithms.

Finally, we need to compare the number of steps it takes to find the fixed point with and without average damping. We'll use the simple average procedure defined in the lecture notes.

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

We'll start with the fixed-point procedure we used in the last exercise. The try sub-procedure looks like a good place to print each guess.

(define tolerance 0.00001)

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

Now we can find the value of x for x^x = 1000 using this procedure. The book warns us not to use 1.0 as a starting value or else the procedure will attempt to divide by log(1) = 0, so let's start with an initial guess of 2.0.

> (fixed-point (lambda (x) (/ (log 1000) (log x))) 2.0)
2.0
9.965784284662087
3.004472209841214
6.279195757507157
3.759850702401539
5.215843784925895
4.182207192401397
4.8277650983445906
4.387593384662677
4.671250085763899
4.481403616895052
4.6053657460929
4.5230849678718865
4.577114682047341
4.541382480151454
4.564903245230833
4.549372679303342
4.559606491913287
4.552853875788271
4.557305529748263
4.554369064436181
4.556305311532999
4.555028263573554
4.555870396702851
4.555315001192079
4.5556812635433275
4.555439715736846
4.555599009998291
4.555493957531389
4.555563237292884
4.555517548417651
4.555547679306398
4.555527808516254
4.555540912917957
4.555532270803653

The procedure took 35 steps to arrive at an answer of 4.555532270803653 without using average damping.

In order to use average damping we just need to modify the input function to average the input value of x with the computed value.

> (fixed-point (lambda (x) (average x (/ (log 1000) (log x)))) 2.0)
2.0
5.9828921423310435
4.922168721308343
4.628224318195455
4.568346513136242
4.5577305909237005
4.555909809045131
4.555599411610624
4.5555465521473675
4.555537551999825

This time it took only 10 steps to come up with approximately the same answer. We can check the two answers by simply raising each result to itself using Scheme's expt primitive.

> (expt 4.555532270803653 4.555532270803653)
999.9913579312362
> (expt 4.555537551999825 4.555537551999825)
1000.0046472054871

So in this case, the procedure using average damping not only arrived at an solution in a fewer number of steps, but the final result happened to be slightly more accurate as well. That won't always be the case, but it's good to know that we're not sacrificing accuracy by using averaging damping.


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.35: Fixed points and the Golden ratio

From SICP section 1.3.3 Procedures as General Methods

Exercise 1.35 asks us to show that the golden ratio ϕ is a fixed point of the transformation x → 1 + 1 / x. We're then asked to use this fact to compute ϕ by means of the fixed-point procedure defined earlier in the chapter.

First we can show that ϕ is one of the roots of x → 1 + 1 / x. We learned the value of ϕ in section 1.2.2 is (1 + √5) / 2, or approximately 1.618.

x = 1 + 1 / x

If we multiply both sides by x we get:

x² = x + 1
x² - x - 1 = 0

Now if we use the quadratic equation (or cheat and use WolframAlpha like I did) we find that the roots of the equation are:

x = 1/2(1 - √5)
x = 1/2(1 + √5)

Computing ϕ by means of the fixed-point procedure is fairly straightforward since the procedure is already given.

(define tolerance 0.00001)

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

The fixed-point procedure takes a function and an initial guess. Since we're trying to find a point where

x = 1 + 1 / x

that's the function we need to pass. The initial guess can really be just about anything, but we already know that we're trying to prove the fixed point is around 1.6, so let's try an initial value close to that.

> (fixed-point (lambda (x) (+ 1 (/ 1 x))) 2.0)
1.6180327868852458

Related:

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

Interesting Miscellany

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

Math
Planck found in "Euler Identity" Crop Circle?!
Ten of the greatest: Math Puzzles
Russian math genius ignores $1 million Millennium Prize
A mathematician's clock...
werewolves and star wars: two exam questions

Science & Technology
Eyeborg bionic eye camera shows winks and all
Update on the diagnosis of the Voyager 2 data system - If you're in IT, this is the ultimate "Turn it off then back on again."

Programming
You don't need anyone's permission to get work experience in software.
Podcast interview with Donald Knuth. - This is a phone interview between Knuth and Larry Felton Johnson. They talk mostly about Literate Programming, but they also touch on TAOCP and a few other topics.


Random Statistics of the week
The LA Times reports that California welfare debit cards are accepted in over half of casinos in the state. Welfare recipients withdrew $1.8 million from casino ATMs between October 2009 and May 2010.


Finally, what I've been working on lately...
StackWrap4J, a Java wrapper for the recently released Stack Exchange API, has finally reached a semi-stable state. I teamed up with another prominent member of the Stack Overflow community, Justin 'jjnguy' Nelson, to work on this over the past couple of months and we've just released version 0.9. You can download it from the StackWrap4J project page on SourceForge. Please give us some feedback.

Now that the API is near stability, I should have more time soon to get back to SICP and my (ir)regular blogging schedule.

SICP Exercise 1.34: Procedures as Arguments

From SICP section 1.3.2 Constructing Procedures Using Lambda

Exercise 1.34 asks us to consider the following procedure:

(define (f g)
  (g 2))

This procedure takes a function as an argument and applies that function to the value 2. We're shown a couple of examples of the procedure in action.

> (f square)
4
> (f (lambda (z) (* z (+ z 1))))
6

We're then asked to consider what would happen if we applied f to itself.

> (f f)
. . procedure application: expected procedure, given: 2; arguments were: 2

We can use the substitution model to explain this failure.

(f f)
(f 2)
(2 2)

In the first substitution the argument f (a procedure) is applied to the value 2, so a recursive call is made. In the second substitution, the argument 2 is applied to 2. Since 2 isn't a procedure, an error is reported.


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.33: Filtered Accumulator

From SICP section 1.3.1 Procedures as Arguments

Exercise 1.33 asks us to us to write an even more general form of the accumulate procedure that we wrote in 1.32. A filtered accumulate procedure should combine only those terms in a specified range that meet a specified condition. The new procedure will take the same arguments as the old one, plus an additional argument that specifies the filtering function.

Once we've written the new procedure we need to test it by using it to write two additional procedures. The first will find the sum of the squares of the prime numbers in a given range, and the second will find the product of all the positive integers less than n that are relatively prime to n.

The filtered-accumulate procedure should be very similar to what we saw in the last exercise. The only difference is that we check each new value to see if it passes through the filter before applying the combiner function.

(define (filtered-accum filter combiner null-value term a next b)
  (if (> a b)
      null-value
      (if (filter a)
          (combiner (term a)
                    (filtered-accum filter combiner null-value term (next a) next b))
          (filtered-accum filter combiner null-value term (next a) next b))))

Note that we're filtering on the value of a, not the value of (term a).

Sum of squared primes

We came up with several different ways to test if a number is prime in section 1.2.6. You can choose any implementation you like, but I'm going to use the fast-prime? procedure from exercises 1.22 and 1.23 (shown here in block structure).

(define (fast-prime? n)
   (define (smallest-divisor n)
      (define (find-divisor n test-divisor)
         (define (next x)
            (if (= x 2) 3 (+ x 2)))
         (define (divides? a b)
            (= (remainder b a) 0))
         (cond ((> (square test-divisor) n) n)
               ((divides? test-divisor n) test-divisor)
               (else (find-divisor n (next test-divisor)))))
      (find-divisor n 2))
   (= n (smallest-divisor n)))

Using fast-prime? as a filtering function, we can implement a procedure to sum the square of primes in a given range.

(define (sum-squared-primes a b)
   (filtered-accum fast-prime? + 0 square a inc b))

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

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

We can use a few small examples to test with.

2² + 3² = 13
2² + 3² + 5² = 38
2² + 3² + 5² + 7² = 87

> (sum-squared-primes 2 3)
13
> (sum-squared-primes 2 6)
38
> (sum-squared-primes 2 10)
87

Procuct of coprimes

The second challenge was to find the product of all the positive integers less than n that are relatively prime (coprime) to n. In other words, multiply together all positive integers i < n such that

gcd(i, n) = 1

As luck would have it, we've already seen a gcd procedure too. SICP section 1.2.5 was all about finding greatest common divisors using Euclid's algorithm.

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

We can use gcd to define a coprime? procedure, and use that as our filter. Note that normally a coprime? procedure should take two parameters. Since filtered-accum expects the filter function to take only one parameter, coprime? is borrowing one of its input parameters, n, from product-of-coprimes.

(define (product-of-coprimes n)
  (define (coprime? i)
     (= 1 (gcd i n)))
  (filtered-accum coprime? * 1 identity 1 inc (- n 1)))

(define (identity x) x)

Again, I'm only going to test with a few small samples, since the results can potentially get very big, very fast. A quick mental check reveals that 3, 7, and 9 are the only values less than 10 that are also coprime to 10.

> (product-of-coprimes 10)
189

All values less than a prime number are coprime to that number, so the product of coprimes to 11 should result in the product of all the values from 2 to 10, or 10!.

> (product-of-coprimes 11)
3628800

You can check that and any other inputs with a calculator.


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.32: Accumulator

From SICP section 1.3.1 Procedures as Arguments

Exercise 1.32 asks us to show how sum and product are special cases of an even more abstract concept called accumulate that combines a collection of terms. We're given the following procedure signature to start with:

(accumulate combiner null-value term a next b)

The arguments term, a, next, and b serve the same purpose as they did in sum and product. The new arguments are combiner, which takes a procedure of two arguments that specifies how the current term should be combined with the accumulation of all the preceding terms, and null-value, which specifies what base value to use when the terms run out.

We need to test accumulate by implementing both sum and product in terms of the new procedure. We also need to implement both recursive and iterative versions of accumulate.

We saw in exercise 1.31 how similar sum and product are. To create a higher-order procedure that can be used to implement both of these ideas, we need to look at where they're different. Let's look at the recursive version first.

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

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

At a casual glance these two functions look almost exactly the same. They're only different in name, the null value used when a > b, and the operator used to combine terms.

(define (<name> term a next b)
(if (> a b)
<null-value>
(<operator> (term a)
(<name> term (next a) next b))))

These are exactly the parameters that need to change to create a higher-order accumulate procedure.

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

Given their similarities, redefining sum and product in terms of accumulate is easy.

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

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

An iterative version can be created using the same technique of substituting what's different about the iterative versions of sum and product from previous exercises.

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

Any of these procedures can be tested using the same tests from the previous exercises and comparing the results.


Related:

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