--
--
-- Differentiation
--
--

∂/∂ $f *x :=
  match f as mathExpr with
    -- symbol
    | #x -> 1
    | ?isSymbol -> 0
    -- function expression
    | func _ $argnames $args _ ->
      sum (map2 (\s r -> (userRefs f [s]) * ∂/∂ r x) argnames args)
    -- function application
    | #`exp $g -> exp g * ∂/∂ g x
    | #`log $g -> 1 / g * ∂/∂ g x
    | #`sqrt $g -> 1 / (2 * sqrt g) * ∂/∂ g x
    | #`(^) $g $h -> f * ∂/∂ (log g * h) x
    | #`cos $g -> (- sin g) * ∂/∂ g x
    | #`sin $g -> cos g * ∂/∂ g x
    | #`arccos $g -> 1 / sqrt (1 - g ^ 2) * ∂/∂ g x
    | apply $g $args ->
      sum (map 2#((capply `(userRefs g [%1]) args) * ∂/∂ %2 x) (zip nats args))
    -- quote
    | quote $g ->
      let g' := ∂/∂ g x
       in if isMonomial g'
            then g'
            else let d := capply gcd (fromPoly g')
                  in d *' '(mapPoly (/' d) g')
    -- 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 1#(∂/∂ %1 x) ts)
    -- quotient
    | $p1 / $p2 ->
      let p1' := ∂/∂ p1 x
          p2' := ∂/∂ p2 x
       in (p1' * p2 - p2' * p1) / p2 ^ 2

d/d := ∂/∂

pd/pd := ∂/∂

∇ := ∂/∂

nabla := ∇

grad := ∇

taylorExpansion $f $x $a := multivariateTaylorExpansion f [|x|] [|a|]

maclaurinExpansion := 2#(taylorExpansion %1 %2 0)

multivariateTaylorExpansion $f %xs %ys :=
  withSymbols [h]
    let hs := generateTensor 1#h_%1 (tensorShape xs)
     in map2
          (*)
          (map 1#(1 / fact %1) nats0)
          (map
             (compose
                1#(V.substitute xs ys %1)
                1#(V.substitute hs (withSymbols [i] xs_i - ys_i) %1))
             (iterate (compose 1#(∇ %1 xs) 1#(V.* hs %1)) f))

multivariateMaclaurinExpansion $f %xs :=
  multivariateTaylorExpansion f xs (tensorMap 1#0 xs)

addUserScript $f $i :=
  let (g, is) := deconsUserScripts f
   in appendUserScripts g (sort (is ++ [i]))