;;; -*- Mode:Common-Lisp;  Base:10 -*-
;; Written by: Richard J. Fateman
;; File: diffrat.lisp
;; (c) 1991 Richard J. Fateman

;; This is a derivative-divides integration program
;; written to use rational forms. Much of the work is
;; in peripheral issues, like computing derivatives of
;; random functions and testing to see if expressions are
;; free of variables.  This program depends on simplification
;; to find out if expressions divide out evenly, and will lose
;; (fail to integrate), if it is required to know identities
;; that are algebraic or transcendental (e.g. sin^2+cos^2 =1).

;; Also, this program produces ANTIDERIVATIVES, not really integrals.
;; That is, the domain of the answer may not be the same as the domain
;; of integration (unstated, but lurking in most applications).
;;(provide 'diffrat)
(in-package :mma)
;;(require "poly")(require "simp1")(require "rat1")

(defun Int(e x)(diffdiv e x)) ;testing program
(defun D(e x)(outof-rat(ratdiff (into-rat e) x)))

;; Some useful utility programs

;; pfreevar returns t if the poly p is free of the variable numbered v.

(defun pfreevar (p v)
  (labels
   ((pf1(x)
	(if (coefp x) t
	  (and (freevar (gethash (svref x 0) revtab) (gethash v revtab))
	       (every #'pf1 x)))))
   (pf1 p)))

;; freevar returns t if  x is free of the variable or kernel v.
;; x and v are in list representation. Only explicit dependencies
;; are discovered here.

(defun freevar (x v)
  (labels ((f1 (x)
	       (cond ((eq x v) nil)
		     ;; here could check for "x depends on v indirectly."
		     ((consp x) (every #'f1 (cdr x)))
		     ;; could check for poly or rat stuff
		     (t t))))
	  (f1 x)))

;; see if a rational form r is free of a variable, in list form, v. v may
;; be inside a kernel as, say (Log v) or a regular variable.
;; As a particular use, we know that if rfreevar (r v) is t, 
;; then r is a constant wrt v, and integrate(r,v) = r*v {+ a const}.

(defun rfreevar(r v) ;; r is a rat form. v is a list-form
  (let ((v (look-up-var v)))
    (flet ((rf1 (x)
		(pfreevar (car x) v)))
	  (and (every #'rf1 (rat-numerator r))
	       (every #'rf1 (rat-denominator r))))))

;; (ratderiv r v) computes the partial derivative
;; of the rational expression r wrt to the list-form expression v,
;; which is presumed to be either an indeterminate, or perhaps
;; a kernel, in which case it should be simplified. The answer is
;; in rational form.

;; In an attempt to do this efficiently in rational form we consider
;; two cases:
;; (a) v is an indeterminate and the only kernel in r that
;; involves v is exactly v.  This is the "fast" case, in the sense
;; that almost all the arithmetic can be done on polynomials.  It is
;; a subset of case b, and if we had only one program to write, we'd
;; have to write case b.
;; (b) r includes kernels that depend on v in other ways.  For example,
;; (Log v) or (Sin v). All the arithmetic must be done in rational form
;; if the derivative is rational (not polynomial) in v. (Log v) is like
;; that.

;; This may seem like a long way to do it, but let's try, anyway.

;; collectvars returns a list of all the variable (numbers) in a
;; rational form r.  E.g. (collectvars #(3 #(1 1 1) 0 1)) -> (1 3)

(defun collectvars (r &aux (cv nil))
  (labels
   ((cv2 (pr) ;; cv2 is applied to each pair: poly . power
	 (cv3 (car pr)))
    (cv3 (p) ;; cv3 is applied to each polynomial
	 (cond ((coefp p) nil)
	       (t (pushnew (svref p 0) cv)		  
		  (do ((i (1- (length p))(1- i)))
		      ((= i 0) nil)
		      (cv3 (svref p i)))))))
	      
   (mapc #'cv2 (rat-numerator r))
   (mapc #'cv2 (rat-denominator r))
   cv))

;; for each variable in a rational form, we want to compute (or remind
;; ourselves that we already know) its derivative wrt *var*


(declaim (special result difftab *var* *varnum* *sign* *r*))


(defun ratdiff(r v)
  (flet
   ((setupderiv (h) ; variable is global, h is a number..
  ;; we store derivative on a special derivative hash-table
		(setf (gethash h difftab) 
		      (into-rat(gendiff (gethash h revtab) *var*)))))
   (let*
      ((c (collectvars r))
       (*var* v)
       (*r* r)
       (*varnum* (gethash v vartab))
       (difftab (make-hash-table :size (+ 2 (length c))))
       (*sign* 1)
       (result (poly2rat 0 1)))

    ;; set up derivatives for all kernels in this rational expression
    (mapc #'setupderiv c)
    ;;   (showhash difftab)   
    ;; for each (poly . pow) do result= pow*r*poly'/poly + result
    (mapc #'ratdiff1 (rat-numerator r))
    (setq *sign* -1)
    (mapc #'ratdiff1 (rat-denominator r))
    result)))

(defun ratdiff1 (pr);

  (let ((poly (car pr))
	(pow (* (cdr pr) *sign*))) ;; takes care of numerator/denominator
    (cond ((coefp poly) nil)
	  ((pfreevar poly *varnum*) nil)
	  (t 
	   ;; for each (poly . pow) do result= pow*r*poly'/poly + result
	   (setq result
		 (rat+ result
		       (rat* *r* (rat* (poly2rat pow 1)
				       (rat* (ratdiff2 poly)
					     (poly2rat poly -1))))))))))

;; ratdiff2 takes a polynomial and returns a rational form
;; which is its derivative

(defun ratdiff2 (p)
  (if (coefp p) (return-from ratdiff2 (poly2rat 0 1)))
  (let ((dp (gethash (svref p 0) difftab))
	(ans (poly2rat 0 1))
	(x (vector (svref p 0) 0 1))) ;;mainvar as poly
 ;; (format t "~%dp=~s,x=~s" dp x)
 (cond ((fpe-coefzero-p (rat-numerator dp))
	;; main var of poly is independent of *var*

	   (do ((i (1- (length p))(1- i)))
	       ((= i 0) ans)
	       (setq ans (rat+ ans (rat* (poly2rat x (1- i))
					 (ratdiff2 (svref p i)))))))
	  (t  ;; main var of poly depends on *var*
	   (do ((i (1- (length p))(1- i)))
	       ((= i 0) ans)
	       (cond ((coefzerop (svref p i))) ; add zero term.. do nuthin.
		     ;; given c*p^n, set
		     ;;ans := ans + n*c*p^(n-1)*dp + dc *p^n
		     (t 
		      (setq ans
			    (rat+ ans
				  (rat+
				   (rat* (ratdiff2 (svref p i))
					 (poly2rat x (1- i)))
				   (rat* (poly2rat(1- i) 1)
					 (rat* (poly2rat (svref p i)1)
					       (rat* (poly2rat x (- i 2))
						     dp)))))))))))))

;;; set up a table of derivatives .. or here we are
;;; using atom property lists.  Is that fair?
;;; maybe a hash table would be better?

;;; this computes derivative wrt 1st and only argument..
;;; how can we store the derivative wrt nth argument?

(setf (get 'Exp 'deriv) '(Exp %)) ;; or should it be (Power E %)???
(setf (get 'Log 'deriv) '(Power % -1))
(setf (get 'Sin 'deriv) '(Cos %))
(setf (get 'Cos 'deriv) '(Times -1 (Sin %)))
(setf (get 'Tan 'deriv) '(Power (Sec %) 2))
(setf (get 'Sec 'deriv) '(Times (Sec %) (Tan %)))
;;... etc rest of Trig functions

(setf (get 'Sinh 'deriv) '(Cosh %))
;;... etc rest of Hyperbolic functions

(setf (get 'ArcSin 'deriv) '(Power (Plus -1 (Power % 2)) -1/2))
(setf (get 'ArcCos 'deriv) '(Times -1 (Power (Plus 1 (Times -1 (Power % 2)))
					     -1/2)))
(setf (get 'ArcTan 'deriv) '(Power (Plus 1 (Power % 2)) -1))
(setf (get 'ArcSec 'deriv) '(Times (Power (Plus 1 (Times -1 (Power % -2)))
					  -1/2)
				   (Power x -2)))
;;... etc rest of ArcTrig functions

(setf (get 'ArcSinh 'deriv) '(Power (Plus 1 (Power % 2)) -1/2))
(setf (get 'ArcCosh 'deriv)
      '(Times (Power (Times (Plus -1 %)(Power (Plus 1 %) -1)) -1/2)
	      (Power (Plus 1 %) -1)))
(setf (get 'ArcTanh 'deriv) '(Power (Plus 1 (Times -1 (Power % 2))) -1))
(setf (get 'ArcSech 'deriv)
      '(Times -1 (Power % -1)
	      (Power (Times (Plus 1 (Times -1 %))
			    (Power (Plus 1 %) -1))
		     -1/2)
	      (Power (Plus 1 x) -1)))

;;... etc could add the rest of ArcHyperbolic functions ..Cot,Csc


;; odds and ends: Erf, ExpIntegralEi, Abs

(setf (get 'Erf 'deriv) '(Times 2
			       (Power
				(Times (Exp (Power % 2)) (Power Pi 1/2)) -1)))

(setf (get 'ExpIntegralEi 'deriv) '(Times (Exp %) (Power % -1)))

;; line below is perhaps a problem when x=0..
	      
(setf (get 'Abs 'deriv) '(Times (Abs %)(Power x -1)))

;;;; integration properties
(setf (get 'Log 'integ) '(Times % (Plus -1 (Log %))))
(setf (get 'Sin 'integ) '(Times -1 (Cos %)))
(setf (get 'Cos 'integ) '(Sin %))
(setf (get 'Tan 'integ) '(Times -1 (Log (Cos %))))
(setf (get 'Sec 'integ) '(Times 2 (ArcTanh(Tan (Times 1/2 %)))))
(setf (get 'Cot 'integ) '(Times -1 (Log (Sin %))))
;;; etc
(setf (get 'ArcSin 'integ)
      '(Plus (Times % (ArcSin %)) 
	     (Power (Plus 1 (Times -1 (Power % 2))) 1/2)))
(setf (get 'ArcCos 'integ)
      '(Plus (Times % (ArcCos %)) 
	     (Times -1 (Power (Plus 1 (Times -1 (Power % 2))) 1/2))))
(setf (get 'ArcTan 'integ)
      '(Plus (Times % (ArcTan %))(Times -1/2 (Log (Plus 1 (Power % 2))))))
;;; etc
(setf (get 'ArcTanh 'integ)
      '(Plus (Times % (ArcTanh %))(Times 1/2 (Log (Plus -1 (Power % 2))))))

(setf (get 'Exp 'integ) '(Exp %))
;;
;;Here's some more odds and ends

(setf (get 'ExpIntegralEi 'integ) '(Plus (Times % (ExpIntegralEi %))
					 (Times -1 (Exp %))))


;;int[ erf(_X) ] := _X * erf(_X) + Pi^(-1/2) * exp(-_X^2)

(setf (get 'Erf 'integ) '(Plus (Times % (Erf %))
			       (Power
				(Times (Exp (Power % 2)) (Power Pi 1/2)) -1)))

(setf (get 'Abs 'integ) '(Times 1/2 % (Abs %)))


;; Well, we have to do some non-rational differentiation, and this
;; program is the one that does it. 
;; Restrictions that remain: it doesn't grok derivatives of
;; symbolic Derivatives.  That is, (gendiff '(((Derivative 1) f) x) 'x)
;; should mean something, namely (((Derivative 2) f) x)

(defun gendiff (h v)
  (cond ((eq h v) 1)
	((freevar h v) 0)
	((not (consp h)) (error "Gendiff of ~s" h))
	((member (car h) '(Plus Times) :test #'eq)
	 (outof-rat (ratdiff (into-rat h) v)))
	((eq (car h) 'Power)
	 (if (integerp (caddr h))
	     (outof-rat (ratdiff (into-rat h) v))
	   (ged h v))) 
	
	;; fake a derivative if necessary
	(t (let(k)
	     (setq k (if (atom (car h))
			 (get (car h) 'deriv)
		       nil)) ;; some kind of operator head like (((Derivative.)))
		   
;;	      (format t "~% deriv=~s h=~s" k h)
	     ;; Note that Mathematica (tm) uses the notation
	     ;; equivalent to (((Derivative 1) f) x) for f'[x].
	     ;; This allows for the handling of
	     ;; f(u(x),v(x)) ...
	     (cond
		   
		  ((null k);;unknown deriv. Use chain rule
		    (do ((i (cdr h)(cdr i))
			 (dlist '(1) (cons 0 dlist))
			 (ans nil; initialize answer
			      ;; all the other times through
			      (let ((thispart (gendiff (car i) v)))
				(if (eql 0 thispart) ans
			      (cons `(Times ,thispart
					    ,(uniq 
					      `(((Derivative ,@dlist),(car h))
						,@(cdr h))))
				    ans)))))
			 ;; termination test
			((null i) (simp (cons 'Plus ans)))
			;;do loop body is empty
			))
		    (t 
		    (simp
		     (uniq (list 'Times (subst (cadr h) '% k)
				 (gendiff (cadr h) v))))))))))

;; generalexponentdiff

(defun ged (e x)
  ;; x is not an integer  and e = (Power a b)
  (let ((a (cadr e))
	(b (caddr e)))
    (simp
     (cond((freevar b x)
	   ;; one form would be b*a^(b-1)*d(a,x)
	   ;; a better form would be b*a^b/a *d(a,x), using same kernels
	   `(Times ,b ,e (Power ,a -1) ,(gendiff a x)))
	  ((eq a 'E);; exponential
	   `(Times ,e ,(gendiff b x)))
	  (t `(Times ,e ,(gendiff `(Times ,b (Log ,a)) x)))))))
	  
      
;; This is the main derivative-divides integration program.
;; Integrate exp-in  wrt var, both in list form.  It returns
;; an answer in list form, as well.

(defun diffdiv (exp-in var) 
  (let*
      ;; factoring exp would be "most" effective, but let's just
      ;; put it into (partially) factored rational form.
      ((exp (into-rat exp-in))
       ;; bind factors to a list of all the (term . power) pairs. reverse
       ;; the sign of the powers in the denominator.
       ;; example,  x^2*y^2/(log(x)+1)^2 comes out as
       ;; ( (1 . 1) (x . 2)  (y . 2) (1 . -1) ((log(x)+1) . -2))
       
       (factors (append (rat-numerator exp)
			(mapcar #'(lambda(x)(cons (car x)(- (cdr x))))
				(rat-denominator exp))))
       ;; initialize constcoef, a rat term that does not 
       ;; involve integration variable
       (constcoef (poly2rat 1 1))
       (vnum (look-up-var var))
       (den nil)
       (ans nil)			; better init val?
       (vfactors nil))
    (do ((f factors (cdr f)))
	((null f) nil)
	;; several possibilities for (car f)  which is a (factor . power)
	(cond
	 ;; perhaps (caar f) is free of the variable var?
	 ((pfreevar (caar f)(look-up-var var))
	  ;;in which case we multiply constcoef by (caar f)^(cdar f)
	  (if (> (cdar f) 0) (setf (rat-numerator constcoef)
				   (fpe-insert (caar f) (cdar f)
					       (rat-numerator constcoef)))
	    (setf (rat-denominator constcoef)
		  (fpe-insert (caar f) (- (cdar f))
			      (rat-denominator constcoef)))))
	 (t (setq vfactors (cons (car f) vfactors)))))
    
  ;;(format t "~%constcoef=~s, vfactors=~s" (outof-rat constcoef) vfactors)
    ;; if all the factors are free of the variable of integration,
    ;; quit now with the answer.
    (if (null vfactors) 
	(return-from diffdiv
		     (outof-rat
		      (rat* constcoef(into-rat var)))))
    ;; vfactors is a list of the (factor . power) pairs containing var.
    (setq exp (rat/ exp constcoef))
    ;; Next, let's do some quick checks.
    ;; if the integrand is just a polynomial, we can do it
    ;; another way. Here's how...
     (if (and
	  ;; only the variable itself depends on var
	  (every
	   #'(lambda(h) (or (eql vnum h) 
			    (freevar (gethash h revtab) var))) 
	 (collectvars exp))
	  ;; and the denominator of exp is free of v entirely
	  (rfreevar (setq den (make-rat :numerator '((1 . 1))
					:denominator (rat-denominator exp)))
		    var))
	 (return-from 
	  diffdiv 
	  (outof-rat 
	   (rat* constcoef
		 (rat* den
		       (pintegrate (fpe-expand (rat-numerator exp))
				   vnum))))))

     ;; ok, we do not have a polynomial in var.  We could check here for
     ;; a simple rational function in var by the same kind of
     ;; check on the denominator (is it a simple poly in var also?)
     ;; But maybe derivative-divides will work on it, anyway..
     (do ((f vfactors (cdr f)))
	;; if we've exhausted vfactors unsuccessfully, return the "input".
	;; this is perhaps not what is wanted -- we could provide, for
	;; example, an error message, or we could call another routine.
	 ((null f) (return (ulist 'integrate exp-in var)))

	 (let* 
	     ((y (caar f))
	      (ypow (cdar f))
	      (yprime (ratdiff (poly2rat y 1) var)))
;;	  (format t "~% y=~s, ypow=~s yprime=~s" y ypow yprime)
;;	  (format t "~% y=~s, ypow=~s yprime=~s" (intol y) ypow (outof-rat yprime))
	  (setq ans (rat/ exp (rat* 
			       (poly2rat y ypow)
			       yprime)))
	  ;;(format t "~%ans=~s" (outof-rat ans))
	  (if (rfreevar ans var);; check if y*y' or y^n*y' divides
	      ;;; AHA -- WE've GOT IT!
	      (cond
	       ;; case of f=power[n], n not -1. we have y^n*y'
	       ((not(eql ypow -1))
		(return-from diffdiv
			     ;; const* (y^(p+1))/(p+1)
			     (outof-rat
			      (rat* (rat* constcoef ans)
				    (rat* (poly2rat y (1+ ypow))
					  (poly2rat (1+ ypow) -1))))))
	       
	       ((= ypow -1);; we have y^-1*y' -> log(y)
		(return-from 
		 diffdiv
		 ;; const* log(y)
		 (outof-rat(rat* constcoef
				 (rat* ans
				       (into-rat (ulist 'Log 
							(outof-rat (poly2rat y 1)))))))

		 ))
)
	    
	    ;; rfreevar test fails on y^n*y'.  Try extracting the head of
	    ;; y, that is, y=f(u), and look for f(u)*u'  (if we know
	    ;; an antiderivative for f, anyway.)
	    (cond 
	     ((not(= ypow 1))nil)
	  (t 
	     (let* ((ylist (intol y));; list form
		    
		    ;; if ylist has as head, a function of 1 variable 
		    ;; e.g. Sin
		    (head nil)
		    (antideriv nil)
		    (arg nil))
	       (cond ((and (consp ylist)
			   (atom (setq head (car ylist))))
		      ;;usual case
		      (setq antideriv (get head 'integ))
		      (setq arg (cadr ylist))))
;;	       (format t "~% y=~s, antideriv=~s yprime=~s" ylist antideriv		       (gendiff arg var))
	       ;; another case is (((Derivative list) fun) x1 ...)
	       ;; We haven't programmed that yet .. How to look up
	       ;; the antiderivative should not be toooo hard.
		     
	       (cond((atom ylist) nil);; if ylist is an atom, deriv is 0 or 1
		    ((or (cddr ylist)(null antideriv)) 
		     ;; we don't know antiderivative 
		     ;; or f has more than one argument
		     nil)
		    ;; check the exact division situation:
		    ((rfreevar
		      (setq ans (rat/ exp (rat* (poly2rat y 1)
						(into-rat (gendiff arg var)))))
		      var)
		     ;; AHA -- WE'VE GOT IT!

	       (return-from diffdiv
			    (outof-rat
			     (rat* (into-rat (subst (cadr ylist) '% antideriv))
				   (rat* ans constcoef))))))))))))))

;;Derivative divides cannot integrate polynomials, in general. This
;; program (Rint) can.  It leaves the answer in factored form, which
;; is sometimes neat, but sometimes distracting.

;;The task: Integrate a polynomial over the integers. 

;;The ploy is as follows:
;; We want to do it as much as possible as a polynomial over the integers,
;; so we can just do x^2 -> x^3/3 which requires rational numbers.
;; But we can accumulate denominators separately.  Then consider, for example,
;; integrating  9 + x + 3*x^2 + 7*x^3 - 8*x^4. 
;; Since the poly is of degree 4 set the denominator to 5!= 120.

;;Let g(k) := 5!/(k+1).  That is g(4) =24, g(3)=30, ... g(0) =120.

;;Then the answer is 1/120 times 
;;   9*g(0)*x + g(1)*x^2 +3*g(2)*x^3+7*g(3)*x^4-8*g(4)*x^5.
;; Now note that we can factor out the common factor of x, so the answer is
;; (x/120)* (9*g(0) +g(1)*x +3*g(2)*x^2 +7*g(3)*x^3-8*g(4)*x^4).
;; or after removing common factors, 
;;                                   2       3       4
;;Out[24] = 1/20 x (180 + 10 x + 20 x  + 35 x  - 32 x )

;; What about  (1+x) + (1+x+x^2)*y , integrated wrt x?  Let the highest
;; degree of x ANYWHERE be the denominator.  The analysis still works.

;;; rint makes the assumption that other kernels in the numerator
;;; are free of the variable of integration.

(defun rint (r v)
  (let ((den (make-rat :numerator '((1 . 1))
		       :denominator (rat-denominator r))))
    (if (rfreevar den v)
	(rat* (pintegrate (fpe-expand (rat-numerator r)) (look-up-var v))
	      den)
      (error "~%Denominator ~s is not free of ~s" (outof-rat den) v) )))

(defun |Rint|(r v)(outof-rat(rint (into-rat r)v)))

;; pintegrate takes a polynomial and a variable but returns a
;; RATIONAL form

(defun pintegrate (p v &aux maxfact (maxdegv 0))
  (labels
   
   ((g (k) (/ maxfact k))
    
    (factorial(k)(if (= k 0) 1 (* k (factorial (1- k)))))
    
    ;; set maxdegv to maximum degree that v occurs in p. return value nil.
    (maxdegree (p)
	       (cond((coefp p))
		    ((var> v (svref p 0)))
		    ((eql (svref p 0) v)
		     (setq maxdegv (max maxdegv (- (length p) 2))))
		    (t (do ((i (1- (length p)) (1- i)))
			   ((= i 0)) ;; returns nil. size effect to maxdegv
			   (maxdegree (svref p i))))))
    (pd1 (p &aux r)
	 (cond ((coefp p) (* p maxfact))
	       ((var> v (svref p 0)) p)
	       ((eql (svref p 0) v)

		(setq r (make-array (length p)))
		(setf (svref r 0) v)
		(do ((i (1- (length p)) (1- i)))
		    ((= i 0)
;;		     (format t "~%integrated p=~s to get r=~s" p r)
		     r)
		    (setf (svref r i)(p* (g i) (svref p i)))))
	       (t (setq r(copy-seq p))
		  (do ((i (1- (length r))(1- i)))
		      ((= i 0) r)
		      (setf (svref r i)(pd1 (svref r i))))))))
   (maxdegree p)
   (setq maxfact (factorial (1+ maxdegv)))
;; (format t "~s, ~s, ~s" (pd1 p) v maxfact)
   (rat/ (rat* (poly2rat (vector v 0 1) 1)
	       (make-rat :numerator (make-fpe (pd1 p) 1)
			 :denominator (list '(1 . 1))))
	(poly2rat maxfact 1) )))


;;; Notes for further heuristics.
;;  for symbolic n, make Int[x^n,x] into (x^(n+1)-1)/(n+1).
;; This has a different constant from the usual, but the nice property that
;; if n-> -1, the answer is Log[x]. (suggested by WK 12/90)

;; similar stuff possible for Cos[n*x]*Sin[m*n] formulas when n=+-m.