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))))