symbolic differentiation in R

library(rlang)

# Functions with names ending in `_impl` take quoted expressions as input.
# This removes the need for constant quoting and unquoting
simplify_sum_impl <- function(e1, e2) {
  if (is_syntactic_literal(e1) & is_syntactic_literal(e2)) {
    return(eval_bare(e1 + e2))
  }
  
  if (e1 == 0) {
    return(e2)
  }
  else if (e2 == 0 ) {
    return(e1)
  }
  else {
    return(expr(!!e1 + !!e2))
  }
}

simplify_difference_impl <- function(e1, e2) {
  if (is_syntactic_literal(e1) & is_syntactic_literal(e2)) {
    return(eval_bare(e1 - e2))
  }
  
  if (e1 == 0) {
    return(expr(-1*!!e2))
  }
  else if (e2 == 0 ) {
    return(e1)
  }
  else {
    return(expr(!!e1 - !!e2))
  }
}


simplify_prod_impl <- function(e1, e2) {
  if (is_syntactic_literal(e1) & is_syntactic_literal(e2)) {
    return(eval_bare(e1 * e2))
  }
  
  if ((e1 == 0) | (e2 == 0)) {
    return(0)
  }
  
  if (e1 == 1) {
    return(e2)
  }
  
  if (e2 == 1) {
    return(e1)
  }
  
  return(expr(!!e1 * !!e2))
}

simplify_pow_impl <- function(e1, e2) {
  if (is_syntactic_literal(e1) & is_syntactic_literal(e2)) {
    return(eval_bare(e1^e2))
  }
  
  if ((e1 != 0) & (e2 == 0)) {
    return(1)
  }
  
  if (e2 == 1) {
    return(e1)
  }
  
  if (e1 == 1) {
    return(1)
  }
  
  return(expr((!!e1)^(!!e2)))
}

d_sum_impl <- function(e1, e2, var) {
  de1 <- d(!!e1, !!var)
  de2 <- d(!!e2, !!var)
  
  simplify_sum_impl(de1, de2)
}

d_difference_impl <- function(e1, e2, var) {
  de1 <- d(!!e1, !!var)
  de2 <- d(!!e2, !!var)
  
  simplify_difference_impl(de1, de2)
}

d_prod_impl <- function(e1, e2, var) {
  de1 <- d(!!e1, !!var)
  de2 <- d(!!e2, !!var)
  
  simplify_sum_impl(simplify_prod_impl(de1, e2), simplify_prod_impl(e1, de2))
}

d_frac_impl <- function(e1, e2, var) {
  de1 <- d(!!e1, !!var)
  de2 <- d(!!e2, !!var)
  
  numerator <- simplify_difference_impl(simplify_prod_impl(de1, e2), simplify_prod_impl(e1, de2))
  denominator <- simplify_pow_impl(e2, 2)
  
  expr(!!numerator/!!denominator)
}

d_pow_impl <- function(e1, e2, var) {
  if (!is_syntactic_literal(e2)) {
    return(expr(d(!!e1^!!e2, !!var))) # cannot differentiate, leave symbolic
  }
  
  de1 <- d(!!e1, !!var)
  newexp <- e2 - 1
  
  simplify_prod_impl(simplify_prod_impl(de1, e2), simplify_pow_impl(e1, newexp))
}

d_parens_impl <- function(e1, var) {
  de1 <- d(!!e1, !!var)
  expr((!!de1))
}

d_exp_impl <- function(e1, var) {
  de1 <- d(!!e1, !!var)
  simplify_prod_impl(de1, expr(exp(!!e1)))
}

d <- function(e, var = x) {
  e <- enexpr(e)
  var <- enexpr(var)
  
  if (is_syntactic_literal(e)) {
    return(0) # derivative of a numeric constant is 0
  }

  if (is_call(e, "+")) {
    return(d_sum_impl(e[[2]], e[[3]], var)) # differentiate sum
  }

  if (is_call(e, "-")) {
    return(d_difference_impl(e[[2]], e[[3]], var)) # differentiate difference
  }
    
  if (is_call(e, "*")) {
    return(d_prod_impl(e[[2]], e[[3]], var)) # differentiate product
  }
  
  if (is_call(e, "/")) {
    return(d_frac_impl(e[[2]], e[[3]], var)) # differentiate fraction
  }
  
  if (is_call(e, "^")) {
    return(d_pow_impl(e[[2]], e[[3]], var)) # differentiate power
  }
  
  if (is_call(e, "exp")) {
    return(d_exp_impl(e[[2]], var)) # differentiate exp function
  }
 
  if (is_call(e, "(")) {
    return(d_parens_impl(e[[2]], var)) # differentiate expression in parentheses
  }
   
  if (is_call(e)) {
    # cannot differentiate, leave symbolic
    return(expr(d(!!e, !!var)))
  }
  
  if (e == var) {
    return(1) # d(x, x) = 1
  }
  else {
    return(0) # d(y, x) = 0
  }
}

d(5*x^2+8*x+5, x)
d(exp(5*y^2+8*x*y+5*x)*(2*y+7), y)
d(exp(5*y^2+8*x*y+5*x)*(2*y+7), x)
d(exp(5*x^2)/(3*x-2), x)
d(exp(2*x*f(x)), x)