SICP Exercise 1.29: Approximating Integrals with Simpson's Rule

Simpson's Rule is a more accurate method of numerical integration than the method illustrated above. Using Simpson's Rule, the integral of a function f between a and b is approximated as

  h
  - [y0 + 4y1 + 2y2 + 4y3 + 2y4 + ... + 2yn-2 + 4yn-1 + yn]
  3

where h = (b - a)/n, for some even integer n, and y_k = f(a + kh). (Increasing n increases the accuracy of the approximation.) Define a procedure that takes as arguments f, a, b, and n and returns the value of the integral, computed using Simpson's Rule. Use your procedure to integrate cube between 0 and 1 (with n = 100 and n = 1000), and compare the results to those of the integral procedure shown above.

Okay, so we want a procedure that takes parameters f, a, b and n and uses sum to perform the sum portion of the equation above. So it's going to have the general form:

(define (simpsons-rule f a b n)
  <definitions of h and kth-term go here>
  (* (sum kth-term 0 inc n)
     (/ h 3.0)))

Producing a definition of h is straightforward. This is a constant value in the context of a particular set of values for a, b and n, and so can be defined simply as:

(define h (/ (- b a) n))

We can also produce a straightforward definition of y_k:

(define (y k) (f (+ a (* k h))))

However to produce the sum we need to know what to multiply y_k by! We need a procedure that multiplies the results of y_k by 1, 2 or 4 depending upon whether k is 0 or n, is even, or is odd respectively:

(define (kth-multiple k)
  (cond ((or (= k 0) (= k n)) 1.0)
        ((even? k) 2.0)
        (else 4.0)))

Using these we can produce kth-term as:

(define (kth-term k)
  (* (kth-multiple k)
     (kth-term k)))

And finally we can put all these together to produce:

(define (simpsons-rule f a b n)
  (define h (/ (- b a) n))
  (define (y k) (f (+ a (* k h))))
  (define (kth-multiple k)
    (cond ((or (= k 0) (= k n)) 1.0)
          ((even? k) 2.0)
          (else 4.0)))
  (define (kth-term k)
    (* (kth-multiple k) (y k)))
  (* (sum kth-term 0 inc n)
     (/ h 3.0)))

Now we could compress this a bit more. We can combine y, kth-multiple and kth-term to produce:

(define (simpsons-rule f a b n)
  (define h (/ (- b a) n))
  (define (kth-term k)
     (* (f (+ a (* k h)))
        (cond ((or (= k 0) (= k n)) 1.0)
              ((even? k) 2.0)
              (else 4.0))))
  (* (sum kth-term 0 inc n)
     (/ h 3.0)))

This is slightly more compact and perhaps a bit more readable.

Right, let's compare Simpson's rule with integral:

> (integral cube 0 1 0.01)
0.24998750000000042
> (simpsons-rule cube 0 1 100)
0.24999999999999992
> (integral cube 0 1 0.001)
0.249999875000001
> (simpsons-rule cube 0 1 1000)
0.2500000000000002

As noted right at the end of section 1.3.1, the exact value of the integral of cube between 0 and 1 is 1/4. As we can see, Simpson's rule provides a more accurate approximation than integral and also increasing n does indeed improve the accuracy of the approximation.

A Quick Note on Progress

As of tomorrow it'll be two weeks since I made my first post stating my intention to write up my solutions to the SICP exercises. In those two weeks I've written up the first 28 exercises, which means I'm averaging 2 per day. Of course some of them are quick, taking 10 minutes tops to write up. Some of them, however, have taken considerably longer and involved significant head-scratching, looking up details of the interpreter I'm using and other niceties.

Now at the moment I'm playing catch-up... Some of my colleagues and I started working our way through SICP nearly 11 weeks ago and, barring one week off for good behaviour, we've been slowly working through the book in a weekly reading group. As of this coming week's meeting, we'll be up to the start of section 2.3. That's a grand total of 97 exercises we've passed on the way, so I'm about a quarter of the way through writing up - but of course it's a moving target!

I should say that we're not dictating that everyone must do every exercise, or write a book report each week. When we get together we chat about any interesting topics raised in the section we've decided to cover or any problems we had with the exercises. I've been a good boy so far and managed to produce a solution to each exercise every week. Something that this blog has reitterated to me though is that writing up work for public consumption is different from doing it for your own benefit.

For example, when I originally tackled exercise 1.14 as part of our weekly reading group homework, I was happy to state that obviously the algorithm had order of growth in steps of Θ(n^m) . However, I realised that in order to turn this into a blog post I was going to have to show why that's the case. I'm still not overly happy with my answer to this particular exercise (I feel like my assertions at the end to get from Θ(n²) to Θ(n^m) are almost cheating), but it's much better than simply stating the growth!

Anyway, if I can carry on at this rate then I should have caught up by Hallowe'en. Of course that's always assuming I can maintain this initial burst of enthusiasm for so long... and that my family continue to let me steal away for an hour or so every so often to do exercises and write-ups! I'm also hoping that at somepoint along the way I'll feel that my technical communication skills have been given a bit of a boost. Which, to be honest, is one of the main reasons for blogging my solutions.

Right. Back to it now. Got another few exercises to finish off before this week's meeting. If I'm really lucky I might also manage to write up another few exercises later on this evening.

SICP Exercise 1.28: The Miller-Rabin Test

One variant of the Fermat test that cannot be fooled is called the Miller-Rabin test (Miller 1976; Rabin 1980). This starts from an alternate form of Fermat's Little Theorem, which states that if n is a prime number and a is any positive integer less than n, then a raised to the (n - 1)^st power is congruent to 1 modulo n. To test the primality of a number n by the Miller-Rabin test, we pick a random number a<n and raise a to the (n - 1)^st power modulo n using the expmod procedure. However, whenever we perform the squaring step in expmod, we check to see if we have discovered a "nontrivial square root of 1 modulo n," that is, a number not equal to 1 or n - 1 whose square is equal to 1 modulo n. It is possible to prove that if such a nontrivial square root of 1 exists, then n is not prime. It is also possible to prove that if n is an odd number that is not prime, then, for at least half the numbers a<n, computing a^n-1 in this way will reveal a nontrivial square root of 1 modulo n. (This is why the Miller-Rabin test cannot be fooled.) Modify the expmod procedure to signal if it discovers a nontrivial square root of 1, and use this to implement the Miller-Rabin test with a procedure analogous to fermat-test. Check your procedure by testing various known primes and non-primes. Hint: One convenient way to make expmod signal is to have it return 0.

Let's start with figuring out how to test whether or not a number is a non-trivial square root of 1 modulo n. We want to produce a predicate that returns true if, and only if, the number passed to it is not 1 or n - 1 and, when you square it and modulo the result by n you get 1. This is pretty straightforward:

(define (non-trivial-sqrt? x n)
  (and (not (= x 1))
       (not (= x (- n 1)))
       (= (remainder (square x) n) 1)))

Now we're given the hint that "one convenient way to make expmod signal [if it's found a non-trival square root of 1 modulo n] is to have [expmod] return 0." How can we do that? First note that we want to apply the non-trivial-sqrt? test to the result of the recursive call, so we need to put our test in between the recursive call and the subsequent squaring of the result of the recursive call and application of modulo m to that. What happens when you square 0? You get 0. What happens when you calculate 0 modulo m? You get 0 for any m. So if the non-trivial-sqrt? test fails then we can return 0 instead of the value of the recursive call, and this value will then propagate back up the chain of calls and become the final result of the expmod call.

To do this we introduce a procedure that performs the test and returns 0 if it's found a non-trivial square root of 1 modulo n, and just returns the number tested otherwise:

(define (zero-if-non-trivial-sqrt x n)
  (if (non-trivial-sqrt? x n) 0 x))

Then we can modify expmod, simply wrapping the recursive call in a call to this. I renamed this version of expmod to mrt-expmod, as it no longer performs a true calculation of base^exp mod m:

(define (mrt-expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp) (remainder
                       (square
                         (zero-if-non-trivial-sqrt
                           (mrt-expmod base (/ exp 2) m) m)) m))
        (else (remainder (* base (mrt-expmod base (- exp 1) m)) m))))

Finally, we can wrap all this up in a test:

(define (miller-rabin-test n)
  (define (try-it a)
    (= (mrt-expmod a (- n 1) n) 1))
  (try-it (+ 1 (random (- n 1)))))

And then try it out:

> (miller-rabin-test 561)
#f
> (miller-rabin-test 1105)
#f
> (miller-rabin-test 1729)
#f
> (miller-rabin-test 2465)
#f
> (miller-rabin-test 2821)
#f
> (miller-rabin-test 6601)
#f

Of course, this is just testing a random value between 1 and n. To be truly sure that the Miller-Rabin test is not defeated by the Charmichael numbers we tested in exercise 1.27, let's check all of the numbers in the range 1 to n using:

(define (miller-rabin-test-all n)
  (define (try-it a)
    (cond ((= a n) #t)
          ((= (mrt-expmod a (- n 1) n) 1) (try-it (+ a 1)))
          (else #f)))
  (try-it 1))

And here are the results of our testing:

> (miller-rabin-test-all 561)
#f
> (miller-rabin-test-all 1105)
#f
> (miller-rabin-test-all 1729)
#f
> (miller-rabin-test-all 2465)
#f
> (miller-rabin-test-all 2821)
#f
> (miller-rabin-test-all 6601)
#f

And just for the sake of sanity:

> (miller-rabin-test-all 1)
#t
> (miller-rabin-test-all 2)
#t
> (miller-rabin-test-all 3)
#t
> (miller-rabin-test-all 4)
#f
> (miller-rabin-test-all 5)
#t
> (miller-rabin-test-all 6)
#f
> (miller-rabin-test-all 7)
#t
> (miller-rabin-test-all 8)
#f
> (miller-rabin-test-all 9)
#f