(define-class Polynomial Object
((= variable immutable:) ; a symbol
(= terms immutable:))) ; a list of terms
(define (Polynomial-variable= var1 var2)
(eq? var1 var2))
;;; a term is a pair (order coeff) (order should really be degree, but ...)
;;; We're going to use a Meroon class for terms to aid debugging.
(define-class term Object
((= order)
(= coeff)))
#|
(define make-term cons)
(define term-order car)
(define term-coeff cdr)
(define term? pair?)
|#
(define (every? p list)
(or (null? list)
(and (p (car list))
(every? p (cdr list)))))
(define (any? p list)
(and (not (null? list))
(or (p (car list))
(any? p (cdr list)))))
(define-method (initialize! (p Polynomial))
(let ((terms (Polynomial-terms p)))
(cond ((not (symbol? (Polynomial-variable p)))
(error "initialize! for Polynomials: The variable is not a symbol: " p))
((not (list? terms))
(error "intialize! for Polynomials: (Polynomial-terms p) is not a list " p))
((not (every? term? terms))
(error "intialize! for Polynomials: Not all entries in (Polynomial-terms p) are terms " p))
((any? (lambda (t) (=zero? (term-coeff t))) terms)
(error "intialize! for Polynomials: Not all term coefficents in are nonzero " p))
(else
(let loop ((terms terms))
(cond ((or (null? terms)
(null? (cdr terms)))
;; We're done checking the orders of terms
)
((> (term-order (car terms))
(term-order (cadr terms)))
(loop (cdr terms)))
(else
(error "intialize! for Polynomials: terms are not in decreasing order: " p)))))))
(call-next-method))
;;; See the Haskell Wiki page
;;; http://www.haskell.org/haskellwiki/Fold
;;; for a good explanation, together with pictures, for how
;;; fold-left and fold-right work.
(define (fold-left operator ;; called with (operator result-so-far (car list))
initial-value
list)
(if (null? list)
initial-value
(fold-left operator
(operator initial-value (car list))
(cdr list))))
(define (fold-right operator ;; called with (operator (car list) result-so-far)
initial-value
list)
(if (null? list)
initial-value
(operator (car list)
(fold-right operator initial-value (cdr list)))))
;;; map is a builtin function that works on lists, but it could
;;; be defined as follows:
(define (my-map f list)
(fold-right (lambda (v l)
(cons (f v) l))
'()
list))
;; You want every term in a Polynomial's term-list to have a nonzero
;; coefficient. To ensure this you always use adjoin-term to add a
;; term to a term-list
(define (adjoin-term term term-list)
(if (=zero? (term-coeff term))
term-list
(cons term term-list)))
;; So when you're building up a term-list by applying a function
;; to another list, always use map-termlist:
(define (map-termlist f list)
(fold-right (lambda (v l)
(adjoin-term (f v) l))
'()
list))
(define-generic (add (x) y)
(if (number? x)
(if (number? y)
;; We know how to add two numbers.
(+ x y)
;; We don't know how to add a number to an arbitrary object,
;; but perhaps there is a method to add the object y to
;; a number, so we swap arguments and try again.
(add y x))
;; If x isn't a number, we don't have a method for it.
(error "add: This generic is not defined on these objects: " x y)))
(define-method (add (p_1 Polynomial) p_2)
(cond ((number? p_2)
(add p_1 (number->Polynomial p_2 (Polynomial-variable p_1))))
((and (Polynomial? p_2)
(Polynomial-variable= (Polynomial-variable p_1)
(Polynomial-variable p_2)))
(instantiate Polynomial
variable: (Polynomial-variable p_1)
terms: (add-terms (Polynomial-terms p_1)
(Polynomial-terms p_2))))
(else
(error "add: p_2 is neither a number nor a polynomial with the same variable as p_1 " p_1 p_2))))
;; We first wrote add-terms, but then realized that the control-flow could
;; be abstracted into merge-lists
;; I don't know whether this is the right abstraction, but ...
;; One can define union on sorted lists (no repetition of elements) by
;; (define (union l1 l2) (merge-lists precedes? (lambda (t1 t2 l) (cons t1 l)) cons l1 l2))
;; One can define merge on sorted lists (repetition of elements) by
;; (define (union l1 l2) (merge-lists precedes? (lambda (t1 t2 l) (cons t1 (cons t2 l))) cons l1 l2))
(define (merge-lists precedes?
combine-and-adjoin
adjoin
l1 l2)
(define (do-it l1 l2)
(cond ((null? l1)
l2)
((null? l2)
l1)
(else
(let ((t1 (car l1))
(t2 (car l2)))
(cond ((precedes? t1 t2)
(adjoin t1 (do-it (cdr l1) l2)))
((precedes? t2 t1)
(adjoin t2 (do-it l1 (cdr l2))))
(else
(combine-and-adjoin t1 t2
(do-it (cdr l1) (cdr l2)))))))))
(do-it l1 l2))
(define (add-terms l1 l2)
(merge-lists (lambda (t1 t2)
(> (term-order t1)
(term-order t2)))
(lambda (t1 t2 rest)
(adjoin-term (make-term (term-order t1)
(add (term-coeff t1)
(term-coeff t2)))
rest))
adjoin-term
l1 l2))
(define (multiply-terms l1 l2)
(fold-right (lambda (t1 result-so-far)
(add-terms (multiply-term-by-all-terms t1 l2)
result-so-far))
'()
l1))
(define (multiply-term-by-all-terms t1 L)
(map-termlist (lambda (t2)
(make-term (+ (term-order t1)
(term-order t2))
(multiply (term-coeff t1)
(term-coeff t2))))
L))
(define-method (show (p Polynomial) . stream)
(let ((port (if (null? stream)
(current-output-port)
(car stream))))
(if (=zero? p)
(display 0)
(show-terms (Polynomial-variable p)
(Polynomial-terms p)
port))
(newline port)))
(define (show-terms variable terms port)
(show-first-term variable (car terms) port)
(for-each (lambda (term)
(show-term variable term port))
(cdr terms)))
(define (show-first-term variable term port)
(let ((coeff (term-coeff term))
(order (term-order term)))
(print port: port
(list (let ((sign (if (negative? coeff) "-" '()))
(abs-coeff (abs coeff)))
(list sign
(if (and (< 0 order)
(= abs-coeff 1))
'()
abs-coeff)))
(cond ((zero? order) '())
((= order 1) variable)
(else
(list variable "^" order)))))))
(define (show-term variable term port)
(let ((coeff (term-coeff term))
(order (term-order term)))
(print port: port
(list (if (negative? coeff)
"-"
"+")
(let ((abs-coeff (abs coeff)))
(if (and (eq? coeff 1)
(< 0 order))
'()
(abs coeff)))
(cond ((zero? order) '())
((= order 1) variable)
(else
(list variable "^" order)))))))
(define-generic (differentiate (f) variable)
(if (number? f)
0
(error "differentiate: argument not of correct type " f)))
(define-generic (integrate (f) variable #!optional (a #f) (b #f))
(if (number? f)
(if (and (number? a)
(number? b))
(multiply f (subtract b a))
(instantiate Polynomial
variable: variable
terms: (adjoin-term (make-term 1 f) '())))
(error "integrate: unknown argument type " f variable)))
(define-method (integrate (p Polynomial) variable #!optional (a #f) (b #f))
(if (Polynomial-variable= (Polynomial-variable p)
variable)
(let ((indefinite-integral
(instantiate Polynomial
variable: variable
terms: (map-termlist integrate-term (Polynomial-terms p)))))
(if (and (number? a)
(number? b))
(subtract (evaluate indefinite-integral b)
(evaluate indefinite-integral a))
indefinite-integral))
(error "integrate: The variable of integration is not the variable of the polynomial " p variable)))
(define (make-inner-product weight variable left right)
(lambda (p q)
(integrate (multiply p (multiply q weight)) ;; weight can be a constant
variable left right)))
;;; The Gauss-Lobatto weight on (-1, 1)
(define (G-L-weight variable)
;; 1-x^2
(let ((X (variable->Polynomial variable)))
(subtract 1 (multiply X X))))
(define (G-L-inner-product variable left right)
(make-inner-product (G-L-weight variable) variable left right))