A Temporary Interruption to Service

Honest, I'm working on SICP exercise 2.97. In fact I've finished it and I'm in the process of writing it up. However, I'm in Seattle for a couple of weeks of work, which is keeping me hectically busy. I've also just started Andrew Ng's Machine Learning course on Coursera. Incidentally, he gave a keynote speech at a conference I was attending the other day. A real smart guy! Add to that my work being busy as ever and the minor additional factor of my beautiful almost-5-month old daughter and I've got no free time at all at the moment!

Normal service will be resumed shortly!

SICP Exercise 2.96: Pseudo-Remainder Terms

1. Implement the procedure pseudoremainder-terms, which is just like remainder-terms except that it multiplies the dividend by the integerizing factor described above before calling div-terms. Modify gcd-terms to use pseudoremainder-terms, and verify that greatest-common-divisor now produces an answer with integer coefficients on the example in exercise 2.95.
2. The GCD now has integer coefficients, but they are larger than those of P₁. Modify gcd-terms so that it removes common factors from the coefficients of the answer by dividing all the coefficients by their (integer) greatest common divisor.

Part A: Introducing the Integerizing Factor

The calculation of the integerizing factor and the implementation of pseudoremainder-terms are described in the book as follows:

[…] if P and Q are polynomials, let O₁ be the order of P (i.e., the order of the largest term of P) and let O₂ be the order of Q. Let c be the leading coefficient of Q. Then it can be shown that, if we multiply P by the integerizing factor c^{1 + O₁ - O₂}, the resulting polynomial can be divided by Q by using the div-terms algorithm without introducing any fractions. The operation of multiplying the dividend by this constant and then dividing is sometimes called the pseudodivision of P by Q. The remainder of the division is called the pseudoremainder.

Note that, in order to calculated the integerizing factor, we need to be able to calculate exponents. Note also that, within our system, the c in c^{1 + O₁ - O₂} is a tagged type. This means that we'll need to add support for calculating exponentiation to our system. On the other hand, O₁ and O₂ are primitive integer values, so we'll need to convert the results of 1 + O₁ - O₂ to a tagged integer value so that it can be used as the exponent.

Let's begin with implementing support for exponentiation.

The book states that the modified GCD algorithm "really works only in the case of polynomials with integer coefficients". This means that our coefficient will always be an integer. As a result we can be pragmatic here and define the exponent operation for the integer package only. While defining it for the rational and real packages would be straightforward, this allows us to avoid having to define a generic logarithm operation in order to support complex numbers, as well as avoiding the complexities of calculating exponents for polynomial types.

We could define the integer implementation using the fast-expt implementation given in section 1.2.4 of the book. However, I chose simply to delegate the calculation to Scheme's built-in expt procedure:

(put 'exponent '(integer integer)
  (lambda (x y) (tag (expt x y))))

With that in place we can define a generic exponent operation in the usual manner:

(define (exponent b e)
  (apply-generic 'exponent b e))

We can now implement pseudoremainder-terms. This will still call through to div-terms and return the cadr of the results as before. However, we need to apply the integerizing factor to the first term-list and then apply div-terms to divide the resulting term-list by second term-list. As noted above, the result of the calculation 1 + O₁ - O₂ needs converted to a tagged integer type using make-integer in order that it can be used as the exponent when applying the generic exponent operation. In the implementation below we then convert this integerizing factor to a zero-order term in order that we can use our existing mul-term-by-all-terms procedure to apply the factor to the first term-list. Here's the code:

(define (pseudoremainder-terms L1 L2)
  (let* ((t1 (first-term L1))
         (t2 (first-term L2))
         (factor (exponent (coeff t2)
                           (make-integer (+ 1 (order t1) (- (order t2))))))
         (pseudoL1 (mul-term-by-all-terms (make-term 0 factor) L1)))
    (cadr (div-terms pseudoL1 L2))))

All that remains in part A is to modify gcd-terms to use pseudoremainder-terms instead of remainder-terms and we're done:

 (define (gcd-terms L1 L2)
  (if (empty-termlist? L2)
      L1
      (gcd-terms L2 (pseudoremainder-terms L1 L2))))

Let's see this in action. If you recall, in exercise 2.95 we defined two polynoials, q1 and q2, and calculating their GCD gave the following result:

> (greatest-common-divisor q1 q2)
(polynomial x dense-terms (rational (integer . 18954) integer . 2197)
                          (rational (integer . -37908) integer . 2197)
                          (rational (integer . 1458) integer . 169))

When we repeat this calculation with pseudoremainder-terms we now get the following result:

> (greatest-common-divisor q1 q2)
(polynomial x dense-terms (integer . 1458) (integer . -2916) (integer . 1458))

As you can see, just like in exercise 2.95, the coefficients hold the same relative ratios as the coefficients of P₁, i.e. 1:-2:1. However, the resulting polynomial now has (large) integer coefficients throughout.

Part B: Removing Common Factors

We now want to remove the common factors from the coefficients of the result of gcd-terms so that we simplify the calculated GCD as much as possible. If we assume we have a procedure, remove-common-factors, which does this then our gcd-terms implementation becomes:

(define (gcd-terms L1 L2)
  (if (empty-termlist? L2)
      (remove-common-factors L1)
      (gcd-terms L2 (pseudoremainder-terms L1 L2))))

So what does remove-common-factors need to do? Well firstly it will need to determine the highest common factor (a.k.a. the GCD) of all of the coefficients and then it will need produce a new term-list by dividing each of the coefficients in the term-list by this common factor.

In order to find the GCD of the coefficients, first observe that the GCD of any two coefficients must itself have the GCD of all of the coefficients as a factor. So in order to determine the GCD of all of the coefficients we can simply start with the first coefficient as our initial GCD candidate and then repeatedly find the GCD of our current candidate and the next coefficient until we've processed the entire term-list. Programmatically this can be expressed as:

(define (find-gcd c L)
  (if (empty-termlist? L)
      c
      (find-gcd (greatest-common-divisor c (coeff (first-term L)))
                (rest-terms L))))

Obviously we can't use this to find the GCD of an entirely empty term-list, as there will be no first term to start with. However, this is just a helper procedure to support the implementation of remove-common-factors. As an empty term-list with its common factors removed is just an empty term-list we won't even need to invoke find-gcd in such a case. This provides a basic skeleton for our implementation of remove-common-factors: we determine whether we're dealing with an empty term-list, returning an empty term-list if we are and actually performing the removal if we're not. Here's the skeleton:

(define (remove-common-factors L)
  (if (empty-termlist? L)
      L
      <RemoveCommonFactors>))

All that remains is to use find-gcd to determine the GCD of the coefficients and then divide all of the coefficients by this GCD. One way of achieving this division would be to write a procedure that iterates through the term-list, applying the division to each coefficient and building a new term-list with the results. However, we already have a method of dividing a term-list by another value: div-terms, so it makes sense to use that instead. Unfortunately div-terms requires a term-list as the divisor, and find-gcd produces a tagged integer value, not a term-list. This means we'll need to first convert our tagged integer value into a zero-order term-list. We can do this by creating an empty termlist (with the-empty-termlist), creating a zero-order term with the GCD of the coefficients as its coefficient, and then adding this term to the empty term-list (with adjoin-term). Note that the division is guaranteed to produce no remainder, so we can simply return the calculated quotient as the simplified GCD.

Here's the full implementation of remove-common-factors:

  (define (remove-common-factors L)
    (if (empty-termlist? L)
        L
        (let* ((gcd-coeff (find-gcd (coeff (first-term L)) (rest-terms L)))
               (divisor (adjoin-term (make-term 0 gcd-coeff) (the-empty-termlist))))
          (car (div-terms L divisor)))))

With that in place we can now retry our calculation of the GCD of q1 and q2 again:

> (greatest-common-divisor q1 q2)
(polynomial x dense-terms (integer . 1) (integer . -2) (integer . 1))

This time we get the expected result. We can now move on to the last exercise in the chapter, where we'll use gcd-poly to reduce rational functions to their simplest terms.

SICP Exercise 2.95: GCD and Non-Integer Values

Define P₁, P₂, and P₃ to be the polynomials

P1: x2 - 2x + 1
P2: 11x2 + 7
P3: 13x + 5

Now define Q₁ to be the product of P₁ and P₂ and Q₂ to be the product of P₁ and P₃, and use greatest-common-divisor (exercise 2.94) to compute the GCD of Q₁ and Q₂. Note that the answer is not the same as P₁. This example introduces noninteger operations into the computation, causing difficulties with the GCD algorithm. To understand what is happening, try tracing gcd-terms while computing the GCD or try performing the division by hand.

First let's define P₁, P₂ and P₃ as p1, p2 and p3 respectively:

> (define p1 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 1)
                                                (make-integer -2)
                                                (make-integer 1))))
> (define p2 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 11)
                                                zero
                                                (make-integer 7))))
> (define p3 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 13)
                                                (make-integer 5))))

We then use these to define Q₁ and Q₂ as q1 and q2 respectively:

> (define q1 (mul p1 p2))
> (define q2 (mul p1 p3))
> q1
(polynomial x dense-terms (integer . 11)
                          (integer . -22)
                          (integer . 18)
                          (integer . -14)
                          (integer . 7))
> q2
(polynomial x dense-terms (integer . 13)
                          (integer . -21)
                          (integer . 3)
                          (integer . 5))

Finally we calculate the GCD of Q₁ and Q₂:

> (greatest-common-divisor q1 q2)
(polynomial x dense-terms (rational (integer . 18954) integer . 2197)
                          (rational (integer . -37908) integer . 2197)
                          (rational (integer . 1458) integer . 169))

As the authors note in the exercise statement, the answer is not the same as P₁, it does introduce non-integer operations into the computation and these then cause difficulties with the GCD algorithm. However, if you look closely at the numbers you might notice a correlation with P₁. Note that 2197 / 13 = 169. We can use this to first express all the coefficients as multiples of 1/169:

  (18954/2197)x2 - (37908/2197)x + 1458/169
= ((18954/13)/(2197/13))x2 - ((37908/13)/(2197/13))x + 1458/169
= (1458/169)x2 - (2916/169)x + 1458/169

Now let's normalize this by multiplying all the coefficients by 169/1458:

  ((1458/169)x2 - (2916/169)x + 1458/169) * (169/1458)
= ((1458*169)/(169*1458))x2 - ((2916*169)/(169*1458))x + (1458*169)/(169*1458)
= ((1458*169)/(169*1458))x2 - ((2916*169)/(169*1458))x + (1458*169)/(169*1458)
= (246402/246402)x2 - (492804/246402)x + 246402/246402
= x2 - 2x + 1
= P1

As you can see, the calculated GCD's term coefficients have the same relative ratios as the coefficients of P₁, i.e. 1:-2:1.

Let's trace the calculation. The Scheme interpreter I'm using has a trace procedure, which takes a procedure as an argument and turns on tracing of all calls to that passed procedure. However, I found this to be limited in use: when a trace output line is "too long" it arbitrarily collapses each argument passed to the procedure to "#". This is obviously no use for what we'd like to achieve here. As a result, I quickly hacked a trace directly into the procedures gcd-terms and div-terms so that I could get a better picture of what's going on.

With this in place my traced call produced the following output, with the point where non-integer values enter the calculation highlighted in yellow:

> (greatest-common-divisor q1 q2)
gcd-terms:
+-> (dense-terms (integer . 11)
|                (integer . -22)
|                (integer . 18)
|                (integer . -14)
|                (integer . 7))
+-> (dense-terms (integer . 13)
|                (integer . -21)
|                (integer . 3)
|                (integer . 5))
| div-terms:
| +-> (dense-terms (integer . 11)
| |                (integer . -22)
| |                (integer . 18)
| |                (integer . -14)
| |                (integer . 7))
| +-> (dense-terms (integer . 13)
| |                (integer . -21)
| |                (integer . 3)
| |                (integer . 5))
| | div-terms:
| | +-> (dense-terms (rational (integer . -55) integer . 13)
| | |                (rational (integer . 201) integer . 13)
| | |                (rational (integer . -237) integer . 13)
| | |                (integer . 7))
| | +-> (dense-terms (integer . 13)
| | |                (integer . -21)
| | |                (integer . 3)
| | |                (integer . 5))
| | | div-terms:
| | | +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | | |                 (term 1 (rational (integer . -37908) integer . 2197))
| | | |                 (term 0 (rational (integer . 1458) integer . 169)))
| | | +-> (dense-terms (integer . 13)
| | | |                (integer . -21)
| | | |                (integer . 3)
| | | |                (integer . 5))
| | | +-= ((sparse-terms)
| | |      (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | |                    (term 1 (rational (integer . -37908) integer . 2197))
| | |                    (term 0 (rational (integer . 1458) integer . 169))))
| | +-= ((sparse-terms (term 0 (rational (integer . -55) integer . 169)))
| |      (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| |                    (term 1 (rational (integer . -37908) integer . 2197))
| |                    (term 0 (rational (integer . 1458) integer . 169))))
| +-= ((sparse-terms (term 1 (rational (integer . 11) integer . 13))
|                    (term 0 (rational (integer . -55) integer . 169)))
|      (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
|                    (term 1 (rational (integer . -37908) integer . 2197))
|                    (term 0 (rational (integer . 1458) integer . 169))))
| gcd-terms:
| +-> (dense-terms (integer . 13)
| |                (integer . -21)
| |                (integer . 3)
| |                (integer . 5))
| +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| |                 (term 1 (rational (integer . -37908) integer . 2197))
| |                 (term 0 (rational (integer . 1458) integer . 169)))
| | div-terms:
| | +-> (dense-terms (integer . 13)
| | |                (integer . -21)
| | |                (integer . 3)
| | |                (integer . 5))
| | +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | |                 (term 1 (rational (integer . -37908) integer . 2197))
| | |                 (term 0 (rational (integer . 1458) integer . 169)))
| | | div-terms:
| | | +-> (dense-terms (rational (integer . 208209690) integer . 41641938)
| | | |                (rational (integer . -32032260) integer . 3203226)
| | | |                (integer . 5))
| | | +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | | |                 (term 1 (rational (integer . -37908) integer . 2197))
| | | |                 (term 0 (rational (integer . 1458) integer . 169)))
| | | | div-terms:
| | | | +-> (sparse-terms)
| | | | +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | | | |                 (term 1 (rational (integer . -37908) integer . 2197))
| | | | |                 (term 0 (rational (integer . 1458) integer . 169)))
| | | | +-= ((sparse-terms)
| | | |      (sparse-terms))
| | | +-= ((sparse-terms (term 0 (rational (integer . 457436688930)
| | |                                       integer . 789281292852)))
| | |      (sparse-terms))
| | +-= ((sparse-terms (term 1 (rational (integer . 28561) integer . 18954))
| |                    (term 0 (rational (integer . 457436688930)
| |                                       integer . 789281292852)))
| |      (sparse-terms))
| | gcd-terms:
| | +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | |                 (term 1 (rational (integer . -37908) integer . 2197))
| | |                 (term 0 (rational (integer . 1458) integer . 169)))
| | +-> (sparse-terms)
| | +-= (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| |                   (term 1 (rational (integer . -37908) integer . 2197))
| |                   (term 0 (rational (integer . 1458) integer . 169)))
| +-= (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
|                   (term 1 (rational (integer . -37908) integer . 2197))
|                   (term 0 (rational (integer . 1458) integer . 169)))
+-= (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
                  (term 1 (rational (integer . -37908) integer . 2197))
                  (term 0 (rational (integer . 1458) integer . 169)))
(polynomial x dense-terms (rational (integer . 18954) integer . 2197)
                          (rational (integer . -37908) integer . 2197)
                          (rational (integer . 1458) integer . 169))

As you can see the non-integer values enter the calculation when the first invocation of div-terms by gcd-terms (via remainder-terms) makes a recursive call to itself. The fraction first enters the calculations when div-terms calculates the coefficient to use for the first term of the result. This is calculated by dividing the coefficient of the highest-order term from the first term-list by the coefficient of the highest-order term from the second term-list. Provided the latter is a factor of the former we will avoid non-integer values entering the computation. However in this case the coefficients are 11 and 13 respectively which, being prime numbers, are relatively prime to each other and so share no factors. Not only does this mean the coefficient of the first term of the result will be non-integer, but it also affects the calculations of both the coefficients for the other terms in the result and the remainder.

You may also have noticed in the trace above that the denominators and numerators of the rational numbers that appear in our calculations get exponentially larger as the algorithm progresses. This is a known limitation of Euclidean division, which is the algorithm we're using to implement div-terms. In the next exercise we'll address this by first implementing a procedure which calculates pseudo-remainders, which guarantees that no non-integer coefficents will result from the calculation, and then incorporating this into our gcd-terms implementation.

SICP Exercise 2.94: A Generic GCD

Using div-terms, implement the procedure remainder-terms and use this to define gcd-terms as above. Now write a procedure gcd-poly that computes the polynomial GCD of two polys. (The procedure should signal an error if the two polys are not in the same variable.) Install in the system a generic operation greatest-common-divisor that reduces to gcd-poly for polynomials and to ordinary gcd for ordinary numbers. As a test, try

(define p1 (make-polynomial 'x '((4 1) (3 -1) (2 -2) (1 2))))
(define p2 (make-polynomial 'x '((3 1) (1 -1))))
(greatest-common-divisor p1 p2)

and check your result by hand.

Calculating the GCD of Polynomials

The procedure div-terms was introduced back in exercise 2.91 and returns a list of the quotient term-list and the remainder term-list. As a result the implementation of remainder-terms is very straightforward. It just needs to extract and return the second item from the list returned by div-terms:

(define (remainder-terms L1 L2)
  (cadr (div-terms L1 L2)))

This can then be used by the implementation of gcd-terms supplied in the book (I've renamed the parameters a and b to L1 and L2 respectively to match the convention used in our rational package):

(define (gcd-terms L1 L2)
  (if (empty-termlist? L2)
      L1
      (gcd-terms L2 (remainder-terms L1 L2))))

We then need to define a gcd-poly for polynomials which raises an error if the polynomials have different indeterminates. This requires that we first check that the indeterminates match, raising the error if they don't, invoking gcd-terms to determine the GCD of the two polynomials' term-lists, and then finally creating a new polynomial with the same indeterminate as the two original polynomials and with the GCD as its term-list.

Note that, with our implementation of the rational package, we can't simply use eq? to compare the indeterminates and then take the indeterminate from the first polynomial as the one to create the resulting polynomial in. Back in exercise 2.90 we introduced the concept of the unbound indeterminate in order to allow us to handle type conversions correctly. This requires that an unbound polynomial is considered to match the indeterminate of any other polynomial. In order to achieve this we also introduced the procedures same-variable? and select-variable, which cope with the unbound indeterminate. As a result, we must implement gcd-poly in terms of these procedures. Other than this the implementation is straightforward:

(define (gcd-poly p1 p2)
  (let ((v1 (variable p1))
        (v2 (variable p2)))
    (if (same-variable? v1 v2)
        (make-poly (select-variable v1 v2)
                   (gcd-terms (term-list p1) (term-list p2)))
        (error (list "Polynomials must have same indeterminate to calculate GCD"
                     p1
                     p2)))))

We can install this in the polynomial package in the usual manner:

(put 'gcd '(polynomial polynomial)
     (lambda (p1 p2) (tag (gcd-poly p1 p2))))

Calculating the GCD of Ordinary Numbers

The exercise statement indicates that our generic greatest-common-divisor operation should reduce to ordinary gcd for "ordinary numbers". Of course we have slowly been building up our type system over the last several exercises. As a result it now supports several different numerical types: integers, rationals, reals and complex numbers. However, the GCD algorithm is only really defined for integers, so is not applicable to any but the first of these. As a result, we can restrict "ordinary numbers" to be integers, and so only need to deal with the integer package here.

We can use the GCD implementation provided in the book, adding this to the integer package, and then install this procedure, tagging its results:

(define (gcd a b)
  (if (= b 0)
      a
      (gcd b (remainder a b))))
(put 'gcd '(integer integer)
     (lambda (x y) (tag (gcd x y))))

Adding a Generic GCD Operation

As we have installed appropriate GCD procedures in our table for the integer and polynomial types, greatest-common-divisor can be implemented as a normal generic operation, invoking apply-generic:

(define (greatest-common-divisor a b)
  (apply-generic 'gcd a b))

Attempts to apply this to types other than integer or polynomial will, of course, fail. In this respect our implementation corresponds pretty close to the Scheme interpreter I'm using:

> (gcd 3.4 2.3)

In standard input:
  10: 0* [gcd {3.4} 2.3]
In .../section-2.5.3/base-definitions.scm:
   9: 1  (if (= b 0) a (gcd b (remainder a b)))
  11: 2  [gcd 2.3 ...
  11: 3*  [remainder {3.4} 2.3]

.../section-2.5.3/base-definitions.scm:11:14:
    In procedure remainder in expression (remainder a b):
.../section-2.5.3/base-definitions.scm:11:14:
    Wrong type argument in position 1: 3.4
ABORT: (wrong-type-arg)
Backtrace:
> (greatest-common-divisor  (make-real 3.4) (make-real 2.3))

In standard input:
  11: 0* [greatest-common-divisor (real . 3.4) (real . 2.3)]
In .../section-2.5.3/sicp-2.94.scm:
 685: 1  [apply-generic gcd (real . 3.4) (real . 2.3)]
In .../section-2.5.3/base-apply-generic.scm:
   ...
  73: 2  (let* ((result #)) (if (in-tower? result) (drop result) result))
  73: 3* [find-and-apply-op]
  59: 4  (let* ((type-tags #) (proc #)) (if proc (apply proc #) (if # #)))
   ...
  70: 5  [error "No method for these types -- APPLY-GENERIC" (gcd (real real))]

.../section-2.5.3/base-apply-generic.scm:70:21:
    In procedure error in expression
    (error "No method for these types -- APPLY-GENERIC"
     list op type-tags)):
.../section-2.5.3/base-apply-generic.scm:70:21:
    No method for these types -- APPLY-GENERIC (gcd (real real))
Backtrace:

Putting It to the Test

We're given a specific example to try in the exercise statement:

> (define p1 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 1)
                                                (make-integer -1)
                                                (make-integer -2)
                                                (make-integer 2)
                                                zero)))
> (define p2 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 1)
                                                zero
                                                (make-integer -1)
                                                zero)))
> (greatest-common-divisor p1 p2)
(polynomial x sparse-terms (term 2 (integer . -1))
                           (term 1 (integer . 1)))

So, according to our system:

GCD(x4 - x3 - 2x2 + 2x, x3 - x) = -x2 + x

One way of checking our results is to first check that this calculated GCD is actually a divisor of both p1 and p2,then compare the calculated GCD to the calculated quotients to ensure that the former is "greater" than the latter in both cases and, finally, to ensure that the calculated quotients are prime relative to each other. We can do the first step by applying div twice, with the calculated GCD as the divisor in both cases, but with p1 as the dividend in the first instance and p2 in the second:

> (div p1 (greatest-common-divisor p1 p2))
((polynomial x sparse-terms (term 2 (rational (integer . 1) integer . -1))
                                    (term 0 (rational (integer . -2) integer . -1)))
 (integer . 0))
> (div p2 (greatest-common-divisor p1 p2))
((polynomial x dense-terms (rational (integer . 1) integer . -1)
                           (integer . -1))
 (integer . 0))

The first thing you'll note is that the calculated quotients have not been reduced to their simplest terms. In both cases the calculated quotient contains at least one rational coefficient that is obviously equivalent to an integer value. This is because, as part of exercise 2.93, we removed the code which reduced rational fractions to their simplest terms. This was a necessary first step in introducing rational functions. Of course we're trying to remedy that removal with this exercise and the remaining exercises in this chapter, but for now we'll just have to live with it.

Having performed these calculations we can see that the calculated GCD is definitely a whole divisor of both p1 and p2 as, in both cases, the remainder is equal to zero. We can also see that, in both cases, the calculated GCD is "greater" than the calculated quotient:

• In the case of p1 the calculated GCD and quotient are both second-order polynomials, but the calculated GCD has a higher coefficient for the first term where the coefficients differ. While both have a coefficient of -1 for the second-order term, the GCD has a coefficient of 1 for the first-order term while the calculated quotient has a coefficient of 0.
• In the case of p2, the calculated quotient is a lower-order polynomial than the GCD. The calculated quotient is first-order, while the calculated GCD is second-order.

All that remains to confirm that the calculated GCD is correct is to check that the quotients are prime relative to each other. We can do this by simply trying to calculate the GCD of the quotients. If this is equal to 1 then they're relatively prime:

> (greatest-common-divisor (car (div p1 (greatest-common-divisor p1 p2)))
                              (car (div p2 (greatest-common-divisor p1 p2))))
(integer . 1)

All looks good! Let's move on to the next exercise.

Addendum: A Simpler gcd-poly

The exercise requires that the procedure should signal an error if the two polys are not in the same variable. However, back in exercise 2.92 we introduced the procedure coerce-and-call, which ensures that two polynomials are expressed in the same variable and then later expanded it so that we ensured that polynomials were expressed in canonical form. This approach allows us to entirely remove the test of the indeterminates from gcd-poly, giving the following implementation:

(define (gcd-poly p1 p2)
  (make-poly (select-variable-from-polys p1 p2)
             (gcd-terms (term-list p1) (term-list p2))))

...which we can use provided we install the procedure as follows:

(put 'gcd '(polynomial polynomial)
     (lambda (p1 p2) (tag (coerce-and-call p1 p2 gcd-poly))))

SICP Exercise 2.93: Rational Functions

Modify the rational-arithmetic package to use generic operations, but change make-rat so that it does not attempt to reduce fractions to lowest terms. Test your system by calling make-rational on two polynomials to produce a rational function

(define p1 (make-polynomial 'x '((2 1)(0 1))))
(define p2 (make-polynomial 'x '((3 1)(0 1))))
(define rf (make-rational p2 p1))

Now add rf to itself, using add. You will observe that this addition procedure does not reduce fractions to lowest terms.

Let's start by looking at our current rational package implementation:

(define (install-rational-package)
  ;; internal procedures
  (define (numer x) (car x))
  (define (denom x) (cdr x))
  (define (make-rat n d)
    (if (and (integer? n) (integer? d))
        (let ((g (gcd n d)))
          (cons (/ n g) (/ d g)))
        (error "non-integer numerator or denominator"
               (list n d))))
  (define (add-rat x y)
    (make-rat (+ (* (numer x) (denom y))
                 (* (numer y) (denom x)))
              (* (denom x) (denom y))))
  (define (sub-rat x y)
    (make-rat (- (* (numer x) (denom y))
                 (* (numer y) (denom x)))
              (* (denom x) (denom y))))
  (define (mul-rat x y)
    (make-rat (* (numer x) (numer y))
              (* (denom x) (denom y))))
  (define (div-rat x y)
    (make-rat (* (numer x) (denom y))
              (* (denom x) (numer y))))
  (define (equ? x y)
    (and (= (numer x) (numer y))
         (= (denom x) (denom y))))
  (define (=zero? x) (= (numer x) 0))
  (define (addd x y z)
    (make-rat (+ (* (numer x) (denom y) (denom z))
                 (* (denom x) (numer y) (denom z))
                 (* (denom x) (denom y) (numer z)))
              (* (denom x) (denom y) (denom z))))
  (define (rational->real r) (make-real (/ (numer r) (denom r))))
  (define (rational->integer r) (make-integer (round (/ (numer r) (denom r)))))
  (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))))
  (define (negate-rat x)
    (make-rat (- (numer x)) (denom x)))
  ;; interface to rest of the system
  (define (tag x) (attach-tag 'rational x))
  (put 'add '(rational rational)
       (lambda (x y) (tag (add-rat x y))))
  (put 'sub '(rational rational)
       (lambda (x y) (tag (sub-rat x y))))
  (put 'mul '(rational rational)
       (lambda (x y) (tag (mul-rat x y))))
  (put 'div '(rational rational)
       (lambda (x y) (tag (div-rat x y))))
  (put 'equ? '(rational rational) equ?)
  (put '=zero? '(rational) =zero?)
  (put 'addd '(rational rational rational)
       (lambda (x y z) (tag (addd x y z))))
  (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))))
  (put 'negate '(rational)
       (lambda (x) (tag (negate-rat x))))
  (put 'make 'rational
       (lambda (n d) (tag (make-rat n d))))
  (put-coercion 'rational 'real rational->real)
  (put-coercion 'rational 'integer rational->integer)
  'done)

(define (make-rational n d)
  ((get 'make 'rational) n d))

We've added a lot in there and made a lot of changes since the rational package was first introduced. As a result it's not simply a case of doing a find-and-replace, mapping the primitive operations of +, -, * and / to their generic equivalents of add, sub, mul and div, although that does form part of the modifications.

Let's go through the package and identify the changes we need to make.

Rational Representation

Changing the values we store for the numerator and denominator from primitive integers to tagged types does not affect the overarching representation; we can still use a pair to hold these. As a result numer and denom remain unchanged. However, we do need to modify the creation of the rational representation provided by make-rat:

• In order to produce mathematically correct rational numbers, make-rat currently restricts us to constructing rational representations using primitive integer numerator and denominator values. Extending the rational package so that it can also be used to represent rational functions requires us to allow polynomials as well. We're doing this by modifying the package so that it uses tagged types for its numerator and denominator. However, while we'll now accept tagged types, in order to ensure mathematical correctness we should ensure that the numerator and denominator are restricted to integer or polynomial types.
• The integer and real packages include procedures for performing type conversions to the rational type (i.e. integer->rational and real->rational respectively). Both invoke make-rational as part of their type conversion, and pass in primitive integer values. If we're going to modify the rational package to use tagged types and generic operations then our changes need to account for these type conversion procedures. One approach would be to modify the integer->rational and real->rational procedures so that they pass tagged integer values to make-rational. However, we can limit the scope of our changes to the rational package alone if, in addition to integer and polynomial tagged types, we also accept primitive integers as the numerator and/or denominator for make-rat and have make-rat automatically convert them to tagged integer types. To do this we simply use a cond in make-rat, with the predicates in the first two clauses testing respectively whether the numerator or the denominator is a primitive integer, creating a tagged integer and recursively calling with the tagged value if the predicate matches.

These changes look as follows:

(define (numer x) (car x))
(define (denom x) (cdr x))
(define (valid-component? c)
  (memq (type-tag c) '(integer polynomial)))
(define (make-rat n d)
  (cond ((integer? n) (make-rat (make-integer n) d))
        ((integer? d) (make-rat n (make-integer d)))
        ((and (valid-component? n) (valid-component? d)) (cons n d))
        (else (error
               "numerator and denominator must both be integer or polynomial types"
               (list n d)))))

We've introduced the procedure valid-component? here which checks whether a value is tagged as an integer or polynomial. This is used by make-rat to perform the final validation of the types of the passed arguments once it's checked that any primitive integers have been wrapped.

Basic Arithmetic Operations

For the four basic operations of addition, subtraction, multiplication and division a find-and-replace will suffice. We replace + with add, - with sub and * with mul (the primitive division operation, /, is unused here):

(define (add-rat x y)
  (make-rat (add (mul (numer x) (denom y))
                 (mul (numer y) (denom x)))
            (mul (denom x) (denom y))))
(define (sub-rat x y)
  (make-rat (sub (mul (numer x) (denom y))
                 (mul (numer y) (denom x)))
            (mul (denom x) (denom y))))
(define (mul-rat x y)
  (make-rat (mul (numer x) (numer y))
            (mul (denom x) (denom y))))
(define (div-rat x y)
  (make-rat (mul (numer x) (denom y))
            (mul (denom x) (numer y))))

The procedure square-rat simply calls through to mul-rat without inspecting its argument, so it can remain unchanged:

(define (square-rat x)
  (mul-rat x x))

There are a few other basic arithmetic operations, equ?, =zero?, addd and negate-rat, which need only minor tweaks:

• We can replace = with a call to the generic equ? operation in the equ? procedure. Similarly, we can replace the comparison of the numerator with 0 in =zero? with a call to the generic =zero? on the numerator. However, as these rational operations have identical names to the generic operations, the globally scoped generic operations are obscured by the rational package's implementations. As a result, the interpreter will not invoke the global generic operations with these names but will instead recursively invoke the rational package's implementations, which will then fail as the passed values will not be of the correct type. To resolve this we'll rename the procedures within the rational package to equ-rat? and =zero-rat? respectively and change the corresponding put invocations which install them appropriately.
• Substituting add and mul for + and * respectively into the procedure addd, which we built as part of exercise 2.82, will not work. The primitive operation * is used as a ternary operation in the current implementation of addd. However, our generic arithmetic operation add is a binary operation. We can address this simply by converting the implementation of addd to use add-rat twice in succession: adding the first two arguments together first, which produces a rational result, and then adding the third argument to the result of that.
• Similarly, simply substituting sub for - into the procedure negate-rat will not work. In this case - is currently used as a unary operation with the semantics that this negates the primitive integer value passed to it. However, our generic arithmetic operation sub is a binary operation. Fixing this is simpler than for addd: we make the substitution, but we pass 0 as the first argument to sub, and x as the second (as 0 - n = -n).

Here are the updated operations:

(define (equ-rat? x y)
  (and (equ? (numer x) (numer y))
       (equ? (denom x) (denom y))))
(define (=zero-rat? x) (=zero? (numer x)))
(define (addd x y z)
  (add-rat (add-rat x y) z))
(define (negate-rat x)
  (make-rat (sub zero (numer x)) (denom x)))
…
(put 'equ? '(rational rational) equ-rat?)
(put '=zero? '(rational) =zero-rat?)

Type Conversion, Square Root and Trigonometric Operations

The remaining operations are:

• The type conversion operations: rational->real and rational->integer.
• The square root operation: sqrt-rat.
• The trigonometric operations: arctan-rat, cosine-rat and sine-rat.

These all share a couple of somewhat problematic properties:

1. Applying these operations to polynomials is either complicated, or does not make sense. We can make our lives easier in this exercise if we state that we'll only support these operations for rational numbers; rational functions are not supported.
2. When applied to integer-tagged values these operations require access to the underlying primitive integer values in order that we can compute their values easily. This breaks the encapsulation of the tagged types but is necessary as all of the operations need to calculate the result of dividing the numerator by the denominator in order to apply a primitive operation to the result. Note that we can't simply use div to calculate this as applying div to the two tagged integer values representing the numerator and denominator will just produce the original tagged rational value we started with! We'll deal with resolving this problem below...

Avoiding the Issue: Using real Operations

Let's deal with sqrt-rat and the trigonometric operations first. Note that, while implementing these directly for rational numbers involves breaking encapsulation, no such problem arises if we first raise a rational number to a real number, and then delegate the calculation to that package! As a result, each of these operations can be achieved by:

1. Validating that the parameter(s) to the operation are rational numbers and not rational functions (i.e. both the numerator and denominator are integers).
2. Raising the rational numbers to real numbers. Of course we'll need to re-tag the rational numbers as such before attempting this, otherwise raise won't know what type of number it's dealing with! Note that, assuming we fix our type conversion operations, this raise will be guaranteed to succeed and result in a real number.
3. Applying the corresponding generic operation (e.g. square-root for sqrt-rat, and so on). The generic operations will then apply the implementations defined for real numbers and so avoid us having to implement them directly!

We can encapsulate the first step of this process in a procedure that checks whether or not all rationals in a list are rational numbers, with integer numerators and denomniators, as follows:

(define (all-rational-numbers? rs)
  (cond ((null? rs) true)
        ((and (eq? 'integer (type-tag (numer (car rs))))
              (eq? 'integer (type-tag (denom (car rs)))))
         (all-rational-numbers? (cdr rs)))
        (else false)))

We can then capture the full process in a procedure which takes a generic operation that is to be applied and a list of rationals to apply it to, tests that they are all rational numbers, tags and raises the rationals and finally applies the generic operation using apply as follows:

(define (apply-as-real-if-all-rational-numbers f rs)
  (if (all-rational-numbers? rs)
      (apply f (map raise (map tag rs)))
      (error (list "numerator and denominator must both be integer to apply " f))))

The final steps are to update sqrt-rat and the trigonometric operations to utilize this procedure, and to change the installation of these procedures via put. The latter is necessary as the results of apply-as-real-if-all-rational-numbers will be a valid tagged value if it succeeds, so there's no need to turn the results into a real number anymore. This allows us to remove the wrapping λ-functions here:

(define (sqrt-rat x)
  (apply-as-real-if-all-rational-numbers square-root (list x)))
(define (arctan-rat x y)
  (apply-as-real-if-all-rational-numbers arctan (list x y)))
(define (cosine-rat x)
  (apply-as-real-if-all-rational-numbers cosine (list x)))
(define (sine-rat x)
  (apply-as-real-if-all-rational-numbers sine (list x)))
…
(put 'sqrt '(rational) sqrt-rat)
(put 'arctan '(rational rational) arctan-rat)
(put 'cosine '(rational) cosine-rat)
(put 'sine '(rational) sine-rat)

Type Conversion and Breaking Encapsulation

All that's left is to deal with the type conversion operations rational->real and rational->integer. While we're only going to support rational numbers here, it's pretty much impossible to avoid the stumbling block that in order to perform these conversions we need to understand the internal representation of the integer type from within the rational package in order that we can perform the primitive division necessary to support the type conversion. I.e. we end up breaking the encapsulation.

In order to contain the scope of the damage I extracted out a procedure, to-primitive, which converts a rational representing a rational number into its equivalent Scheme primitive numerical value. This restricts the breakage to this procedure alone. The two type conversion operations invoke this, blissfully unaware of the issue, and do not need to access the internals of the numerator and denominator in order to perform the division themselves. Note that to-primitive returns 0 if either the numerator or denominator is a polynomial. This is to allow the existing drop behaviour to continue function correctly - if it's not possible to perform a conversion then the value remains unchanged.

Here's the code:

(define (to-primitive x)
  (let ((n (numer x))
        (d (denom x)))
    (if (and (eq? 'integer (type-tag n)) (eq? 'integer (type-tag d)))
        (/ (contents n) (contents d))
        0)))
(define (rational->real r)
      (make-real (to-primitive r)))
(define (rational->integer r)
      (make-integer (round (to-primitive r))))

Yuck! I do apologise!

Of course, if anyone's got a better solution here then please shout!

Rational Functions In Action

Okay, having completed our implementation, let's see it in action... Here's the tests from the book:

> (define p1 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 1)
                                                zero
                                                (make-integer 1))))
> (define p2 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 1)
                                                zero
                                                zero
                                                (make-integer 1))))
> (define rf (make-rational p2 p1))
> (add rf rf)
(rational (polynomial x sparse-terms (term 5 (integer . 2))
                                     (term 3 (integer . 2))
                                     (term 2 (integer . 2))
                                     (term 0 (integer . 2)))
           polynomial x sparse-terms (term 4 (integer . 1))
                                     (term 2 (integer . 2))
                                     (term 0 (integer . 1)))

As you can see this adds the rational function to itself correctly but, as noted in the exercise statement, fails to reduce the result to lowest terms.

Now let's exercise the package with some rational fractions and check we haven't broken anything obvious:

> (define r1 (make-rational 10 3))
> (define r2 (make-rational 3 4))
> (add r1 r2)
(rational (integer . 49) integer . 12)
> (sub r1 r2)
(rational (integer . 31) integer . 12)
> (mul r1 r2)
(rational (integer . 30) integer . 12)
> (div r1 r2)
(rational (integer . 40) integer . 9)
> (equ? r1 r2)
#f
> (equ? r1 r1)
#t
> (=zero? r1)
#f
> (=zero? (make-rational zero 4))
#t
> (addd r1 r1 r1)
(rational (integer . 270) integer . 27)
> (square r2)
(rational (integer . 9) integer . 16)
> (square-root (square r2))
(rational (integer . 3) integer . 4)
> (negate r1)
(rational (integer . -10) integer . 3)

For the trigonometric operations we'll compare the results from our package with the results of Scheme's built-in trigonometric operations. These should be identical as the former ultimately delegates to the latter...

> (arctan r1 r2)
(rational (integer . 3038763055969609) integer . 2251799813685248)
> (exact->inexact (/ 3038763055969609 2251799813685248))
1.34948188444711
> (atan 10/3 3/4)
1.34948188444711
> (cosine r2)
(rational (integer . 6590467434422559) integer . 9007199254740992)
> (exact->inexact (/ 6590467434422559 9007199254740992))
0.731688868873821
> (cos 3/4)
0.731688868873821
> (sine r2)
(rational (integer . 6139656131284749) integer . 9007199254740992)
> (exact->inexact (/ 6139656131284749 9007199254740992))
0.681638760023334
> (sin 3/4)
0.681638760023334

All looks good from our brief test... Now let's move on to writing GCD for polynomials in exercise 2.94!

Guess What?

Yep, I stalled again.

I've got a couple of excuses. No really I have.

First it was holiday season. The combination of Thanksgiving, Christmas and so on make things somewhat busy at my work. The spare time I normally spend doing this was diverted away to more work-related stuff. It was challenging-but-enjoyable stuff, but meant I couldn't make much progress on SICP.

Secondly my second child was born in December. If you've kids of your own you'll know just how much time a tiny baby can absorb!

Anyway, I'm starting to get back on track now. I've got a working implementation of exercise 2.93 which I'll try to write up and post over the next week or so. I'll also go back and make a minor correction to some of the recent work we've been doing as I realized I'd left in the support for the 'scheme-number tagged type (i.e. Scheme's numerical primitive types). This was somewhat messing up the work needed for 2.93 so I've stripped it out.