;;;from comm
;;; We rewrote sdiff, which was originally presented as
;;; a big cond -- a dispatch on the caar of e.
;;; e.g. (cond ...
;;;           ((eq (caar e) 'mrat) (ratdx e x)).....)

;;; Here we do the lookup via 
;;;  (setf fun (get (caar e) 'sdiff_function))  and apply that function if it exist.
;;; to set up the functions we do this..

;;; (setf (get 'mrat) 'sdiff_function #'ratdx)
;;; and so on.

;;; This design is used in the simplifier (circa 1963.)
;;; and for gradients of functions like sin, cos in sdiffgrad.

;;; Why was mtimes, mplus etc written in a big cond? perhaps it was
;;; thought to be faster?  Maybe at one time, but sdiff has gotten crowded,
;;; and before it gets to look for %sin, it looks for such things as
;;; mnctimes and hypergeometric functions. That can't be fast.



(defmfun sdiff (e x &aux fun) ; The args to SDIFF are assumed to be simplified.
  (cond ((alike1 e x) 1)
	((mnump e) 0)
	((or (atom e) (memq 'array (cdar e))) (chainrule e x))
	((mbagp e) (cons (car e) (sdiffmap (cdr e) x)))
	;;((not (depends e x)) 0) ;; depends does not work on CRE exprs? optimize?
	((setf fun (get (caar e) 'sdiff-fun))
	 (if fun (funcall fun e x)))
	;; no special diff function; look for gradient
	(t (sdiffgrad e x))))



(setf (get 'mrat 'sdiff-fun) 'ratdx)
(setf (get 'mplus 'sdiff-fun) #'(lambda(e x)(addn (sdiffmap (cdr e) x) t)))
(setf (get '%sum 'sdiff-fun) 'diffsumprod)
(setf (get '%product 'sdiff-fun) 'diffsumprod)
(setf (get '%lsum 'sdiff-fun) 'difflsum)
(setf (get '%at 'sdiff-fun) 'diff-%at)
(setf (get 'mtimes 'sdiff-fun) #'(lambda(e x) (addn (sdifftimes (cdr e) x) t)))
(setf (get 'mexpt 'sdiff-fun) 'diffexpt)
(setf (get 'mnctimes 'sdiff-fun) #'(lambda(e x)
	 (let (($dotdistrib t))
	   (add2 (ncmuln (cons (sdiff (cadr e) x) (cddr e)) t)
		 (ncmul2 (cadr e) (sdiff (cons '(mnctimes) (cddr e)) x))))))
(setf (get '|$~| 'sdiff-fun) #'(lambda(e x)
				 (if $vect_cross
				     (add2* `((|$~|) ,(cadr e) ,(sdiff (caddr e) x))
					    `((|$~|) ,(sdiff (cadr e) x) ,(caddr e))))))




(setf (get 'mncexpt 'sdiff-fun) 'diffncexpt)
(setf (get '%log 'sdiff-fun) 
  #'(lambda(e x)
      (sdiffgrad (cond ((and (not (atom (cadr e))) (eq (caaadr e) 'mabs))
			(cons (car e) (cdadr e)))
		       (t e))
		 x)))
(setf (get '%plog 'sdiff-fun) (get '%log 'sdiff-fun))
(setf (get '%derivative 'sdiff-fun) 
  #'(lambda(e x)
      (cond ((or (atom (cadr e)) (memq 'array (cdaadr e))) (chainrule e x))
	    ((freel (cddr e) x) (diff%deriv (cons (sdiff (cadr e) x) (cddr e))))
	    (t (diff%deriv (list e x 1))))))
(setf (get '%binomial 'sdiff-fun) 
  #'(lambda (e x)
      (let ((efact ($makefact e)))
	(mul2 (factor (sdiff efact x)) (div e efact)))))

(setf (get '$beta 'sdiff-fun) (get '%binomial 'sdiff-fun))
(setf (get '%integrate 'sdiff-fun) 'diffint)
(setf (get '%laplace 'sdiff-fun) 'difflaplace)
(setf (get '%realpart 'sdiff-fun) 
  #'(lambda(e x) (list (cons (caar e) nil) (sdiff (cadr e) x))))
(setf (get '%imagpart 'sdiff-fun) (get '%realpart 'sdiff-fun))
(setf (get 'mqapply 'sdiff-fun)
  #'(lambda(e x)
  (if  (eq (caaadr e) '$%f)
	 ;; Handle %f, hypergeometric function
	 ;;
	 ;; The derivative of %f[p,q]([a1,...,ap],[b1,...,bq],z) is
	 ;;
	 ;; a1*a2*...*ap/(b1*b2*...*bq)
	 ;;   *%f[p,q]([a1+1,a2+1,...,ap+1],[b1+1,b2+1,...,bq+1],z)
	 (let* ((arg1 (cdr (third e)))
		(arg2 (cdr (fourth e)))
		(v (fifth e)))
	   (mul (sdiff v x)
		(div (mull arg1) (mull arg2))
		`((mqapply) (($%f array) ,(length arg1) ,(length arg2))
		  ((mlist) ,@(incr1 arg1))
		  ((mlist) ,@(incr1 arg2))
		  ,v)))
    (sdiffgrad e x))))