;;;;;
;;;;;
;;;;; Differentiation
;;;;;
;;;;;

(define $∂/∂
  (lambda [$f $x]
    (match f math-expr
      {; symbol
       [,x 1]
       [?symbol? 0]
       ; function application
       [(,exp $g) (* (exp g) (∂/∂ g x))]
       [(,log $g) (* (/ 1 g) (∂/∂ g x))]
       [(,cos $g) (* (* -1 (sin g)) (∂/∂ g x))]
       [(,sin $g) (* (cos g) (∂/∂ g x))]
       [(,sqrt $g) (* (/ 1 (* 2 (sqrt g))) (∂/∂ g x))]
       [(,** $g $h) (* f (∂/∂ (* (log g) h) x))]
       [<apply $g $args>
        (sum (map 2#(* (capply `g|%1 args) (∂/∂ %2 x))
                  (zip nats args)))]
       ; quote
       [<quote $g> (∂/∂ g x)]
       ; term (constant)
       [,0 0]
       [(* _ ,1) 0]
       ; term (multiplication)
       [(* ,1 $fx^$n) (* n (** fx (- n 1)) (∂/∂ fx x))]
       [(* $a $fx^$n $r)
        (+ (* a (∂/∂ (**' fx n) x) r)
           (* a (**' fx n) (∂/∂ r x)))]
       ; polynomial
       [<poly $ts> (sum (map (∂/∂ $ x) ts))]
       ; quotient
       [(/ $p1 $p2)
        (let {[$p1' (∂/∂ p1 x)]
              [$p2' (∂/∂ p2 x)]}
          (/ (- (* p1' p2) (* p2' p1)) (** p2 2)))]
       })))

(define $d/d ∂/∂)
(define $pd/pd ∂/∂)

(define $∇ ∂/∂)
(define $nabla ∇)

(define $grad ∇)

(define $taylor-expansion
  (lambda [%f %xs %as]
    (with-symbols {h}
      (let {[$hs (generate-tensor 1#h_%1 (tensor-size xs))]}
        (map2 *
              (map 1#(/ 1 (fact %1)) nats0)
              (map (compose (V.substitute xs as $)
                            (V.substitute hs (with-symbols {i} (- xs_i as_i)) $))
                   (iterate (compose 1#(∇ %1 xs) 1#(V.* hs %1)) f)))))))

(define $maclaurin-expansion
  (lambda [%f %xs]
    (multivariate-taylor-expansion f xs (tensor-map 1#0 xs))))