SICP Exercise 2.90: Supporting Dense and Sparse Polynomials - Part 3

Suppose we want to have a polynomial system that is efficient for both sparse and dense polynomials. One way to do this is to allow both kinds of term-list representations in our system. The situation is analogous to the complex-number example of section 2.4, where we allowed both rectangular and polar representations. To do this we must distinguish different types of term lists and make the operations on term lists generic. Redesign the polynomial system to implement this generalization. This is a major effort, not a local change.

Okay, we've talked about it at length and updated the system so that polynomials form a part of our overall generic arithmetic system. Let's put it to the test and check it all works.

First we should define a few polynomials to work with:

(define dense
        (make-polynomial-from-coeffs 'x
                                     (list (make-integer 4)
                                           (make-integer 3)
                                           (make-integer 2)
                                           (make-integer 1)
                                           zero)))

(define dense-with-many-zeros
        (make-polynomial-from-coeffs 'x
                                     (list (make-integer 42)
                                           zero
                                           zero
                                           zero
                                           zero
                                           zero
                                           (make-integer -1))))

(define sparse
        (make-polynomial-from-terms 'x 
                                    (list (make-term 5 (make-integer 5))
                                          (make-term 3 (make-integer 3))
                                          (make-term 1 (make-integer 1)))))

(define another-sparse
        (make-polynomial-from-terms 'x
                                    (list (make-term 5 (make-integer 5))
                                          (make-term 3 (make-integer 3))
                                          (make-term 1 (make-integer 1))
                                          (make-term 0 (make-integer 3)))))

(define very-sparse
        (make-polynomial-from-terms 'x
                                    (list (make-term 50 (make-integer 150))
                                          (make-term 10 (make-integer 11))
                                          (make-term 0 (make-integer 1)))))

(define polypoly
        (make-polynomial-from-coeffs
         'x
         (list (make-polynomial-from-coeffs 'y
                                            (list (make-integer 2)
                                                  (make-integer 1))))))

Now let's try out a bunch of operations and see what happens:

> (add polypoly dense)
(polynomial x
            dense-terms
            (integer . 4)
            (integer . 3)
            (integer . 2)
            (integer . 1)
            (polynomial y
                        dense-terms
                        (integer . 2)
                        (integer . 1)))
> (add polypoly polypoly)
(polynomial x
            dense-terms
            (polynomial y
                        dense-terms
                        (integer . 4)
                        (integer . 2)))
> (add (add polypoly polypoly) (make-integer 3))
(polynomial x
            dense-terms
            (polynomial y
                        dense-terms
                        (integer . 4)
                        (integer . 5)))
> (add dense dense-with-many-zeros)
(polynomial x
            dense-terms
            (integer . 42)
            (integer . 0)
            (integer . 4)
            (integer . 3)
            (integer . 2)
            (integer . 1)
            (integer . -1))
> (add dense-with-many-zeros dense-with-many-zeros)
(polynomial x
            sparse-terms
            (term 6 (integer . 84))
            (term 0 (integer . -2)))
> (add sparse sparse)
(polynomial x
            sparse-terms
            (term 5 (integer . 10))
            (term 3 (integer . 6))
            (term 1 (integer . 2)))
> (add sparse another-sparse)
(polynomial x
            sparse-terms
            (term 5 (integer . 10))
            (term 3 (integer . 6))
            (term 1 (integer . 2))
            (term 0 (integer . 3)))
> (add very-sparse sparse)
(polynomial x
            sparse-terms
            (term 50 (integer . 150))
            (term 10 (integer . 11))
            (term 5 (integer . 5))
            (term 3 (integer . 3))
            (term 1 (integer . 1))
            (term 0 (integer . 1)))
> (mul sparse dense)
(polynomial x
            sparse-terms
            (term 9 (integer . 20))
            (term 8 (integer . 15))
            (term 7 (integer . 22))
            (term 6 (integer . 14))
            (term 5 (integer . 10))
            (term 4 (integer . 6))
            (term 3 (integer . 2))
            (term 2 (integer . 1)))
> (add dense sparse)
(polynomial x
            dense-terms
            (integer . 5)
            (integer . 4)
            (integer . 6)
            (integer . 2)
            (integer . 2)
            (integer . 0))
> (sub sparse dense)
(polynomial x
            sparse-terms
            (term 5 (integer . 5))
            (term 4 (integer . -4))
            (term 2 (integer . -2)))
> (negate very-sparse)
(polynomial x
            sparse-terms
            (term 50 (integer . -150))
            (term 10 (integer . -11))
            (term 0 (integer . -1)))
> (sub (add dense (make-integer 1)) dense)
(integer . 1)

SICP Exercise 2.90: Supporting Dense and Sparse Polynomials - Part 2

Suppose we want to have a polynomial system that is efficient for both sparse and dense polynomials. One way to do this is to allow both kinds of term-list representations in our system. The situation is analogous to the complex-number example of section 2.4, where we allowed both rectangular and polar representations. To do this we must distinguish different types of term lists and make the operations on term lists generic. Redesign the polynomial system to implement this generalization. This is a major effort, not a local change.

We've already discussed how we're going to approach this, so let's get onto the code. For a change I'm going to post all of the code involved in the system so far, rather than just the changes, so I'll include all the code for the integer, rational, real and complex packages too, along with all the supporting framework. We've been working on the system for quite a while now, so this'll help to provide a check point for where we've got to.

I'll intersperse the code with commentary where I feel it's appropriate... Oh, and it's worth noting that I've tweaked and tidied code from previous exercises as I've gone along, so some bits won't correspond exactly to what's gone before.

Here we go...

; If you're using DrRacket you might need to uncomment these...
;(require rnrs/base-6)
;(require rnrs/mutable-pairs-6)


;;;
;;; Basic definitions
;;;

(define false #f)
(define true #t)


;;;
;;; Look-up table from Creating local tables in section 3.3.3
;;;

(define (assoc key records)
  (cond ((null? records) false)
        ((equal? key (caar records)) (car records))
        (else (assoc key (cdr records)))))

(define (make-table)
  (let ((local-table (list '*table*)))
    (define (lookup key-1 key-2)
      (let ((subtable (assoc key-1 (cdr local-table))))
        (if subtable
            (let ((record (assoc key-2 (cdr subtable))))
              (if record
                  (cdr record)
                  false))
            false)))
    (define dump local-table)
    (define (insert! key-1 key-2 value)
      (let ((subtable (assoc key-1 (cdr local-table))))
        (if subtable
            (let ((record (assoc key-2 (cdr subtable))))
              (if record
                  (set-cdr! record value)
                  (set-cdr! subtable
                            (cons (cons key-2 value)
                                  (cdr subtable)))))
            (set-cdr! local-table
                      (cons (list key-1
                                  (cons key-2 value))
                            (cdr local-table)))))
      'ok)    
    (define (dispatch m)
      (cond ((eq? m 'lookup-proc) lookup)
            ((eq? m 'insert-proc!) insert!)
            ((eq? m 'dump) dump)
            (else (error "Unknown operation -- TABLE" m))))
    dispatch))

(define operation-table (make-table))
(define get (operation-table 'lookup-proc))
(define put (operation-table 'insert-proc!))
(define dump-table (operation-table 'dump))

(define (put-coercion source-type target-type proc)
  (put 'coercion (list source-type target-type) proc))

(define (get-coercion source-type target-type)
  (get 'coercion (list source-type target-type)))


;;;
;;; Tagging utils
;;;

(define (attach-tag type-tag contents)
  (cons type-tag contents))
(define (type-tag datum)
  (if (pair? datum)
      (car datum)
      (error "Bad tagged datum -- TYPE-TAG" datum)))
(define (contents datum)
  (if (pair? datum)
      (cdr datum)
      (error "Bad tagged datum -- CONTENTS" datum)))


;;;
;;; Generic Arithmetic operations
;;;

(define (add x y) (apply-generic 'add x y))
(define (sub x y) (apply-generic 'sub x y))
(define (mul x y) (apply-generic 'mul x y))
(define (div x y) (apply-generic 'div x y))
(define (equ? x y) (apply-generic 'equ? x y))
(define (=zero? x) (apply-generic '=zero? x))
(define (addd x y z) (apply-generic 'addd x y z))
(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))
(define (negate x) (apply-generic 'negate x))

First item of note: we update the tower-of-types to include the polynomial package. We also add in the dense-terms and sparse-terms packages to allow coercion between these two types, but we'll not install coercion procedures for polynomial to dense-terms or vice versa, effectively splitting the tower-of-types into two towers.

;;;
;;; Tower-of-types
;;;

(define tower-of-types '(integer rational real complex polynomial
                         dense-terms sparse-terms))


;;;
;;; Generic Operation Framework
;;;

(define (find-highest-type list-of-types)
  (define (filter-type type list-to-filter)
    (cond ((null? list-to-filter) '())
          ((eq? (car list-to-filter) type) (filter-type type (cdr list-to-filter)))
          (else (cons (car list-to-filter) (filter-type type (cdr list-to-filter))))))
  (define (find-highest highest remaining-tower remaining-list)
    (cond ((null? remaining-list) highest)
          ((null? remaining-tower)
           (error "Cannot find highest type from non-tower types -- FIND-HIGHEST-TYPE"
                  remaining-list))
          (else (find-highest (car remaining-tower)
                              (cdr remaining-tower)
                              (filter-type (car remaining-tower) remaining-list)))))
  (find-highest #f tower-of-types list-of-types))
  
(define (raise-to type value)
  (cond ((eq? type (type-tag value)) value)
        ((memq type tower-of-types) (raise-to type (raise value)))
        (else (error "Cannot raise to non-tower type -- RAISE-TO"
                     (list type tower-of-types)))))

(define (raise-all-to type values)
  (if (null? values)
      '()
      (cons (raise-to type (car values)) (raise-all-to type (cdr values)))))

As noted in the preceding discussion, in order to allow a missing coercion operation to indicate a point at which a type should not be lowered any further (and so to allow multiple towers-of-types to be represented by a single tower) we modify apply-raise so that it returns the value passed to it, rather than raising an error, if it can't find an appropriate coercion operation.

(define (apply-raise x types)
  (cond ((null? types)
         (error "Type not found in the tower-of-types"
                (list (type-tag x) tower-of-types)))
        ((eq? (type-tag x) (car types))
         (if (null? (cdr types))
             x
             (let ((raiser (get-coercion (type-tag x) (cadr types))))
               (if raiser
                   (raiser (contents x))
                   x))))
        (else (apply-raise x (cdr types)))))

(define (raise x)
  (apply-raise x tower-of-types))

(define (project x)
  (apply-raise x (reverse tower-of-types)))

(define (drop x)
  (let ((dropped (project x)))
    (if (eq? (type-tag x) (type-tag dropped))
        x
        (let ((raised (raise dropped)))
          (if (not (equ? x raised))
              x
              (drop dropped))))))

(define (raise-to type value)
  (cond ((eq? type (type-tag value)) value)
        ((memq type tower-of-types) (raise-to type (raise value)))
        (else (error "Cannot raise to non-tower type -- RAISE-TO"
                     (list type tower-of-types)))))

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

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

(define (apply-generic op . args)
  (define (find-and-apply-op)
    (let* ((type-tags (map type-tag args))
           (proc (get op type-tags)))
      (if proc
          (apply proc (map contents args))
          (if (> (length args) 1)
              (let* ((highest-type (find-highest-type type-tags))
                     (mapped-args (raise-all-to highest-type args))
                     (mapped-types (map type-tag mapped-args))
                     (mapped-proc (get op mapped-types)))
                (if mapped-proc
                    (apply mapped-proc (map contents mapped-args))
                    (error
                     "No method for these types -- APPLY-GENERIC"
                     (list op type-tags))))))))
  (let ((result (find-and-apply-op)))
    (if (in-tower? result)
        (drop result)
        result)))

The standard numerical type representations follow. They're mostly unchanged apart from the addition of support for the negate operation, and the addition of a coercion operation to raise a complex type to a polynomial.

;;;
;;; Integer package
;;;

(define (install-integer-package)
  (define (tag x)
    (attach-tag 'integer x))    
  (define (integer->rational i) (make-rational i 1))
  (put 'add '(integer integer)
       (lambda (x y) (tag (+ x y))))
  (put 'sub '(integer integer)
       (lambda (x y) (tag (- x y))))
  (put 'mul '(integer integer)
       (lambda (x y) (tag (* x y))))
  (put 'div '(integer integer)
       (lambda (x y) (make-rational x y)))
  (put 'equ? '(integer integer) =)
  (put '=zero? '(integer)
       (lambda (x) (= 0 x)))
  (put 'addd '(integer integer integer)
       (lambda (x y z) (tag (+ x y z))))
  (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))))
  (put 'negate '(integer)
       (lambda (x) (tag (- x))))
  (put 'make 'integer
       (lambda (x) (if (integer? x)
                       (tag x)
                       (error "non-integer value" x))))
  (put-coercion 'integer 'rational integer->rational)
  'done)

(define (make-integer n)
  ((get 'make 'integer) n))


;;;
;;; Rational numbers package
;;;

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


;;;
;;; Real package
;;;

(define (install-real-package)
  (define (tag x)
    (attach-tag 'real x))    
  (define (real->complex r) (make-complex-from-real-imag (tag r) (tag 0)))
  (define (real->rational r) (make-rational (inexact->exact (numerator r))
                                            (inexact->exact (denominator r))))
  (put 'add '(real real)
       (lambda (x y) (tag (+ x y))))
  (put 'sub '(real real)
       (lambda (x y) (tag (- x y))))
  (put 'mul '(real real)
       (lambda (x y) (tag (* x y))))
  (put 'div '(real real)
       (lambda (x y) (tag (/ x y))))
  (put 'equ? '(real real) =)
  (put '=zero? '(real)
       (lambda (x) (= 0 x)))
  (put 'addd '(real real real)
       (lambda (x y z) (tag (+ x y z))))
  (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))))
  (put 'negate '(real)
       (lambda (x) (tag (- x))))
  (put 'make 'real
       (lambda (x) (if (real? x)
                       (tag x)
                       (error "non-real value" x))))
  (put-coercion 'real 'complex real->complex)
  (put-coercion 'real 'rational real->rational)
  'done)

(define (make-real n)
  ((get 'make 'real) n))


;;;
;;; Complex (rectangular) number package
;;;

(define (install-rectangular-package)
  ;; internal procedures
  (define (real-part z) (car z))
  (define (imag-part z) (cdr z))
  (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
  (define (tag x) (attach-tag 'rectangular x))
  (put 'real-part '(rectangular) real-part)
  (put 'imag-part '(rectangular) imag-part)
  (put 'magnitude '(rectangular) magnitude)
  (put 'angle '(rectangular) angle)
  (put 'make-from-real-imag 'rectangular 
       (lambda (x y) (tag (make-from-real-imag x y))))
  (put 'make-from-mag-ang 'rectangular 
       (lambda (r a) (tag (make-from-mag-ang r a))))
  'done)

;;;
;;; Complex (polar) number package
;;;

(define (install-polar-package)
  ;; internal procedures
  (define (magnitude z) (car z))
  (define (angle z) (cdr z))
  (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
  (define (tag x) (attach-tag 'polar x))
  (put 'real-part '(polar) real-part)
  (put 'imag-part '(polar) imag-part)
  (put 'magnitude '(polar) magnitude)
  (put 'angle '(polar) angle)
  (put 'make-from-real-imag 'polar
       (lambda (x y) (tag (make-from-real-imag x y))))
  (put 'make-from-mag-ang 'polar 
       (lambda (r a) (tag (make-from-mag-ang r a))))
  'done)


;;;
;;; Generic complex number procedures
;;;

;;; N.B. Name-clash with complex number support in rnrs/base-6
;;; - prefixing complex ops with complex-

(define (complex-real-part z) (apply-generic 'real-part z))
(define (complex-imag-part z) (apply-generic 'imag-part z))
(define (complex-magnitude z) (apply-generic 'magnitude z))
(define (complex-angle z) (apply-generic 'angle z))

Note the body of the complex->polynomial coercion operation below: (make-zero-order-polynomial-from-coeff (drop (tag z))). We re-tag the complex number as such, as it will have no tag at this point, and then drop it as far as possible before using a special creation operation. The former is necessary in order to simplify the type as much as possible before creating the polynomial (and so to undo any coercion that may have occurred getting to this level in the tree). The latter allows us to encapsulate the special 'unbound marker within the polynomial package.

;;;
;;; Complex number package
;;;

(define (install-complex-package)
  ;; imported procedures from rectangular and polar packages
  (define (make-from-real-imag x y)
    ((get 'make-from-real-imag 'rectangular) x y))
  (define (make-from-mag-ang r a)
    ((get 'make-from-mag-ang 'polar) r a))
  ;; 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))))
  (define (complex->real z) (raise-to 'real (complex-real-part z)))
  (define (complex->polynomial z) (make-zero-order-polynomial-from-coeff (drop (tag z))))
  (define (negate-complex z)
    (make-from-real-imag (negate (complex-real-part z))
                         (negate (complex-imag-part z))))
  ;; interface to rest of the system
  (define (tag z) (attach-tag 'complex z))
  (put 'add '(complex complex)
       (lambda (z1 z2) (tag (add-complex z1 z2))))
  (put 'sub '(complex complex)
       (lambda (z1 z2) (tag (sub-complex z1 z2))))
  (put 'mul '(complex complex)
       (lambda (z1 z2) (tag (mul-complex z1 z2))))
  (put 'div '(complex complex)
       (lambda (z1 z2) (tag (div-complex z1 z2))))
  (put 'equ? '(complex complex) equ-complex?)
  (put '=zero? '(complex) =zero-complex?)
  (put 'addd '(complex complex complex)
       (lambda (z1 z2 z3) (tag (addd-complex z1 z2 z3))))
  ;;; N.B. Name-clash with complex number support in rnrs/base-6
  ;;; - prefixing complex ops with complex-
  (put 'real-part '(complex) complex-real-part)
  (put 'imag-part '(complex) complex-imag-part)
  (put 'magnitude '(complex) complex-magnitude)
  (put 'angle '(complex) complex-angle)
  (put 'negate '(complex)
       (lambda (z) (tag (negate-complex z))))
  (put 'make-from-real-imag 'complex
       (lambda (x y) (tag (make-from-real-imag x y))))
  (put 'make-from-mag-ang 'complex
       (lambda (r a) (tag (make-from-mag-ang r a))))
  (put-coercion 'complex 'real complex->real)
  (put-coercion 'complex 'polynomial complex->polynomial)
  'done)

(define (make-complex-from-real-imag x y)
  ((get 'make-from-real-imag 'complex) x y))
(define (make-complex-from-mag-ang r a)
  ((get 'make-from-mag-ang 'complex) r a))


;;;
;;; Install numerical type packages
;;;

(install-rectangular-package)
(install-polar-package)

(install-integer-package)
(install-rational-package)
(install-real-package)
(install-complex-package)

Zero is a useful concept... We used it in exercise 2.89 to avoid constantly recreating zero coefficients. We'll use it again here!

;;;
;;; Zero
;;;
(define zero (make-integer 0))

We're going to be able to construct both sparse and dense term-lists using terms, so we pull out the concept of a term into its own package. While it uses generic operations to encapsulate its internal representation, it doesn't exist within the tower-of-types.

;;;
;;; Terms
;;;

(define (install-term-package)
  ;; internal procedures
  ;; representation of terms
  (define (make-term order coeff) (list order coeff))
  (define (order-term term) (car term))
  (define (coeff-term term) (cadr term))
  (define (equ-term? t1 t2)
    (and (= (order-term t1) (order-term t2))
         (equ? (coeff-term t1) (coeff-term t2))))

  ;; interface to rest of the system
  (define (tag t) (attach-tag 'term t))
  (put 'order '(term) order-term)
  (put 'coeff '(term) coeff-term)
  (put 'equ? '(term term) equ-term?)
  (put 'make 'term
       (lambda (order coeff) (tag (make-term order coeff))))
  'done)

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

(define (order term)
  (apply-generic 'order term))

(define (coeff term)
  (apply-generic 'coeff term))

Before we get into the specific implementations of the term-lists, we need to know when to use which representation. As there are only two available representations in our system at the moment this is a binary choice, so we produce a predicate that can tell us when to use the sparse term-list representation, and then fall back to using the dense term-list representation in all other cases. Here's the predicate:

;;;
;;; Sparse/dense representation selector
;;;
(define (store-as-sparse? highest-order zero-terms)
  (if (>= highest-order 10)
      (> (/ zero-terms highest-order) 0.1)
      (> zero-terms (/ highest-order 5))))

The sparse term-list package is mostly derived from the implementation presented in section 2.5.3 and our work in exercises 2.87 and 2.88, with operations related to polynomial manipulation stripped out. There are some slight tweaks to the behaviour, of course...

• We ensure that we always have at least one term in any polynomial, even if its a single order 0 term with a coefficient of 0. This makes some manipulations easier, especially coercing from polynomial to complex, as we don't have to check especially for an empty term list. We use this same trick with the dense term-list package too.
• We add a creation operation that will create sparse term-lists from coefficients. This operation is mostly implemented by convert-to-term-list, which assigns an order to each coefficient based upon its position in the list.
• There are various additions in order to allow coercing a sparse term-list to a dense term-list. These include operations to support testing whether or not it's valid to represent a sparse term-list as a dense term-list as well as the coercion procedure itself, which simply calls through to the appropriate dense term-list creation operation if it's a viable coercion. An implementation of the equ? operation is also added, as drop coercions need to be able to test whether projecting and then re-raising a sparse polynomial gives a value equal to the original before they can allow drops to succeed.
• There are also various operations added to support the coercion of polynomial types to complex. The implementation of equ? is also used here, and a new operation, get-constant, is added, which returns the coefficient from the order 0 term of the term-list (or returns zero if there is no such term in the representation).
• The interface to the rest of the system changes, so that it exposes operations on term-lists rather than polynomials, and we also expose a few more operations.

;;;
;;; Sparse term-lists
;;; 
(define (install-sparse-terms-package)
  ;; internal procedures
  ;; 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))

  ;; Ensures a term list always has at least one term
  (define (ensure-valid-term-list terms)
    (if (empty-termlist? terms)
        (list (make-term 0 zero))
        terms))
               
  ;; creation
  (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-terms terms result)
    (if (null? terms)
        result
        (build-terms (cdr terms) (insert-term (car terms) result))))
  (define (make-from-terms terms)
    (build-terms terms (the-empty-termlist)))

  (define (convert-to-term-list coeffs)
    (if (null? coeffs)
        (the-empty-termlist)
        (adjoin-term (make-term (- (length coeffs) 1) (car coeffs))
                     (convert-to-term-list (cdr coeffs)))))
  (define (make-from-coeffs coeffs)
    (convert-to-term-list coeffs))

  ;; procedures used by add-poly
  (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-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))))))

  ;; procedures used by div-poly
  (define (div-terms L1 L2)
    (if (empty-termlist? L1)
        (list (the-empty-termlist) (the-empty-termlist))
        (let ((t1 (first-term L1))
              (t2 (first-term L2)))
          (if (> (order t2) (order t1))
              (list (the-empty-termlist) L1)
              (let* ((new-c (div (coeff t1) (coeff t2)))
                     (new-o (- (order t1) (order t2)))
                     (new-t (make-term new-o new-c))
                     (rest-of-result
                      (div-terms
                       (add-terms L1
                                  (negate-terms (mul-term-by-all-terms new-t L2)))
                        L2)))
                  (list (adjoin-term new-t (car rest-of-result))
                        (cadr rest-of-result)))))))

  ;; zero test
  (define (=zero-all-terms? L)
    (cond ((empty-termlist? L) #t)
          ((not (=zero? (coeff (first-term L)))) #f)
          (else (=zero-all-terms? (rest-terms L)))))
        
  ;; Negation
  (define (negate-terms L)
    (if (empty-termlist? L)
        (the-empty-termlist)
        (let ((term (first-term L)))
          (adjoin-term (make-term (order term)
                                  (negate (coeff term)))
                       (negate-terms (rest-terms L))))))

  ;; Equality
  (define (equ-terms? L1 L2)
    (cond ((empty-termlist? L1) (empty-termlist? L2))
          ((empty-termlist? L2) #f)
          (else (and (equ? (first-term L1) (first-term L2))
                     (equ-terms? (rest-terms L1) (rest-terms L2))))))

  ;; Constant
  (define (get-sparse-constant L)
    (cond ((empty-termlist? L) zero)
          ((= (order (first-term L)) 0) (coeff (first-term L)))
          (else (get-sparse-constant (rest-terms L)))))
  
  ;; Coercion
  (define (calculate-zero-terms first rest)
    (if (empty-termlist? rest)
        (order first)
        (let ((next (first-term rest)))
          (+ (- (order first) (order next) 1)
             (calculate-zero-terms next (rest-terms rest))))))
         
  (define (keep-as-sparse? L)
    (if (empty-termlist? L)
        #f
        (let ((highest-order (order (first-term L)))
              (zero-terms (calculate-zero-terms (first-term L) (rest-terms L))))
          (store-as-sparse? highest-order zero-terms))))

  (define (sparse-terms->dense-terms L)
    (if (keep-as-sparse? L)
        (tag L)
        ((get 'make-from-terms 'dense-terms) L)))

  ;; interface to rest of the system
  (define (tag t) (attach-tag 'sparse-terms (ensure-valid-term-list t)))
  (put 'add '(sparse-terms sparse-terms) 
       (lambda (t1 t2) (tag (add-terms t1 t2))))
  (put 'mul '(sparse-terms sparse-terms) 
       (lambda (t1 t2) (tag (mul-terms t1 t2))))
  (put 'div '(sparse-terms sparse-terms) 
       (lambda (t1 t2)
         (let ((result (div-terms t1 t2)))
           (list (sparse-terms->dense-terms (car result))
                 (sparse-terms->dense-terms (cadr result))))))
  (put 'equ? '(sparse-terms sparse-terms) equ-terms?)
  (put '=zero? '(sparse-terms) =zero-all-terms?)
  (put 'negate '(sparse-terms)
       (lambda (t) (tag (negate-terms t))))
  (put 'get-constant '(sparse-terms) get-sparse-constant)
  (put 'make-from-terms 'sparse-terms
       (lambda (terms) (tag (make-from-terms terms))))
  (put 'make-from-coeffs 'sparse-terms
       (lambda (coeffs) (tag (make-from-coeffs coeffs))))
  (put-coercion 'sparse-terms 'dense-terms sparse-terms->dense-terms)
  'done)

Similarly, the dense term-list package is mostly derived from the our work in exercise 2.89, with operations related to polynomial manipulation stripped out and a similar bunch of slight tweaks to the behaviour...

• We ensure that we always have at least one term in any polynomial.
• We add a creation operation that will create dense term-lists from terms. This operation is derived from the sparse term-list's equivalent creation operation.
• Coercing a dense term-list to a sparse term-list is much simpler than the reverse - we add a coercion operation which simply calls through to the appropriate sparse term-list creation. However, we won't get automatic raising to sparse term-lists when a dense term-list representation is inappropriate, so we add an operation, to-best-representation, which takes the result of an arithmetic operation and either tags it as a dense term-list or raises it to a sparse term-list as appropriate. We then modify all of the arithmetic operations to use this instead of tag.
• The same set of operations are added to support the coercion of polynomial types to complex as were added for sparse term-lists. I.e. implementations of equ? and get-constant.
• The interface to the rest of the system changes in an identical manner to sparse term-lists.

;;;
;;; Dense term-lists
;;; 
(define (install-dense-terms-package)
  ;; internal procedures
  ;; representation of term lists
  (define (adjoin-term term-order term-coeff term-list)
    (cond ((=zero? term-coeff) term-list)
          ((= term-order (+ 1 (term-list-order term-list)))
           (cons term-coeff term-list))
          ((> term-order (term-list-order term-list))
           (adjoin-term term-order
                        term-coeff
                        (cons zero term-list)))
          (else (error "Cannot adjoin term of lower order than term list -- ADJOIN-TERM"
                       (list term-order term-coeff 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 (term-list-order term-list)
    (- (length term-list) 1))

  ;; Ensures a term list always has at least one term
  (define (ensure-valid-term-list terms)
    (if (empty-termlist? terms)
        (list zero)
        terms))

  ;; Creation
  (define (strip-leading-zeros coeffs)
    (cond ((empty-termlist? coeffs) (the-empty-termlist))
          ((not (=zero? (first-term coeffs))) coeffs)
          (else (make-from-coeffs (rest-terms coeffs)))))
  (define (make-from-coeffs coeffs) coeffs)
  
  (define (insert-term term terms)
    (if (empty-termlist? terms)
        (adjoin-term (order term) (coeff term) (the-empty-termlist))
        (let ((head-order (term-list-order terms))
              (term-order (order term)))
          (cond ((> term-order head-order)
                 (adjoin-term term-order (coeff term) terms))
                ((= term-order head-order)
                 (adjoin-term term-order
                              (add (coeff term) (car terms))
                              (rest-terms terms)))
                (else (adjoin-term head-order
                                   (car terms)
                                   (insert-term term (rest-terms terms))))))))
  (define (build-terms terms result)
    (if (null? terms)
        result
        (build-terms (cdr terms) (insert-term (car terms) result))))
  (define (make-from-terms terms)
    (build-terms terms (the-empty-termlist)))
  
  ;; procedures used by add-poly
  (define (add-terms L1 L2)
    (cond ((empty-termlist? L1) L2)
          ((empty-termlist? L2) L1)
          (else
           (let* ((t1 (first-term L1))
                  (t2 (first-term L2))
                  (o1 (term-list-order L1))
                  (o2 (term-list-order L2)))
             (cond ((> o1 o2)
                    (adjoin-term
                     o1 t1 (add-terms (rest-terms L1) L2)))
                   ((< o1 o2)
                    (adjoin-term
                     o2 t2 (add-terms L1 (rest-terms L2))))
                   (else
                    (adjoin-term
                     o1
                     (add t1 t2)
                     (add-terms (rest-terms L1)
                                (rest-terms L2)))))))))

  ;; procedures used by mul-poly
  (define (mul-terms L1 L2)
    (if (empty-termlist? L1)
        (the-empty-termlist)
        (add-terms (mul-term-by-all-terms (term-list-order L1) (first-term L1) L2)
                   (mul-terms (rest-terms L1) L2))))
  (define (mul-term-by-all-terms o1 c1 L)
    (if (empty-termlist? L)
        (the-empty-termlist)
        (let* ((t2 (first-term L))
               (new-order (+ o1 (term-list-order L))))
          (adjoin-term
           new-order
           (mul c1 t2)
           (mul-term-by-all-terms o1 c1 (rest-terms L))))))

  ;; zero test
  (define (=zero-all-terms? L)
    (cond ((empty-termlist? L) #t)
          ((not (=zero? (first-term L))) #f)
          (else (=zero-all-terms? (rest-terms L)))))
        
  ;; Negation
  (define (negate-terms L)
    (if (empty-termlist? L)
        (the-empty-termlist)
        (let ((term (first-term L)))
          (adjoin-term (term-list-order L)
                       (negate term)
                       (negate-terms (rest-terms L))))))
  
  ;; Equality
  (define (equ-terms? L1 L2)
    (cond ((empty-termlist? L1) (empty-termlist? L2))
          ((empty-termlist? L2) #f)
          (else (and (equ? (first-term L1) (first-term L2))
                     (equ-terms? (rest-terms L1) (rest-terms L2))))))
  
  ;; Constant
  (define (get-dense-constant L)
    (cond ((empty-termlist? L) zero)
          ((= (term-list-order L) 0) (first-term L))
          (else (get-dense-constant (rest-terms L)))))
  
  ;; Coercion
  (define (dense-terms->sparse-terms L)
    ((get 'make-from-coeffs 'sparse-terms) L))
  
  (define (count-zero-terms L)
    (if (empty-termlist? L)
        0
        (+ (if (=zero? (first-term L)) 1 0)
           (count-zero-terms (rest-terms L)))))
  
  (define (to-best-representation L)
    (if (store-as-sparse? (term-list-order L) (count-zero-terms L))
        (dense-terms->sparse-terms L)
        (tag L)))
  
  ;; interface to rest of the system
  (define (tag t) (attach-tag 'dense-terms (ensure-valid-term-list t)))
  (put 'add '(dense-terms dense-terms) 
       (lambda (t1 t2) (to-best-representation (add-terms t1 t2))))
  (put 'mul '(dense-terms dense-terms) 
       (lambda (t1 t2) (to-best-representation (mul-terms t1 t2))))
  (put 'equ? '(dense-terms dense-terms) equ-terms?)
  (put '=zero? '(dense-terms) =zero-all-terms?)
  (put 'negate '(dense-terms)
       (lambda (t) (to-best-representation (negate-terms t))))
  (put 'get-constant '(dense-terms) get-dense-constant)
  (put 'make-from-terms 'dense-terms
       (lambda (terms) (tag (make-from-terms terms))))
  (put 'make-from-coeffs 'dense-terms
       (lambda (coeffs) (tag (make-from-coeffs coeffs))))
  (put-coercion 'dense-terms 'sparse-terms dense-terms->sparse-terms)
  'done)

(define (get-constant term-list)
  (apply-generic 'get-constant term-list))

Finally we have the polynomial package that wraps the sparse and dense term-list representations. This basically takes the polynomial operations from the preceding exercises, strips out the term-list operations and replaces them with calls to the appropriate generic operations (i.e. add-terms becomes add, mul-terms becomes mul, and so on. And then we add a few tweaks...

• There's the previously noted changes to same-variable? and the addition of select-variable to support zero-order "unbound" polynomials
• We add an implementation of equ? to support coercion from polynomial to complex.
• We add the coercion operation itself. Note that this takes the coefficient of the order 0 term from the underlying term-list and raises it to a complex so that it's of the correct type. Note that the coefficient may already be a complex (or may in fact be a polynomial itself), so we only try to raise-to a complex if it's lower than a complex.

One thing I didn't do in the polynomial package was to automatically select the term-list representation on construction. I decided to allow the clients of the system to determine the initial representation for their polynomials (or I was being lazy - you decide!).

;;;
;;; Polynomial wrapper package
;;; 
(define (install-polynomial-package)
  ;; internal procedures
  ;; representation of poly
  (define (make-poly variable term-list)
    (cons variable term-list))
  (define (make-from-coeffs variable coeffs)
    (make-poly variable
               ((get 'make-from-coeffs 'dense-terms) coeffs)))
  (define (make-from-terms variable terms)
    (make-poly variable
               ((get 'make-from-terms 'sparse-terms) terms)))
  
  (define (variable p) (car p))
  (define (term-list p) (cdr p))
  
  ;; variable tests and selection
  (define (variable? x) (symbol? x))
  (define (same-variable? v1 v2)
    (and (variable? v1)
         (variable? v2)
         (or (eq? v1 v2)
             (eq? v1 'unbound)
             (eq? v2 'unbound))))

  (define (select-variable p1 p2)
    (let ((v1 (variable p1)))
      (if (eq? v1 'unbound)
          (variable p2)
          v1)))

  ;; procedures used by add-poly
  (define (add-poly p1 p2)
    (if (same-variable? (variable p1) (variable p2))
        (make-poly (select-variable p1 p2)
                   (add (term-list p1)
                        (term-list p2)))
        (error "Polys not in same var -- ADD-POLY"
               (list p1 p2))))

  ;; procedures used by mul-poly
  (define (mul-poly p1 p2)
    (if (same-variable? (variable p1) (variable p2))
        (make-poly (select-variable p1 p2)
                   (mul (term-list p1)
                        (term-list p2)))
        (error "Polys not in same var -- MUL-POLY"
               (list p1 p2))))

  ;; Subtraction
  (define (sub-poly p1 p2)
    (add-poly p1 (negate-poly p2)))
  
  ;; zero test
  (define (=zero-poly? p)
    (=zero? (term-list p)))
        
  ;; Negation
  (define (negate-poly p)
    (make-poly (variable p)
               (negate (term-list p))))
  
  ;; Equality
  (define (equ-poly? p1 p2)
    (and (same-variable? (variable p1) (variable p2))
         (equ? (term-list p1) (term-list p2))))
  
  ;; Coercion
  (define (polynomial->complex p)
    (let ((constant (get-constant (term-list p))))
      (if (is-lower? constant 'complex)
          (raise-to 'complex constant)
          constant)))

  ;; 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 'equ? '(polynomial polynomial) equ-poly?)
  (put '=zero? '(polynomial) =zero-poly?)
  (put 'negate '(polynomial)
       (lambda (p) (tag (negate-poly p))))
  (put 'sub '(polynomial polynomial)
       (lambda (p1 p2) (tag (sub-poly p1 p2))))
  (put 'make 'polynomial
       (lambda (var terms) (tag (make-poly var terms))))
  (put 'make-from-terms 'polynomial
       (lambda (variable terms) (tag (make-from-terms variable terms))))
  (put 'make-from-coeffs 'polynomial
       (lambda (variable coeffs) (tag (make-from-coeffs variable coeffs))))
  (put-coercion 'polynomial 'complex polynomial->complex)
  'done)

(define (make-polynomial-from-coeffs variable coeffs)
  ((get 'make-from-coeffs 'polynomial) variable coeffs))

(define (make-polynomial-from-terms variable terms)
  ((get 'make-from-terms 'polynomial) variable terms))

We have a special constructor for zero-order polynomials. As noted in the previous discussion, we're using a special marker, 'unbound, for zero-order polynomials that can participate in arithmetic operations with any other polynomial, regardless of their variable. This encapsulates the marker to support coercion from complex to polynomial types without complex having to know about 'unbound.

(define (make-zero-order-polynomial-from-coeff coeff)
  ((get 'make-from-coeffs 'polynomial) 'unbound (list coeff)))


;;;
;;; Install the Polynomial Packages
;;;
(install-term-package)
(install-sparse-terms-package)
(install-dense-terms-package)
(install-polynomial-package)

Okay, that's a big enough post for now... We'll take it for a spin in the next one!

SICP Exercise 2.90: Supporting Dense and Sparse Polynomials - Part 1

Suppose we want to have a polynomial system that is efficient for both sparse and dense polynomials. One way to do this is to allow both kinds of term-list representations in our system. The situation is analogous to the complex-number example of section 2.4, where we allowed both rectangular and polar representations. To do this we must distinguish different types of term lists and make the operations on term lists generic. Redesign the polynomial system to implement this generalization. This is a major effort, not a local change.

They're not kidding - it's not a local change. I've split the solution into two posts. In this post I'll discuss the changes required, and only give a few snippets of code. In the second post I'll present the code itself and talk about various interesting bits of it.

Using the Analogy

We're told that the situation is analogous to the complex-number example, so let's first examine the complex-number example to see how it handles its different representations. Then we can think about how we can apply that approach to this exercise.

In the complex-number example the rectangular and polar representations are extracted to their own packages, named rectangular and polar respectively. These two packages expose identical sets of primitive operations: constructors and accessors for the components of the complex number. They do not expose, or implement in any way, any complex operations such as arithmetic operations. Nor are data objects of these types directly constructed. Instead an overarching complex type is added to support these.

This complex package exposes two different construction operations that allow for the creation of complex number representations appropriate to the semantics of the passed parameters. I.e. rectangular representations are used for {real, complex} tuples and polar representations for {magnitude, angle} tuples. These delegate much of the construction of the representation to the rectangular and polar packages, simply tagging the results as 'complex types. It also provides representation-agnostic implementations of arithmetic operations that utilize the set of primitive operations exposed by the two representations and that select the best representation for the result based upon the operation.

The authors show how this fits into the type system in figure 2.23 at the start of section 2.5. Here it is again:

Programs that use numbers
				add sub mul div
Generic arithmetic package
	add-rat sub-rat mul-rat div-rat			add-complex sub-complex mul-complex div-complex				+ - * /
Rational arithmetic			Complex arithmetic				Ordinary arithmetic
			Rectangular		Polar

By analogy this suggest that we should identify a set of primitive operations that are common to both the sparse and dense polynomial representations that will allow for the implementation of the higher-level operations to be defined in terms of these operations only. Doing so will allow us to provide an over-arching polynomial package that can perform operations on either (or both!) of the two representations without bypassing this abstraction layer.

We can show how this will fit into the system by extending the figure above:

Programs that use numbers

add sub mul div

Generic arithmetic package

add-rat sub-rat mul-rat div-rat

add-complex sub-complex mul-complex div-complex

add-poly sub-poly mul-poly div-poly

+ - * /

Rational arithmetic

Complex arithmetic

Polynomial arithmetic

Ordinary arithmetic

Rectangular

Polar

Dense

Sparse

Identifying the Interface

If we look at the original polynomial package implementation (which you can find in its entirety at the top of my solution to exercise 2.87), then you can see that the authors have carefully split out the term-list-specific portions of the operations in the package. add-terms is split out from add-poly, mul-terms is split out from mul-poly and so on.

We've retained this split throughout the subsequent exercises. As a result, if we compare the sparse and dense polynomial packages we've produced then we can see that they differ only from the level of these operations down. This suggests that we should treat these operations as the primitive operations that we should extract. In other words, we'll use different representations for the term-lists, but the surrounding polynomial representation (i.e. a tuple consisting of a variable and a term-list) will remain unchanged.

If we do choose these as our primitive operations then we'll quickly hit a stumbling block: we can't add two polynomials together which have different term-list representations. The implementation of add-terms for sparse term-lists will only work with, and will only be installed for, pairs of sparse term-lists. Similarly the implementation of add-terms for dense term-lists will only work with, and will only be installed for, pairs of dense term-lists. So we'll get an error if we try to add a sparse and a dense term-list together.

Obviously we can't work around this by simply e.g. installing the dense term-list operations for all permutations of term-list representations. The internal representations are different and so errors will occur. Instead we'll need to coerce the types of parameters to some common type (probably sparse term-lists, as we'd not want to expand out a sparse term-list with, say, an order 1,000,000 term to a dense term-list with all the zero-padding that would involve!).

Coercing Dense and Sparse Term-Lists

In section 2.5.2 and the related exercises we developed a mechanism for automatic coercion of types based upon a tower-of-types. Surely all we need to do then is to add the sparse and dense term-list types into the tower-of-types and all will work, won't it?

Well, not quite.

In our solution to this exercise sparse and dense term-lists will not be "primary" types. I.e like the rectangular and polar representations of complex numbers, they are not directly constructable and do not expose the full set of arithmetic operations. Instead they will be "secondary" (or sub-) types of the polynomial type. As a result they do not fit into the tower-of-types we have been using so far. What we need is an entirely separate tower-of-types for term-list representations and have apply-generic support this.

Thankfully this isn't as difficult as it sounds. Our current implementation for raising and dropping types is dependant upon there being an unbroken chain of coercion operations that can raise a type from any point in the tower-of-types right the way up to the top of the tower. This requirement is made explicit in apply-raise, where it tries to obtain an appropriate coercion operation to raise the type one level in the tower. If no such operation can be found it raises an error:

…
(let ((raiser (get-coercion (type-tag x) (cadr types))))
  (if raiser
      (raiser (contents x))
      (error "No coercion procedure found for types"
             (list (type-tag x) (cadr types)))))
…

And as apply-raise is also used by project to lower types, (by passing in the reversed tower-of-types to apply-raise), the current implementation also depends upon their being a similar chain of coercion operations that can lower any type all the way to the bottom of the tower.

If we loosen this restriction and allow for breaks in the chain of coercion operations up and down the tower then this effectively allows multiple towers-of-types to be represented with a single tower. All we need is for the implementation to assume that a missing coercion operation indicates a point at which a type should not be lowered any further. To do this all we need to do is to return the value passed in instead of raising an error. Here's the updated apply-raise:

(define (apply-raise x types)
  (cond ((null? types)
         (error "Type not found in the tower-of-types"
                (list (type-tag x) tower-of-types)))
        ((eq? (type-tag x) (car types))
         (if (null? (cdr types))
             x
             (let ((raiser (get-coercion (type-tag x) (cadr types))))
               (if raiser
                   (raiser (contents x))
                   x))))
        (else (apply-raise x (cdr types)))))

Having made this change we can just add the sparse and dense term-list types to the tower, and add the appropriate coercion procedures, and it'll all work. We'll leave the coercion procedures for now, as they'll fall out of the work we do to create the term-lists.

Coercing Polynomials

While we're talking about coercions it's worth noting that we can also fit our overarching polynomial package into the tower-of-types too.

Let's start with raising a complex type to a polynomial. What we need to produce is a polynomial of order 0 whose coefficient of the order 0 term is the complex value we're raising. We need to consider, however, what variable we will create the polynomial with. If we choose e.g. 'x as our variable then we'll be unable to use it in any calculations with polynomials whose variable is 'y and vice versa - despite the fact that the polynomial has no actual terms in which the variable is important.

One way of addressing this is to introduce a special marker variable that is only used for polynomials of order 0 that indicates that it can be used in calculations with any other polynomial. Effectively it late-binds the variable at the point the calculation is performed. I chose to call this marker 'unbound. Supporting this first requires a simple extension to same-variable? so that it matches 'unbound with any other variable:

(define (same-variable? v1 v2)
  (and (variable? v1)
       (variable? v2)
       (or (eq? v1 v2)
           (eq? v1 'unbound)
           (eq? v2 'unbound))))

It also requires a modification to the arithmetic operations so that they select the appropriate variable for any polynomial resulting from a calculation. At the moment our arithmetic operations just use the variable of the first polynomial involved in any calculation. All we need to do is to modify these so that they identify whether or not one this variable is 'unbound and, if so, use the other polynomial's variable instead. This selection can be encapsulated in a separate operation as follows:

(define (select-variable p1 p2)
  (let ((v1 (variable p1)))
    (if (eq? v1 'unbound)
        (variable p2)
        v1)))

Once we can raise complex types to polynomial types in this manner, it's a trivial matter to drop them again. We know that the order 0 term of a polynomial is equivalent to a constant numerical value. As a result, we can successfully drop a polynomial of order 0 to a complex by simply taking the coefficient of the order 0 term. Now due to the manner in drop works (i.e. it projects down the tree, which is allowed to drop information, then raises the result and checks that this is the same), this means that our coercion procedure can follow this approach, regardless of the order of the polynomial. I.e. it takes the coefficient of the order 0 term, and discards the rest of the polynomial.

After we incorporate these changes we'll have a new tower-of-types that effectively looks like this:

    sparse-terms
         ↑
    dense-terms

     polynomial
         ↑
      complex
         ↑
      rational
         ↑
      integer

Creating Polynomials

I noted above that coercing between sparse and dense term-lists would fall out of the work on creating polynomials, so let's look at that next.

The implementations of the polynomial package that we've dealt with in the preceding exercises for sparse and dense polynomials have produced two ways in which we can construct polynomials. For sparse polynomials we create a polynomial using a list of terms, with each term being a tuple of the term's order and its coefficient. For dense polynomials we create a polynomial using a list of the term coefficients, with the order of the term corresponding to a particular coefficient being implicit in the position of the coefficient in that list.

Now let's suppose we implement our system so that we can create either term-list representation in either manner. I.e. we can create a term list with either a list of terms or a list of coefficients, regardless of whether we're using the sparse or dense representations. In such a system we would be able to easily convert (i.e. coerce) a sparse term-list into a dense term-list by simply creating a dense term-list using the sparse term-list's terms. Similarly, we would be able to easily convert a dense term-list into a sparse term-list by creating a sparse term-list using the dense term-list's coefficients.

So, as you can see, by implementing our system in this manner the coercion operations will fall out of the work on creating polynomials.

The creation operations themselves are pretty straightforward. We already have appropriate operations for creating sparse polynomials from term lists and dense polynomials from coefficient lists. In order to create a sparse polynomial from a coefficient list we simply run through the list and add a term of the appropriate order for each coefficient (i.e. the term order for a coefficient is dictated by its position in the list). To create a dense polynomial from a term list we perform a similar operation to the one of building a sparse polynomial from a term list, except that the list we build consists of coefficients only, and so we determine whether or not we're inserting at the right location based upon the length of the list.

Dense or Sparse?

One side-effect of supporting coercion between dense and sparse term-list representations is that, thanks to the work we did in exercise 2.85, apply-generic will automatically drop a sparse term-list to a dense term-list. Now we don't necessarily want to do this, as the result of an arithmetic operation on two sparse term-lists may itself be a term-list that's best represented sparsely. We'll want to block the coercion in this case. One way of doing this is to simply have the coercion procedure decide whether or not it should coerce and, if not, just return the term-list unchanged (well, apart from re-tagging it). Our implementation of drop already checks for coercion procedures doing this and uses it to indicate that the value could not be dropped any further.

To make things a bit more complex, I'd like to generalize my statement above... The result of an arithmetic operation on any pair of term-lists, regardless of their representations, may be a term-list that's best represented sparsely. We should also make sure that the result of dense term-list arithmetic operations is provided in the best representation. The raising won't happen automatically for us here, so we'll need to do it ourselves.

So when should we use a dense term-list and when a sparse term-list? The book defines the difference as follows: "A polynomial is said to be dense if it has nonzero coefficients in terms of most orders. If it has many zero terms it is said to be sparse." There's no hard and fast rule about what having "nonzero coefficients in terms of most orders" or "many zero terms" mean, so I opted for the following: we'll represent a term-list sparsely when more than 10% of its terms have a zero coefficient, with a special rule to deal with term-lists whose highest order term is order 10 or lower.

Right, I think that's enough discussion of it. Let's move on to the code!