;; Symbolic Differentiator
;; Adapted from Abelson & Sussman, Structure and Interpretion of Computer Programs
;; Chapter 2
(require test-engine/scheme-gui)

;; d/dx( x + 3 ) = 1
(check-expect (deriv '(+ x 3) 'x) 1)

;; d/dx( xy ) = y
(check-expect (deriv '(* x y) 'x) 'y)

;; d/dy( 3x+2x ) = 5
(check-expect (deriv '(+ (* 3 x) (* 2 x)) 'x) 5)

;; d/dx( 3x^2 ) = 6x
;(check-expect (deriv '(* 3 (* x x)) 'x) '(* 6 x))
(check-expect (deriv '(* 3 (* x x)) 'x) '(* 3 (+ x x)))

;; d/dx( xy(x+3) ) = xy+y(x+3)
(check-expect (deriv '(* (* x y) (+ x 3)) 'x) 
	      '(+ (* x y) (* y (+ x 3))))

;; deriv : exp var -> exp
;; apply the following differentiation rules to differentiate
;; given expression containing sums and products wrt to given variable
;;     - dx/dx = 1 , dc/dx = 0
;;     - d(u+v)/dx = du/dx + dv/dx
;;     - d(uv)/dx = u(dv/dx) + v(du/dx)
(define (deriv exp var)
  (cond [(number? exp) 0]
        [(symbol? exp)
         (cond [(symbol=? exp var) 1] 
               [else 0])]
        [(sum? exp) 
         (make-sum (deriv (sum-arg1 exp) var)
                   (deriv (sum-arg2 exp) var))]
        [(product? exp) 
         (make-sum
           (make-product (product-arg1 exp)
                         (deriv (product-arg2 exp) var))
           (make-product (deriv (product-arg1 exp) var)
                         (product-arg2 exp)))]
        [else
         (error "unknown expression type -- DERIV" exp)]))

;; an exp is either
;;   - (make-sum exp exp)
;;   - (make-product exp exp)
;;   - a symbol
;;   - a number 

(define (make-sum-raw a1 a2) (list '+ a1 a2))

(define (sum? x)
  (and (cons? x) (symbol=? (first x) '+)))

(define (sum-arg1 s) (second s))
(define (sum-arg2 s) (third s))

(define (make-product-raw m1 m2) (list '* m1 m2))

(define (product? x)
  (and (cons? x) (symbol=? (first x) '*)))

(define (product-arg1 p) (second p))
(define (product-arg2 p) (third p))

;; simplification

(define (make-sum a1 a2)
  (cond [(=number? a1 0) a2]
        [(=number? a2 0) a1]
        [(and (number? a1) (number? a2)) (+ a1 a2)]
        [else (make-sum-raw a1 a2)]))

(define (make-product m1 m2)
  (cond [(or (=number? m1 0) (=number? m2 0)) 0]
        [(=number? m1 1) m2]
        [(=number? m2 1) m1]
        [(and (number? m1) (number? m2)) (* m1 m2)]
        [else (make-product-raw m1 m2)]))

(define (=number? exp num)
  (and (number? exp) (= exp num)))

;; Numerical Differentiator
;; Adapted from Felleisen et al, How to Design Programs
;; Section 23

;; prime : (number -> number) -> (number -> number)
;; return the derivative of given function
(define (prime f)
  (let ([epsilon 0.001])
    (lambda (x)
      (/ (- (f (+ x epsilon)) (f (- x epsilon)))
         (* 2 epsilon)))))

;; line with slope 2
(check-within ((prime (lambda (x) (* 2 x))) 0) 2 0.001)

(define (step-fun x) (cond [(< x 0) 0] [else 1]))

(check-within ((prime step-fun) -1) 0 0.001)
(check-within ((prime step-fun) 1) 0 0.001)
(check-within ((prime step-fun) 0) 500 0.001)

;; compare : exp variable number -> void
;; print the results of the derivative of given function expression
;; computed symbolically and numerically
(define (compare exp var x)
  (printf "symbolic: ~a~n" 
          (eval (list (list 'lambda (list var) (deriv exp var))
                      x)))
  (printf "numerical: ~a~n" 
          ((prime (eval (list 'lambda (list var) exp)))
           x)))

(test)

;; parabola at origin
(compare '(* 3 (* x x)) 'x 0)
(compare '(* 3 (* x x)) 'x 100)