2012-05-31

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!

2 comments:

  1. FYI, made a slight change to the way the ensure-valid-term-list operations are used within the sparse and dense term-lists. This is now called by tag prior to attaching the tag. As everything is tagged in one way or another this ensures that the result of all term-list operations will return a non-empty term-list.

    ReplyDelete
  2. FYI, reverted the attach-tag, type-tag and contents procedures to their pre-exercise 2.78 (http://jots-jottings.blogspot.com/2012/02/sicp-exercise-278-using-primitive-types.html) state. The 'scheme-number type is no longer supported, so primitive numerical type support is not necessary (and gets in the way of exercise 2.93).

    ReplyDelete