SICP Exercise 2.87 Addendum: Making Polynomials

If you were paying close attention at the end of exercise 2.87 you would have noticed that I created polynomials for testing =zero? as follows:

(make-polynomial 'x '())
(make-polynomial 'x (list (list 4 (make-integer 3))
                          (list 2 (make-integer 1))
                          (list 0 (make-real 2.3))))
(make-polynomial 'x (list (list 3 (make-real 0))
                          (list 2 (make-rational 0 4))
                          (list 1 (make-integer 0))))

Do you notice anything about the term list?

Yep, I had to create the term list representation directly as the operations within the polynomial package for doing so are not exposed externally. This means that I, and anyone using the system, needs to have knowledge of the internal representation used for terms and term lists in order to create polynomials. They would need to know that a term is a list consisting of the order of the indeterminate that the term is for and the coefficient to apply. They would also need to know that a term list is a list of such terms, ordered in decreasing value of the order of the terms.

Exposing the internal representation in this manner has a further implication. It means that once our system is in use we will be unable to change the internal representation without having to find all clients of the system and update them to use the new representation.

Yuck.

We can address these issues by encapsulating the internal representation of terms. To do this we'll need a way for clients to create valid term representations for passing to make-polynomial. We'll also want to update make-polynomial so that clients can pass us in any list of terms and we'll reorganize them into our desired internal representation.

The first of these steps is straightforward. The polynomial package already has an internal operation, make-term for constructing polynomial terms. We just need to expose this externally in the normal manner.

As for accepting any list of polynomials and converting them into a valid representation, let's consider what we'll have to do. We'll need to go through the passed-in list of terms and build up the valid internal representation as we go. To do this we'll need to start with a valid internal representation of an empty term list (which we can get by evaluating (the-empty-termlist)). We then need to iterate through the passed-in terms and, for each one, insert it into the correct position in the term list we're building up. We should also cope with the situation where there are two or more terms passed in with the same order. We could either raise an error in this case, or add the terms together. I chose the latter.

Here's the updates required to do all this:

(define (install-polynomial-package)
  …
  (define (make-poly variable term-list)
    (cons variable (build-term-list term-list (the-empty-termlist))))
  …
  (define (insert-term term terms)
    (if (empty-termlist? terms)
        (adjoin-term term (the-empty-termlist))
        (let* ((head (first-term terms))
               (head-order (order head))
               (term-order (order term)))
          (cond ((> term-order head-order) (adjoin-term term terms))
                ((= term-order head-order)
                 (adjoin-term (make-term term-order (add (coeff term) (coeff head)))
                              (rest-terms terms)))
                (else (adjoin-term head (insert-term term (rest-terms terms))))))))
  (define (build-term-list terms result)
    (if (null? terms)
        result
        (build-term-list (cdr terms) (insert-term (car terms) result))))
  …
  (put 'make 'polynomial-term make-term)
  …
  'done)

(define (make-polynomial-term order coeff)
  ((get 'make 'polynomial-term) order coeff))

We can then make polynomials as follows:

> (make-polynomial 'x
                   (list (make-polynomial-term 4 (make-integer 3))
                         (make-polynomial-term 2 (make-integer 1))
                         (make-polynomial-term 0 (make-real 2.3))))
(polynomial x (4 (integer . 3)) (2 (integer . 1)) (0 (real . 2.3)))
> (make-polynomial 'y
                   (list (make-polynomial-term 2 (make-integer 1))
                         (make-polynomial-term 0 (make-real 2.3))
                         (make-polynomial-term 4 (make-integer 3))))
(polynomial y (4 (integer . 3)) (2 (integer . 1)) (0 (real . 2.3)))
> (make-polynomial 'z
                   (list (make-polynomial-term 0 (make-integer 0))
                         (make-polynomial-term 3 (make-real 0.0))
                         (make-polynomial-term 2 (make-rational 1 2))))
(polynomial z (2 (rational 1 . 2)))
> (make-polynomial 'a
                   (list (make-polynomial-term 2 (make-integer 1))
                         (make-polynomial-term 3 (make-real 3.5))
                         (make-polynomial-term 2 (make-rational 1 2))))
(polynomial a (3 (real . 3.5)) (2 (rational 3 . 2)))

SICP Exercise 2.87: Testing Zero for Polynomials

Install =zero? for polynomials in the generic arithmetic package. This will allow adjoin-term to work for polynomials with coefficients that are themselves polynomials.

The polynomial Package

In section 2.5.3 the authors provide us with various operations for manipulating polynomials and their terms, along with the skeleton of a polynomial package. Before we can start into the exercise we need to put all this together so that we have a package we can install into the system we've been building from section 2.5.1.

To save you the hassle of doing it yourself, here's the basic polynomial package:

(define (install-polynomial-package)
  ;; internal procedures
  ;; representation of poly
  (define (make-poly variable term-list)
    (cons variable term-list))
  (define (variable p) (car p))
  (define (term-list p) (cdr p))
  
  ;; procedures same-variable? and variable? from section 2.3.2
  (define (variable? x) (symbol? x))
  (define (same-variable? v1 v2)
    (and (variable? v1) (variable? v2) (eq? v1 v2)))

  ;; representation of terms and term lists
  (define (adjoin-term term term-list)
    (if (=zero? (coeff term))
        term-list
        (cons term term-list)))
  (define (the-empty-termlist) '())
  (define (first-term term-list) (car term-list))
  (define (rest-terms term-list) (cdr term-list))
  (define (empty-termlist? term-list) (null? term-list))
  (define (make-term order coeff) (list order coeff))
  (define (order term) (car term))
  (define (coeff term) (cadr term))

  ;; procedures used by add-poly
  (define (add-poly p1 p2)
    (if (same-variable? (variable p1) (variable p2))
        (make-poly (variable p1)
                   (add-terms (term-list p1)
                              (term-list p2)))
        (error "Polys not in same var -- ADD-POLY"
               (list p1 p2))))
  (define (add-terms L1 L2)
    (cond ((empty-termlist? L1) L2)
          ((empty-termlist? L2) L1)
          (else
           (let ((t1 (first-term L1)) (t2 (first-term L2)))
             (cond ((> (order t1) (order t2))
                    (adjoin-term
                     t1 (add-terms (rest-terms L1) L2)))
                   ((< (order t1) (order t2))
                    (adjoin-term
                     t2 (add-terms L1 (rest-terms L2))))
                   (else
                    (adjoin-term
                     (make-term (order t1)
                                (add (coeff t1) (coeff t2)))
                     (add-terms (rest-terms L1)
                                (rest-terms L2)))))))))

  ;; procedures used by mul-poly
  (define (mul-poly p1 p2)
    (if (same-variable? (variable p1) (variable p2))
        (make-poly (variable p1)
                   (mul-terms (term-list p1)
                              (term-list p2)))
        (error "Polys not in same var -- MUL-POLY"
               (list p1 p2))))
  (define (mul-terms L1 L2)
    (if (empty-termlist? L1)
        (the-empty-termlist)
        (add-terms (mul-term-by-all-terms (first-term L1) L2)
                   (mul-terms (rest-terms L1) L2))))
  (define (mul-term-by-all-terms t1 L)
    (if (empty-termlist? L)
        (the-empty-termlist)
        (let ((t2 (first-term L)))
          (adjoin-term
           (make-term (+ (order t1) (order t2))
                      (mul (coeff t1) (coeff t2)))
           (mul-term-by-all-terms t1 (rest-terms L))))))
  
  ;; interface to rest of the system
  (define (tag p) (attach-tag 'polynomial p))
  (put 'add '(polynomial polynomial) 
       (lambda (p1 p2) (tag (add-poly p1 p2))))
  (put 'mul '(polynomial polynomial) 
       (lambda (p1 p2) (tag (mul-poly p1 p2))))
  (put 'make 'polynomial
       (lambda (var terms) (tag (make-poly var terms))))
  'done)

(define (make-polynomial var terms)
  ((get 'make 'polynomial) var terms))

Don't forget to install this by evaluating (install-polynomial-package) of course, or you won't get very far!

Zero Testing for Polynomials

Now onto the exercise itself!

In order to implement =zero? for the polynomial package we need to ask ourselves how we can test whether a polynomial is zero or not. We can't determine equality with zero by looking at the indeterminate itself. However, we can determine it from the coefficients.

If there are no coefficients then the polynomial is equal to zero. On the other hand, when there are coefficients, we can only say that the polynomial is equal to zero iff all of the coefficients are equal to zero. As soon as one of them is non-zero we're unable to tell whether or not the polynomial is equal to zero until we bind the indeterminate to a particular value. Of course, as the coefficients are types from our generic arithmetic system we can test these for equality with zero by recursively evaluating =zero? on each of them.

So we have a recursive approach. If the list of coefficients is empty then the polynomial is equal to zero. Otherwise, we test whether or not the first coefficient is equal to zero. If it is then the polynomial may be equal to zero, but we need to test the remaining coefficients for equality with zero before we can confirm this; if it isn't then the polynomial isn't equal to zero.

Let's put this into code:

(define (install-polynomial-package)
  …
  (define (=zero-all-terms? L)
    (cond ((empty-termlist? L) #t)
          ((not (=zero? (coeff (first-term L)))) #f)
          (else (=zero-all-terms? (rest-terms L)))))
  (define (=zero-poly? p)
    (=zero-all-terms? (term-list p)))
  …
  (put '=zero? '(polynomial) =zero-poly?)
  …
  'done)

...and let's give it a spin:

> (=zero? (make-polynomial 'x '()))
#t
> (=zero? (make-polynomial 'x (list (list 4 (make-integer 3))
                                    (list 2 (make-integer 1))
                                    (list 0 (make-real 2.3)))))
#f
> (=zero? (make-polynomial 'x (list (list 3 (make-real 0))
                                    (list 2 (make-rational 0 4))
                                    (list 1 (make-integer 0)))))
#t

SICP Exercise 2.86: Complex Complex Numbers

Suppose we want to handle complex numbers whose real parts, imaginary parts, magnitudes, and angles can be either ordinary numbers, rational numbers, or other numbers we might wish to add to the system. Describe and implement the changes to the system needed to accommodate this. You will have to define operations such as sine and cosine that are generic over ordinary numbers and rational numbers.

Type Validation

If you recall, back in exercise 2.83 we augmented all of the procedures installed under the 'make key to validate the (Scheme built-in) types of the numeric operands passed to them to ensure that we could only construct valid representations of numbers. In the case of the complex package we restricted the real and imaginary parts, magnitudes and angles to Scheme real numbers by using real? to test each operand passed to the make-from-real-imag and make-from-mag-ang procedures defined in the rectangular and polar packages.

I'm assuming that, as we're going to allow the use of the integer, rational and real number representations provided by the integer, rational and real packages as the numeric operands for constructing complex representations, we're now going to disallow the use of the built-in Scheme real numbers. This means we need to replace the existing tests with ones that test that the operands are using numeric representations from within our tower-of-types. To ensure that we're not trying to build complex numbers with components that are themselves complex numbers we'll go one step further and ensure that the type we're using is lower in the tower than 'complex.

In the preceding exercise we devised a simple test for the first part of this restriction, which we used in apply-generic to determine whether or not we should drop the result down the tower-of-types: we check that the value is a pair and, if it is, get its type-tag and check that this is in the tower-of-types using memq. We can extract this to its own procedure:

(define (in-tower? value)
  (and (pair? value) (memq (type-tag value) tower-of-types)))

To perform the second part of the test, i.e. to check that the value's type is lower than 'complex, we can use the fact that the tower-of-types is a list ordered from lowest to highest in the tower. Recall that if the procedure memq finds the requested element in the provided list then it returns the tail of that list from that element onwards. If we use memq with the type we want to ensure is higher in the tower ('complex in this case, but let's call it T) then, provided T is present in the tower-of-types, memq will return a list consisting of T and all higher types. We can then check whether or not the type of the value we want to check is lower than T by testing to see whether it's present in the list of types that are T or higher. Iff it's not present in that list then it's a lower type.

Here's our procedure:

(define (is-lower? value type)
  (let ((type-and-higher (memq type tower-of-types)))
    (if (and type-and-higher
             (in-tower? value))
        (not (memq (type-tag value) type-and-higher))
        (error "Either value's type or type is not in tower-of-types"
               (list value type)))))

Let's see it in action:

> (is-lower? (make-integer 4) 'complex)
#t
> (is-lower? (make-rational 3 4) 'real)
#t
> (is-lower? (make-rational 3 4) 'rational)
#f
> (is-lower? (make-rational 3 4) 'integer)
#f

We can then update the rectangular and polar packages to use this test instead of real?:

(define (install-rectangular-package)
  ;; internal procedures
  …
  (define (make-from-real-imag x y)
    (if (and (is-lower? x 'complex) (is-lower? y 'complex))
        (cons x y)
        (error "non-real real or imaginary value" (list x y))))
  …
  (define (make-from-mag-ang r a) 
    (if (and (is-lower? r 'complex) (is-lower? a 'complex))
        (cons (* r (cos a)) (* r (sin a)))
        (error "non-real magnitude or angle" (list r a))))
  …
  ;; interface to the rest of the system
  …
  'done)

(define (install-polar-package)
  ;; internal procedures
  …
  (define (make-from-mag-ang r a)
    (if (and (is-lower? r 'complex) (is-lower? a 'complex))
        (cons r a)
        (error "non-real magnitude or angle" (list r a))))
  …
  (define (make-from-real-imag x y) 
    (if (and (is-lower? x 'complex) (is-lower? y 'complex))
        (cons (sqrt (add (square x) (square y)))
              (atan y x))
        (error "non-real real or imaginary value" (list x y))))
  …
  ;; interface to the rest of the system
  …
  'done)

Generic Arithmetic Operations

Next let's think about the further changes we need to make to the rectangular and polar packages, as well as the changes we'll need to make to the complex package, in order that the operations defined within these packages can work with our system's numeric representations.

The exercise says "you will have to define operations such as sine and cosine that are generic over ordinary numbers and rational numbers." We'll need to go further than just defining generic equivalents for the trigonometric functions though. We'll need to define generic equivalents for all mathematical functions used within the complex packages. Once we've done that we can then replace all of the built-in mathematical functions within the complex packages with their generic equivalents.

Looking at the complex package and the two sub-packages, rectangular and polar we can see that the built-in mathematical functions that we'll need generic equivalents for are: +, -, *, /, =, sqrt, square, atan, cos and sin. We've already got generic equivalents for the first five: add, sub, mul and divide respectively for the first four, and equ?/=zero? for = depending upon the value we're comparing.

Let's move on to the remaining functions.

We can produce the generic operations in the usual manner:

(define (square-root x) (apply-generic 'sqrt x))
(define (square x) (apply-generic 'square x))
(define (arctan x y) (apply-generic 'arctan x y))
(define (cosine x) (apply-generic 'cosine x))
(define (sine x) (apply-generic 'sine x))

All that remains is to install appropriate implementations of these in each of the packages lower in the tower-of-types than 'complex.

The procedures we need to install for 'square are straightforward. For 'integer and 'real values we can just use the built-in * operation and multiply the underlying numeric value by itself. For 'rational numbers we can simply use mul-rat to multiply the rational number by itself. The result can then be tagged as normal.

As for the square root and trigonometric functions, while we could calculate these using the techniques we learned back in section 1.3.3 and its associated exercises (e.g. such as exercise 1.39), we'll go easy on ourselves and just delegate to the corresponding built-in operations appropriately for each type. However, we do need to treat the results of these operations carefully:

• Given Scheme integer, rational or real values, the results of the built-in trigonometric operations may themselves be Scheme real numbers so, apart from in the real package, we can't just tag the results of the delegated operation normally. Instead we'll have to create real representations from them.
• Worse still, the result of the built-in sqrt operation may be a Scheme complex number. In this case we'll need to extract out the real and imaginary parts of the result and use these to create a complex representation... And as we're updating the system so that complex numbers can only be constructed using our system's own number representations, we'll have to turn the real and imaginary parts into real representations in order for us to construct the complex representation.

Once we've constructed the appropriate representation for the results we can then rely upon the work we did in exercise 2.85 to drop the number's type to the simplest representation available.

It's possibly worth noting that the complexity of dealing with the results of sqrt points to another piece of work we need to do. As complex representations can only be created using representations from within the tower-of-types we'll need to update the coercion procedure real->complex to that it too turns the components into real representations.

Given all that, here are the changes required to the integer, rational and real packages in order to install the procedures and support the complex package restrictions:

(define (install-integer-package)
  …
  (put 'sqrt '(integer)
       (lambda (x)
         (let ((root (sqrt x)))
           (make-complex-from-real-imag (make-real (real-part root))
                                        (make-real (imag-part root))))))
  (put 'square '(integer)
       (lambda (x) (tag (* x x))))
  (put 'arctan '(integer integer)
       (lambda (x y) (make-real (atan x y))))
  (put 'cosine '(integer)
       (lambda (x) (make-real (cos x))))
  (put 'sine '(integer)
       (lambda (x) (make-real (sin x))))
  …
  'done)

(define (install-rational-package)
  ;; internal procedures
  …
  (define (sqrt-rat x)
    (let ((root (sqrt (/ (numer x) (denom x)))))
      (make-complex-from-real-imag (make-real (real-part root))
                                   (make-real (imag-part root)))))
  (define (square-rat x)
    (mul-rat x x))
  (define (arctan-rat x y)
    (atan (/ (numer x) (denom x))
          (/ (numer y) (denom y))))
  (define (cosine-rat x)
    (cos (/ (numer x) (denom x))))
  (define (sine-rat x)
    (sin (/ (numer x) (denom x))))
  …
  ;; interface to rest of the system
  …
  (put 'sqrt '(rational)
       (lambda (x) (make-real (sqrt-rat x))))
  (put 'square '(rational)
       (lambda (x) (tag (square-rat x))))
  (put 'arctan '(rational rational)
       (lambda (x y) (make-real (arctan-rat x y))))
  (put 'cosine '(rational)
       (lambda (x) (make-real (cosine-rat x))))
  (put 'sine '(rational)
       (lambda (x) (make-real (sine-rat x))))
  …
  'done)

(define (install-real-package)
  …
  (define (real->complex r) (make-complex-from-real-imag (tag r) (tag 0)))
  …
  (put 'sqrt '(real)
       (lambda (x)
         (let ((root (sqrt x)))
           (make-complex-from-real-imag (tag (real-part root))
                                        (tag (imag-part root))))))
  (put 'square '(real)
       (lambda (x) (tag (* x x))))
  (put 'arctan '(real real)
       (lambda (x y) (tag (atan x y))))
  (put 'cosine '(real)
       (lambda (x) (tag (cos x))))
  (put 'sine '(real)
       (lambda (x) (tag (sin x))))
  …
  'done)

Before we go any further, let's give this a spin:

> (square-root (make-real 16.0))
(integer . 4)
> (square (make-rational 3 5))
(rational 9 . 25)
> (cosine (make-integer 0))
(integer . 1)
> (cos 0)
1.0
> (sine (make-integer 3))
(rational 5084384125703515 . 36028797018963968)
> (/ 5084384125703515.0 36028797018963968.0)
0.141120008059867
> (sin 3)
0.141120008059867
> (arctan (make-integer 3) (make-integer 4))
(rational 5796142707547873 . 9007199254740992)
> (/ 5796142707547873.0 9007199254740992.0)
0.643501108793284
> (atan 3 4)
0.643501108793284

Updating the complex Packages

All that remains is to go through the rectangular, polar and complex packages and replace all the built-in operations with our generic equivalents.

Note that the changes we'll make to make-from-mag-ang in the rectangular package and to make-from-real-imag in the polar package will have the effect of dropping the calculated components of the created complex representation. To keep everything consistent, let's explicitly drop the components passed to make-from-real-imag in the rectangular package and to make-from-mag-ang in the polar package.

Here are the changes to the packages:

(define (install-rectangular-package)
  ;; internal procedures
  …
  (define (make-from-real-imag x y)
    (if (and (is-lower? x 'complex) (is-lower? y 'complex))
        (cons (drop x) (drop y))
        (error "non-real real or imaginary value" (list x y))))
  (define (magnitude z)
    (square-root (add (square (real-part z))
                 (square (imag-part z)))))
  (define (angle z)
    (arctan (imag-part z) (real-part z)))
  (define (make-from-mag-ang r a) 
    (if (and (is-lower? r 'complex) (is-lower? a 'complex))
        (cons (mul r (cosine a)) (mul r (sine a)))
        (error "non-real magnitude or angle" (list r a))))
  …
  ;; interface to the rest of the system
  …
  'done)

(define (install-polar-package)
  ;; internal procedures
  …
  (define (make-from-mag-ang r a)
    (if (and (is-lower? r 'complex) (is-lower? a 'complex))
        (cons (drop r) (drop a))
        (error "non-real magnitude or angle" (list r a))))
  (define (real-part z)
    (mul (magnitude z) (cosine (angle z))))
  (define (imag-part z)
    (mul (magnitude z) (sine (angle z))))
  (define (make-from-real-imag x y) 
    (if (and (is-lower? x 'complex) (is-lower? y 'complex))
        (cons (square-root (add (square x) (square y)))
              (arctan y x))
        (error "non-real real or imaginary value" (list x y))))
  …
  ;; interface to the rest of the system
  …
  'done)

(define (install-complex-package)
  ;; imported procedures from rectangular and polar packages
  …
  ;; internal procedures
  …
  (define (add-complex z1 z2)
    (make-from-real-imag (add (complex-real-part z1) (complex-real-part z2))
                         (add (complex-imag-part z1) (complex-imag-part z2))))
  (define (sub-complex z1 z2)
    (make-from-real-imag (sub (complex-real-part z1) (complex-real-part z2))
                         (sub (complex-imag-part z1) (complex-imag-part z2))))
  (define (mul-complex z1 z2)
    (make-from-mag-ang (mul (complex-magnitude z1) (complex-magnitude z2))
                       (add (complex-angle z1) (complex-angle z2))))
  (define (div-complex z1 z2)
    (make-from-mag-ang (div (complex-magnitude z1) (complex-magnitude z2))
                       (sub (complex-angle z1) (complex-angle z2))))
  (define (equ-complex? z1 z2)
    (and (equ? (complex-real-part z1) (complex-real-part z2))
         (equ? (complex-imag-part z1) (complex-imag-part z2))))
  (define (=zero-complex? x) (=zero? (complex-magnitude x)))
  (define (addd-complex z1 z2 z3)
    (make-from-real-imag (addd (complex-real-part z1)
                               (complex-real-part z2)
                               (complex-real-part z3))
                         (addd (complex-imag-part z1)
                               (complex-imag-part z2)
                               (complex-imag-part z3))))
  …
  ;; interface to rest of the system
  …
  'done)

I've included the addd procedure we added in exercise 2.82 as well for completeness.

Okay, let's see this all in action:

> (square-root (make-integer -1))
(complex rectangular (integer . 0) integer . 1)
> (make-complex-from-real-imag (make-real 3.0) (make-rational 1 2))
(complex rectangular (integer . 3) rational 1 . 2)
> (add (make-complex-from-real-imag (make-integer 4) (make-rational 2 4))
       (make-integer 4))
(complex rectangular (integer . 8) rational 1 . 2)
> (sub (make-complex-from-real-imag (make-real 1.5) (make-integer 3))
       (make-complex-from-real-imag (make-rational 1 2) (make-integer 3)))
(integer . 1)