SICP Exercise 1.38: Calculating e as a k-Term Continued Fraction

In 1737, the Swiss mathematician Leonhard Euler published a memoir De Fractionibus Continuis, which included a continued fraction expansion for e - 2, where e is the base of the natural logarithms. In this fraction, the N_i are all 1, and the D_i are successively 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, .... Write a program that uses your cont-frac procedure from exercise 1.37 to approximate e, based on Euler's expansion.

We know from the previous exercise, where we used cont-frac to approximate 1/ϕ, that if N_i are all 1 then we can express this using (lambda (i) 1.0). So all we have to do to approximate e is to figure out how to calculate the sequence for D_i (namely 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, ...), produce a lambda expression for this, plug them into cont-frac and add 2 to the result.

Let's have a look at the sequence: 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, .... Notice how the numbers come in groups of 3: {1, 2, 1}, {1, 4, 1}, {1, 6, 1}, {1, 8, 1}, .... In each group of three the first and third elements are 1, while the middle element is 2g, where g is the index of the group. How can we calculate this from the sequence 1, 2, 3, ..., k?

Well, first of all, note that we can determine whether or not term i is the middle of a group of 3 by looking at i modulo 3: if it's in the middle of a group of 3 then i modulo 3 will be 2, and we need to calculate 2g, otherwise it's the first or third element of the group, and so D_i will be 1. If it's in the middle of a group of 3 then we can calculate the group number by adding 1 to it (which will make it a multiple of 3) and then dividing by three.

Using this we can produce a procedure for approximating e:

(define (approximate-e k)
  (+ (cont-frac (lambda (i) 1.0)
                (lambda (i) (if (= (remainder i 3) 2)
                                (* 2 (/ (+ i 1) 3))
                                1))
                k)
     2))

And then we can put it into practice. e is approximately 2.7183 to four decimal places, so let's see how quickly our procedure produces a number of similar accuracy:

> (approximate-e 1)
3.0
> (approximate-e 2)
2.6666666666666665
> (approximate-e 3)
2.75
> (approximate-e 4)
2.7142857142857144
> (approximate-e 5)
2.71875
> (approximate-e 6)
2.717948717948718
> (approximate-e 7)
2.7183098591549295

SICP Exercise 1.37: Calculating ϕ as a k-Term Continued Fraction

1. An infinite continued fraction is an expression of the form

   f =          N1         
       D1 +       N2      
            D2 +    N3   
                 D3 + …
   

   As an example, one can show that the infinite continued fraction expansion with the N_i and the D_i all equal to 1 produces 1/ϕ, where ϕ is the golden ratio (described in section 1.2.2). One way to approximate an infinite continued fraction is to truncate the expansion after a given number of terms. Such a truncation -- a so-called k-term finite continued fraction -- has the form

         N1      
   D1 +    N2   
        ⋱ + NK
             DK
   

   Suppose that n and d are procedures of one argument (the term index i) that return the N_i and D_i of the terms of the continued fraction. Define a procedure cont-frac such that evaluating (cont-frac n d k) computes the value of the k-term finite continued fraction. Check your procedure by approximating 1/ϕ using

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

   for successive values of k. How large must you make k in order to get an approximation that is accurate to 4 decimal places?
2. If your cont-frac procedure generates a recursive process, write one that generates an iterative process. If it generates an iterative process, write one that generates a recursive process.

Let's start with the iterative approach for a change...

If we look at the formula for a k-term continued fraction we can see that it looks "rather difficult" to start from the 1^st term and accumulate the result as we go until we reach the k^th term. Using the 1^st terms we could compute N₁/D₁. To then iterate on this value and accumulate the 2^nd terms into this we would need to know both N₁ and D₁ so we could separate out the numerator and denominator and then add N₂/D₂ to the denominator. To accumulate the third term we'd need to know N₁, N₂, D₁ and D₂, and so on. In other words, for any i, 0 < i ≤ k, there is not a separably calculable portion of the k-term continued fraction that contains only the terms 1…i.

What about if we flip it around. We know we can calculate the k^th term without knowing any of the preceeding terms. It's just N_k/D_k. Using this value and the values of N_k-1/D_k-1 we can then calculate the portion of the total sum that is made up of the k^th and (k - 1)^th terms alone. In other words we can calculate:

   Nk-1  
Dk-1 + Nk
      DK

We can then iterate all the way down to the 1^st term, accumulating as we go and produce the result we're after. So:

• We have an iterative step: for any value i, 0 < i ≤ k, divide N_i by the sum of D_i and the accumulated results of the terms i+1 … k.
• We know that we're only accumulating for values in the range 1 … k, so for our initial step (i.e. i = k), we can use 0 as the accumulated results of the terms i+1 … k
• Finally, we know when to terminate: when we get to i = 0 we've accumulated all the terms we need to.

Putting these together we get:

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

Let's now use this to calculate approximations to 1/ϕ. We're asked to try to find the value of k that produces an approximation to 1/ϕ that is accurate to 4 decimal places. 1/ϕ = 2 / (1 + √5) ≈ 0.6180 to four decimal places, so let's start at k = 1 and stop when we find our answer:

> (cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 1)
1.0
> (cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 2)
0.5
> (cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 3)
0.6666666666666666
> (cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 4)
0.6000000000000001
> (cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 5)
0.625
> (cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 6)
0.6153846153846154
> (cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 7)
0.6190476190476191
> (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
> (cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 11)
0.6180555555555556
> (cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 12)
0.6180257510729613

So we need k = 12 to get an approximation to 1/ϕ that is accurate to 4 decimal places.

For part (b) we're asked to produce the corresponding iterative or recursive procedure depending upon whether we originally produced a recursive or iterative procedure respectively. Well we produced an iterative procedure to begin with, so we need to produce a recursive procedure here.

Note that we can't directly translate the iterative procedure to a recursive one (i.e. going through the terms in the order k…1) for a similar reason to why the iterative procedure has to go through the terms in that order. For the recursive procedure we want to calculate the i^th term by performing some operation on the next term in the sequence in the order we are traversing it. If we were to follow the same traversal order as the iterative process that would mean that we would want to calculate the i^th term by performing some operation on the (i-1)^th term. But, by the definition of the formula, to calculate the (i-1)^th term we have to first calculate the i^th term, and we would have an infinite regress.

Instead simply notice that the formula for a k-term continued fraction can be described easily as a recursive process: for any i, 0 < i ≤ k, the portion of the formula produced by the i^th to k^th terms can be calculated by dividing N_i by the sum of D_i and the portion of the formula produced by the (i+1)^th to k^th terms. All we need to do is ensure that we terminate the process after we've calculated the k^th term and we're done. We can do this by having the recursion return 0 when i>k:

(define (cont-frac n d k)
  (define (recur i)
    (if (> i k)
        0
        (/ (n i) (+ (d i) (recur (+ i 1))))))
  (recur 1))

SICP Exercise 1.36: Average Damping

Modify fixed-point so that it prints the sequence of approximations it generates, using the newline and display primitives shown in exercise 1.22. Then find a solution to x^x = 1000 by finding a fixed point of x → log(1000)/log(x). (Use Scheme's primitive log procedure, which computes natural logarithms.) Compare the number of steps this takes with and without average damping. (Note that you cannot start fixed-point with a guess of 1, as this would cause division by log(1) = 0.)

Here's the original fixed-point:

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

Let's report a couple of things: the initial guess, and the value of next calculated at each call to try:

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

We can now use this to find a fixed point of x → log(1000)/log(x) with and without average damping. We're reminded that we can't start with a guess of 1, so let's pick 1.1 as our initial guess.

To find the fixed point without average damping we simply translate the transformation directly:

> (fixed-point (lambda (x) (/ (log 1000) (log x))) 1.1)
first-guess: 1.1
next: 72.47657378429035
next: 1.6127318474109593
next: 14.45350138636525
next: 2.5862669415385087
next: 7.269672273367045
next: 3.4822383620848467
next: 5.536500810236703
next: 4.036406406288111
next: 4.95053682041456
next: 4.318707390180805
next: 4.721778787145103
next: 4.450341068884912
next: 4.626821434106115
next: 4.509360945293209
next: 4.586349500915509
next: 4.535372639594589
next: 4.568901484845316
next: 4.546751100777536
next: 4.561341971741742
next: 4.551712230641226
next: 4.558059671677587
next: 4.55387226495538
next: 4.556633177654167
next: 4.554812144696459
next: 4.556012967736543
next: 4.555220997683307
next: 4.555743265552239
next: 4.555398830243649
next: 4.555625974816275
next: 4.555476175432173
next: 4.555574964557791
next: 4.555509814636753
next: 4.555552779647764
next: 4.555524444961165
next: 4.555543131130589
next: 4.555530807938518
next: 4.555538934848503
4.555538934848503

That took 37 calculations of next to find a fixed point that's within the set tolerance.

Now let's do it with average damping. To do this we need to average the preceding guess, x, with the current guess, log(1000)/log(x). Assuming we have the procedure average, which is defined earlier on in the book, we can do this as follows:

> (fixed-point (lambda (x) (average x (/ (log 1000) (log x)))) 1.1)
first-guess: 1.1
next: 36.78828689214517
next: 19.352175531882512
next: 10.84183367957568
next: 6.870048352141772
next: 5.227224961967156
next: 4.701960195159289
next: 4.582196773201124
next: 4.560134229703681
next: 4.5563204194309606
next: 4.555669361784037
next: 4.555558462975639
next: 4.55553957996306
next: 4.555536364911781
4.555536364911781

Here we have 13 calculations of next to find a fixed point that's within the set tolerance. So we see that average damping causes the function to converge quicker.