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.

SICP Exercise 1.31: Product of a Series

From SICP section 1.3.1 Procedures as Arguments

Exercise 1.31 asks us to write a product procedure (analogous to sum) that computes the product of a function at points over a given range. We need to show how to define factorial in terms of the new product procedure, and use product to compute approximations to π using the following formula:


Finally, if our product procedure is recursive, we need to write an iterative version, and vice versa.

Since summing a series and multiplying a series are extremely similar ideas, it's not surprising to find that it's very simple to modify sum to produce product. Let's use the recursive version from exercise 1.29 first.

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

The product procedure really only needs to perform a different operation and end on a different value than sum. The last step (when a > b) should be to multiply by 1 instead of adding 0.

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

We can test this out by implementing factorial in terms of product using the identity and inc procedures that we've used before.

(define (identity x) x)

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

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

> (factorial 3)
6
> (factorial 4)
24
> (factorial 5)
120
> (factorial 10)
3628800
> (factorial 20)
2432902008176640000

The next test is to approximate π using something called Wallis' product or the Wallis Formula. (Lucky for us that mathematicians have already expressed this as the product of a series, saving us the work of coming up with a function for generating the terms.)


We'll multiply the product on the right hand side by the denominator on the left hand side so our procedure will just approximate π directly.

(define (wallis-pi n)
   (define (term x)
      (/ (* 4.0 (square x))
         (- (* 4.0 (square x)) 1)))
   (* 2.0 (product term 1 inc n)))

> (wallis-pi 10)
3.0677038066434994
> (wallis-pi 100)
3.133787490628163
> (wallis-pi 1000)
3.1408077460304042
> (wallis-pi 10000)
3.1415141186819313

Since our product procedure is implemented recursively, we need to show an iterative version. We can go back to the iterative sum procedure for inspiration.

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

Again, there are very few changes required to transform this procedure.

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

You can plug product-iter into the factorial and wallis-pi procedures above to verify that they give the same results.

The similarities between our two versions of sum and product is no accident. This section of SICP is all about higher-order procedures, so in the next exercise we're going to see how pull the similarities out of these procedures into something even more abstract.


Related:

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