Symbolic Differentiation
; ; This is essentially the program given by Ravi Sethi in "Programming ; Languages: Concepts and Constructs," which performs symbolic ; differentiation of a polynomial. ; ; First, establish an error code for bad parsing. (errcode ERR_BADEXPR) ; ; The polynomials are represented in Lisp notation, so sums (+ ...) ; and produces are (* ...), etc. The top-level function decides what ; the sort of thing it's working on, and applies an appropriate worker ; function as below. (define (d x E) (cond ((constant? E) (diff-constant x E)) ((variable? E) (diff-variable x E)) ((sum? E) (diff-sum x E)) ((product? E) (diff-product x E)) (#t (error ERR_BADEXPR "Cannot parse expression.")) ) ) ; ; A constant, for our implementation, is an integer, and the differential ; of a constant is zero. (setq constant? int?) (define (diff-constant x E) 0) ; ; A variable is just a id atom, and the dx x = 1, dx y = 0. (setq variable? id?) (define (diff-variable x E) (if (equal? x E) 1 0) ) ; ; Utilities for sums and products. (define (terms E) (cdr E)) (define (factors E) (cdr E)) (define (sum? E) (and (pair? E) (equal? '+ (car E))) ) (define (product? E) (and (pair? E) (equal? '* (car E))) ) (define (make-sum x) (cons '+ x)) (define (make-product x) (cons '* x)) ; ; The differential of a sum is the sum of the differentials: ; d (E1 + E2) = d E1 + d E2 (define (diff-sum x E) (make-sum (map (lambda (expr) (d x expr)) (terms E)) ) ) ; This beast simply takes the terms of E, runs d x on each, then makes the ; sum of the result. ; ; The differential of a produce uses the chain rule: ; d E1 E2 = E1 d E2 + E2 d E1 ; This is actually implemented by chain-rule; the diff-product function ; filters out the cases where the product has fewer than two factors. (define (diff-product x E) (let ; See how many things we're got. ((nfact (length (factors E)))) (cond ; No factors. Just a long-winded way to write the constant 1, ; whose differential is zero. ((equal? nfact 0) 0) ; One factor. Not really much of a product. Loose the multiply ; and take the differential of the single factor. ((equal? nfact 1) (d x (car (factors E)))) ; Real case. (#t (chain-rule x E)) ) ) ) ; ; Here's the actual chain rule function. ; d E1 E2 = E1 d E2 + E2 d E1 (define (chain-rule x E) (let ( (E1 (car (factors E))) ; First factor. (E2 (make-product (cdr (factors E)))) ; (* Other factors) (dE1 (d x E1)) ; d E1 (dE2 (d x E2)) ; d E2 ) (make-sum (list (make-product (list E1 dE2)) (make-product (list E2 dE1)) ) ) ) ) ; ; Some basic simplification. Not easy to do "completely," but this helps ; a lot. ; ; Top-level simplify. (define (simplify E) (cond ((sum? E) (simplify-sum E)) ((product? E) (simplify-product E)) (#t E) ) ) ; ; The sum and product simplifiers are mainly calls to simpl, with some ; appropriate control parameters. The parameters are the corresponding ; identifier and make- function, and the identity for that operation. (define (simplify-sum E) (simpl sum? make-sum 0 E) ) (define (simplify-product E) (simpl product? make-product 1 E) ) ; ; Here's the simplifier. (define (simpl isit? addop ident E) (let ( (parts (cdr E)) ; Terms or factors. (sparts (map simplify parts)) ; Terms or factors simplified. (fparts (flat isit? sparts)) ; Simp (* x (* y z)) to (* x y z) (zout (replace-zero fparts)) ; Reduce (* ... 0 ...) to 0. (unid (select (lambda (x) (not (equal? x ident))) zout)) ; Remove identity (0 for + 1 for *) ) (proper addop ident unid) ; Cleanup; see below. ) ) ; ; The flat function looks for subexpressions of the same operator and merges ; them in. For instance, change (+ x y (+ z w) (+ q 4) g) to ; (+ x y z w q 4 g). (define (flat isit args) (cond ((null? args) ()) ; Empty is empty. ((not (pair? args)) (list args)); I don't see how this happens, but ok. ((isit (car args)) ; If first arg same op, combine. (append (flat isit (cdar args)) (flat isit (cdr args))) ) (#t ; Default: go on to the next. (cons (car args) (flat isit (cdr args))) ) ) ) ; ; This simply adds the operator back to a list of terms or factors, but it ; avoids turning the empty list into (+) or the singleton list int (* 17). (define (proper addop ident args) (cond ((null? args) ident) ; () becomes 0 or 1. ((null? (cdr args)) (car args)) ; (x) becomes x (#t (addop args)) ; (x y z) to (+/* x y z) ) ) ; ; See if the expression is a multiplication containing zero (hence equal ; to zero). (define (is-zero-mult? E) (and (product? E) (not (equal? (select (lambda (x) (equal? x 0)) (factors E)) ())) ) ) ; ; Replace zero-valued multiply items with zero. This works only on ; the top level. (define (replace-zero list) (if (null? list) () (cons (if (is-zero-mult? (car list)) 0 (car list)) (replace-zero (cdr list)) ) ) )