Saturday, May 1, 2010

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.

3 comments:

Guilherme said...

Hi Bill,

I've made another solution. There you go:

(define (pi4 n)
(define (term k)
(if (even? k)
(/ (+ 2 k) (+ 1 k))
(/ (+ 1 k) (+ 2 k))
))
(product term 1 inc n))

Anonymous said...

This is actually the correct implementation to the SICP provided formula.

Unknown said...

Here's an alternative approach:
This is recognizing that if we factor out the first 2 in the numerator we are left with a series of square terms that are offset by 1 (4^2/3^2)*(6^2/5^2). One thing to be mindful of, is that this will add one too many terms in the numerator x in this case) so we need to divide by it (e.g. Notice in the book the 8 in the numerator is not squared). The last term definition just handles the cases where x is odd where we end up with the same multiplication but need to divide by one plus x.

(define (pi-prod x)
(define last-term
(if (even? x)
x
(+ x 1)))
(define (pi-term a)
(/ (sqr (+ a 1))
(sqr a)))
(define (pi-next a)
(+ a 2))
(/ (* 2 (product pi-term 3.0 pi-next x)) last-term))


Evaluating (* 4 (pi-prod 10000)) gives us 3.1417497371492993