#lang racket

;; This is the complete symbolic algebra program from section 2.5.3, including
;; polynomials in more than one variable, a hierarchy of variables and coercion,
;; and rational functions, as outlined in the exercises.

;;;;;;;;;;;;;;;;;;;;;;
;; GLOBAL VARIABLES ;;
;;;;;;;;;;;;;;;;;;;;;;

;; Specify the output representation for polynomial term lists here.
(define TERMLIST-OUTPUT-TYPE 'dense) ; 'dense or 'sparse

;; Only allow these variables in polynomials.
(define VARIABLE-LIST '(x y z))
(define HIGHEST-VARIABLE (first VARIABLE-LIST))
(define LOWEST-VARIABLE (last VARIABLE-LIST))

;; Only simplify the results of these operations.
(define SIMPLIFY-LIST '(add sub mul div expo
                          sqroot square
                          arctan cosine sine
                          re-part im-part mag-part ang-part
                          first-term constant-term
                          ))
(define (can-simplify? op) (member op SIMPLIFY-LIST))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; PROCEDURE TABLES AND VARIABLE DICTIONARIES ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; procedure table for generic operations
(define op-table (make-hash))
(define (put op tag-list proc)
  (hash-set! op-table (cons op tag-list) proc))
(define (can-apply? op tag-list)
  (member (cons op tag-list) (hash-keys op-table)))
(define (get op tag-list) (hash-ref op-table (cons op tag-list)))

;; procedure table for type coercion operations
(define coercion-table (make-hash))
(define (put-project type proc)
  (hash-set! coercion-table (cons 'project type) proc))
(define (can-project? arg)
  (member (cons 'project (type-tag arg)) (hash-keys coercion-table)))
(define (project arg)
  ((hash-ref coercion-table (cons 'project (type-tag arg))) arg))
;; raising the type
(define (put-raise type proc)
  (hash-set! coercion-table (cons 'raise type) proc))
(define (can-raise? arg)
  (member (cons 'raise (type-tag arg)) (hash-keys coercion-table)))
(define (raise arg)
  ((hash-ref coercion-table (cons 'raise (type-tag arg))) arg))

;; variable dictionaries
(define next-higher-variable (make-hash))
(define next-lower-variable (make-hash))
(define (install-variable-dictionaries)
  (define (loop var)
    (cond [(eq? var LOWEST-VARIABLE) 'done]
          [else
            (define next-lower (second (member var VARIABLE-LIST)))
            (hash-set! next-lower-variable var next-lower)
            (hash-set! next-higher-variable next-lower var)
            (loop next-lower)]))
  (loop HIGHEST-VARIABLE))
(install-variable-dictionaries)

(define (next-higher var) (hash-ref next-higher-variable var))
(define (next-lower var) (hash-ref next-lower-variable var))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; APPLY-GENERIC AND HELPER FUNCTIONS ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; procedures for handling type tags
(define (attach-tag type-tag datum)
  (cond [(exact-integer? datum) datum]
        [(inexact-real? datum) datum]
        [else (cons type-tag datum)]))
(define (type-tag datum)
  (cond [(exact-integer? datum) 'integer]
        [(inexact-real? datum) 'real]
        [(pair? datum) (first datum)]
        [else (error "Bad tagged datum -- TYPE-TAG:" datum)]))
(define (contents datum)
  (cond [(exact-integer? datum) datum]
        [(inexact-real? datum) datum]
        [(pair? datum) (cdr datum)]
        [else (error "Bad tagged datum -- CONTENTS:" datum)]))

;; type coercion and simplification
(define (simplify arg) 
  ; Simplify arg to the lowest possible correct type.
  (cond [(can-project? arg)
         (define projected-arg (project arg))
         (if (equ? arg (raise projected-arg))
           (simplify projected-arg)
           arg)]
        [else arg]))
(define (highest-type tag-list)
  ; Return the highest type in tag-list.
  (define (loop type arg)
    (cond [(can-raise? arg)
           (define raised-arg (raise arg))
           (define raised-type (type-tag raised-arg))
           (if (member raised-type tag-list)
             (loop raised-type raised-arg)
             (loop type raised-arg))]
          [else type]))
  (loop 'integer 1))

(define (coerce arg target-type)
  ; Raise arg until it has type target-type. Assume type of arg <= target-type.
  (if (equal? (type-tag arg) target-type)
    arg
    (coerce (raise arg) target-type)))
(define (coerce-all args target-type)
  ; Coerce all arguments to the target type.
  (map (lambda (arg) (coerce arg target-type)) args))

;; the generic application procedure
(define (apply-generic op . args)
  ; Dispatch by type.
  (define (simplify-result datum)
    (if (can-simplify? op)
      (simplify datum)
      datum))
  (define type-tags (map type-tag args))
  (cond [(can-apply? op type-tags)
         (simplify-result (apply (get op type-tags) (map contents args)))]
        [(andmap (lambda (type) (equal? type (first type-tags))) type-tags)
         (error "No method available -- APPLY-GENERIC:" op type-tags)]
        [else
          (define coerced-args (coerce-all args (highest-type type-tags)))
          (simplify-result (apply (get op (map type-tag coerced-args))
                                  (map contents coerced-args)))]))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; THE GENERIC OPERATIONS ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; constructors and selectors
(define (make-integer a) ((get 'make 'integer) a))
(define (make-rational-number n d) ((get 'make 'rational-number) n d))
(define (make-rational-function n d) ((get 'make 'rational-function) n d))
(define (make-real x) ((get 'make 'real) x))
;; complex numbers
(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))
(define (re-part z) (apply-generic 're-part z))
(define (im-part z) (apply-generic 'im-part z))
(define (mag-part z) (apply-generic 'mag-part z))
(define (ang-part z) (apply-generic 'ang-part z))
;; polynomials
(define (make-poly-from-dense-termlist var term-list)
  ((get 'make-from-dense-termlist `(polynomial ,var)) term-list))
(define (make-poly-from-sparse-termlist var term-list)
  ((get 'make-from-sparse-termlist `(polynomial ,var)) term-list))

;; generic operations
(define (=zero? x) (apply-generic '=zero? x))
(define (absolute x) (apply-generic 'absolute x))
(define (add x y) (apply-generic 'add x y))
(define (arctan y x) (apply-generic 'arctan y x))
(define (constant-term p) (apply-generic 'constant-term p))
(define (cosine x) (apply-generic 'cosine x))
(define (div x y) (apply-generic 'div x y))
(define (equ? x y) (apply-generic 'equ? x y))
(define (expo x y) (apply-generic 'expo x y))
(define (first-term p) (apply-generic 'first-term p))
(define (greatest-common-divisor x y)
  (apply-generic 'greatest-common-divisor x y))
(define (mul x y) (apply-generic 'mul x y))
(define (polynomial-division p1 p2)
  (apply-generic 'polynomial-division p1 p2))
(define (reduce x y) (apply-generic 'reduce x y))
(define (sine x) (apply-generic 'sine x))
(define (sqroot x) (apply-generic 'sqroot x))
(define (square x) (apply-generic 'square x))
(define (sub x y) (apply-generic 'sub x y))

;;;;;;;;;;;;;;;
;; INTEGERS  ;;
;;;;;;;;;;;;;;;

(define (install-integer-package)
  ;; internal procedures
  (define (make-integer a)
    (define int (floor a))
    (if (exact? int)
      int
      (inexact->exact int)))
  (define (=zero-integer? a) (= a 0))
  (define (gcd-integer a b)
    (if (= b 0) a (gcd-integer b (remainder a b))))
  (define (reduce-integer n d)
    (define g (gcd-integer n d))
    (list (/ n g) (/ d g)))
  (define arctan (compose make-real atan))
  (define cosine (compose make-real cos))
  (define sine (compose make-real sin))
  (define sqroot (compose make-real sqrt))
  ;; coercion procedures
  (define (raise a) (make-rational-number a 1))
  (put-raise 'integer raise)
  ;; interface to the rest of the system
  (put 'make 'integer make-integer)
  (put '=zero? '(integer) =zero-integer?)
  (put 'absolute '(integer) abs)
  (put 'add '(integer integer) +)
  (put 'arctan '(integer integer) arctan)
  (put 'cosine '(integer) cosine)
  (put 'div '(integer integer) make-rational-number)
  (put 'equ? '(integer integer) =)
  (put 'expo '(integer integer) expt)
  (put 'greatest-common-divisor '(integer integer) gcd-integer)
  (put 'mul '(integer integer) *)
  (put 'reduce '(integer integer) reduce-integer)
  (put 'sine '(integer) sine)
  (put 'sqroot '(integer) sqroot)
  (put 'square '(integer) sqr)
  (put 'sub '(integer integer) -)
  )
(install-integer-package)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; RATIONAL NUMBERS AND RATIONAL FUNCTIONS ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(define (install-rational-package)

  ;; internal procedures
  (define (make-rat n d) (reduce n d))
  (define (numer x) (first x))
  (define (denom x) (second x))
  (define (=zero-rat? x) (=zero? (numer x)))
  (define (equ-rat? x y) (equ? (mul (numer x) (denom y))
                               (mul (denom x) (numer y))))
  (define (add-rat x y)
    (make-rat (add (mul (numer x) (denom y))
                   (mul (numer y) (denom x)))
              (mul (denom x) (denom y))))
  (define (sub-rat x y)
    (make-rat (sub (mul (numer x) (denom y))
                   (mul (numer y) (denom x)))
              (mul (denom x) (denom y))))
  (define (mul-rat x y)
    (make-rat (mul (numer x) (numer y))
              (mul (denom x) (denom y))))
  (define (div-rat x y)
    (make-rat (mul (numer x) (denom y))
              (mul (denom x) (numer y))))

  ;; procedures for rational numbers only
  (define (absolute x) (make-real (abs (/ (numer x) (denom x)))))
  (define (arctan y x) (make-real (atan (/ (numer y) (denom y)))))
  (define (cosine x) (make-real (cos (/ (numer x) (denom x)))))
  (define (expo-ratnum x y) (make-real (expt (/ (numer x) (denom x))
                                             (/ (numer y) (denom y)))))
  (define (sine x) (make-real (sin (/ (numer x) (denom x)))))
  (define (square x) (make-rat (sqr (numer x)) (sqr (denom x))))
  (define (sqroot x) (make-real (sqrt (/ (numer x) (denom x)))))

  ;; procedures for rational functions only
  (define (expo-ratfunc x n)
    (cond [(and (integer? n) (positive? n))
           (if (= n 1)
             x
             (mul-rat x (expo-ratfunc x (sub1 n))))]
          [else
            (error "Exponent must be a positive integer -- EXPO-RATFUNC:" x n)]))

  ;; coercion procedures for rational numbers
  (define (project-rat x) (make-integer (floor (/ (numer x) (denom x)))))
  (define (raise-rat x) (make-real (/ (numer x) (denom x))))
  (put-project 'rational-number (compose project-rat contents))
  (put-raise 'rational-number (compose raise-rat contents))

  ;; interface to the rest of the system for rational numbers
  (define (tag-as-number x) (attach-tag 'rational-number x))
  (put 'make 'rational-number
       (compose tag-as-number make-rat))
  (put '=zero? '(rational-number) =zero-rat?)
  (put 'absolute '(rational-number) absolute)
  (put 'add '(rational-number rational-number)
       (compose tag-as-number add-rat))
  (put 'arctan '(rational-number rational-number) arctan)
  (put 'cosine '(rational-number) cosine)
  (put 'div '(rational-number rational-number)
       (compose tag-as-number div-rat))
  (put 'equ? '(rational-number rational-number) equ-rat?)
  (put 'expo '(rational-number rational-number) expo-ratnum)
  (put 'mul '(rational-number rational-number)
       (compose tag-as-number mul-rat))
  (put 'sine '(rational-number) sine)
  (put 'sqroot '(rational-number) sqroot)
  (put 'square '(rational-number) (compose tag-as-number square))
  (put 'sub '(rational-number rational-number)
       (compose tag-as-number sub-rat))

  ;; interface to the rest of the system for rational functions
  (define (tag-as-function x) (attach-tag 'rational-function x))
  (put 'make 'rational-function
       (compose tag-as-function make-rat))
  (put '=zero? '(rational-function) =zero-rat?)
  (put 'add '(rational-function rational-function)
       (compose tag-as-function add-rat))
  (put 'div '(rational-function rational-function)
       (compose tag-as-function div-rat))
  (put 'equ? '(rational-function rational-function) equ-rat?)
  (put 'expo '(rational-function integer)
       (compose tag-as-function expo-ratfunc))
  (put 'mul '(rational-function rational-function)
       (compose tag-as-function mul-rat))
  (put 'sub '(rational-function rational-function)
       (compose tag-as-function sub-rat))
  )
(install-rational-package)

;;;;;;;;;;;;;;;;;;
;; REAL NUMBERS ;;
;;;;;;;;;;;;;;;;;;

(define (install-real-package)
  ;; internal procedures
  (define (make-real x) ; x is a racket number or has lower type
    (cond [(pair? x) (raise x)]
          [(exact? x) (exact->inexact x)]
          [else x]))
  ; Account for floating-point error.
  (define (=zero-real? x) (< (abs x) 0.000000001))
  (define (equ-real? x y) (< (abs (- x y)) 0.000000001))
  ;; coercion procedures
  (define (project-real x)
    (make-rational-number (make-integer (floor (* x 1000000000))) 1000000000))
  (define (raise-real x) (make-complex-from-real-imag x 0))
  (put-project 'real project-real)
  (put-raise 'real raise-real)
  ;; interface to the rest of the system
  (put 'make 'real make-real)
  (put '=zero? '(real) =zero-real?)
  (put 'absolute '(real) abs)
  (put 'add '(real real) +)
  (put 'arctan '(real real) atan)
  (put 'cosine '(real) cos)
  (put 'div '(real real) /)
  (put 'equ? '(real real) =)
  (put 'expo '(real real) expt)
  (put 'mul '(real real) *)
  (put 'sine '(real) sin)
  (put 'sqroot '(real) sqrt)
  (put 'square '(real) sqr)
  (put 'sub '(real real) -)
  )
(install-real-package)

;;;;;;;;;;;;;;;;;;;;;
;; COMPLEX NUMBERS ;;
;;;;;;;;;;;;;;;;;;;;;

(define (install-complex-package)

  ;; rectangular coordinates
  (define (install-rectangular-package)
    ;; internal procedures
    (define (make-from-real-imag x y) (list x y))
    (define (re-part z) (first z))
    (define (im-part z) (second z))
    (define (make-from-mag-ang r a)
      (list (mul r (cosine a))
            (mul r (sine a))))
    (define (mag-part z)
      (sqroot (add (square (re-part z))
                   (square (im-part z)))))
    (define (ang-part z) (arctan (im-part z) (re-part z)))
    ;; interface to the rest of the system
    (define (tag z) (attach-tag 'rectangular z))
    (put 'make-from-mag-ang 'rectangular (compose tag make-from-mag-ang))
    (put 'make-from-real-imag 'rectangular (compose tag make-from-real-imag))
    (put 're-part '(rectangular) re-part)
    (put 'im-part '(rectangular) im-part)
    (put 'mag-part '(rectangular) mag-part)
    (put 'ang-part '(rectangular) ang-part)
    )
  (install-rectangular-package)

  ;; polar coordinates
  (define (install-polar-package)
    ;; internal procedures
    (define (make-from-mag-ang r a) (list r a))
    (define (mag-part z) (first z))
    (define (ang-part z) (second z))
    (define (make-from-real-imag x y)
      (list (sqroot (add (square x) (square y)))
            (arctan y x)))
    (define (re-part z) (mul (mag-part z) (cosine (ang-part z))))
    (define (im-part z) (mul (mag-part z) (sine (ang-part z))))
    ;; interface to the rest of the system
    (define (tag x) (attach-tag 'polar x))
    (put 'make-from-mag-ang 'polar (compose tag make-from-mag-ang))
    (put 'make-from-real-imag 'polar (compose tag make-from-real-imag))
    (put 're-part '(polar) re-part)
    (put 'im-part '(polar) im-part)
    (put 'mag-part '(polar) mag-part)
    (put 'ang-part '(polar) ang-part)
    )
  (install-polar-package)

  ;; imported procedures from rectangular and polar subpackages
  (define make-from-real-imag (get 'make-from-real-imag 'rectangular))
  (define make-from-mag-ang (get 'make-from-mag-ang 'polar))

  ;; internal complex number procedures
  (define (=zero-complex? z) ; account for floating-point error
    (and (< (absolute (re-part z)) 0.000000001)
         (< (absolute (im-part z)) 0.000000001)))
  (define (equ? z1 z2) ; account for floating-point error
    (and (< (absolute (sub (re-part z1) (re-part z2))) 0.000000001)
         (< (absolute (sub (im-part z1) (im-part z2))) 0.000000001)))
  (define (add-complex z1 z2)
    (make-from-real-imag (add (re-part z1) (re-part z2))
                         (add (im-part z1) (im-part z2))))
  (define (sub-complex z1 z2)
    (make-from-real-imag (sub (re-part z1) (re-part z2))
                         (sub (im-part z1) (im-part z2))))
  (define (mul-complex z1 z2)
    (make-from-mag-ang (mul (mag-part z1) (mag-part z2))
                       (add (ang-part z1) (ang-part z2))))
  (define (div-complex z1 z2)
    (make-from-mag-ang (div (mag-part z1) (mag-part z2))
                       (sub (ang-part z1) (ang-part z2))))
  (define (expo-complex z1 z2)
    (define w (expt (+ (make-real (re-part z1))
                       (* (make-real (im-part z1)) 0+i))
                    (+ (make-real (re-part z2))
                       (* (make-real (im-part z2)) 0+i))))
    (make-from-real-imag (real-part w) (imag-part w)))
  ;; coercion procedures
  (define (project-complex z) (make-real (re-part z)))
  (define (raise-complex z)
    (make-poly-from-dense-termlist LOWEST-VARIABLE (list (simplify z))))
  (put-project 'complex project-complex)
  (put-raise 'complex raise-complex)
  ;; interface to the rest of the system
  (define (tag z) (attach-tag 'complex z))
  (put 'make-from-mag-ang 'complex
       (compose tag make-from-mag-ang))
  (put 'make-from-real-imag 'complex
       (compose tag make-from-real-imag))
  (put '=zero? '(complex) =zero-complex?)
  (put 'add '(complex complex) (compose tag add-complex))
  (put 'ang-part '(complex) ang-part)
  (put 'div '(complex complex) (compose tag div-complex))
  (put 'equ? '(complex complex) equ?)
  (put 'expo '(complex complex) (compose tag expo-complex))
  (put 'im-part '(complex) im-part)
  (put 'mag-part '(complex) mag-part)
  (put 'mul '(complex complex) (compose tag mul-complex))
  (put 're-part '(complex) re-part)
  (put 'sub '(complex complex) (compose tag sub-complex))
  )
(install-complex-package)

;;;;;;;;;;;;;;;;;
;; POLYNOMIALS ;;
;;;;;;;;;;;;;;;;;

(define (install-polynomial-package)

  ;; term list representations

  ;; sparse term lists
  (define (install-sparse-termlist-package)
    ;; constructors and selectors
    (define (make-term order coeff) (list order coeff))
    (define (order term) (first term))
    (define (coeff term) (second term))
    (define (the-empty-termlist) empty)
    (define empty-termlist? null?)
    (define (first-term L)
      (if (empty-termlist? L)
        (make-term 0 0)
        (first L)))
    (define (rest-terms L) (cdr L))
    (define (dense->sparse DL)
      (cond [(null? DL) (the-empty-termlist)]
            [(=zero? (first DL)) (dense->sparse (cdr DL))]
            [else (cons (make-term (- (length DL) 1)
                                   (first DL))
                        (dense->sparse (cdr DL)))]))
    ;; procedures for generic operations
    (define (adjoin-term term L)
      (if (=zero? (coeff term))
        L
        (cons term L)))
    (define (constant-term L)
      (cond [(empty-termlist? L)
             (make-term 0 0)]
            [else
              (define last-term (last L))
              (if (zero? (order last-term))
                last-term
                (make-term 0 0))]))
    (define (negate-terms L)
      (map (lambda (term)
             (make-term (order term)
                        (mul -1 (coeff term))))
           L))
    (define (reduce-coefficients L) ; divide all coeffs by their gcd
      (define (list-of-coefficients term-list)
        ; Return a list of coefficients in no particular order.
        (foldr (lambda (term a-list)
                 (cons (second term) a-list))
               empty
               term-list))
      (define gcf (foldr greatest-common-divisor
                         0
                         (list-of-coefficients L)))
      (map (lambda (term)
             (make-term (order term)
                        (div (coeff term) gcf)))
           L))
    ;; interface to the rest of the system
    (define (tag-as-sparse x) (attach-tag 'sparse x))
    (define (tag-as-term x) (attach-tag 'term x))
    (put 'make-from-dense-termlist 'sparse
         (compose tag-as-sparse dense->sparse))
    (put 'make-from-sparse-termlist 'sparse tag-as-sparse)
    (put 'adjoin-term '(term sparse)
         (compose tag-as-sparse adjoin-term))
    (put 'constant-term '(sparse)
         (compose tag-as-term constant-term))
    (put 'first-term '(sparse)
         (compose tag-as-term first-term))
    (put 'negate-terms '(sparse)
         (compose tag-as-sparse negate-terms))
    (put 'reduce-coefficients '(sparse)
         (compose tag-as-sparse reduce-coefficients))
    (put 'rest-terms '(sparse)
         (compose tag-as-sparse rest-terms))
    )
  (install-sparse-termlist-package)

  ;; dense term lists
  (define (install-dense-termlist-package)
    ;; internal procedures
    (define (make-term order coeff) (list order coeff))
    (define (order term) (first term))
    (define (coeff term) (second term))
    (define (the-empty-termlist) empty)
    (define empty-termlist? null?)
    (define (first-term L)
      (if (empty-termlist? L)
        (make-term 0 0)
        (make-term (- (length L) 1)
                   (first L))))
    (define (rest-terms L) (cdr L))
    (define (sparse->dense SL)
      (define (list-of-zeroes n) (for/list ([i (in-range n)]) 0))
      (cond [(empty? SL)
             (the-empty-termlist)]
            [else
              (define first-term (first SL))
              (cond [(=zero? (coeff first-term))
                     (sparse->dense (rest SL))]
                    [(empty? (rest SL)) ; only one term in SL
                     (cons (coeff first-term)
                           (list-of-zeroes (order first-term)))]
                    [else
                      (define second-term (second SL))
                      (append (cons (coeff first-term)
                                    (list-of-zeroes (- (order first-term)
                                                       (order second-term)
                                                       1)))
                              (sparse->dense (rest SL)))])]))
    (define (adjoin-term term L) ; order >= length L
      (define (loop term-list)
        (if (> (order term) (length term-list))
          (loop (cons 0 term-list))
          (cons (coeff term) term-list)))
      (if (=zero? (coeff term))
        L
        (loop L)))
    (define (constant-term L)
      (if (empty-termlist? L)
        (make-term 0 0)
        (make-term 0 (last L))))
    (define (negate-terms L)
      (map (lambda (coefficient)
             (mul -1 coefficient))
           L))
    (define (reduce-coefficients L) ; divide all coeffs by their gcd
      (define gcf (foldr greatest-common-divisor 0 L))
      (map (lambda (coeff)
             (div coeff gcf))
           L))
    ;; interface to the rest of the system
    (define (tag-as-dense x) (attach-tag 'dense x))
    (define (tag-as-term x) (attach-tag 'term x))
    (put 'make-from-dense-termlist 'dense tag-as-dense)
    (put 'make-from-sparse-termlist 'dense
         (compose tag-as-dense sparse->dense))
    (put 'adjoin-term '(term dense)
         (compose tag-as-dense adjoin-term))
    (put 'constant-term '(dense)
         (compose tag-as-term constant-term))
    (put 'first-term '(dense)
         (compose tag-as-term first-term))
    (put 'negate-terms '(dense)
         (compose tag-as-dense negate-terms))
    (put 'reduce-coefficients '(dense)
         (compose tag-as-dense reduce-coefficients))
    (put 'rest-terms '(dense)
         (compose tag-as-dense rest-terms))
    )
  (install-dense-termlist-package)

  ;; imported procedures from the term list subpackages
  (define make-from-dense-termlist
    (get 'make-from-dense-termlist TERMLIST-OUTPUT-TYPE))
  (define make-from-sparse-termlist
    (get 'make-from-sparse-termlist TERMLIST-OUTPUT-TYPE))

  ;; term and term list operations
  (define (make-term order coeff) (list 'term order coeff))
  (define (order term) (second term))
  (define (coeff term) (third term))
  (define (rest-terms L) (apply-generic 'rest-terms L))
  (define (the-empty-termlist) (list TERMLIST-OUTPUT-TYPE))
  (define (empty-termlist? L) (null? (contents L)))
  (define (=zero-termlist? L)
    (or (empty-termlist? L)
        (and (=zero? (coeff (first-term L)))
             (=zero-termlist? (rest-terms L)))))
  (define (reduce-coefficients L)
    (apply-generic 'reduce-coefficients L))
  (define (adjoin-term term L) (apply-generic 'adjoin-term term L))
  (define (negate-terms L) (apply-generic 'negate-terms L))
  (define (add-terms L1 L2)
    (cond [(empty-termlist? L1) L2]
          [(empty-termlist? L2) L1]
          [else
            (define t1 (first-term L1))
            (define 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)))])]))
  (define (sub-terms L1 L2) (add-terms L1 (negate-terms L2)))
  (define (equ-termlist? L1 L2)
    (=zero-termlist? (sub-terms L1 L2)))
  (define (mul-term-by-all-terms t1 L)
    (cond [(empty-termlist? L)
           (the-empty-termlist)]
          [else
            (define 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)))]))
  (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 (div-terms L1 L2)
    (cond [(empty-termlist? L1)
           (list (the-empty-termlist) (the-empty-termlist))]
          [else
            (define t1 (first-term L1))
            (define t2 (first-term L2))
            (cond [(> (order t2) (order t1))
                   (list (the-empty-termlist) L1)]
                  [else
                    (define divisor-term
                      (make-term (- (order t1) (order t2))
                                 (div (coeff t1) (coeff t2))))
                    (define rest-of-result
                      (div-terms (sub-terms L1
                                            (mul-term-by-all-terms divisor-term
                                                                    L2))
                                 L2))
                    (list (adjoin-term divisor-term (first rest-of-result))
                          (second rest-of-result))])]))
  (define (divisor-termlist L1 L2) (first (div-terms L1 L2)))
  (define (remainder-termlist L1 L2) (second (div-terms L1 L2)))
  (define (expo-terms L n)
    (cond [(positive? n)
           (if (= n 1)
             L
             (mul-terms L (expo-terms L (sub1 n))))]
          [else
            (error "Exponent must be a positive integer -- EXPO-TERMS:" L n)]))
  (define (gcd-terms a b)
    (if (empty-termlist? b)
      a
      (gcd-terms b (remainder-termlist a b))))
  (define (pseudo-remainder-termlist L1 L2)
    ; Clear denominators from remainder-termlist.
    (remainder-termlist (mul-term-by-all-terms
                          (make-term 0 (expt (coeff (first-term L2))
                                             (+ 1
                                                (order (first-term L1))
                                                (- (order (first-term L2))))))
                          L1)
                        L2))
  (define (pseudo-gcd-terms a b)
    ; Return a gcd term list with integer coefficients.
    (if (empty-termlist? b)
      a
      (pseudo-gcd-terms b (pseudo-remainder-termlist a b))))
  (define (reduce-terms a b)
    (define (same-sign? L1 L2)
      (positive? (make-real
                   (re-part
                     (coerce (mul (coeff (first-term L1))
                                  (coeff (first-term L2)))
                             'complex)))))
    (define (fix-signs L1 L2)
      ; Ensure the leading coefficient of L1/L2 is the same as that of a/b.
      (cond [(and (same-sign? a L1) (not (same-sign? b L2)))
             (list L1 (negate-terms L2))]
            [(and (not (same-sign? a L1)) (same-sign? b L2))
             (list (negate-terms L1) L2)]
            [(and (not (same-sign? a L1)) (not (same-sign? b L2)))
             (list (negate-terms L1) (negate-terms L2))]
            [else (list L1 L2)]))
    (define pseudo-gcf (pseudo-gcd-terms a b))
    (define (integerizing-factor L1 L2)
      (expt (coeff (first-term pseudo-gcf))
            (+ 1 (max (order (first-term L1)) (order (first-term L2)))
               (- (order (first-term pseudo-gcf))))))
    (define (clear-denominators L)
      (mul-term-by-all-terms (make-term 0 (integerizing-factor a b))
                             L))
    (apply fix-signs
           (map (lambda (term-list) (reduce-coefficients
                                      (divisor-termlist
                                        (clear-denominators term-list)
                                        pseudo-gcf)))
                (list a b))))
  ;; polynomial operations
  (define (make-poly var term-list)
    (attach-tag `(polynomial ,var) term-list))
  (define (variable poly) (second (first poly)))
  (define (term-list poly) (rest poly))

  ;; coercion operations
  ; Coercion procedures act on tagged polynomials.
  (define (project-poly-to-complex p)
    (coerce (coeff (constant-term (term-list p))) 'complex))
  (define (project-poly-to-poly p)
    (coerce (coeff (constant-term (term-list p)))
            `(polynomial ,(next-lower (variable p)))))
  (define (raise-poly p)
    (make-poly-from-dense-termlist (next-higher (variable p)) (list p)))
  (define (install-coercion-procedures)
    (define (loop var) ; assume there is more than one variable
      (cond [(eq? var HIGHEST-VARIABLE)       ; start of the loop
             (put-project `(polynomial ,var) project-poly-to-poly)
             (loop (next-lower var))]
            [(eq? var LOWEST-VARIABLE)        ; end of the loop
             (put-project `(polynomial ,var) project-poly-to-complex)
             (put-raise `(polynomial ,var) raise-poly)
             'done]
            [else                               ; middle of the loop
              (put-project `(polynomial ,var) project-poly-to-poly)
              (put-raise `(polynomial ,var) raise-poly)
              (loop (next-lower var))]))
    (if (eq? HIGHEST-VARIABLE LOWEST-VARIABLE) ; only one variable
      (put-project `(polynomial ,LOWEST-VARIABLE) project-poly-to-complex)
      (loop HIGHEST-VARIABLE)))
  (install-coercion-procedures)

  ;; interface to the rest of the system
  ; Generic procedures act on the contents of polynomials, ie term lists.
  (define (install-generic-procedures)
    (define (loop var)
      (put 'make-from-dense-termlist `(polynomial ,var)
           (compose (curry make-poly var) make-from-dense-termlist))
      (put 'make-from-sparse-termlist `(polynomial ,var)
           (compose (curry make-poly var) make-from-sparse-termlist))
      (put '=zero? `((polynomial ,var)) =zero-termlist?)
      (put 'equ? `((polynomial ,var) (polynomial ,var)) equ-termlist?)
      (put 'add `((polynomial ,var) (polynomial ,var))
           (compose (curry make-poly var) add-terms))
      (put 'sub `((polynomial ,var) (polynomial ,var))
           (compose (curry make-poly var) sub-terms))
      (put 'mul `((polynomial ,var) (polynomial ,var))
           (compose (curry make-poly var) mul-terms))
      (put 'expo `((polynomial ,var) integer)
           (compose (curry make-poly var) expo-terms))
      (put 'first-term `((polynomial ,var)) first-term)
      (put 'constant-term `((polynomial ,var)) constant-term)
      (if (eq? var LOWEST-VARIABLE)
        (begin
          (put 'polynomial-division
               `((polynomial ,LOWEST-VARIABLE)
                 (polynomial ,LOWEST-VARIABLE))
               (lambda (L1 L2)
                 (map (curry make-poly LOWEST-VARIABLE)
                      (div-terms L1 L2))))
          (put 'greatest-common-divisor
               `((polynomial ,LOWEST-VARIABLE)
                 (polynomial ,LOWEST-VARIABLE))
               ; return a gcd with reduced integer coefficients
               (compose (curry make-poly LOWEST-VARIABLE)
                        reduce-coefficients
                        pseudo-gcd-terms))
          (put 'reduce
               `((polynomial ,LOWEST-VARIABLE)
                 (polynomial ,LOWEST-VARIABLE))
               (lambda (L1 L2) (map (curry make-poly LOWEST-VARIABLE)
                                    (reduce-terms L1 L2))))
          )
        (loop (next-lower var))))
    (loop HIGHEST-VARIABLE))
  (install-generic-procedures)
  )
(install-polynomial-package)

;;;;;;;;;;;
;; TESTS ;;
;;;;;;;;;;;

;; NUMBERS

(=zero? (sub (make-complex-from-real-imag (make-rational-number 1 2)
                                          0.75)
             (make-complex-from-real-imag 0.5
                                          (make-rational-number 3 4))))
;; #t
(equ? (expo (make-complex-from-mag-ang 1 1) pi)
      -1) ; e ^ (i * pi) = -1
;; #t
(mul 0.25 (make-complex-from-real-imag 4 0))
;; 1

;; POLYNOMIALS IN DIFFERENT VARIABLES

(define p (make-poly-from-dense-termlist 'x '(1 2 3)))
(define q (make-poly-from-dense-termlist 'y '(4 5 6)))
(add p q)
;; '((polynomial x) dense 1 2 ((polynomial y) dense 4 5 9))
(sub p q)
;; '((polynomial x) dense 1 2 ((polynomial y) dense -4 -5 -3))
(mul p q)
;; '((polynomial x)
  ;; dense
  ;; ((polynomial y) dense 4 5 6)
  ;; ((polynomial y) dense 8 10 12)
  ;; ((polynomial y) dense 12 15 18))

;; POLYNOMIALS IN THE SAME VARIABLE

(define p1 (make-poly-from-dense-termlist 'z '(1 1)))
(define p2 (make-poly-from-dense-termlist 'z '(1 0 1)))
(mul p1 p2)
;; '((polynomial z) dense 1 1 1 1)
(expo p1 3)
;; '((polynomial z) dense 1 3 3 1)
(reduce (mul p1 p2) (expo p1 3))
;; '(((polynomial z) dense 1 0 1) ((polynomial z) dense 1 2 1))
(polynomial-division
  (make-poly-from-dense-termlist 'z '(32 0 0 0 0 -243))
  (make-poly-from-dense-termlist 'z '(2 -3)))
;; '(((polynomial z) dense 16 24 36 54 81) ((polynomial z) dense))

;; RATIONAL FUNCTIONS

(define r1 (make-rational-function
             (make-poly-from-dense-termlist 'z '(1 3 3 1))
             (make-poly-from-dense-termlist 'z '(1 -2 1)))) ; (x+1)^3/(x-1)^2
r1
;; '(rational-function
  ;; ((polynomial z) dense 1 3 3 1)
  ;; ((polynomial z) dense 1 -2 1))
(define r2 (make-rational-function
             (make-poly-from-dense-termlist 'z '(1 -3 3 -1))
             (make-poly-from-dense-termlist 'z '(1 2 1)))) ; (x-1)^3/(x+1)^2
r2
;; '(rational-function
  ;; ((polynomial z) dense 1 -3 3 -1)
  ;; ((polynomial z) dense 1 2 1))
(add r1 r2)
;; '(rational-function
  ;; ((polynomial z) dense 1 0 10 0 5 0)
  ;; ((polynomial z) dense 1 0 -2 0 1))
(mul r1 r2)
;; '(rational-function ((polynomial z) dense 1 0 -1) ((polynomial z) dense 1))
(div r1 r2)
;; '(rational-function
  ;; ((polynomial z) dense 1 5 10 10 5 1)
  ;; ((polynomial z) dense 1 -5 10 -10 5 -1))

;; WITH SPARSE OUTPUT

(define s (make-poly-from-dense-termlist
            'y
            `(1
              2.0
              ,(make-rational-number 6 2)
              ,(make-complex-from-real-imag 4 0)
              ,(make-complex-from-mag-ang -5 pi))))
(expo s 6) ; (y^4 + 2*y^3 + 3*y^2 + 4*y + 5) ^ 6
;; '((polynomial y)
  ;; sparse
  ;; (24 1)
  ;; (23 12)
  ;; (22 78)
  ;; (21 364)
  ;; (20 1365)
  ;; (19 4332)
  ;; (18 11974)
  ;; (17 29376)
  ;; (16 64818)
  ;; (15 129740)
  ;; (14 236958)
  ;; (13 396516)
  ;; (12 609389)
  ;; (11 860772)
  ;; (10 1117050)
  ;; (9 1329584)
  ;; (8 1446498)
  ;; (7 1430532)
  ;; (6 1276546)
  ;; (5 1016220)
  ;; (4 709125)
  ;; (3 422500)
  ;; (2 206250)
  ;; (1 75000)
  ;; (0 15625))