# derivative



#define F p3
#define X p4
#define N p5

Eval_derivative = ->

  # evaluate 1st arg to get function F

  i = 0
  p1 = cdr(p1)
  push(car(p1))
  Eval()

  # evaluate 2nd arg and then...

  # example  result of 2nd arg  what to do
  #
  # d(f)    nil      guess X, N = nil
  # d(f,2)  2      guess X, N = 2
  # d(f,x)  x      X = x, N = nil
  # d(f,x,2)  x      X = x, N = 2
  # d(f,x,y)  x      X = x, N = y

  p1 = cdr(p1)
  push(car(p1))
  Eval()

  p2 = pop()
  if (p2 == symbol(NIL))
    guess()
    push(symbol(NIL))
  else if (isNumericAtom(p2))
    guess()
    push(p2)
  else
    push(p2)
    p1 = cdr(p1)
    push(car(p1))
    Eval()

  p5 = pop(); # p5 is N
  p4 = pop(); # p4 is X
  p3 = pop(); # p3 is F

  while (1)

    # p5 (N) might be a symbol instead of a number

    if (isNumericAtom(p5)) # p5 is N
      push(p5); # p5 is N
      n = pop_integer()
      if (isNaN(n))
        stop("nth derivative: check n")
    else
      n = 1

    push(p3); # p3 is F

    if (n >= 0)
      for i in [0...n]
        push(p4); # p4 is X
        derivative()
    else
      n = -n
      for i in [0...n]
        push(p4); # p4 is X
        integral()

    p3 = pop(); # p3 is F

    # if p5 (N) is nil then arglist is exhausted

    if (p5 == symbol(NIL)) # p5 is N
      break

    # otherwise...

    # N    arg1    what to do
    #
    # number  nil    break
    # number  number    N = arg1, continue
    # number  symbol    X = arg1, N = arg2, continue
    #
    # symbol  nil    X = N, N = nil, continue
    # symbol  number    X = N, N = arg1, continue
    # symbol  symbol    X = N, N = arg1, continue

    if (isNumericAtom(p5)) # p5 is N
      p1 = cdr(p1)
      push(car(p1))
      Eval()
      p5 = pop(); # p5 is N
      if (p5 == symbol(NIL)) # p5 is N
        break;    # arglist exhausted
      if (isNumericAtom(p5)) # p5 is N
        doNothing = 1    # N = arg1
      else
        p4 = p5; # p5 is N   # p4 is X  # X = arg1
        p1 = cdr(p1)
        push(car(p1))
        Eval()
        p5 = pop(); # p5 is N  # N = arg2
    else
      p4 = p5;   # p5 is N   # p4 is X  # X = N
      p1 = cdr(p1)
      push(car(p1))
      Eval()
      p5 = pop();   # p5 is N  # N = arg1

  push(p3);  # p3 is F # final result

derivative = ->
  save()
  p2 = pop()
  p1 = pop()
  if (isNumericAtom(p2))
    stop("undefined function")
  if (istensor(p1))
    if (istensor(p2))
      d_tensor_tensor()
    else
      d_tensor_scalar()
  else
    if (istensor(p2))
      d_scalar_tensor()
    else
      d_scalar_scalar()
  restore()

d_scalar_scalar = ->
  if (issymbol(p2))
    d_scalar_scalar_1()
  else
    # Example: d(sin(cos(x)),cos(x))
    # Replace cos(x) <- X, find derivative, then do X <- cos(x)
    push(p1);    # sin(cos(x))
    push(p2);    # cos(x)
    push(symbol(SECRETX));  # X
    subst();    # sin(cos(x)) -> sin(X)
    push(symbol(SECRETX));  # X
    derivative()
    push(symbol(SECRETX));  # X
    push(p2);    # cos(x)
    subst();    # cos(X) -> cos(cos(x))

d_scalar_scalar_1 = ->
  # d(x,x)?

  if (equal(p1, p2))
    push(one)
    return

  # d(a,x)?

  if (!iscons(p1))
    push(zero)
    return

  if (isadd(p1))
    dsum()
    return

  if (car(p1) == symbol(MULTIPLY))
    dproduct()
    return

  if (car(p1) == symbol(POWER))
    dpower()
    return

  if (car(p1) == symbol(DERIVATIVE))
    dd()
    return

  if (car(p1) == symbol(LOG))
    dlog()
    return

  if (car(p1) == symbol(SIN))
    dsin()
    return

  if (car(p1) == symbol(COS))
    dcos()
    return

  if (car(p1) == symbol(TAN))
    dtan()
    return

  if (car(p1) == symbol(ARCSIN))
    darcsin()
    return

  if (car(p1) == symbol(ARCCOS))
    darccos()
    return

  if (car(p1) == symbol(ARCTAN))
    darctan()
    return

  if (car(p1) == symbol(SINH))
    dsinh()
    return

  if (car(p1) == symbol(COSH))
    dcosh()
    return

  if (car(p1) == symbol(TANH))
    dtanh()
    return

  if (car(p1) == symbol(ARCSINH))
    darcsinh()
    return

  if (car(p1) == symbol(ARCCOSH))
    darccosh()
    return

  if (car(p1) == symbol(ARCTANH))
    darctanh()
    return

  if (car(p1) == symbol(ABS))
    dabs()
    return

  if (car(p1) == symbol(SGN))
    dsgn()
    return

  if (car(p1) == symbol(HERMITE))
    dhermite()
    return

  if (car(p1) == symbol(ERF))
    derf()
    return

  if (car(p1) == symbol(ERFC))
    derfc()
    return

  if (car(p1) == symbol(BESSELJ))
    if (isZeroAtomOrTensor(caddr(p1)))
      dbesselj0()
    else
      dbesseljn()
    return

  if (car(p1) == symbol(BESSELY))
    if (isZeroAtomOrTensor(caddr(p1)))
      dbessely0()
    else
      dbesselyn()
    return

  if (car(p1) == symbol(INTEGRAL) && caddr(p1) == p2)
    derivative_of_integral()
    return

  dfunction()

dsum = ->
  h = tos
  p1 = cdr(p1)
  while (iscons(p1))
    push(car(p1))
    push(p2)
    derivative()
    p1 = cdr(p1)
  add_all(tos - h)

dproduct = ->
  i = 0
  j = 0
  n = 0
  n = length(p1) - 1
  for i in [0...n]
    p3 = cdr(p1)
    for j in [0...n]
      push(car(p3))
      if (i == j)
        push(p2)
        derivative()
      p3 = cdr(p3)
    multiply_all(n)
  add_all(n)

#-----------------------------------------------------------------------------
#
#       v
#  y = u
#
#  log y = v log u
#
#  1 dy   v du           dv
#  - -- = - -- + (log u) --
#  y dx   u dx           dx
#
#  dy    v  v du           dv
#  -- = u  (- -- + (log u) --)
#  dx       u dx           dx
#
#-----------------------------------------------------------------------------

dpower = ->
  push(caddr(p1));  # v/u
  push(cadr(p1))
  divide()

  push(cadr(p1));    # du/dx
  push(p2)
  derivative()

  multiply()

  push(cadr(p1));    # log u
  logarithm()

  push(caddr(p1));  # dv/dx
  push(p2)
  derivative()

  multiply()

  add()

  push(p1);    # u^v

  multiply()

dlog = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  divide()

#  derivative of derivative
#
#  example: d(d(f(x,y),y),x)
#
#  p1 = d(f(x,y),y)
#
#  p2 = x
#
#  cadr(p1) = f(x,y)
#
#  caddr(p1) = y

dd = ->
  # d(f(x,y),x)

  push(cadr(p1))
  push(p2)
  derivative()

  p3 = pop()

  if (car(p3) == symbol(DERIVATIVE))

    # sort dx terms

    push_symbol(DERIVATIVE)
    push_symbol(DERIVATIVE)
    push(cadr(p3))

    if (lessp(caddr(p3), caddr(p1)))
      push(caddr(p3))
      list(3)
      push(caddr(p1))
    else
      push(caddr(p1))
      list(3)
      push(caddr(p3))

    list(3)

  else
    push(p3)
    push(caddr(p1))
    derivative()

# derivative of a generic function

dfunction = ->
  p3 = cdr(p1);  # p3 is the argument list for the function

  if (p3 == symbol(NIL) || Find(p3, p2))
    push_symbol(DERIVATIVE)
    push(p1)
    push(p2)
    list(3)
  else
    push(zero)

dsin = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  cosine()
  multiply()

dcos = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  sine()
  multiply()
  negate()

dtan = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  cosine()
  push_integer(-2)
  power()
  multiply()

darcsin = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push_integer(1)
  push(cadr(p1))
  push_integer(2)
  power()
  subtract()
  push_rational(-1, 2)
  power()
  multiply()

darccos = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push_integer(1)
  push(cadr(p1))
  push_integer(2)
  power()
  subtract()
  push_rational(-1, 2)
  power()
  multiply()
  negate()

#        Without simplify  With simplify
#
#  d(arctan(y/x),x)  -y/(x^2*(y^2/x^2+1))  -y/(x^2+y^2)
#
#  d(arctan(y/x),y)  1/(x*(y^2/x^2+1))  x/(x^2+y^2)

darctan = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push_integer(1)
  push(cadr(p1))
  push_integer(2)
  power()
  add()
  inverse()
  multiply()
  simplify()

dsinh = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  ycosh()
  multiply()

dcosh = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  ysinh()
  multiply()

dtanh = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  ycosh()
  push_integer(-2)
  power()
  multiply()

darcsinh = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  push_integer(2)
  power()
  push_integer(1)
  add()
  push_rational(-1, 2)
  power()
  multiply()

darccosh = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  push_integer(2)
  power()
  push_integer(-1)
  add()
  push_rational(-1, 2)
  power()
  multiply()

darctanh = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push_integer(1)
  push(cadr(p1))
  push_integer(2)
  power()
  subtract()
  inverse()
  multiply()

dabs = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  sgn()
  multiply()

dsgn = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  dirac()
  multiply()
  push_integer(2)
  multiply()

dhermite = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push_integer(2)
  push(caddr(p1))
  multiply()
  multiply()
  push(cadr(p1))
  push(caddr(p1))
  push_integer(-1)
  add()
  hermite()
  multiply()

derf = ->
  push(cadr(p1))
  push_integer(2)
  power()
  push_integer(-1)
  multiply()
  exponential()
  if evaluatingAsFloats
    push_double(Math.PI)
  else
    push_symbol(PI)
  push_rational(-1,2)
  power()
  multiply()
  push_integer(2)
  multiply()
  push(cadr(p1))
  push(p2)
  derivative()
  multiply()

derfc = ->
  push(cadr(p1))
  push_integer(2)
  power()
  push_integer(-1)
  multiply()
  exponential()
  if evaluatingAsFloats
    push_double(Math.PI)
  else
    push_symbol(PI)
  push_rational(-1,2)
  power()
  multiply()
  push_integer(-2)
  multiply()
  push(cadr(p1))
  push(p2)
  derivative()
  multiply()

dbesselj0 = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  push_integer(1)
  besselj()
  multiply()
  push_integer(-1)
  multiply()

dbesseljn = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  push(caddr(p1))
  push_integer(-1)
  add()
  besselj()
  push(caddr(p1))
  push_integer(-1)
  multiply()
  push(cadr(p1))
  divide()
  push(cadr(p1))
  push(caddr(p1))
  besselj()
  multiply()
  add()
  multiply()

dbessely0 = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  push_integer(1)
  besselj()
  multiply()
  push_integer(-1)
  multiply()

dbesselyn = ->
  push(cadr(p1))
  push(p2)
  derivative()
  push(cadr(p1))
  push(caddr(p1))
  push_integer(-1)
  add()
  bessely()
  push(caddr(p1))
  push_integer(-1)
  multiply()
  push(cadr(p1))
  divide()
  push(cadr(p1))
  push(caddr(p1))
  bessely()
  multiply()
  add()
  multiply()

derivative_of_integral = ->
  push(cadr(p1))