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