;;;; Symbolic Integrator, First Stage
;;;   based on SIN (Joel Moses, MIT PhD Thesis, 1967).

;;; Joel Moses: "Symbolic Integration, The Stormy Decade"
;;;  In CACM, Vol.14, No.8, August 1971.

;;; This code must be run under Scmutils (Mechanics system).

;;; More to do.  Only first stage is implemented.
;;; Second stage, integration by parts

;;; in Scmutils
;;; convert sqrt->expt $ 1/2

(define (integrate integrand variable)
  (let ((result
	 (let intlp ((expr (new-simplify integrand)))
	   (cond ((indep-of? expr variable)
		  (symb:* expr variable))
		 ((sum? expr)
		  (let ((subexprs (map intlp (operands expr))))
		    (and (&and subexprs)
			 (symb:add:n subexprs))))
		 ((quotient? expr)
		  (cond ((sum? (symb:numerator expr))
			 (let ((d (symb:denominator expr)))
			   (let ((subexprs
				  (map (lambda (n)
					 (intlp
					  (new-simplify
					   (symb:/ num d))))
				       (operands
					(symb:numerator expr)))))
			     (and (&and subexprs)
				  (symb:add:n subexprs)))))
			(else (sin-stage-1 expr variable))))
		 (else
		  (sin-stage-1 expr variable))))))
    (if (and result *debugging-integrate*)
	;; Paranoid Programming
	(let ((check (derivative result variable)))
	  (if (not (zero? (new-simplify (symb:- integrand check)))) 
	      (pp `(derivative-differs ,check ,integrand)))))
    (and result (symb:+ result (new-constant)))))

(define *debugging-integrate* #f)


(define (sin-stage-1 integrand variable)
  (derivative-divides integrand variable))

(define (derivative-divides expression x)
  (define (elementary-finish results)
    (and *debugging-derivative-divides*
	 (pp `(elementary ,@results)))
    (and results
	 (let ((factors (new-simplify (car results)))
	       (op (cadr results))
	       (u (caddr results)))
	   (let ((du/dx (derivative u x)))
	     (let ((c (new-simplify (symb:/ factors du/dx))))
	       (if (indep-of? c x)
		   (let ((intv (elementary? op)))
		     (cond (intv
			    (new-simplify
			     (symb:* c (substitute u '$ (cadr intv)))))
			   (else (error "Should not get here"))))
		   #f))))))
  (define (common-finish results)
    (and *debugging-derivative-divides*
	 (pp `(common ,@results)))
    (and results
	 (let ((factors (new-simplify (car results)))
	       (essence (cadr results))
	       (u (caddr results)))
	   (let ((du/dx (derivative u x)))
	     (let ((c (new-simplify (symb:/ factors du/dx))))
	       (if (indep-of? c x)
		   (new-simplify (symb:* c essence))
		   #f))))))
  (or (elementary-match `(,expression ,x) elementary-finish)
      (common-patterns-match `(,expression ,x) common-finish)
      ))

(define *debugging-derivative-divides* #f)

;;; Tables... numbers refer to entries in Gradshteyn and Ryzhik

(define (elementary? op)
  (assq op integrals-of-elementary-functions))

(define elementary-match
  (make-match
   (list				; 3,5,6
    (rule ((// (** ((? op elementary?) (? u))
		   #!rest
		   (? n))
	       (? d))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (? op)
	   (? u)))

    (rule ((// (** (? u) #!rest (? n)) (? d))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   identity
	   (? u)))
    )))

(define integrals-of-elementary-functions
  '(
    (identity (/ (expt $ 2) 2))
    (exp (exp $))
    (sin (- (cos $)))
    (cos (sin $))
    (tan (- (log (cos $))))
    (cot (log (sin $)))
    (sec (- (log (+ (cos (/ $ 2)) (sin (/ $ 2))))
	    (log (- (cos (/ $ 2)) (sin (/ $ 2))))))
    (csc (- (log (sin (/ $ 2))) (log (cos (/ $ 2)))))
    (asin (+ (* $ (asin $)) (sqrt (- 1 (expt $ 2)))))
    (acos (- (* $ (acos $)) (sqrt (- 1 (expt $ 2)))))
    (atan (- (* $ (atan $)) (* 1/2 (log (+ 1 (expt $ 2))))))
    (acot (+ (* $ (acot $)) (* 1/2 (log (+ 1 (expt $ 2))))))
    (asec (- (* $ (asec $))
	     (log (* $ (+ 1 (sqrt (/ (- (expt $ 2) 1) (expt $ 2))))))))
    (acsc (+ (* $ (acsc $))
	     (log (* $ (+ 1 (sqrt (/ (- (expt $ 2) 1) (expt $ 2))))))))
    (log (* $ (- (log $) 1)))
    (cosh (sinh $))
    (sinh (cosh $))
    (tanh (log (cosh $)))
    (asinh (- (* $ (asinh $)) (sqrt (+ 1 (expt $ 2)))))
    (acosh (- (* $ (acosh $)) (* (sqrt (+ $ 1)) (sqrt (- $ 1)))))
    (atanh (+ (* $ (atanh $)) (* 1/2 (log (- (expt $ 2) 1)))))
    ))

(define common-patterns-match
  (make-match
   (list
    (rule ((// (? n)			;2
	       (** (? u) #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (log (? u))
	   (? u)))

    (rule ((// (** (^^ (? u) (? n))	;1
		   #!rest (? f))
	       (? d))
	   (? x))
	  (and (depends-on? u x)
	       (indep-of? n x)
	       (not (equal? n -1)))
	  ((/ (? f) (? d))
	   (/ (expt (? u) (? (+ n 1)))
	      (? (+ n 1)))
	   (? u)))

    (rule ((// (** (^^ (? b) (? u))	;4
		   #!rest (? n))
	       (? d))
	   (? x))
	  (and (depends-on? u x) (indep-of? b x))
	  ((/ (? n) (? d))
	   (/ (expt (? b) (? u)) (log (? b)))
	   (? u)))

    (rule ((// (? n)			;7
	      (** (^^ (sin (? u)) 2)
		  #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (- (/ (cos (? u)) (sin (? u))))
	   (? u)))

    (rule ((// (? n)			;8
	       (** (expt (cos (? u)) 2)
		   #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (/ (sin (? u)) (cos (? u)))
	   (? u)))

    (rule ((// (** (sin (? u)) 		;9
		   #!rest (? n))
	       (** (expt (cos (? u)) 2)
		   #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (/ 1 (cos (? u)))
	   (? u)))

    (rule ((// (** (cos (? u))		;10
		   #!rest (? n))
	       (** (expt (sin (? u)) 2)
		   #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (/ -1 (sin (? u)))
	   (? u)))

    (rule ((// (** (sin (? u))		;11
		   #!rest (? n))
	       (** (cos (? u))
		   #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (- (log (cos (? u))))
	   (? u)))

    (rule ((// (** (cos (? u))		;12
		   #!rest (? n))
	       (** (sin (? u))
		   #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (log (sin (? u)))
	   (? u)))

    (rule ((// (? n)			;13
	      (** (sin (? u))
		  #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (log (/ (sin (/ (? u) 2)) (cos (/ (? u) 2))))
	   (? u)))

    (rule ((// (? n)			;14
	      (** (cos (? u))
		  #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (log
	    (+ (/ (sin (/ (? u) 2)) (cos (/ (? u) 2)))
	       (/ 1 (cos (/ (? u) 2)))))
	   (? u)))

    (rule ((// (? n)			;15
	       (** (++ 1 (^^ (? u) 2))
		   #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (atan (? u))
	   (? u)))

    (rule ((// (? n)			;16
	      (** (++ 1 (** -1 (^^ (? u) 2)))
		  #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (* 1/2 (log (/ (+ 1 (? u)) (- 1 (? u)))))
	   (? u)))

    (rule ((// (? n)			;17
	       (** (^^ (++ 1 (** -1 (^^ (? u) 2))) 1/2)
		   #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (asin (? u))
	   (? u)))

    (rule ((// (? n)			;18
	      (** (^^ (++ 1 (^^ (? u) 2)) 1/2)
		  #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (log (+ (? u) (sqrt (+ 1 (expt (? u) 2)))))
	   (? u)))

    (rule ((// (? n)			;19
	       (** (^^ (++ (^^ (? u) 2) -1) 1/2)
		   #!rest (? d)))
	   (? x))
	  (depends-on? u x)
	  ((/ (? n) (? d))
	   (log (+ (? u) (sqrt (- (expt (? u) 2) 1))))
	   (? u)))

    )))

#|
(pp (integrate '(+ (* a (expt x 2)) (* b x) c) 'x))
(+ (* 1/3 a (expt x 3)) (* 1/2 b (expt x 2)) (* c x) C30)

(pe (integrate '(* x (exp (expt x 2))) 'x))
(+ C7 (* 1/2 (exp (expt x 2))))

(pe (integrate '(* 4 (cos (+ (* 2 x) 3))) 'x))
(+ C17 (* 2 (sin (+ 3 (* 2 x)))))

(pe (integrate '(/ (exp x) (+ 1 (exp x))) 'x))
(+ C20 (log (+ 1 (exp x))))

(pe (integrate '(* (sin x) (cos x)) 'x))
(+ (* -1/2 (expt (cos x) 2)) C8)

(pe (integrate '(* (cos x) (sin x)) 'x))
(+ (* -1/2 (expt (cos x) 2)) C11)

(pe (integrate '(* x (expt (+ (expt x 2) 1) 1/2)) 'x))
(+ C43 (* 1/3 (expt (+ 1 (expt x 2)) 3/2)))

(pe (integrate
     '(* (expt (cos (exp x)) 2) (sin (exp x)) (exp x))
     'x))
(+ (* -1/3 (expt (cos (exp x)) 3)) C61)

|#

#|
;;; A problem
(pe (integrate '(expt (exp x) 2) 'x))
#f

;;; But!
(pe (derivative-divides '(* (exp x) (exp x)) 'x))
(* 1/2 (expt (exp x) 2))

;;; and
(pe (derivative-divides '(expt e (* 2 x)) 'x))
(/ (* 1/2 (expt e (* 2 x))) (log e))


(pe (derivative-divides '(exp (* 2 x)) 'x))
(* 1/2 (expt (exp x) 2))

(pe (integrate '(exp (* 2 x)) 'x))
#f

;;; Ahhh.  Simplifier makes this hard to recognize
(pe '(exp (* 2 x)))
(expt (exp x) 2)

;;; Not clear how to handle this...


;;; For second stage...  by parts...  This can blow up.
(pe (derivative-divides '(* 2 x (exp x)) 'x))
#f
|#