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.