view lisp/calc/calcalg2.el @ 56905:661d52db56de

(isearch-toggle-regexp): Set `isearch-success' and `isearch-adjusted' to `t'. (isearch-toggle-case-fold): Set `isearch-success' to `t'. (isearch-message-prefix): Add "pending" for isearch-adjusted. (isearch-other-meta-char): Restore isearch-point unconditionally. (isearch-query-replace): Add new arg `regexp-flag' and use it. Set point to start of match if region is not active in transient mark mode (to include the current match to region boundaries). Push the search string to `query-replace-from-history-variable'. Add prompt "Query replace regexp" for isearch-regexp. Add region beginning/end as last arguments of `perform-replace.' (isearch-query-replace-regexp): Replace code by the call to `isearch-query-replace' with arg `t'.

;;; calcalg2.el --- more algebraic functions for Calc

;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.

;; Author: David Gillespie <[e-mail]>
;; Maintainers: D. Goel <[e-mail]>
;;              Colin Walters <[e-mail]>

;; This file is part of GNU Emacs.

;; GNU Emacs is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY.  No author or distributor
;; accepts responsibility to anyone for the consequences of using it
;; or for whether it serves any particular purpose or works at all,
;; unless he says so in writing.  Refer to the GNU Emacs General Public
;; License for full details.

;; Everyone is granted permission to copy, modify and redistribute
;; GNU Emacs, but only under the conditions described in the
;; GNU Emacs General Public License.   A copy of this license is
;; supposed to have been given to you along with GNU Emacs so you
;; can know your rights and responsibilities.  It should be in a
;; file named COPYING.  Among other things, the copyright notice
;; and this notice must be preserved on all copies.

;;; Commentary:

;;; Code:

;; This file is autoloaded from calc-ext.el.
(require 'calc-ext)

(require 'calc-macs)

(defun calc-Need-calc-alg-2 () nil)


(defun calc-derivative (var num)
  (interactive "sDifferentiate with respect to: \np")
  (calc-slow-wrapper
   (when (< num 0)
     (error "Order of derivative must be positive"))
   (let ((func (if (calc-is-hyperbolic) 'calcFunc-tderiv 'calcFunc-deriv))
	 n expr)
     (if (or (equal var "") (equal var "$"))
	 (setq n 2
	       expr (calc-top-n 2)
	       var (calc-top-n 1))
       (setq var (math-read-expr var))
       (when (eq (car-safe var) 'error)
	 (error "Bad format in expression: %s" (nth 1 var)))
       (setq n 1
	     expr (calc-top-n 1)))
     (while (>= (setq num (1- num)) 0)
       (setq expr (list func expr var)))
     (calc-enter-result n "derv" expr))))

(defun calc-integral (var)
  (interactive "sIntegration variable: ")
  (calc-slow-wrapper
   (if (or (equal var "") (equal var "$"))
       (calc-enter-result 2 "intg" (list 'calcFunc-integ
					 (calc-top-n 2)
					 (calc-top-n 1)))
     (let ((var (math-read-expr var)))
       (if (eq (car-safe var) 'error)
	   (error "Bad format in expression: %s" (nth 1 var)))
       (calc-enter-result 1 "intg" (list 'calcFunc-integ
					 (calc-top-n 1)
					 var))))))

(defun calc-num-integral (&optional varname lowname highname)
  (interactive "sIntegration variable: ")
  (calc-tabular-command 'calcFunc-ninteg "Integration" "nint"
			nil varname lowname highname))

(defun calc-summation (arg &optional varname lowname highname)
  (interactive "P\nsSummation variable: ")
  (calc-tabular-command 'calcFunc-sum "Summation" "sum"
			arg varname lowname highname))

(defun calc-alt-summation (arg &optional varname lowname highname)
  (interactive "P\nsSummation variable: ")
  (calc-tabular-command 'calcFunc-asum "Summation" "asum"
			arg varname lowname highname))

(defun calc-product (arg &optional varname lowname highname)
  (interactive "P\nsIndex variable: ")
  (calc-tabular-command 'calcFunc-prod "Index" "prod"
			arg varname lowname highname))

(defun calc-tabulate (arg &optional varname lowname highname)
  (interactive "P\nsIndex variable: ")
  (calc-tabular-command 'calcFunc-table "Index" "tabl"
			arg varname lowname highname))

(defun calc-tabular-command (func prompt prefix arg varname lowname highname)
  (calc-slow-wrapper
   (let (var (low nil) (high nil) (step nil) stepname stepnum (num 1) expr)
     (if (consp arg)
	 (setq stepnum 1)
       (setq stepnum 0))
     (if (or (equal varname "") (equal varname "$") (null varname))
	 (setq high (calc-top-n (+ stepnum 1))
	       low (calc-top-n (+ stepnum 2))
	       var (calc-top-n (+ stepnum 3))
	       num (+ stepnum 4))
       (setq var (if (stringp varname) (math-read-expr varname) varname))
       (if (eq (car-safe var) 'error)
	   (error "Bad format in expression: %s" (nth 1 var)))
       (or lowname
	   (setq lowname (read-string (concat prompt " variable: " varname
					      ", from: "))))
       (if (or (equal lowname "") (equal lowname "$"))
	   (setq high (calc-top-n (+ stepnum 1))
		 low (calc-top-n (+ stepnum 2))
		 num (+ stepnum 3))
	 (setq low (if (stringp lowname) (math-read-expr lowname) lowname))
	 (if (eq (car-safe low) 'error)
	     (error "Bad format in expression: %s" (nth 1 low)))
	 (or highname
	     (setq highname (read-string (concat prompt " variable: " varname
						 ", from: " lowname
						 ", to: "))))
	 (if (or (equal highname "") (equal highname "$"))
	     (setq high (calc-top-n (+ stepnum 1))
		   num (+ stepnum 2))
	   (setq high (if (stringp highname) (math-read-expr highname)
			highname))
	   (if (eq (car-safe high) 'error)
	       (error "Bad format in expression: %s" (nth 1 high)))
	   (if (consp arg)
	       (progn
		 (setq stepname (read-string (concat prompt " variable: "
						     varname
						     ", from: " lowname
						     ", to: " highname
						     ", step: ")))
		 (if (or (equal stepname "") (equal stepname "$"))
		     (setq step (calc-top-n 1)
			   num 2)
		   (setq step (math-read-expr stepname))
		   (if (eq (car-safe step) 'error)
		       (error "Bad format in expression: %s"
			      (nth 1 step)))))))))
     (or step
	 (if (consp arg)
	     (setq step (calc-top-n 1))
	   (if arg
	       (setq step (prefix-numeric-value arg)))))
     (setq expr (calc-top-n num))
     (calc-enter-result num prefix (append (list func expr var low high)
					   (and step (list step)))))))

(defun calc-solve-for (var)
  (interactive "sVariable to solve for: ")
  (calc-slow-wrapper
   (let ((func (if (calc-is-inverse)
		   (if (calc-is-hyperbolic) 'calcFunc-ffinv 'calcFunc-finv)
		 (if (calc-is-hyperbolic) 'calcFunc-fsolve 'calcFunc-solve))))
     (if (or (equal var "") (equal var "$"))
	 (calc-enter-result 2 "solv" (list func
					   (calc-top-n 2)
					   (calc-top-n 1)))
       (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
			   (not (string-match "\\[" var)))
		      (math-read-expr (concat "[" var "]"))
		    (math-read-expr var))))
	 (if (eq (car-safe var) 'error)
	     (error "Bad format in expression: %s" (nth 1 var)))
	 (calc-enter-result 1 "solv" (list func
					   (calc-top-n 1)
					   var)))))))

(defun calc-poly-roots (var)
  (interactive "sVariable to solve for: ")
  (calc-slow-wrapper
   (if (or (equal var "") (equal var "$"))
       (calc-enter-result 2 "prts" (list 'calcFunc-roots
					 (calc-top-n 2)
					 (calc-top-n 1)))
     (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
			 (not (string-match "\\[" var)))
		    (math-read-expr (concat "[" var "]"))
		  (math-read-expr var))))
       (if (eq (car-safe var) 'error)
	   (error "Bad format in expression: %s" (nth 1 var)))
       (calc-enter-result 1 "prts" (list 'calcFunc-roots
					 (calc-top-n 1)
					 var))))))

(defun calc-taylor (var nterms)
  (interactive "sTaylor expansion variable: \nNNumber of terms: ")
  (calc-slow-wrapper
   (let ((var (math-read-expr var)))
     (if (eq (car-safe var) 'error)
	 (error "Bad format in expression: %s" (nth 1 var)))
     (calc-enter-result 1 "tylr" (list 'calcFunc-taylor
				       (calc-top-n 1)
				       var
				       (prefix-numeric-value nterms))))))


(defun math-derivative (expr)   ; uses global values: deriv-var, deriv-total.
  (cond ((equal expr deriv-var)
	 1)
	((or (Math-scalarp expr)
	     (eq (car expr) 'sdev)
	     (and (eq (car expr) 'var)
		  (or (not deriv-total)
		      (math-const-var expr)
		      (progn
			(math-setup-declarations)
			(memq 'const (nth 1 (or (assq (nth 2 expr)
						      math-decls-cache)
						math-decls-all)))))))
	 0)
	((eq (car expr) '+)
	 (math-add (math-derivative (nth 1 expr))
		   (math-derivative (nth 2 expr))))
	((eq (car expr) '-)
	 (math-sub (math-derivative (nth 1 expr))
		   (math-derivative (nth 2 expr))))
	((memq (car expr) '(calcFunc-eq calcFunc-neq calcFunc-lt
					calcFunc-gt calcFunc-leq calcFunc-geq))
	 (list (car expr)
	       (math-derivative (nth 1 expr))
	       (math-derivative (nth 2 expr))))
	((eq (car expr) 'neg)
	 (math-neg (math-derivative (nth 1 expr))))
	((eq (car expr) '*)
	 (math-add (math-mul (nth 2 expr)
			     (math-derivative (nth 1 expr)))
		   (math-mul (nth 1 expr)
			     (math-derivative (nth 2 expr)))))
	((eq (car expr) '/)
	 (math-sub (math-div (math-derivative (nth 1 expr))
			     (nth 2 expr))
		   (math-div (math-mul (nth 1 expr)
				       (math-derivative (nth 2 expr)))
			     (math-sqr (nth 2 expr)))))
	((eq (car expr) '^)
	 (let ((du (math-derivative (nth 1 expr)))
	       (dv (math-derivative (nth 2 expr))))
	   (or (Math-zerop du)
	       (setq du (math-mul (nth 2 expr)
				  (math-mul (math-normalize
					     (list '^
						   (nth 1 expr)
						   (math-add (nth 2 expr) -1)))
					    du))))
	   (or (Math-zerop dv)
	       (setq dv (math-mul (math-normalize
				   (list 'calcFunc-ln (nth 1 expr)))
				  (math-mul expr dv))))
	   (math-add du dv)))
	((eq (car expr) '%)
	 (math-derivative (nth 1 expr)))   ; a reasonable definition
	((eq (car expr) 'vec)
	 (math-map-vec 'math-derivative expr))
	((and (memq (car expr) '(calcFunc-conj calcFunc-re calcFunc-im))
	      (= (length expr) 2))
	 (list (car expr) (math-derivative (nth 1 expr))))
	((and (memq (car expr) '(calcFunc-subscr calcFunc-mrow calcFunc-mcol))
	      (= (length expr) 3))
	 (let ((d (math-derivative (nth 1 expr))))
	   (if (math-numberp d)
	       0    ; assume x and x_1 are independent vars
	     (list (car expr) d (nth 2 expr)))))
	(t (or (and (symbolp (car expr))
		    (if (= (length expr) 2)
			(let ((handler (get (car expr) 'math-derivative)))
			  (and handler
			       (let ((deriv (math-derivative (nth 1 expr))))
				 (if (Math-zerop deriv)
				     deriv
				   (math-mul (funcall handler (nth 1 expr))
					     deriv)))))
		      (let ((handler (get (car expr) 'math-derivative-n)))
			(and handler
			     (funcall handler expr)))))
	       (and (not (eq deriv-symb 'pre-expand))
		    (let ((exp (math-expand-formula expr)))
		      (and exp
			   (or (let ((deriv-symb 'pre-expand))
				 (catch 'math-deriv (math-derivative expr)))
			       (math-derivative exp)))))
	       (if (or (Math-objvecp expr)
		       (eq (car expr) 'var)
		       (not (symbolp (car expr))))
		   (if deriv-symb
		       (throw 'math-deriv nil)
		     (list (if deriv-total 'calcFunc-tderiv 'calcFunc-deriv)
			   expr
			   deriv-var))
		 (let ((accum 0)
		       (arg expr)
		       (n 1)
		       derv)
		   (while (setq arg (cdr arg))
		     (or (Math-zerop (setq derv (math-derivative (car arg))))
			 (let ((func (intern (concat (symbol-name (car expr))
						     "'"
						     (if (> n 1)
							 (int-to-string n)
						       ""))))
			       (prop (cond ((= (length expr) 2)
					    'math-derivative-1)
					   ((= (length expr) 3)
					    'math-derivative-2)
					   ((= (length expr) 4)
					    'math-derivative-3)
					   ((= (length expr) 5)
					    'math-derivative-4)
					   ((= (length expr) 6)
					    'math-derivative-5))))
			   (setq accum
				 (math-add
				  accum
				  (math-mul
				   derv
				   (let ((handler (get func prop)))
				     (or (and prop handler
					      (apply handler (cdr expr)))
					 (if (and deriv-symb
						  (not (get func
							    'calc-user-defn)))
					     (throw 'math-deriv nil)
					   (cons func (cdr expr))))))))))
		     (setq n (1+ n)))
		   accum))))))

(defun calcFunc-deriv (expr deriv-var &optional deriv-value deriv-symb)
  (let* ((deriv-total nil)
	 (res (catch 'math-deriv (math-derivative expr))))
    (or (eq (car-safe res) 'calcFunc-deriv)
	(null res)
	(setq res (math-normalize res)))
    (and res
	 (if deriv-value
	     (math-expr-subst res deriv-var deriv-value)
	   res))))

(defun calcFunc-tderiv (expr deriv-var &optional deriv-value deriv-symb)
  (math-setup-declarations)
  (let* ((deriv-total t)
	 (res (catch 'math-deriv (math-derivative expr))))
    (or (eq (car-safe res) 'calcFunc-tderiv)
	(null res)
	(setq res (math-normalize res)))
    (and res
	 (if deriv-value
	     (math-expr-subst res deriv-var deriv-value)
	   res))))

(put 'calcFunc-inv\' 'math-derivative-1
     (function (lambda (u) (math-neg (math-div 1 (math-sqr u))))))

(put 'calcFunc-sqrt\' 'math-derivative-1
     (function (lambda (u) (math-div 1 (math-mul 2 (list 'calcFunc-sqrt u))))))

(put 'calcFunc-deg\' 'math-derivative-1
     (function (lambda (u) (math-div-float '(float 18 1) (math-pi)))))

(put 'calcFunc-rad\' 'math-derivative-1
     (function (lambda (u) (math-pi-over-180))))

(put 'calcFunc-ln\' 'math-derivative-1
     (function (lambda (u) (math-div 1 u))))

(put 'calcFunc-log10\' 'math-derivative-1
     (function (lambda (u)
		 (math-div (math-div 1 (math-normalize '(calcFunc-ln 10)))
			   u))))

(put 'calcFunc-lnp1\' 'math-derivative-1
     (function (lambda (u) (math-div 1 (math-add u 1)))))

(put 'calcFunc-log\' 'math-derivative-2
     (function (lambda (x b)
		 (and (not (Math-zerop b))
		      (let ((lnv (math-normalize
				  (list 'calcFunc-ln b))))
			(math-div 1 (math-mul lnv x)))))))

(put 'calcFunc-log\'2 'math-derivative-2
     (function (lambda (x b)
		 (let ((lnv (list 'calcFunc-ln b)))
		   (math-neg (math-div (list 'calcFunc-log x b)
				       (math-mul lnv b)))))))

(put 'calcFunc-exp\' 'math-derivative-1
     (function (lambda (u) (math-normalize (list 'calcFunc-exp u)))))

(put 'calcFunc-expm1\' 'math-derivative-1
     (function (lambda (u) (math-normalize (list 'calcFunc-expm1 u)))))

(put 'calcFunc-sin\' 'math-derivative-1
     (function (lambda (u) (math-to-radians-2 (math-normalize
					       (list 'calcFunc-cos u))))))

(put 'calcFunc-cos\' 'math-derivative-1
     (function (lambda (u) (math-neg (math-to-radians-2
				      (math-normalize
				       (list 'calcFunc-sin u)))))))

(put 'calcFunc-tan\' 'math-derivative-1
     (function (lambda (u) (math-to-radians-2
			    (math-div 1 (math-sqr
					 (math-normalize
					  (list 'calcFunc-cos u))))))))

(put 'calcFunc-arcsin\' 'math-derivative-1
     (function (lambda (u)
		 (math-from-radians-2
		  (math-div 1 (math-normalize
			       (list 'calcFunc-sqrt
				     (math-sub 1 (math-sqr u)))))))))

(put 'calcFunc-arccos\' 'math-derivative-1
     (function (lambda (u)
		 (math-from-radians-2
		  (math-div -1 (math-normalize
				(list 'calcFunc-sqrt
				      (math-sub 1 (math-sqr u)))))))))

(put 'calcFunc-arctan\' 'math-derivative-1
     (function (lambda (u) (math-from-radians-2
			    (math-div 1 (math-add 1 (math-sqr u)))))))

(put 'calcFunc-sinh\' 'math-derivative-1
     (function (lambda (u) (math-normalize (list 'calcFunc-cosh u)))))

(put 'calcFunc-cosh\' 'math-derivative-1
     (function (lambda (u) (math-normalize (list 'calcFunc-sinh u)))))

(put 'calcFunc-tanh\' 'math-derivative-1
     (function (lambda (u) (math-div 1 (math-sqr
					(math-normalize
					 (list 'calcFunc-cosh u)))))))

(put 'calcFunc-arcsinh\' 'math-derivative-1
     (function (lambda (u)
		 (math-div 1 (math-normalize
			      (list 'calcFunc-sqrt
				    (math-add (math-sqr u) 1)))))))

(put 'calcFunc-arccosh\' 'math-derivative-1
     (function (lambda (u)
		  (math-div 1 (math-normalize
			       (list 'calcFunc-sqrt
				     (math-add (math-sqr u) -1)))))))

(put 'calcFunc-arctanh\' 'math-derivative-1
     (function (lambda (u) (math-div 1 (math-sub 1 (math-sqr u))))))

(put 'calcFunc-bern\'2 'math-derivative-2
     (function (lambda (n x)
		 (math-mul n (list 'calcFunc-bern (math-add n -1) x)))))

(put 'calcFunc-euler\'2 'math-derivative-2
     (function (lambda (n x)
		 (math-mul n (list 'calcFunc-euler (math-add n -1) x)))))

(put 'calcFunc-gammag\'2 'math-derivative-2
     (function (lambda (a x) (math-deriv-gamma a x 1))))

(put 'calcFunc-gammaG\'2 'math-derivative-2
     (function (lambda (a x) (math-deriv-gamma a x -1))))

(put 'calcFunc-gammaP\'2 'math-derivative-2
     (function (lambda (a x) (math-deriv-gamma a x
					       (math-div
						1 (math-normalize
						   (list 'calcFunc-gamma
							 a)))))))

(put 'calcFunc-gammaQ\'2 'math-derivative-2
     (function (lambda (a x) (math-deriv-gamma a x
					       (math-div
						-1 (math-normalize
						    (list 'calcFunc-gamma
							  a)))))))

(defun math-deriv-gamma (a x scale)
  (math-mul scale
	    (math-mul (math-pow x (math-add a -1))
		      (list 'calcFunc-exp (math-neg x)))))

(put 'calcFunc-betaB\' 'math-derivative-3
     (function (lambda (x a b) (math-deriv-beta x a b 1))))

(put 'calcFunc-betaI\' 'math-derivative-3
     (function (lambda (x a b) (math-deriv-beta x a b
						(math-div
						 1 (list 'calcFunc-beta
							 a b))))))

(defun math-deriv-beta (x a b scale)
  (math-mul (math-mul (math-pow x (math-add a -1))
		      (math-pow (math-sub 1 x) (math-add b -1)))
	    scale))

(put 'calcFunc-erf\' 'math-derivative-1
     (function (lambda (x) (math-div 2
				     (math-mul (list 'calcFunc-exp
						     (math-sqr x))
					       (if calc-symbolic-mode
						   '(calcFunc-sqrt
						     (var pi var-pi))
						 (math-sqrt-pi)))))))

(put 'calcFunc-erfc\' 'math-derivative-1
     (function (lambda (x) (math-div -2
				     (math-mul (list 'calcFunc-exp
						     (math-sqr x))
					       (if calc-symbolic-mode
						   '(calcFunc-sqrt
						     (var pi var-pi))
						 (math-sqrt-pi)))))))

(put 'calcFunc-besJ\'2 'math-derivative-2
     (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besJ
						       (math-add v -1)
						       z)
						 (list 'calcFunc-besJ
						       (math-add v 1)
						       z))
				       2))))

(put 'calcFunc-besY\'2 'math-derivative-2
     (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besY
						       (math-add v -1)
						       z)
						 (list 'calcFunc-besY
						       (math-add v 1)
						       z))
				       2))))

(put 'calcFunc-sum 'math-derivative-n
     (function
      (lambda (expr)
	(if (math-expr-contains (cons 'vec (cdr (cdr expr))) deriv-var)
	    (throw 'math-deriv nil)
	  (cons 'calcFunc-sum
		(cons (math-derivative (nth 1 expr))
		      (cdr (cdr expr))))))))

(put 'calcFunc-prod 'math-derivative-n
     (function
      (lambda (expr)
	(if (math-expr-contains (cons 'vec (cdr (cdr expr))) deriv-var)
	    (throw 'math-deriv nil)
	  (math-mul expr
		    (cons 'calcFunc-sum
			  (cons (math-div (math-derivative (nth 1 expr))
					  (nth 1 expr))
				(cdr (cdr expr)))))))))

(put 'calcFunc-integ 'math-derivative-n
     (function
      (lambda (expr)
	(if (= (length expr) 3)
	    (if (equal (nth 2 expr) deriv-var)
		(nth 1 expr)
	      (math-normalize
	       (list 'calcFunc-integ
		     (math-derivative (nth 1 expr))
		     (nth 2 expr))))
	  (if (= (length expr) 5)
	      (let ((lower (math-expr-subst (nth 1 expr) (nth 2 expr)
					    (nth 3 expr)))
		    (upper (math-expr-subst (nth 1 expr) (nth 2 expr)
					    (nth 4 expr))))
		(math-add (math-sub (math-mul upper
					      (math-derivative (nth 4 expr)))
				    (math-mul lower
					      (math-derivative (nth 3 expr))))
			  (if (equal (nth 2 expr) deriv-var)
			      0
			    (math-normalize
			     (list 'calcFunc-integ
				   (math-derivative (nth 1 expr)) (nth 2 expr)
				   (nth 3 expr) (nth 4 expr)))))))))))

(put 'calcFunc-if 'math-derivative-n
     (function
      (lambda (expr)
	(and (= (length expr) 4)
	     (list 'calcFunc-if (nth 1 expr)
		   (math-derivative (nth 2 expr))
		   (math-derivative (nth 3 expr)))))))

(put 'calcFunc-subscr 'math-derivative-n
     (function
      (lambda (expr)
	(and (= (length expr) 3)
	     (list 'calcFunc-subscr (nth 1 expr)
		   (math-derivative (nth 2 expr)))))))


(defvar math-integ-var '(var X ---))
(defvar math-integ-var-2 '(var Y ---))
(defvar math-integ-vars (list 'f math-integ-var math-integ-var-2))
(defvar math-integ-var-list (list math-integ-var))
(defvar math-integ-var-list-list (list math-integ-var-list))

(defmacro math-tracing-integral (&rest parts)
  (list 'and
	'trace-buffer
	(list 'save-excursion
	      '(set-buffer trace-buffer)
	      '(goto-char (point-max))
	      (list 'and
		    '(bolp)
		    '(insert (make-string (- math-integral-limit
					     math-integ-level) 32)
			     (format "%2d " math-integ-depth)
			     (make-string math-integ-level 32)))
	      ;;(list 'condition-case 'err
		    (cons 'insert parts)
		;;    '(error (insert (prin1-to-string err))))
	      '(sit-for 0))))

;;; The following wrapper caches results and avoids infinite recursion.
;;; Each cache entry is: ( A B )          Integral of A is B;
;;;			 ( A N )          Integral of A failed at level N;
;;;			 ( A busy )	  Currently working on integral of A;
;;;			 ( A parts )	  Currently working, integ-by-parts;
;;;			 ( A parts2 )	  Currently working, integ-by-parts;
;;;			 ( A cancelled )  Ignore this cache entry;
;;;			 ( A [B] )        Same result as for cur-record = B.
(defun math-integral (expr &optional simplify same-as-above)
  (let* ((simp cur-record)
	 (cur-record (assoc expr math-integral-cache))
	 (math-integ-depth (1+ math-integ-depth))
	 (val 'cancelled))
    (math-tracing-integral "Integrating "
			   (math-format-value expr 1000)
			   "...\n")
    (and cur-record
	 (progn
	   (math-tracing-integral "Found "
				  (math-format-value (nth 1 cur-record) 1000))
	   (and (consp (nth 1 cur-record))
		(math-replace-integral-parts cur-record))
	   (math-tracing-integral " => "
				  (math-format-value (nth 1 cur-record) 1000)
				  "\n")))
    (or (and cur-record
	     (not (eq (nth 1 cur-record) 'cancelled))
	     (or (not (integerp (nth 1 cur-record)))
		 (>= (nth 1 cur-record) math-integ-level)))
	(and (math-integral-contains-parts expr)
	     (progn
	       (setq val nil)
	       t))
	(unwind-protect
	    (progn
	      (let (math-integ-msg)
		(if (eq calc-display-working-message 'lots)
		    (progn
		      (calc-set-command-flag 'clear-message)
		      (setq math-integ-msg (format
					    "Working... Integrating %s"
					    (math-format-flat-expr expr 0)))
		      (message math-integ-msg)))
		(if cur-record
		    (setcar (cdr cur-record)
			    (if same-as-above (vector simp) 'busy))
		  (setq cur-record
			(list expr (if same-as-above (vector simp) 'busy))
			math-integral-cache (cons cur-record
						  math-integral-cache)))
		(if (eq simplify 'yes)
		    (progn
		      (math-tracing-integral "Simplifying...")
		      (setq simp (math-simplify expr))
		      (setq val (if (equal simp expr)
				    (progn
				      (math-tracing-integral " no change\n")
				      (math-do-integral expr))
				  (math-tracing-integral " simplified\n")
				  (math-integral simp 'no t))))
		  (or (setq val (math-do-integral expr))
		      (eq simplify 'no)
		      (let ((simp (math-simplify expr)))
			(or (equal simp expr)
			    (progn
			      (math-tracing-integral "Trying again after "
						     "simplification...\n")
			      (setq val (math-integral simp 'no t))))))))
	      (if (eq calc-display-working-message 'lots)
		  (message math-integ-msg)))
	  (setcar (cdr cur-record) (or val
				       (if (or math-enable-subst
					       (not math-any-substs))
					   math-integ-level
					 'cancelled)))))
    (setq val cur-record)
    (while (vectorp (nth 1 val))
      (setq val (aref (nth 1 val) 0)))
    (setq val (if (memq (nth 1 val) '(parts parts2))
		  (progn
		    (setcar (cdr val) 'parts2)
		    (list 'var 'PARTS val))
		(and (consp (nth 1 val))
		     (nth 1 val))))
    (math-tracing-integral "Integral of "
			   (math-format-value expr 1000)
			   "  is  "
			   (math-format-value val 1000)
			   "\n")
    val))
(defvar math-integral-cache nil)
(defvar math-integral-cache-state nil)

(defun math-integral-contains-parts (expr)
  (if (Math-primp expr)
      (and (eq (car-safe expr) 'var)
	   (eq (nth 1 expr) 'PARTS)
	   (listp (nth 2 expr)))
    (while (and (setq expr (cdr expr))
		(not (math-integral-contains-parts (car expr)))))
    expr))

(defun math-replace-integral-parts (expr)
  (or (Math-primp expr)
      (while (setq expr (cdr expr))
	(and (consp (car expr))
	     (if (eq (car (car expr)) 'var)
		 (and (eq (nth 1 (car expr)) 'PARTS)
		      (consp (nth 2 (car expr)))
		      (if (listp (nth 1 (nth 2 (car expr))))
			  (progn
			    (setcar expr (nth 1 (nth 2 (car expr))))
			    (math-replace-integral-parts (cons 'foo expr)))
			(setcar (cdr cur-record) 'cancelled)))
	       (math-replace-integral-parts (car expr)))))))

(defun math-do-integral (expr)
  (let (t1 t2)
    (or (cond ((not (math-expr-contains expr math-integ-var))
	       (math-mul expr math-integ-var))
	      ((equal expr math-integ-var)
	       (math-div (math-sqr expr) 2))
	      ((eq (car expr) '+)
	       (and (setq t1 (math-integral (nth 1 expr)))
		    (setq t2 (math-integral (nth 2 expr)))
		    (math-add t1 t2)))
	      ((eq (car expr) '-)
	       (and (setq t1 (math-integral (nth 1 expr)))
		    (setq t2 (math-integral (nth 2 expr)))
		    (math-sub t1 t2)))
	      ((eq (car expr) 'neg)
	       (and (setq t1 (math-integral (nth 1 expr)))
		    (math-neg t1)))
	      ((eq (car expr) '*)
	       (cond ((not (math-expr-contains (nth 1 expr) math-integ-var))
		      (and (setq t1 (math-integral (nth 2 expr)))
			   (math-mul (nth 1 expr) t1)))
		     ((not (math-expr-contains (nth 2 expr) math-integ-var))
		      (and (setq t1 (math-integral (nth 1 expr)))
			   (math-mul t1 (nth 2 expr))))
		     ((memq (car-safe (nth 1 expr)) '(+ -))
		      (math-integral (list (car (nth 1 expr))
					   (math-mul (nth 1 (nth 1 expr))
						     (nth 2 expr))
					   (math-mul (nth 2 (nth 1 expr))
						     (nth 2 expr)))
				     'yes t))
		     ((memq (car-safe (nth 2 expr)) '(+ -))
		      (math-integral (list (car (nth 2 expr))
					   (math-mul (nth 1 (nth 2 expr))
						     (nth 1 expr))
					   (math-mul (nth 2 (nth 2 expr))
						     (nth 1 expr)))
				     'yes t))))
	      ((eq (car expr) '/)
	       (cond ((and (not (math-expr-contains (nth 1 expr)
						    math-integ-var))
			   (not (math-equal-int (nth 1 expr) 1)))
		      (and (setq t1 (math-integral (math-div 1 (nth 2 expr))))
			   (math-mul (nth 1 expr) t1)))
		     ((not (math-expr-contains (nth 2 expr) math-integ-var))
		      (and (setq t1 (math-integral (nth 1 expr)))
			   (math-div t1 (nth 2 expr))))
		     ((and (eq (car-safe (nth 1 expr)) '*)
			   (not (math-expr-contains (nth 1 (nth 1 expr))
						    math-integ-var)))
		      (and (setq t1 (math-integral
				     (math-div (nth 2 (nth 1 expr))
					       (nth 2 expr))))
			   (math-mul t1 (nth 1 (nth 1 expr)))))
		     ((and (eq (car-safe (nth 1 expr)) '*)
			   (not (math-expr-contains (nth 2 (nth 1 expr))
						    math-integ-var)))
		      (and (setq t1 (math-integral
				     (math-div (nth 1 (nth 1 expr))
					       (nth 2 expr))))
			   (math-mul t1 (nth 2 (nth 1 expr)))))
		     ((and (eq (car-safe (nth 2 expr)) '*)
			   (not (math-expr-contains (nth 1 (nth 2 expr))
						    math-integ-var)))
		      (and (setq t1 (math-integral
				     (math-div (nth 1 expr)
					       (nth 2 (nth 2 expr)))))
			   (math-div t1 (nth 1 (nth 2 expr)))))
		     ((and (eq (car-safe (nth 2 expr)) '*)
			   (not (math-expr-contains (nth 2 (nth 2 expr))
						    math-integ-var)))
		      (and (setq t1 (math-integral
				     (math-div (nth 1 expr)
					       (nth 1 (nth 2 expr)))))
			   (math-div t1 (nth 2 (nth 2 expr)))))
		     ((eq (car-safe (nth 2 expr)) 'calcFunc-exp)
		      (math-integral
		       (math-mul (nth 1 expr)
				 (list 'calcFunc-exp
				       (math-neg (nth 1 (nth 2 expr)))))))))
	      ((eq (car expr) '^)
	       (cond ((not (math-expr-contains (nth 1 expr) math-integ-var))
		      (or (and (setq t1 (math-is-polynomial (nth 2 expr)
							    math-integ-var 1))
			       (math-div expr
					 (math-mul (nth 1 t1)
						   (math-normalize
						    (list 'calcFunc-ln
							  (nth 1 expr))))))
			  (math-integral
			   (list 'calcFunc-exp
				 (math-mul (nth 2 expr)
					   (math-normalize
					    (list 'calcFunc-ln
						  (nth 1 expr)))))
			   'yes t)))
		     ((not (math-expr-contains (nth 2 expr) math-integ-var))
		      (if (and (integerp (nth 2 expr)) (< (nth 2 expr) 0))
			  (math-integral
			   (list '/ 1 (math-pow (nth 1 expr) (- (nth 2 expr))))
			   nil t)
			(or (and (setq t1 (math-is-polynomial (nth 1 expr)
							      math-integ-var
							      1))
				 (setq t2 (math-add (nth 2 expr) 1))
				 (math-div (math-pow (nth 1 expr) t2)
					   (math-mul t2 (nth 1 t1))))
			    (and (Math-negp (nth 2 expr))
				 (math-integral
				  (math-div 1
					    (math-pow (nth 1 expr)
						      (math-neg
						       (nth 2 expr))))
				  nil t))
			    nil))))))

	;; Integral of a polynomial.
	(and (setq t1 (math-is-polynomial expr math-integ-var 20))
	     (let ((accum 0)
		   (n 1))
	       (while t1
		 (if (setq accum (math-add accum
					   (math-div (math-mul (car t1)
							       (math-pow
								math-integ-var
								n))
						     n))
			   t1 (cdr t1))
		     (setq n (1+ n))))
	       accum))

	;; Try looking it up!
	(cond ((= (length expr) 2)
	       (and (symbolp (car expr))
		    (setq t1 (get (car expr) 'math-integral))
		    (progn
		      (while (and t1
				  (not (setq t2 (funcall (car t1)
							 (nth 1 expr)))))
			(setq t1 (cdr t1)))
		      (and t2 (math-normalize t2)))))
	      ((= (length expr) 3)
	       (and (symbolp (car expr))
		    (setq t1 (get (car expr) 'math-integral-2))
		    (progn
		      (while (and t1
				  (not (setq t2 (funcall (car t1)
							 (nth 1 expr)
							 (nth 2 expr)))))
			(setq t1 (cdr t1)))
		      (and t2 (math-normalize t2))))))

	;; Integral of a rational function.
	(and (math-ratpoly-p expr math-integ-var)
	     (setq t1 (calcFunc-apart expr math-integ-var))
	     (not (equal t1 expr))
	     (math-integral t1))

	;; Try user-defined integration rules.
	(and has-rules
	     (let ((math-old-integ (symbol-function 'calcFunc-integ))
		   (input (list 'calcFunc-integtry expr math-integ-var))
		   res part)
	       (unwind-protect
		   (progn
		     (fset 'calcFunc-integ 'math-sub-integration)
		     (setq res (math-rewrite input
					     '(var IntegRules var-IntegRules)
					     1))
		     (fset 'calcFunc-integ math-old-integ)
		     (and (not (equal res input))
			  (if (setq part (math-expr-calls
					  res '(calcFunc-integsubst)))
			      (and (memq (length part) '(3 4 5))
				   (let ((parts (mapcar
						 (function
						  (lambda (x)
						    (math-expr-subst
						     x (nth 2 part)
						     math-integ-var)))
						 (cdr part))))
				     (math-integrate-by-substitution
				      expr (car parts) t
				      (or (nth 2 parts)
					  (list 'calcFunc-integfailed
						math-integ-var))
				      (nth 3 parts))))
			    (if (not (math-expr-calls res
						      '(calcFunc-integtry
							calcFunc-integfailed)))
				res))))
		 (fset 'calcFunc-integ math-old-integ))))

	;; See if the function is a symbolic derivative.
	(and (string-match "'" (symbol-name (car expr)))
	     (let ((name (symbol-name (car expr)))
		   (p expr) (n 0) (which nil) (bad nil))
	       (while (setq n (1+ n) p (cdr p))
		 (if (equal (car p) math-integ-var)
		     (if which (setq bad t) (setq which n))
		   (if (math-expr-contains (car p) math-integ-var)
		       (setq bad t))))
	       (and which (not bad)
		    (let ((prime (if (= which 1) "'" (format "'%d" which))))
		      (and (string-match (concat prime "\\('['0-9]*\\|$\\)")
					 name)
			   (cons (intern
				  (concat
				   (substring name 0 (match-beginning 0))
				   (substring name (+ (match-beginning 0)
						      (length prime)))))
				 (cdr expr)))))))

	;; Try transformation methods (parts, substitutions).
	(and (> math-integ-level 0)
	     (math-do-integral-methods expr))

	;; Try expanding the function's definition.
	(let ((res (math-expand-formula expr)))
	  (and res
	       (math-integral res))))))

(defun math-sub-integration (expr &rest rest)
  (or (if (or (not rest)
	      (and (< math-integ-level math-integral-limit)
		   (eq (car rest) math-integ-var)))
	  (math-integral expr)
	(let ((res (apply math-old-integ expr rest)))
	  (and (or (= math-integ-level math-integral-limit)
		   (not (math-expr-calls res 'calcFunc-integ)))
	       res)))
      (list 'calcFunc-integfailed expr)))

(defun math-do-integral-methods (expr)
  (let ((so-far math-integ-var-list-list)
	rat-in)

    ;; Integration by substitution, for various likely sub-expressions.
    ;; (In first pass, we look only for sub-exprs that are linear in X.)
    (or (if math-enable-subst
	    (math-integ-try-substitutions expr)
	  (math-integ-try-linear-substitutions expr))

	;; If function has sines and cosines, try tan(x/2) substitution.
	(and (let ((p (setq rat-in (math-expr-rational-in expr))))
	       (while (and p
			   (memq (car (car p)) '(calcFunc-sin
						 calcFunc-cos
						 calcFunc-tan))
			   (equal (nth 1 (car p)) math-integ-var))
		 (setq p (cdr p)))
	       (null p))
	     (or (and (math-integ-parts-easy expr)
		      (math-integ-try-parts expr t))
		 (math-integrate-by-good-substitution
		  expr (list 'calcFunc-tan (math-div math-integ-var 2)))))

	;; If function has sinh and cosh, try tanh(x/2) substitution.
	(and (let ((p rat-in))
	       (while (and p
			   (memq (car (car p)) '(calcFunc-sinh
						 calcFunc-cosh
						 calcFunc-tanh
						 calcFunc-exp))
			   (equal (nth 1 (car p)) math-integ-var))
		 (setq p (cdr p)))
	       (null p))
	     (or (and (math-integ-parts-easy expr)
		      (math-integ-try-parts expr t))
		 (math-integrate-by-good-substitution
		  expr (list 'calcFunc-tanh (math-div math-integ-var 2)))))

	;; If function has square roots, try sin, tan, or sec substitution.
	(and (let ((p rat-in))
	       (setq t1 nil)
	       (while (and p
			   (or (equal (car p) math-integ-var)
			       (and (eq (car (car p)) 'calcFunc-sqrt)
				    (setq t1 (math-is-polynomial
					      (nth 1 (setq t2 (car p)))
					      math-integ-var 2)))))
		 (setq p (cdr p)))
	       (and (null p) t1))
	     (if (cdr (cdr t1))
		 (if (math-guess-if-neg (nth 2 t1))
		     (let* ((c (math-sqrt (math-neg (nth 2 t1))))
			    (d (math-div (nth 1 t1) (math-mul -2 c)))
			    (a (math-sqrt (math-add (car t1) (math-sqr d)))))
		       (math-integrate-by-good-substitution
			expr (list 'calcFunc-arcsin
				   (math-div-thru
				    (math-add (math-mul c math-integ-var) d)
				    a))))
		   (let* ((c (math-sqrt (nth 2 t1)))
			  (d (math-div (nth 1 t1) (math-mul 2 c)))
			  (aa (math-sub (car t1) (math-sqr d))))
		     (if (and nil (not (and (eq d 0) (eq c 1))))
			 (math-integrate-by-good-substitution
			  expr (math-add (math-mul c math-integ-var) d))
		       (if (math-guess-if-neg aa)
			   (math-integrate-by-good-substitution
			    expr (list 'calcFunc-arccosh
				       (math-div-thru
					(math-add (math-mul c math-integ-var)
						  d)
					(math-sqrt (math-neg aa)))))
			 (math-integrate-by-good-substitution
			  expr (list 'calcFunc-arcsinh
				     (math-div-thru
				      (math-add (math-mul c math-integ-var)
						d)
				      (math-sqrt aa))))))))
	       (math-integrate-by-good-substitution expr t2)) )

	;; Try integration by parts.
	(math-integ-try-parts expr)

	;; Give up.
	nil)))

(defun math-integ-parts-easy (expr)
  (cond ((Math-primp expr) t)
	((memq (car expr) '(+ - *))
	 (and (math-integ-parts-easy (nth 1 expr))
	      (math-integ-parts-easy (nth 2 expr))))
	((eq (car expr) '/)
	 (and (math-integ-parts-easy (nth 1 expr))
	      (math-atomic-factorp (nth 2 expr))))
	((eq (car expr) '^)
	 (and (natnump (nth 2 expr))
	      (math-integ-parts-easy (nth 1 expr))))
	((eq (car expr) 'neg)
	 (math-integ-parts-easy (nth 1 expr)))
	(t t)))

(defun math-integ-try-parts (expr &optional math-good-parts)
  ;; Integration by parts:
  ;;   integ(f(x) g(x),x) = f(x) h(x) - integ(h(x) f'(x),x)
  ;;     where h(x) = integ(g(x),x).
  (or (let ((exp (calcFunc-expand expr)))
	(and (not (equal exp expr))
	     (math-integral exp)))
      (and (eq (car expr) '*)
	   (let ((first-bad (or (math-polynomial-p (nth 1 expr)
						   math-integ-var)
				(equal (nth 2 expr) math-prev-parts-v))))
	     (or (and first-bad   ; so try this one first
		      (math-integrate-by-parts (nth 1 expr) (nth 2 expr)))
		 (math-integrate-by-parts (nth 2 expr) (nth 1 expr))
		 (and (not first-bad)
		      (math-integrate-by-parts (nth 1 expr) (nth 2 expr))))))
      (and (eq (car expr) '/)
	   (math-expr-contains (nth 1 expr) math-integ-var)
	   (let ((recip (math-div 1 (nth 2 expr))))
	     (or (math-integrate-by-parts (nth 1 expr) recip)
		 (math-integrate-by-parts recip (nth 1 expr)))))
      (and (eq (car expr) '^)
	   (math-integrate-by-parts (math-pow (nth 1 expr)
					      (math-sub (nth 2 expr) 1))
				    (nth 1 expr)))))

(defun math-integrate-by-parts (u vprime)
  (let ((math-integ-level (if (or math-good-parts
				  (math-polynomial-p u math-integ-var))
			      math-integ-level
			    (1- math-integ-level)))
	(math-doing-parts t)
	v temp)
    (and (>= math-integ-level 0)
	 (unwind-protect
	     (progn
	       (setcar (cdr cur-record) 'parts)
	       (math-tracing-integral "Integrating by parts, u = "
				      (math-format-value u 1000)
				      ", v' = "
				      (math-format-value vprime 1000)
				      "\n")
	       (and (setq v (math-integral vprime))
		    (setq temp (calcFunc-deriv u math-integ-var nil t))
		    (setq temp (let ((math-prev-parts-v v))
				 (math-integral (math-mul v temp) 'yes)))
		    (setq temp (math-sub (math-mul u v) temp))
		    (if (eq (nth 1 cur-record) 'parts)
			(calcFunc-expand temp)
		      (setq v (list 'var 'PARTS cur-record)
			    var-thing (list 'vec (math-sub v temp) v)
			    temp (let (calc-next-why)
				   (math-solve-for (math-sub v temp) 0 v nil)))
		      (and temp (not (integerp temp))
			   (math-simplify-extended temp)))))
	   (setcar (cdr cur-record) 'busy)))))

;;; This tries two different formulations, hoping the algebraic simplifier
;;; will be strong enough to handle at least one.
(defun math-integrate-by-substitution (expr u &optional user uinv uinvprime)
  (and (> math-integ-level 0)
       (let ((math-integ-level (max (- math-integ-level 2) 0)))
	 (math-integrate-by-good-substitution expr u user uinv uinvprime))))

(defun math-integrate-by-good-substitution (expr u &optional user
						 uinv uinvprime)
  (let ((math-living-dangerously t)
	deriv temp)
    (and (setq uinv (if uinv
			(math-expr-subst uinv math-integ-var
					 math-integ-var-2)
		      (let (calc-next-why)
			(math-solve-for u
					math-integ-var-2
					math-integ-var nil))))
	 (progn
	   (math-tracing-integral "Integrating by substitution, u = "
				  (math-format-value u 1000)
				  "\n")
	   (or (and (setq deriv (calcFunc-deriv u
						math-integ-var nil
						(not user)))
		    (setq temp (math-integral (math-expr-subst
					       (math-expr-subst
						(math-expr-subst
						 (math-div expr deriv)
						 u
						 math-integ-var-2)
						math-integ-var
						uinv)
					       math-integ-var-2
					       math-integ-var)
					      'yes)))
	       (and (setq deriv (or uinvprime
				    (calcFunc-deriv uinv
						    math-integ-var-2
						    math-integ-var
						    (not user))))
		    (setq temp (math-integral (math-mul
					       (math-expr-subst
						(math-expr-subst
						 (math-expr-subst
						  expr
						  u
						  math-integ-var-2)
						 math-integ-var
						 uinv)
						math-integ-var-2
						math-integ-var)
					       deriv)
					      'yes)))))
	 (math-simplify-extended
	  (math-expr-subst temp math-integ-var u)))))

;;; Look for substitutions of the form u = a x + b.
(defun math-integ-try-linear-substitutions (sub-expr)
  (and (not (Math-primp sub-expr))
       (or (and (not (memq (car sub-expr) '(+ - * / neg)))
		(not (and (eq (car sub-expr) '^)
			  (integerp (nth 2 sub-expr))))
		(math-expr-contains sub-expr math-integ-var)
		(let ((res nil))
		  (while (and (setq sub-expr (cdr sub-expr))
			      (or (not (math-linear-in (car sub-expr)
						       math-integ-var))
				  (assoc (car sub-expr) so-far)
				  (progn
				    (setq so-far (cons (list (car sub-expr))
						       so-far))
				    (not (setq res
					       (math-integrate-by-substitution
						expr (car sub-expr))))))))
		  res))
	   (let ((res nil))
	     (while (and (setq sub-expr (cdr sub-expr))
			 (not (setq res (math-integ-try-linear-substitutions
					 (car sub-expr))))))
	     res))))

;;; Recursively try different substitutions based on various sub-expressions.
(defun math-integ-try-substitutions (sub-expr &optional allow-rat)
  (and (not (Math-primp sub-expr))
       (not (assoc sub-expr so-far))
       (math-expr-contains sub-expr math-integ-var)
       (or (and (if (and (not (memq (car sub-expr) '(+ - * / neg)))
			 (not (and (eq (car sub-expr) '^)
				   (integerp (nth 2 sub-expr)))))
		    (setq allow-rat t)
		  (prog1 allow-rat (setq allow-rat nil)))
		(not (eq sub-expr expr))
		(or (math-integrate-by-substitution expr sub-expr)
		    (and (eq (car sub-expr) '^)
			 (integerp (nth 2 sub-expr))
			 (< (nth 2 sub-expr) 0)
			 (math-integ-try-substitutions
			  (math-pow (nth 1 sub-expr) (- (nth 2 sub-expr)))
			  t))))
	   (let ((res nil))
	     (setq so-far (cons (list sub-expr) so-far))
	     (while (and (setq sub-expr (cdr sub-expr))
			 (not (setq res (math-integ-try-substitutions
					 (car sub-expr) allow-rat)))))
	     res))))

(defun math-expr-rational-in (expr)
  (let ((parts nil))
    (math-expr-rational-in-rec expr)
    (mapcar 'car parts)))

(defun math-expr-rational-in-rec (expr)
  (cond ((Math-primp expr)
	 (and (equal expr math-integ-var)
	      (not (assoc expr parts))
	      (setq parts (cons (list expr) parts))))
	((or (memq (car expr) '(+ - * / neg))
	     (and (eq (car expr) '^) (integerp (nth 2 expr))))
	 (math-expr-rational-in-rec (nth 1 expr))
	 (and (nth 2 expr) (math-expr-rational-in-rec (nth 2 expr))))
	((and (eq (car expr) '^)
	      (eq (math-quarter-integer (nth 2 expr)) 2))
	 (math-expr-rational-in-rec (list 'calcFunc-sqrt (nth 1 expr))))
	(t
	 (and (not (assoc expr parts))
	      (math-expr-contains expr math-integ-var)
	      (setq parts (cons (list expr) parts))))))

(defun math-expr-calls (expr funcs &optional arg-contains)
  (if (consp expr)
      (if (or (memq (car expr) funcs)
	      (and (eq (car expr) '^) (eq (car funcs) 'calcFunc-sqrt)
		   (eq (math-quarter-integer (nth 2 expr)) 2)))
	  (and (or (not arg-contains)
		   (math-expr-contains expr arg-contains))
	       expr)
	(and (not (Math-primp expr))
	     (let ((res nil))
	       (while (and (setq expr (cdr expr))
			   (not (setq res (math-expr-calls
					   (car expr) funcs arg-contains)))))
	       res)))))

(defun math-fix-const-terms (expr except-vars)
  (cond ((not (math-expr-depends expr except-vars)) 0)
	((Math-primp expr) expr)
	((eq (car expr) '+)
	 (math-add (math-fix-const-terms (nth 1 expr) except-vars)
		   (math-fix-const-terms (nth 2 expr) except-vars)))
	((eq (car expr) '-)
	 (math-sub (math-fix-const-terms (nth 1 expr) except-vars)
		   (math-fix-const-terms (nth 2 expr) except-vars)))
	(t expr)))

;; Command for debugging the Calculator's symbolic integrator.
(defun calc-dump-integral-cache (&optional arg)
  (interactive "P")
  (let ((buf (current-buffer)))
    (unwind-protect
	(let ((p math-integral-cache)
	      cur-record)
	  (display-buffer (get-buffer-create "*Integral Cache*"))
	  (set-buffer (get-buffer "*Integral Cache*"))
	  (erase-buffer)
	  (while p
	    (setq cur-record (car p))
	    (or arg (math-replace-integral-parts cur-record))
	    (insert (math-format-flat-expr (car cur-record) 0)
		    " --> "
		    (if (symbolp (nth 1 cur-record))
			(concat "(" (symbol-name (nth 1 cur-record)) ")")
		      (math-format-flat-expr (nth 1 cur-record) 0))
		    "\n")
	    (setq p (cdr p)))
	  (goto-char (point-min)))
      (set-buffer buf))))

(defun math-try-integral (expr)
  (let ((math-integ-level math-integral-limit)
	(math-integ-depth 0)
	(math-integ-msg "Working...done")
	(cur-record nil)   ; a technicality
	(math-integrating t)
	(calc-prefer-frac t)
	(calc-symbolic-mode t)
	(has-rules (calc-has-rules 'var-IntegRules)))
    (or (math-integral expr 'yes)
	(and math-any-substs
	     (setq math-enable-subst t)
	     (math-integral expr 'yes))
	(and (> math-max-integral-limit math-integral-limit)
	     (setq math-integral-limit math-max-integral-limit
		   math-integ-level math-integral-limit)
	     (math-integral expr 'yes)))))

(defun calcFunc-integ (expr var &optional low high)
  (cond
   ;; Do these even if the parts turn out not to be integrable.
   ((eq (car-safe expr) '+)
    (math-add (calcFunc-integ (nth 1 expr) var low high)
	      (calcFunc-integ (nth 2 expr) var low high)))
   ((eq (car-safe expr) '-)
    (math-sub (calcFunc-integ (nth 1 expr) var low high)
	      (calcFunc-integ (nth 2 expr) var low high)))
   ((eq (car-safe expr) 'neg)
    (math-neg (calcFunc-integ (nth 1 expr) var low high)))
   ((and (eq (car-safe expr) '*)
	 (not (math-expr-contains (nth 1 expr) var)))
    (math-mul (nth 1 expr) (calcFunc-integ (nth 2 expr) var low high)))
   ((and (eq (car-safe expr) '*)
	 (not (math-expr-contains (nth 2 expr) var)))
    (math-mul (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr)))
   ((and (eq (car-safe expr) '/)
	 (not (math-expr-contains (nth 1 expr) var))
	 (not (math-equal-int (nth 1 expr) 1)))
    (math-mul (nth 1 expr)
	      (calcFunc-integ (math-div 1 (nth 2 expr)) var low high)))
   ((and (eq (car-safe expr) '/)
	 (not (math-expr-contains (nth 2 expr) var)))
    (math-div (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr)))
   ((and (eq (car-safe expr) '/)
	 (eq (car-safe (nth 1 expr)) '*)
	 (not (math-expr-contains (nth 1 (nth 1 expr)) var)))
    (math-mul (nth 1 (nth 1 expr))
	      (calcFunc-integ (math-div (nth 2 (nth 1 expr)) (nth 2 expr))
			      var low high)))
   ((and (eq (car-safe expr) '/)
	 (eq (car-safe (nth 1 expr)) '*)
	 (not (math-expr-contains (nth 2 (nth 1 expr)) var)))
    (math-mul (nth 2 (nth 1 expr))
	      (calcFunc-integ (math-div (nth 1 (nth 1 expr)) (nth 2 expr))
			      var low high)))
   ((and (eq (car-safe expr) '/)
	 (eq (car-safe (nth 2 expr)) '*)
	 (not (math-expr-contains (nth 1 (nth 2 expr)) var)))
    (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 2 (nth 2 expr)))
			      var low high)
	      (nth 1 (nth 2 expr))))
   ((and (eq (car-safe expr) '/)
	 (eq (car-safe (nth 2 expr)) '*)
	 (not (math-expr-contains (nth 2 (nth 2 expr)) var)))
    (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 1 (nth 2 expr)))
			      var low high)
	      (nth 2 (nth 2 expr))))
   ((eq (car-safe expr) 'vec)
    (cons 'vec (mapcar (function (lambda (x) (calcFunc-integ x var low high)))
		       (cdr expr))))
   (t
    (let ((state (list calc-angle-mode
		       ;;calc-symbolic-mode
		       ;;calc-prefer-frac
		       calc-internal-prec
		       (calc-var-value 'var-IntegRules)
		       (calc-var-value 'var-IntegSimpRules))))
      (or (equal state math-integral-cache-state)
	  (setq math-integral-cache-state state
		math-integral-cache nil)))
    (let* ((math-max-integral-limit (or (and (boundp 'var-IntegLimit)
					     (natnump var-IntegLimit)
					     var-IntegLimit)
					3))
	   (math-integral-limit 1)
	   (sexpr (math-expr-subst expr var math-integ-var))
	   (trace-buffer (get-buffer "*Trace*"))
	   (calc-language (if (eq calc-language 'big) nil calc-language))
	   (math-any-substs t)
	   (math-enable-subst nil)
	   (math-prev-parts-v nil)
	   (math-doing-parts nil)
	   (math-good-parts nil)
	   (res
	    (if trace-buffer
		(let ((calcbuf (current-buffer))
		      (calcwin (selected-window)))
		  (unwind-protect
		      (progn
			(if (get-buffer-window trace-buffer)
			    (select-window (get-buffer-window trace-buffer)))
			(set-buffer trace-buffer)
			(goto-char (point-max))
			(or (assq 'scroll-stop (buffer-local-variables))
			    (progn
			      (make-local-variable 'scroll-step)
			      (setq scroll-step 3)))
			(insert "\n\n\n")
			(set-buffer calcbuf)
			(math-try-integral sexpr))
		    (select-window calcwin)
		      (set-buffer calcbuf)))
	      (math-try-integral sexpr))))
      (if res
	  (progn
	    (if (calc-has-rules 'var-IntegAfterRules)
		(setq res (math-rewrite res '(var IntegAfterRules
						  var-IntegAfterRules))))
	    (math-simplify
	     (if (and low high)
		 (math-sub (math-expr-subst res math-integ-var high)
			   (math-expr-subst res math-integ-var low))
	       (setq res (math-fix-const-terms res math-integ-vars))
	       (if low
		   (math-expr-subst res math-integ-var low)
		 (math-expr-subst res math-integ-var var)))))
	(append (list 'calcFunc-integ expr var)
		(and low (list low))
		(and high (list high))))))))


(math-defintegral calcFunc-inv
  (math-integral (math-div 1 u)))

(math-defintegral calcFunc-conj
  (let ((int (math-integral u)))
    (and int
	 (list 'calcFunc-conj int))))

(math-defintegral calcFunc-deg
  (let ((int (math-integral u)))
    (and int
	 (list 'calcFunc-deg int))))

(math-defintegral calcFunc-rad
  (let ((int (math-integral u)))
    (and int
	 (list 'calcFunc-rad int))))

(math-defintegral calcFunc-re
  (let ((int (math-integral u)))
    (and int
	 (list 'calcFunc-re int))))

(math-defintegral calcFunc-im
  (let ((int (math-integral u)))
    (and int
	 (list 'calcFunc-im int))))

(math-defintegral calcFunc-sqrt
  (and (equal u math-integ-var)
       (math-mul '(frac 2 3)
		 (list 'calcFunc-sqrt (math-pow u 3)))))

(math-defintegral calcFunc-exp
  (or (and (equal u math-integ-var)
	   (list 'calcFunc-exp u))
      (let ((p (math-is-polynomial u math-integ-var 2)))
	(and (nth 2 p)
	     (let ((sqa (math-sqrt (math-neg (nth 2 p)))))
	       (math-div
		(math-mul
		 (math-mul (math-div (list 'calcFunc-sqrt '(var pi var-pi))
				     sqa)
			   (math-normalize
			    (list 'calcFunc-exp
				  (math-div (math-sub (math-mul (car p)
								(nth 2 p))
						      (math-div
						       (math-sqr (nth 1 p))
						       4))
					    (nth 2 p)))))
		 (list 'calcFunc-erf
		       (math-sub (math-mul sqa math-integ-var)
				 (math-div (nth 1 p) (math-mul 2 sqa)))))
		2))))))

(math-defintegral calcFunc-ln
  (or (and (equal u math-integ-var)
	   (math-sub (math-mul u (list 'calcFunc-ln u)) u))
      (and (eq (car u) '*)
	   (math-integral (math-add (list 'calcFunc-ln (nth 1 u))
				    (list 'calcFunc-ln (nth 2 u)))))
      (and (eq (car u) '/)
	   (math-integral (math-sub (list 'calcFunc-ln (nth 1 u))
				    (list 'calcFunc-ln (nth 2 u)))))
      (and (eq (car u) '^)
	   (math-integral (math-mul (nth 2 u)
				    (list 'calcFunc-ln (nth 1 u)))))))

(math-defintegral calcFunc-log10
  (and (equal u math-integ-var)
       (math-sub (math-mul u (list 'calcFunc-ln u))
		 (math-div u (list 'calcFunc-ln 10)))))

(math-defintegral-2 calcFunc-log
  (math-integral (math-div (list 'calcFunc-ln u)
			   (list 'calcFunc-ln v))))

(math-defintegral calcFunc-sin
  (or (and (equal u math-integ-var)
	   (math-neg (math-from-radians-2 (list 'calcFunc-cos u))))
      (and (nth 2 (math-is-polynomial u math-integ-var 2))
	   (math-integral (math-to-exponentials (list 'calcFunc-sin u))))))

(math-defintegral calcFunc-cos
  (or (and (equal u math-integ-var)
	   (math-from-radians-2 (list 'calcFunc-sin u)))
      (and (nth 2 (math-is-polynomial u math-integ-var 2))
	   (math-integral (math-to-exponentials (list 'calcFunc-cos u))))))

(math-defintegral calcFunc-tan
  (and (equal u math-integ-var)
       (math-neg (math-from-radians-2
		  (list 'calcFunc-ln (list 'calcFunc-cos u))))))

(math-defintegral calcFunc-arcsin
  (and (equal u math-integ-var)
       (math-add (math-mul u (list 'calcFunc-arcsin u))
		 (math-from-radians-2
		  (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))))

(math-defintegral calcFunc-arccos
  (and (equal u math-integ-var)
       (math-sub (math-mul u (list 'calcFunc-arccos u))
		 (math-from-radians-2
		  (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))))

(math-defintegral calcFunc-arctan
  (and (equal u math-integ-var)
       (math-sub (math-mul u (list 'calcFunc-arctan u))
		 (math-from-radians-2
		  (math-div (list 'calcFunc-ln (math-add 1 (math-sqr u)))
			    2)))))

(math-defintegral calcFunc-sinh
  (and (equal u math-integ-var)
       (list 'calcFunc-cosh u)))

(math-defintegral calcFunc-cosh
  (and (equal u math-integ-var)
       (list 'calcFunc-sinh u)))

(math-defintegral calcFunc-tanh
  (and (equal u math-integ-var)
       (list 'calcFunc-ln (list 'calcFunc-cosh u))))

(math-defintegral calcFunc-arcsinh
  (and (equal u math-integ-var)
       (math-sub (math-mul u (list 'calcFunc-arcsinh u))
		 (list 'calcFunc-sqrt (math-add (math-sqr u) 1)))))

(math-defintegral calcFunc-arccosh
  (and (equal u math-integ-var)
       (math-sub (math-mul u (list 'calcFunc-arccosh u))
		 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))

(math-defintegral calcFunc-arctanh
  (and (equal u math-integ-var)
       (math-sub (math-mul u (list 'calcFunc-arctan u))
		 (math-div (list 'calcFunc-ln
				 (math-add 1 (math-sqr u)))
			   2))))

;;; (Ax + B) / (ax^2 + bx + c)^n forms.
(math-defintegral-2 /
  (math-integral-rational-funcs u v))

(defun math-integral-rational-funcs (u v)
  (let ((pu (math-is-polynomial u math-integ-var 1))
	(vpow 1) pv)
    (and pu
	 (catch 'int-rat
	   (if (and (eq (car-safe v) '^) (natnump (nth 2 v)))
	       (setq vpow (nth 2 v)
		     v (nth 1 v)))
	   (and (setq pv (math-is-polynomial v math-integ-var 2))
		(let ((int (math-mul-thru
			    (car pu)
			    (math-integral-q02 (car pv) (nth 1 pv)
					       (nth 2 pv) v vpow))))
		  (if (cdr pu)
		      (setq int (math-add int
					  (math-mul-thru
					   (nth 1 pu)
					   (math-integral-q12
					    (car pv) (nth 1 pv)
					    (nth 2 pv) v vpow)))))
		  int))))))

(defun math-integral-q12 (a b c v vpow)
  (let (q)
    (cond ((not c)
	   (cond ((= vpow 1)
		  (math-sub (math-div math-integ-var b)
			    (math-mul (math-div a (math-sqr b))
				      (list 'calcFunc-ln v))))
		 ((= vpow 2)
		  (math-div (math-add (list 'calcFunc-ln v)
				      (math-div a v))
			    (math-sqr b)))
		 (t
		  (let ((nm1 (math-sub vpow 1))
			(nm2 (math-sub vpow 2)))
		    (math-div (math-sub
			       (math-div a (math-mul nm1 (math-pow v nm1)))
			       (math-div 1 (math-mul nm2 (math-pow v nm2))))
			      (math-sqr b))))))
	  ((math-zerop
	    (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b))))
	   (let ((part (math-div b (math-mul 2 c))))
	     (math-mul-thru (math-pow c vpow)
			    (math-integral-q12 part 1 nil
					       (math-add math-integ-var part)
					       (* vpow 2)))))
	  ((= vpow 1)
	   (and (math-ratp q) (math-negp q)
		(let ((calc-symbolic-mode t))
		  (math-ratp (math-sqrt (math-neg q))))
		(throw 'int-rat nil))  ; should have used calcFunc-apart first
	   (math-sub (math-div (list 'calcFunc-ln v) (math-mul 2 c))
		     (math-mul-thru (math-div b (math-mul 2 c))
				    (math-integral-q02 a b c v 1))))
	  (t
	   (let ((n (1- vpow)))
	     (math-sub (math-neg (math-div
				  (math-add (math-mul b math-integ-var)
					    (math-mul 2 a))
				  (math-mul n (math-mul q (math-pow v n)))))
		       (math-mul-thru (math-div (math-mul b (1- (* 2 n)))
						(math-mul n q))
				      (math-integral-q02 a b c v n))))))))

(defun math-integral-q02 (a b c v vpow)
  (let (q rq part)
    (cond ((not c)
	   (cond ((= vpow 1)
		  (math-div (list 'calcFunc-ln v) b))
		 (t
		  (math-div (math-pow v (- 1 vpow))
			    (math-mul (- 1 vpow) b)))))
	  ((math-zerop
	    (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b))))
	   (let ((part (math-div b (math-mul 2 c))))
	     (math-mul-thru (math-pow c vpow)
			    (math-integral-q02 part 1 nil
					       (math-add math-integ-var part)
					       (* vpow 2)))))
	  ((progn
	     (setq part (math-add (math-mul 2 (math-mul c math-integ-var)) b))
	     (> vpow 1))
	   (let ((n (1- vpow)))
	     (math-add (math-div part (math-mul n (math-mul q (math-pow v n))))
		       (math-mul-thru (math-div (math-mul (- (* 4 n) 2) c)
						(math-mul n q))
				      (math-integral-q02 a b c v n)))))
	  ((math-guess-if-neg q)
	   (setq rq (list 'calcFunc-sqrt (math-neg q)))
	   ;;(math-div-thru (list 'calcFunc-ln
	   ;;			(math-div (math-sub part rq)
	   ;;				  (math-add part rq)))
	   ;;		  rq)
	   (math-div (math-mul -2 (list 'calcFunc-arctanh
					(math-div part rq)))
		     rq))
	  (t
	   (setq rq (list 'calcFunc-sqrt q))
	   (math-div (math-mul 2 (math-to-radians-2
				  (list 'calcFunc-arctan
					(math-div part rq))))
		     rq)))))


(math-defintegral calcFunc-erf
  (and (equal u math-integ-var)
       (math-add (math-mul u (list 'calcFunc-erf u))
		 (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u))
				       (list 'calcFunc-sqrt
					     '(var pi var-pi)))))))

(math-defintegral calcFunc-erfc
  (and (equal u math-integ-var)
       (math-sub (math-mul u (list 'calcFunc-erfc u))
		 (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u))
				       (list 'calcFunc-sqrt
					     '(var pi var-pi)))))))




(defvar math-tabulate-initial nil)
(defvar math-tabulate-function nil)
(defun calcFunc-table (expr var &optional low high step)
  (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf)))
  (or high (setq high low low 1))
  (and (or (math-infinitep low) (math-infinitep high))
       (not step)
       (math-scan-for-limits expr))
  (and step (math-zerop step) (math-reject-arg step 'nonzerop))
  (let ((known (+ (if (Math-objectp low) 1 0)
		  (if (Math-objectp high) 1 0)
		  (if (or (null step) (Math-objectp step)) 1 0)))
	(count '(var inf var-inf))
	vec)
    (or (= known 2)   ; handy optimization
	(equal high '(var inf var-inf))
	(progn
	  (setq count (math-div (math-sub high low) (or step 1)))
	  (or (Math-objectp count)
	      (setq count (math-simplify count)))
	  (if (Math-messy-integerp count)
	      (setq count (math-trunc count)))))
    (if (Math-negp count)
	(setq count -1))
    (if (integerp count)
	(let ((var-DUMMY nil)
	      (vec math-tabulate-initial)
	      (math-working-step-2 (1+ count))
	      (math-working-step 0))
	  (setq expr (math-evaluate-expr
		      (math-expr-subst expr var '(var DUMMY var-DUMMY))))
	  (while (>= count 0)
	    (setq math-working-step (1+ math-working-step)
		  var-DUMMY low
		  vec (cond ((eq math-tabulate-function 'calcFunc-sum)
			     (math-add vec (math-evaluate-expr expr)))
			    ((eq math-tabulate-function 'calcFunc-prod)
			     (math-mul vec (math-evaluate-expr expr)))
			    (t
			     (cons (math-evaluate-expr expr) vec)))
		  low (math-add low (or step 1))
		  count (1- count)))
	  (if math-tabulate-function
	      vec
	    (cons 'vec (nreverse vec))))
      (if (Math-integerp count)
	  (calc-record-why 'fixnump high)
	(if (Math-num-integerp low)
	    (if (Math-num-integerp high)
		(calc-record-why 'integerp step)
	      (calc-record-why 'integerp high))
	  (calc-record-why 'integerp low)))
      (append (list (or math-tabulate-function 'calcFunc-table)
		    expr var)
	      (and (not (and (equal low '(neg (var inf var-inf)))
			     (equal high '(var inf var-inf))))
		   (list low high))
	      (and step (list step))))))

(defun math-scan-for-limits (x)
  (cond ((Math-primp x))
	((and (eq (car x) 'calcFunc-subscr)
	      (Math-vectorp (nth 1 x))
	      (math-expr-contains (nth 2 x) var))
	 (let* ((calc-next-why nil)
		(low-val (math-solve-for (nth 2 x) 1 var nil))
		(high-val (math-solve-for (nth 2 x) (1- (length (nth 1 x)))
					  var nil))
		temp)
	   (and low-val (math-realp low-val)
		high-val (math-realp high-val))
	   (and (Math-lessp high-val low-val)
		(setq temp low-val low-val high-val high-val temp))
	   (setq low (math-max low (math-ceiling low-val))
		 high (math-min high (math-floor high-val)))))
	(t
	 (while (setq x (cdr x))
	   (math-scan-for-limits (car x))))))


(defvar math-disable-sums nil)
(defun calcFunc-sum (expr var &optional low high step)
  (if math-disable-sums (math-reject-arg))
  (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2)))
		(math-sum-rec expr var low high step)))
	 (math-disable-sums t))
    (math-normalize res)))

(defun math-sum-rec (expr var &optional low high step)
  (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf)))
  (and low (not high) (setq high low low 1))
  (let (t1 t2 val)
    (setq val
	  (cond
	   ((not (math-expr-contains expr var))
	    (math-mul expr (math-add (math-div (math-sub high low) (or step 1))
				     1)))
	   ((and step (not (math-equal-int step 1)))
	    (if (math-negp step)
		(math-sum-rec expr var high low (math-neg step))
	      (let ((lo (math-simplify (math-div low step))))
		(if (math-known-num-integerp lo)
		    (math-sum-rec (math-normalize
				   (math-expr-subst expr var
						    (math-mul step var)))
				  var lo (math-simplify (math-div high step)))
		  (math-sum-rec (math-normalize
				 (math-expr-subst expr var
						  (math-add (math-mul step var)
							    low)))
				var 0
				(math-simplify (math-div (math-sub high low)
							 step)))))))
	   ((memq (setq t1 (math-compare low high)) '(0 1))
	    (if (eq t1 0)
		(math-expr-subst expr var low)
	      0))
	   ((setq t1 (math-is-polynomial expr var 20))
	    (let ((poly nil)
		  (n 0))
	      (while t1
		(setq poly (math-poly-mix poly 1
					  (math-sum-integer-power n) (car t1))
		      n (1+ n)
		      t1 (cdr t1)))
	      (setq n (math-build-polynomial-expr poly high))
	      (if (memq low '(0 1))
		  n
		(math-sub n (math-build-polynomial-expr poly
							(math-sub low 1))))))
	   ((and (memq (car expr) '(+ -))
		 (setq t1 (math-sum-rec (nth 1 expr) var low high)
		       t2 (math-sum-rec (nth 2 expr) var low high))
		 (not (and (math-expr-calls t1 '(calcFunc-sum))
			   (math-expr-calls t2 '(calcFunc-sum)))))
	    (list (car expr) t1 t2))
	   ((and (eq (car expr) '*)
		 (setq t1 (math-sum-const-factors expr var)))
	    (math-mul (car t1) (math-sum-rec (cdr t1) var low high)))
	   ((and (eq (car expr) '*) (memq (car-safe (nth 1 expr)) '(+ -)))
	    (math-sum-rec (math-add-or-sub (math-mul (nth 1 (nth 1 expr))
						     (nth 2 expr))
					   (math-mul (nth 2 (nth 1 expr))
						     (nth 2 expr))
					   nil (eq (car (nth 1 expr)) '-))
			  var low high))
	   ((and (eq (car expr) '*) (memq (car-safe (nth 2 expr)) '(+ -)))
	    (math-sum-rec (math-add-or-sub (math-mul (nth 1 expr)
						     (nth 1 (nth 2 expr)))
					   (math-mul (nth 1 expr)
						     (nth 2 (nth 2 expr)))
					   nil (eq (car (nth 2 expr)) '-))
			  var low high))
	   ((and (eq (car expr) '/)
		 (not (math-primp (nth 1 expr)))
		 (setq t1 (math-sum-const-factors (nth 1 expr) var)))
	    (math-mul (car t1)
		      (math-sum-rec (math-div (cdr t1) (nth 2 expr))
				    var low high)))
	   ((and (eq (car expr) '/)
		 (setq t1 (math-sum-const-factors (nth 2 expr) var)))
	    (math-div (math-sum-rec (math-div (nth 1 expr) (cdr t1))
				    var low high)
		      (car t1)))
	   ((eq (car expr) 'neg)
	    (math-neg (math-sum-rec (nth 1 expr) var low high)))
	   ((and (eq (car expr) '^)
		 (not (math-expr-contains (nth 1 expr) var))
		 (setq t1 (math-is-polynomial (nth 2 expr) var 1)))
	    (let ((x (math-pow (nth 1 expr) (nth 1 t1))))
	      (math-div (math-mul (math-sub (math-pow x (math-add 1 high))
					    (math-pow x low))
				  (math-pow (nth 1 expr) (car t1)))
			(math-sub x 1))))
	   ((and (setq t1 (math-to-exponentials expr))
		 (setq t1 (math-sum-rec t1 var low high))
		 (not (math-expr-calls t1 '(calcFunc-sum))))
	    (math-to-exps t1))
	   ((memq (car expr) '(calcFunc-ln calcFunc-log10))
	    (list (car expr) (calcFunc-prod (nth 1 expr) var low high)))
	   ((and (eq (car expr) 'calcFunc-log)
		 (= (length expr) 3)
		 (not (math-expr-contains (nth 2 expr) var)))
	    (list 'calcFunc-log
		  (calcFunc-prod (nth 1 expr) var low high)
		  (nth 2 expr)))))
    (if (equal val '(var nan var-nan)) (setq val nil))
    (or val
	(let* ((math-tabulate-initial 0)
	       (math-tabulate-function 'calcFunc-sum))
	  (calcFunc-table expr var low high)))))

(defun calcFunc-asum (expr var low &optional high step no-mul-flag)
  (or high (setq high low low 1))
  (if (and step (not (math-equal-int step 1)))
      (if (math-negp step)
	  (math-mul (math-pow -1 low)
		    (calcFunc-asum expr var high low (math-neg step) t))
	(let ((lo (math-simplify (math-div low step))))
	  (if (math-num-integerp lo)
	      (calcFunc-asum (math-normalize
			      (math-expr-subst expr var
					       (math-mul step var)))
			     var lo (math-simplify (math-div high step)))
	    (calcFunc-asum (math-normalize
			    (math-expr-subst expr var
					     (math-add (math-mul step var)
						       low)))
			   var 0
			   (math-simplify (math-div (math-sub high low)
						    step))))))
    (math-mul (if no-mul-flag 1 (math-pow -1 low))
	      (calcFunc-sum (math-mul (math-pow -1 var) expr) var low high))))

(defun math-sum-const-factors (expr var)
  (let ((const nil)
	(not-const nil)
	(p expr))
    (while (eq (car-safe p) '*)
      (if (math-expr-contains (nth 1 p) var)
	  (setq not-const (cons (nth 1 p) not-const))
	(setq const (cons (nth 1 p) const)))
      (setq p (nth 2 p)))
    (if (math-expr-contains p var)
	(setq not-const (cons p not-const))
      (setq const (cons p const)))
    (and const
	 (cons (let ((temp (car const)))
		 (while (setq const (cdr const))
		   (setq temp (list '* (car const) temp)))
		 temp)
	       (let ((temp (or (car not-const) 1)))
		 (while (setq not-const (cdr not-const))
		   (setq temp (list '* (car not-const) temp)))
		 temp)))))

(defvar math-sum-int-pow-cache (list '(0 1)))
;; Following is from CRC Math Tables, 27th ed, pp. 52-53.
(defun math-sum-integer-power (pow)
  (let ((calc-prefer-frac t)
	(n (length math-sum-int-pow-cache)))
    (while (<= n pow)
      (let* ((new (list 0 0))
	     (lin new)
	     (pp (cdr (nth (1- n) math-sum-int-pow-cache)))
	     (p 2)
	     (sum 0)
	     q)
	(while pp
	  (setq q (math-div (car pp) p)
		new (cons (math-mul q n) new)
		sum (math-add sum q)
		p (1+ p)
		pp (cdr pp)))
	(setcar lin (math-sub 1 (math-mul n sum)))
	(setq math-sum-int-pow-cache
	      (nconc math-sum-int-pow-cache (list (nreverse new)))
	      n (1+ n))))
    (nth pow math-sum-int-pow-cache)))

(defun math-to-exponentials (expr)
  (and (consp expr)
       (= (length expr) 2)
       (let ((x (nth 1 expr))
	     (pi (if calc-symbolic-mode '(var pi var-pi) (math-pi)))
	     (i (if calc-symbolic-mode '(var i var-i) '(cplx 0 1))))
	 (cond ((eq (car expr) 'calcFunc-exp)
		(list '^ '(var e var-e) x))
	       ((eq (car expr) 'calcFunc-sin)
		(or (eq calc-angle-mode 'rad)
		    (setq x (list '/ (list '* x pi) 180)))
		(list '/ (list '-
			       (list '^ '(var e var-e) (list '* x i))
			       (list '^ '(var e var-e)
				     (list 'neg (list '* x i))))
		      (list '* 2 i)))
	       ((eq (car expr) 'calcFunc-cos)
		(or (eq calc-angle-mode 'rad)
		    (setq x (list '/ (list '* x pi) 180)))
		(list '/ (list '+
			       (list '^ '(var e var-e)
				     (list '* x i))
			       (list '^ '(var e var-e)
				     (list 'neg (list '* x i))))
		      2))
	       ((eq (car expr) 'calcFunc-sinh)
		(list '/ (list '-
			       (list '^ '(var e var-e) x)
			       (list '^ '(var e var-e) (list 'neg x)))
		      2))
	       ((eq (car expr) 'calcFunc-cosh)
		(list '/ (list '+
			       (list '^ '(var e var-e) x)
			       (list '^ '(var e var-e) (list 'neg x)))
		      2))
	       (t nil)))))

(defun math-to-exps (expr)
  (cond (calc-symbolic-mode expr)
	((Math-primp expr)
	 (if (equal expr '(var e var-e)) (math-e) expr))
	((and (eq (car expr) '^)
	      (equal (nth 1 expr) '(var e var-e)))
	 (list 'calcFunc-exp (nth 2 expr)))
	(t
	 (cons (car expr) (mapcar 'math-to-exps (cdr expr))))))


(defvar math-disable-prods nil)
(defun calcFunc-prod (expr var &optional low high step)
  (if math-disable-prods (math-reject-arg))
  (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2)))
		(math-prod-rec expr var low high step)))
	 (math-disable-prods t))
    (math-normalize res)))

(defun math-prod-rec (expr var &optional low high step)
  (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf)))
  (and low (not high) (setq high '(var inf var-inf)))
  (let (t1 t2 t3 val)
    (setq val
	  (cond
	   ((not (math-expr-contains expr var))
	    (math-pow expr (math-add (math-div (math-sub high low) (or step 1))
				     1)))
	   ((and step (not (math-equal-int step 1)))
	    (if (math-negp step)
		(math-prod-rec expr var high low (math-neg step))
	      (let ((lo (math-simplify (math-div low step))))
		(if (math-known-num-integerp lo)
		    (math-prod-rec (math-normalize
				    (math-expr-subst expr var
						     (math-mul step var)))
				   var lo (math-simplify (math-div high step)))
		  (math-prod-rec (math-normalize
				  (math-expr-subst expr var
						   (math-add (math-mul step
								       var)
							     low)))
				 var 0
				 (math-simplify (math-div (math-sub high low)
							  step)))))))
	   ((and (memq (car expr) '(* /))
		 (setq t1 (math-prod-rec (nth 1 expr) var low high)
		       t2 (math-prod-rec (nth 2 expr) var low high))
		 (not (and (math-expr-calls t1 '(calcFunc-prod))
			   (math-expr-calls t2 '(calcFunc-prod)))))
	    (list (car expr) t1 t2))
	   ((and (eq (car expr) '^)
		 (not (math-expr-contains (nth 2 expr) var)))
	    (math-pow (math-prod-rec (nth 1 expr) var low high)
		      (nth 2 expr)))
	   ((and (eq (car expr) '^)
		 (not (math-expr-contains (nth 1 expr) var)))
	    (math-pow (nth 1 expr)
		      (calcFunc-sum (nth 2 expr) var low high)))
	   ((eq (car expr) 'sqrt)
	    (math-normalize (list 'calcFunc-sqrt
				  (list 'calcFunc-prod (nth 1 expr)
					var low high))))
	   ((eq (car expr) 'neg)
	    (math-mul (math-pow -1 (math-add (math-sub high low) 1))
		      (math-prod-rec (nth 1 expr) var low high)))
	   ((eq (car expr) 'calcFunc-exp)
	    (list 'calcFunc-exp (calcFunc-sum (nth 1 expr) var low high)))
	   ((and (setq t1 (math-is-polynomial expr var 1))
		 (setq t2
		       (cond
			((or (and (math-equal-int (nth 1 t1) 1)
				  (setq low (math-simplify
					     (math-add low (car t1)))
					high (math-simplify
					      (math-add high (car t1)))))
			     (and (math-equal-int (nth 1 t1) -1)
				  (setq t2 low
					low (math-simplify
					     (math-sub (car t1) high))
					high (math-simplify
					      (math-sub (car t1) t2)))))
			 (if (or (math-zerop low) (math-zerop high))
			     0
			   (if (and (or (math-negp low) (math-negp high))
				    (or (math-num-integerp low)
					(math-num-integerp high)))
			       (if (math-posp high)
				   0
				 (math-mul (math-pow -1
						     (math-add
						      (math-add low high) 1))
					   (list '/
						 (list 'calcFunc-fact
						       (math-neg low))
						 (list 'calcFunc-fact
						       (math-sub -1 high)))))
			     (list '/
				   (list 'calcFunc-fact high)
				   (list 'calcFunc-fact (math-sub low 1))))))
			((and (or (and (math-equal-int (nth 1 t1) 2)
				       (setq t2 (math-simplify
						 (math-add (math-mul low 2)
							   (car t1)))
					     t3 (math-simplify
						 (math-add (math-mul high 2)
							   (car t1)))))
				  (and (math-equal-int (nth 1 t1) -2)
				       (setq t2 (math-simplify
						 (math-sub (car t1)
							   (math-mul high 2)))
					     t3 (math-simplify
						 (math-sub (car t1)
							   (math-mul low
								     2))))))
			      (or (math-integerp t2)
				  (and (math-messy-integerp t2)
				       (setq t2 (math-trunc t2)))
				  (math-integerp t3)
				  (and (math-messy-integerp t3)
				       (setq t3 (math-trunc t3)))))
			 (if (or (math-zerop t2) (math-zerop t3))
			     0
			   (if (or (math-evenp t2) (math-evenp t3))
			       (if (or (math-negp t2) (math-negp t3))
				   (if (math-posp high)
				       0
				     (list '/
					   (list 'calcFunc-dfact
						 (math-neg t2))
					   (list 'calcFunc-dfact
						 (math-sub -2 t3))))
				 (list '/
				       (list 'calcFunc-dfact t3)
				       (list 'calcFunc-dfact
					     (math-sub t2 2))))
			     (if (math-negp t3)
				 (list '*
				       (list '^ -1
					     (list '/ (list '- (list '- t2 t3)
							    2)
						   2))
				       (list '/
					     (list 'calcFunc-dfact
						   (math-neg t2))
					     (list 'calcFunc-dfact
						   (math-sub -2 t3))))
			       (if (math-posp t2)
				   (list '/
					 (list 'calcFunc-dfact t3)
					 (list 'calcFunc-dfact
					       (math-sub t2 2)))
				 nil))))))))
	    t2)))
    (if (equal val '(var nan var-nan)) (setq val nil))
    (or val
	(let* ((math-tabulate-initial 1)
	       (math-tabulate-function 'calcFunc-prod))
	  (calcFunc-table expr var low high)))))




(defvar math-solve-ranges nil)
;;; Attempt to reduce lhs = rhs to solve-var = rhs', where solve-var appears
;;; in lhs but not in rhs or rhs'; return rhs'.
;;; Uses global values: solve-*.
(defun math-try-solve-for (lhs rhs &optional sign no-poly)
  (let (t1 t2 t3)
    (cond ((equal lhs solve-var)
	   (setq math-solve-sign sign)
	   (if (eq solve-full 'all)
	       (let ((vec (list 'vec (math-evaluate-expr rhs)))
		     newvec var p)
		 (while math-solve-ranges
		   (setq p (car math-solve-ranges)
			 var (car p)
			 newvec (list 'vec))
		   (while (setq p (cdr p))
		     (setq newvec (nconc newvec
					 (cdr (math-expr-subst
					       vec var (car p))))))
		   (setq vec newvec
			 math-solve-ranges (cdr math-solve-ranges)))
		 (math-normalize vec))
	     rhs))
	  ((Math-primp lhs)
	   nil)
	  ((and (eq (car lhs) '-)
		(eq (car-safe (nth 1 lhs)) (car-safe (nth 2 lhs)))
		(Math-zerop rhs)
		(= (length (nth 1 lhs)) 2)
		(= (length (nth 2 lhs)) 2)
		(setq t1 (get (car (nth 1 lhs)) 'math-inverse))
		(setq t2 (funcall t1 '(var SOLVEDUM SOLVEDUM)))
		(eq (math-expr-contains-count t2 '(var SOLVEDUM SOLVEDUM)) 1)
		(setq t3 (math-solve-above-dummy t2))
		(setq t1 (math-try-solve-for (math-sub (nth 1 (nth 1 lhs))
						       (math-expr-subst
							t2 t3
							(nth 1 (nth 2 lhs))))
					     0)))
	   t1)
	  ((eq (car lhs) 'neg)
	   (math-try-solve-for (nth 1 lhs) (math-neg rhs)
			       (and sign (- sign))))
	  ((and (not (eq solve-full 't)) (math-try-solve-prod)))
	  ((and (not no-poly)
		(setq t2 (math-decompose-poly lhs solve-var 15 rhs)))
	   (setq t1 (cdr (nth 1 t2))
		 t1 (let ((math-solve-ranges math-solve-ranges))
		      (cond ((= (length t1) 5)
			     (apply 'math-solve-quartic (car t2) t1))
			    ((= (length t1) 4)
			     (apply 'math-solve-cubic (car t2) t1))
			    ((= (length t1) 3)
			     (apply 'math-solve-quadratic (car t2) t1))
			    ((= (length t1) 2)
			     (apply 'math-solve-linear (car t2) sign t1))
			    (solve-full
			     (math-poly-all-roots (car t2) t1))
			    (calc-symbolic-mode nil)
			    (t
			     (math-try-solve-for
			      (car t2)
			      (math-poly-any-root (reverse t1) 0 t)
			      nil t)))))
	   (if t1
	       (if (eq (nth 2 t2) 1)
		   t1
		 (math-solve-prod t1 (math-try-solve-for (nth 2 t2) 0 nil t)))
	     (calc-record-why "*Unable to find a symbolic solution")
	     nil))
	  ((and (math-solve-find-root-term lhs nil)
		(eq (math-expr-contains-count lhs t1) 1))   ; just in case
	   (math-try-solve-for (math-simplify
				(math-sub (if (or t3 (math-evenp t2))
					      (math-pow t1 t2)
					    (math-neg (math-pow t1 t2)))
					  (math-expand-power
					   (math-sub (math-normalize
						      (math-expr-subst
						       lhs t1 0))
						     rhs)
					   t2 solve-var)))
			       0))
	  ((eq (car lhs) '+)
	   (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
		  (math-try-solve-for (nth 2 lhs)
				      (math-sub rhs (nth 1 lhs))
				      sign))
		 ((not (math-expr-contains (nth 2 lhs) solve-var))
		  (math-try-solve-for (nth 1 lhs)
				      (math-sub rhs (nth 2 lhs))
				      sign))))
	  ((eq (car lhs) 'calcFunc-eq)
	   (math-try-solve-for (math-sub (nth 1 lhs) (nth 2 lhs))
			       rhs sign no-poly))
	  ((eq (car lhs) '-)
	   (cond ((or (and (eq (car-safe (nth 1 lhs)) 'calcFunc-sin)
			   (eq (car-safe (nth 2 lhs)) 'calcFunc-cos))
		      (and (eq (car-safe (nth 1 lhs)) 'calcFunc-cos)
			   (eq (car-safe (nth 2 lhs)) 'calcFunc-sin)))
		  (math-try-solve-for (math-sub (nth 1 lhs)
						(list (car (nth 1 lhs))
						      (math-sub
						       (math-quarter-circle t)
						       (nth 1 (nth 2 lhs)))))
				      rhs))
		 ((not (math-expr-contains (nth 1 lhs) solve-var))
		  (math-try-solve-for (nth 2 lhs)
				      (math-sub (nth 1 lhs) rhs)
				      (and sign (- sign))))
		 ((not (math-expr-contains (nth 2 lhs) solve-var))
		  (math-try-solve-for (nth 1 lhs)
				      (math-add rhs (nth 2 lhs))
				      sign))))
	  ((and (eq solve-full 't) (math-try-solve-prod)))
	  ((and (eq (car lhs) '%)
		(not (math-expr-contains (nth 2 lhs) solve-var)))
	   (math-try-solve-for (nth 1 lhs) (math-add rhs
						     (math-solve-get-int
						      (nth 2 lhs)))))
	  ((eq (car lhs) 'calcFunc-log)
	   (cond ((not (math-expr-contains (nth 2 lhs) solve-var))
		  (math-try-solve-for (nth 1 lhs) (math-pow (nth 2 lhs) rhs)))
		 ((not (math-expr-contains (nth 1 lhs) solve-var))
		  (math-try-solve-for (nth 2 lhs) (math-pow
						   (nth 1 lhs)
						   (math-div 1 rhs))))))
	  ((and (= (length lhs) 2)
		(symbolp (car lhs))
		(setq t1 (get (car lhs) 'math-inverse))
		(setq t2 (funcall t1 rhs)))
	   (setq t1 (get (car lhs) 'math-inverse-sign))
	   (math-try-solve-for (nth 1 lhs) (math-normalize t2)
			       (and sign t1
				    (if (integerp t1)
					(* t1 sign)
				      (funcall t1 lhs sign)))))
	  ((and (symbolp (car lhs))
		(setq t1 (get (car lhs) 'math-inverse-n))
		(setq t2 (funcall t1 lhs rhs)))
	   t2)
	  ((setq t1 (math-expand-formula lhs))
	   (math-try-solve-for t1 rhs sign))
	  (t
	   (calc-record-why "*No inverse known" lhs)
	   nil))))


(defun math-try-solve-prod ()
  (cond ((eq (car lhs) '*)
	 (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
		(math-try-solve-for (nth 2 lhs)
				    (math-div rhs (nth 1 lhs))
				    (math-solve-sign sign (nth 1 lhs))))
	       ((not (math-expr-contains (nth 2 lhs) solve-var))
		(math-try-solve-for (nth 1 lhs)
				    (math-div rhs (nth 2 lhs))
				    (math-solve-sign sign (nth 2 lhs))))
	       ((Math-zerop rhs)
		(math-solve-prod (let ((math-solve-ranges math-solve-ranges))
				   (math-try-solve-for (nth 2 lhs) 0))
				 (math-try-solve-for (nth 1 lhs) 0)))))
	((eq (car lhs) '/)
	 (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
		(math-try-solve-for (nth 2 lhs)
				    (math-div (nth 1 lhs) rhs)
				    (math-solve-sign sign (nth 1 lhs))))
	       ((not (math-expr-contains (nth 2 lhs) solve-var))
		(math-try-solve-for (nth 1 lhs)
				    (math-mul rhs (nth 2 lhs))
				    (math-solve-sign sign (nth 2 lhs))))
	       ((setq t1 (math-try-solve-for (math-sub (nth 1 lhs)
						       (math-mul (nth 2 lhs)
								 rhs))
					     0))
		t1)))
	((eq (car lhs) '^)
	 (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
		(math-try-solve-for
		 (nth 2 lhs)
		 (math-add (math-normalize
			    (list 'calcFunc-log rhs (nth 1 lhs)))
			   (math-div
			    (math-mul 2
				      (math-mul '(var pi var-pi)
						(math-solve-get-int
						 '(var i var-i))))
			    (math-normalize
			     (list 'calcFunc-ln (nth 1 lhs)))))))
	       ((not (math-expr-contains (nth 2 lhs) solve-var))
		(cond ((and (integerp (nth 2 lhs))
			    (>= (nth 2 lhs) 2)
			    (setq t1 (math-integer-log2 (nth 2 lhs))))
		       (setq t2 rhs)
		       (if (and (eq solve-full t)
				(math-known-realp (nth 1 lhs)))
			   (progn
			     (while (>= (setq t1 (1- t1)) 0)
			       (setq t2 (list 'calcFunc-sqrt t2)))
			     (setq t2 (math-solve-get-sign t2)))
			 (while (>= (setq t1 (1- t1)) 0)
			   (setq t2 (math-solve-get-sign
				     (math-normalize
				      (list 'calcFunc-sqrt t2))))))
		       (math-try-solve-for
			(nth 1 lhs)
			(math-normalize t2)))
		      ((math-looks-negp (nth 2 lhs))
		       (math-try-solve-for
			(list '^ (nth 1 lhs) (math-neg (nth 2 lhs)))
			(math-div 1 rhs)))
		      ((and (eq solve-full t)
			    (Math-integerp (nth 2 lhs))
			    (math-known-realp (nth 1 lhs)))
		       (setq t1 (math-normalize
				 (list 'calcFunc-nroot rhs (nth 2 lhs))))
		       (if (math-evenp (nth 2 lhs))
			   (setq t1 (math-solve-get-sign t1)))
		       (math-try-solve-for
			(nth 1 lhs) t1
			(and sign
			     (math-oddp (nth 2 lhs))
			     (math-solve-sign sign (nth 2 lhs)))))
		      (t (math-try-solve-for
			  (nth 1 lhs)
			  (math-mul
			   (math-normalize
			    (list 'calcFunc-exp
				  (if (Math-realp (nth 2 lhs))
				      (math-div (math-mul
						 '(var pi var-pi)
						 (math-solve-get-int
						  '(var i var-i)
						  (and (integerp (nth 2 lhs))
						       (math-abs
							(nth 2 lhs)))))
						(math-div (nth 2 lhs) 2))
				    (math-div (math-mul
					       2
					       (math-mul
						'(var pi var-pi)
						(math-solve-get-int
						 '(var i var-i)
						 (and (integerp (nth 2 lhs))
						      (math-abs
						       (nth 2 lhs))))))
					      (nth 2 lhs)))))
			   (math-normalize
			    (list 'calcFunc-nroot
				  rhs
				  (nth 2 lhs))))
			  (and sign
			       (math-oddp (nth 2 lhs))
			       (math-solve-sign sign (nth 2 lhs)))))))))
	(t nil)))

(defun math-solve-prod (lsoln rsoln)
  (cond ((null lsoln)
	 rsoln)
	((null rsoln)
	 lsoln)
	((eq solve-full 'all)
	 (cons 'vec (append (cdr lsoln) (cdr rsoln))))
	(solve-full
	 (list 'calcFunc-if
	       (list 'calcFunc-gt (math-solve-get-sign 1) 0)
	       lsoln
	       rsoln))
	(t lsoln)))

;;; This deals with negative, fractional, and symbolic powers of "x".
(defun math-solve-poly-funny-powers (sub-rhs)    ; uses "t1", "t2"
  (setq t1 lhs)
  (let ((pp math-poly-neg-powers)
	fac)
    (while pp
      (setq fac (math-pow (car pp) (or math-poly-mult-powers 1))
	    t1 (math-mul t1 fac)
	    rhs (math-mul rhs fac)
	    pp (cdr pp))))
  (if sub-rhs (setq t1 (math-sub t1 rhs)))
  (let ((math-poly-neg-powers nil))
    (setq t2 (math-mul (or math-poly-mult-powers 1)
		       (let ((calc-prefer-frac t))
			 (math-div 1 math-poly-frac-powers)))
	  t1 (math-is-polynomial (math-simplify (calcFunc-expand t1)) b 50))))

;;; This converts "a x^8 + b x^5 + c x^2" to "(a (x^3)^2 + b (x^3) + c) * x^2".
(defun math-solve-crunch-poly (max-degree)   ; uses "t1", "t3"
  (let ((count 0))
    (while (and t1 (Math-zerop (car t1)))
      (setq t1 (cdr t1)
	    count (1+ count)))
    (and t1
	 (let* ((degree (1- (length t1)))
		(scale degree))
	   (while (and (> scale 1) (= (car t3) 1))
	     (and (= (% degree scale) 0)
		  (let ((p t1)
			(n 0)
			(new-t1 nil)
			(okay t))
		    (while (and p okay)
		      (if (= (% n scale) 0)
			  (setq new-t1 (nconc new-t1 (list (car p))))
			(or (Math-zerop (car p))
			    (setq okay nil)))
		      (setq p (cdr p)
			    n (1+ n)))
		    (if okay
			(setq t3 (cons scale (cdr t3))
			      t1 new-t1))))
	     (setq scale (1- scale)))
	   (setq t3 (list (math-mul (car t3) t2) (math-mul count t2)))
	   (<= (1- (length t1)) max-degree)))))

(defun calcFunc-poly (expr var &optional degree)
  (if degree
      (or (natnump degree) (math-reject-arg degree 'fixnatnump))
    (setq degree 50))
  (let ((p (math-is-polynomial expr var degree 'gen)))
    (if p
	(if (equal p '(0))
	    (list 'vec)
	  (cons 'vec p))
      (math-reject-arg expr "Expected a polynomial"))))

(defun calcFunc-gpoly (expr var &optional degree)
  (if degree
      (or (natnump degree) (math-reject-arg degree 'fixnatnump))
    (setq degree 50))
  (let* ((math-poly-base-variable var)
	 (d (math-decompose-poly expr var degree nil)))
    (if d
	(cons 'vec d)
      (math-reject-arg expr "Expected a polynomial"))))

(defun math-decompose-poly (lhs solve-var degree sub-rhs)
  (let ((rhs (or sub-rhs 1))
	t1 t2 t3)
    (setq t2 (math-polynomial-base
	      lhs
	      (function
	       (lambda (b)
		 (let ((math-poly-neg-powers '(1))
		       (math-poly-mult-powers nil)
		       (math-poly-frac-powers 1)
		       (math-poly-exp-base t))
		   (and (not (equal b lhs))
			(or (not (memq (car-safe b) '(+ -))) sub-rhs)
			(setq t3 '(1 0) t2 1
			      t1 (math-is-polynomial lhs b 50))
			(if (and (equal math-poly-neg-powers '(1))
				 (memq math-poly-mult-powers '(nil 1))
				 (eq math-poly-frac-powers 1)
				 sub-rhs)
			    (setq t1 (cons (math-sub (car t1) rhs)
					   (cdr t1)))
			  (math-solve-poly-funny-powers sub-rhs))
			(math-solve-crunch-poly degree)
			(or (math-expr-contains b solve-var)
			    (math-expr-contains (car t3) solve-var))))))))
    (if t2
	(list (math-pow t2 (car t3))
	      (cons 'vec t1)
	      (if sub-rhs
		  (math-pow t2 (nth 1 t3))
		(math-div (math-pow t2 (nth 1 t3)) rhs))))))

(defun math-solve-linear (var sign b a)
  (math-try-solve-for var
		      (math-div (math-neg b) a)
		      (math-solve-sign sign a)
		      t))

(defun math-solve-quadratic (var c b a)
  (math-try-solve-for
   var
   (if (math-looks-evenp b)
       (let ((halfb (math-div b 2)))
	 (math-div
	  (math-add
	   (math-neg halfb)
	   (math-solve-get-sign
	    (math-normalize
	     (list 'calcFunc-sqrt
		   (math-add (math-sqr halfb)
			     (math-mul (math-neg c) a))))))
	  a))
     (math-div
      (math-add
       (math-neg b)
       (math-solve-get-sign
	(math-normalize
	 (list 'calcFunc-sqrt
	       (math-add (math-sqr b)
			 (math-mul 4 (math-mul (math-neg c) a)))))))
      (math-mul 2 a)))
   nil t))

(defun math-solve-cubic (var d c b a)
  (let* ((p (math-div b a))
	 (q (math-div c a))
	 (r (math-div d a))
	 (psqr (math-sqr p))
	 (aa (math-sub q (math-div psqr 3)))
	 (bb (math-add r
		       (math-div (math-sub (math-mul 2 (math-mul psqr p))
					   (math-mul 9 (math-mul p q)))
				 27)))
	 m)
    (if (Math-zerop aa)
	(math-try-solve-for (math-pow (math-add var (math-div p 3)) 3)
			    (math-neg bb) nil t)
      (if (Math-zerop bb)
	  (math-try-solve-for
	   (math-mul (math-add var (math-div p 3))
		     (math-add (math-sqr (math-add var (math-div p 3)))
			       aa))
	   0 nil t)
	(setq m (math-mul 2 (list 'calcFunc-sqrt (math-div aa -3))))
	(math-try-solve-for
	 var
	 (math-sub
	  (math-normalize
	   (math-mul
	    m
	    (list 'calcFunc-cos
		  (math-div
		   (math-sub (list 'calcFunc-arccos
				   (math-div (math-mul 3 bb)
					     (math-mul aa m)))
			     (math-mul 2
				       (math-mul
					(math-add 1 (math-solve-get-int
						     1 3))
					(math-half-circle
					 calc-symbolic-mode))))
		   3))))
	  (math-div p 3))
	 nil t)))))

(defun math-solve-quartic (var d c b a aa)
  (setq a (math-div a aa))
  (setq b (math-div b aa))
  (setq c (math-div c aa))
  (setq d (math-div d aa))
  (math-try-solve-for
   var
   (let* ((asqr (math-sqr a))
	  (asqr4 (math-div asqr 4))
	  (y (let ((solve-full nil)
		   calc-next-why)
	       (math-solve-cubic solve-var
				 (math-sub (math-sub
					    (math-mul 4 (math-mul b d))
					    (math-mul asqr d))
					   (math-sqr c))
				 (math-sub (math-mul a c)
					   (math-mul 4 d))
				 (math-neg b)
				 1)))
	  (rsqr (math-add (math-sub asqr4 b) y))
	  (r (list 'calcFunc-sqrt rsqr))
	  (sign1 (math-solve-get-sign 1))
	  (de (list 'calcFunc-sqrt
		    (math-add
		     (math-sub (math-mul 3 asqr4)
			       (math-mul 2 b))
		     (if (Math-zerop rsqr)
			 (math-mul
			  2
			  (math-mul sign1
				    (list 'calcFunc-sqrt
					  (math-sub (math-sqr y)
						    (math-mul 4 d)))))
		       (math-sub
			(math-mul sign1
				  (math-div
				   (math-sub (math-sub
					      (math-mul 4 (math-mul a b))
					      (math-mul 8 c))
					     (math-mul asqr a))
				   (math-mul 4 r)))
			rsqr))))))
     (math-normalize
      (math-sub (math-add (math-mul sign1 (math-div r 2))
			  (math-solve-get-sign (math-div de 2)))
		(math-div a 4))))
   nil t))

(defvar math-symbolic-solve nil)
(defvar math-int-coefs nil)
(defun math-poly-all-roots (var p &optional math-factoring)
  (catch 'ouch
    (let* ((math-symbolic-solve calc-symbolic-mode)
	   (roots nil)
	   (deg (1- (length p)))
	   (orig-p (reverse p))
	   (math-int-coefs nil)
	   (math-int-scale nil)
	   (math-double-roots nil)
	   (math-int-factors nil)
	   (math-int-threshold nil)
	   (pp p))
      ;; If rational coefficients, look for exact rational factors.
      (while (and pp (Math-ratp (car pp)))
	(setq pp (cdr pp)))
      (if pp
	  (if (or math-factoring math-symbolic-solve)
	      (throw 'ouch nil))
	(let ((lead (car orig-p))
	      (calc-prefer-frac t)
	      (scale (apply 'math-lcm-denoms p)))
	  (setq math-int-scale (math-abs (math-mul scale lead))
		math-int-threshold (math-div '(float 5 -2) math-int-scale)
		math-int-coefs (cdr (math-div (cons 'vec orig-p) lead)))))
      (if (> deg 4)
	  (let ((calc-prefer-frac nil)
		(calc-symbolic-mode nil)
		(pp p)
		(def-p (copy-sequence orig-p)))
	    (while pp
	      (if (Math-numberp (car pp))
		  (setq pp (cdr pp))
		(throw 'ouch nil)))
	    (while (> deg (if math-symbolic-solve 2 4))
	      (let* ((x (math-poly-any-root def-p '(float 0 0) nil))
		     b c pp)
		(if (and (eq (car-safe x) 'cplx)
			 (math-nearly-zerop (nth 2 x) (nth 1 x)))
		    (setq x (calcFunc-re x)))
		(or math-factoring
		    (setq roots (cons x roots)))
		(or (math-numberp x)
		    (setq x (math-evaluate-expr x)))
		(setq pp def-p
		      b (car def-p))
		(while (setq pp (cdr pp))
		  (setq c (car pp))
		  (setcar pp b)
		  (setq b (math-add (math-mul x b) c)))
		(setq def-p (cdr def-p)
		      deg (1- deg))))
	    (setq p (reverse def-p))))
      (if (> deg 1)
	  (let ((solve-var '(var DUMMY var-DUMMY))
		(math-solve-sign nil)
		(math-solve-ranges nil)
		(solve-full 'all))
	    (if (= (length p) (length math-int-coefs))
		(setq p (reverse math-int-coefs)))
	    (setq roots (append (cdr (apply (cond ((= deg 2)
						   'math-solve-quadratic)
						  ((= deg 3)
						   'math-solve-cubic)
						  (t
						   'math-solve-quartic))
					    solve-var p))
				roots)))
	(if (> deg 0)
	    (setq roots (cons (math-div (math-neg (car p)) (nth 1 p))
			      roots))))
      (if math-factoring
	  (progn
	    (while roots
	      (math-poly-integer-root (car roots))
	      (setq roots (cdr roots)))
	    (list math-int-factors (nreverse math-int-coefs) math-int-scale))
	(let ((vec nil) res)
	  (while roots
	    (let ((root (car roots))
		  (solve-full (and solve-full 'all)))
	      (if (math-floatp root)
		  (setq root (math-poly-any-root orig-p root t)))
	      (setq vec (append vec
				(cdr (or (math-try-solve-for var root nil t)
					 (throw 'ouch nil))))))
	    (setq roots (cdr roots)))
	  (setq vec (cons 'vec (nreverse vec)))
	  (if math-symbolic-solve
	      (setq vec (math-normalize vec)))
	  (if (eq solve-full t)
	      (list 'calcFunc-subscr
		    vec
		    (math-solve-get-int 1 (1- (length orig-p)) 1))
	    vec))))))

(defun math-lcm-denoms (&rest fracs)
  (let ((den 1))
    (while fracs
      (if (eq (car-safe (car fracs)) 'frac)
	  (setq den (calcFunc-lcm den (nth 2 (car fracs)))))
      (setq fracs (cdr fracs)))
    den))

(defun math-poly-any-root (p x polish)    ; p is a reverse poly coeff list
  (let* ((newt (if (math-zerop x)
		   (math-poly-newton-root
		    p '(cplx (float 123 -6) (float 1 -4)) 4)
		 (math-poly-newton-root p x 4)))
	 (res (if (math-zerop (cdr newt))
		  (car newt)
		(if (and (math-lessp (cdr newt) '(float 1 -3)) (not polish))
		    (setq newt (math-poly-newton-root p (car newt) 30)))
		(if (math-zerop (cdr newt))
		    (car newt)
		  (math-poly-laguerre-root p x polish)))))
    (and math-symbolic-solve (math-floatp res)
	 (throw 'ouch nil))
    res))

(defun math-poly-newton-root (p x iters)
  (let* ((calc-prefer-frac nil)
	 (calc-symbolic-mode nil)
	 (try-integer math-int-coefs)
	 (dx x) b d)
    (while (and (> (setq iters (1- iters)) 0)
		(let ((pp p))
		  (math-working "newton" x)
		  (setq b (car p)
			d 0)
		  (while (setq pp (cdr pp))
		    (setq d (math-add (math-mul x d) b)
			  b (math-add (math-mul x b) (car pp))))
		  (not (math-zerop d)))
		(progn
		  (setq dx (math-div b d)
			x (math-sub x dx))
		  (if try-integer
		      (let ((adx (math-abs-approx dx)))
			(and (math-lessp adx math-int-threshold)
			     (let ((iroot (math-poly-integer-root x)))
			       (if iroot
				   (setq x iroot dx 0)
				 (setq try-integer nil))))))
		  (or (not (or (eq dx 0)
			       (math-nearly-zerop dx (math-abs-approx x))))
		      (progn (setq dx 0) nil)))))
    (cons x (if (math-zerop x)
		1 (math-div (math-abs-approx dx) (math-abs-approx x))))))

(defun math-poly-integer-root (x)
  (and (math-lessp (calcFunc-xpon (math-abs-approx x)) calc-internal-prec)
       math-int-coefs
       (let* ((calc-prefer-frac t)
	      (xre (calcFunc-re x))
	      (xim (calcFunc-im x))
	      (xresq (math-sqr xre))
	      (ximsq (math-sqr xim)))
	 (if (math-lessp ximsq (calcFunc-scf xresq -1))
	     ;; Look for linear factor
	     (let* ((rnd (math-div (math-round (math-mul xre math-int-scale))
				   math-int-scale))
		    (icp math-int-coefs)
		    (rem (car icp))
		    (newcoef nil))
	       (while (setq icp (cdr icp))
		 (setq newcoef (cons rem newcoef)
		       rem (math-add (car icp)
				     (math-mul rem rnd))))
	       (and (math-zerop rem)
		    (progn
		      (setq math-int-coefs (nreverse newcoef)
			    math-int-factors (cons (list (math-neg rnd))
						   math-int-factors))
		      rnd)))
	   ;; Look for irreducible quadratic factor
	   (let* ((rnd1 (math-div (math-round
				   (math-mul xre (math-mul -2 math-int-scale)))
				  math-int-scale))
		  (sqscale (math-sqr math-int-scale))
		  (rnd0 (math-div (math-round (math-mul (math-add xresq ximsq)
							sqscale))
				  sqscale))
		  (rem1 (car math-int-coefs))
		  (icp (cdr math-int-coefs))
		  (rem0 (car icp))
		  (newcoef nil)
		  (found (assoc (list rnd0 rnd1 (math-posp xim))
				math-double-roots))
		  this)
	     (if found
		 (setq math-double-roots (delq found math-double-roots)
		       rem0 0 rem1 0)
	       (while (setq icp (cdr icp))
		 (setq this rem1
		       newcoef (cons rem1 newcoef)
		       rem1 (math-sub rem0 (math-mul this rnd1))
		       rem0 (math-sub (car icp) (math-mul this rnd0)))))
	     (and (math-zerop rem0)
		  (math-zerop rem1)
		  (let ((aa (math-div rnd1 -2)))
		    (or found (setq math-int-coefs (reverse newcoef)
				    math-double-roots (cons (list
							     (list
							      rnd0 rnd1
							      (math-negp xim)))
							    math-double-roots)
				    math-int-factors (cons (cons rnd0 rnd1)
							   math-int-factors)))
		    (math-add aa
			      (let ((calc-symbolic-mode math-symbolic-solve))
				(math-mul (math-sqrt (math-sub (math-sqr aa)
							       rnd0))
					  (if (math-negp xim) -1 1)))))))))))

;;; The following routine is from Numerical Recipes, section 9.5.
(defun math-poly-laguerre-root (p x polish)
  (let* ((calc-prefer-frac nil)
	 (calc-symbolic-mode nil)
	 (iters 0)
	 (m (1- (length p)))
	 (try-newt (not polish))
	 (tried-newt nil)
	 b d f x1 dx dxold)
    (while
	(and (or (< (setq iters (1+ iters)) 50)
		 (math-reject-arg x "*Laguerre's method failed to converge"))
	     (let ((err (math-abs-approx (car p)))
		   (abx (math-abs-approx x))
		   (pp p))
	       (setq b (car p)
		     d 0 f 0)
	       (while (setq pp (cdr pp))
		 (setq f (math-add (math-mul x f) d)
		       d (math-add (math-mul x d) b)
		       b (math-add (math-mul x b) (car pp))
		       err (math-add (math-abs-approx b) (math-mul abx err))))
	       (math-lessp (calcFunc-scf err (- -2 calc-internal-prec))
			   (math-abs-approx b)))
	     (or (not (math-zerop d))
		 (not (math-zerop f))
		 (progn
		   (setq x (math-pow (math-neg b) (list 'frac 1 m)))
		   nil))
	     (let* ((g (math-div d b))
		    (g2 (math-sqr g))
		    (h (math-sub g2 (math-mul 2 (math-div f b))))
		    (sq (math-sqrt
			 (math-mul (1- m) (math-sub (math-mul m h) g2))))
		    (gp (math-add g sq))
		    (gm (math-sub g sq)))
	       (if (math-lessp (calcFunc-abssqr gp) (calcFunc-abssqr gm))
		   (setq gp gm))
	       (setq dx (math-div m gp)
		     x1 (math-sub x dx))
	       (if (and try-newt
			(math-lessp (math-abs-approx dx)
				    (calcFunc-scf (math-abs-approx x) -3)))
		   (let ((newt (math-poly-newton-root p x1 7)))
		     (setq tried-newt t
			   try-newt nil)
		     (if (math-zerop (cdr newt))
			 (setq x (car newt) x1 x)
		       (if (math-lessp (cdr newt) '(float 1 -6))
			   (let ((newt2 (math-poly-newton-root
					 p (car newt) 20)))
			     (if (math-zerop (cdr newt2))
				 (setq x (car newt2) x1 x)
			       (setq x (car newt))))))))
	       (not (or (eq x x1)
			(math-nearly-equal x x1))))
	     (let ((cdx (math-abs-approx dx)))
	       (setq x x1
		     tried-newt nil)
	       (prog1
		   (or (<= iters 6)
		       (math-lessp cdx dxold)
		       (progn
			 (if polish
			     (let ((digs (calcFunc-xpon
					  (math-div (math-abs-approx x) cdx))))
			       (calc-record-why
				"*Could not attain full precision")
			       (if (natnump digs)
				   (let ((calc-internal-prec (max 3 digs)))
				     (setq x (math-normalize x))))))
			 nil))
		 (setq dxold cdx)))
	     (or polish
		 (math-lessp (calcFunc-scf (math-abs-approx x)
					   (- calc-internal-prec))
			     dxold))))
    (or (and (math-floatp x)
	     (math-poly-integer-root x))
	x)))

(defun math-solve-above-dummy (x)
  (and (not (Math-primp x))
       (if (and (equal (nth 1 x) '(var SOLVEDUM SOLVEDUM))
		(= (length x) 2))
	   x
	 (let ((res nil))
	   (while (and (setq x (cdr x))
		       (not (setq res (math-solve-above-dummy (car x))))))
	   res))))

(defun math-solve-find-root-term (x neg)    ; sets "t2", "t3"
  (if (math-solve-find-root-in-prod x)
      (setq t3 neg
	    t1 x)
    (and (memq (car-safe x) '(+ -))
	 (or (math-solve-find-root-term (nth 1 x) neg)
	     (math-solve-find-root-term (nth 2 x)
					(if (eq (car x) '-) (not neg) neg))))))

(defun math-solve-find-root-in-prod (x)
  (and (consp x)
       (math-expr-contains x solve-var)
       (or (and (eq (car x) 'calcFunc-sqrt)
		(setq t2 2))
	   (and (eq (car x) '^)
		(or (and (memq (math-quarter-integer (nth 2 x)) '(1 2 3))
			 (setq t2 2))
		    (and (eq (car-safe (nth 2 x)) 'frac)
			 (eq (nth 2 (nth 2 x)) 3)
			 (setq t2 3))))
	   (and (memq (car x) '(* /))
		(or (and (not (math-expr-contains (nth 1 x) solve-var))
			 (math-solve-find-root-in-prod (nth 2 x)))
		    (and (not (math-expr-contains (nth 2 x) solve-var))
			 (math-solve-find-root-in-prod (nth 1 x))))))))


(defun math-solve-system (exprs solve-vars solve-full)
  (setq exprs (mapcar 'list (if (Math-vectorp exprs)
				(cdr exprs)
			      (list exprs)))
	solve-vars (if (Math-vectorp solve-vars)
		       (cdr solve-vars)
		     (list solve-vars)))
  (or (let ((math-solve-simplifying nil))
	(math-solve-system-rec exprs solve-vars nil))
      (let ((math-solve-simplifying t))
	(math-solve-system-rec exprs solve-vars nil))))

;;; The following backtracking solver works by choosing a variable
;;; and equation, and trying to solve the equation for the variable.
;;; If it succeeds it calls itself recursively with that variable and
;;; equation removed from their respective lists, and with the solution
;;; added to solns as well as being substituted into all existing
;;; equations.  The algorithm terminates when any solution path
;;; manages to remove all the variables from var-list.

;;; To support calcFunc-roots, entries in eqn-list and solns are
;;; actually lists of equations.

(defun math-solve-system-rec (eqn-list var-list solns)
  (if var-list
      (let ((v var-list)
	    (res nil))

	;; Try each variable in turn.
	(while
	    (and
	     v
	     (let* ((vv (car v))
		    (e eqn-list)
		    (elim (eq (car-safe vv) 'calcFunc-elim)))
	       (if elim
		   (setq vv (nth 1 vv)))

	       ;; Try each equation in turn.
	       (while
		   (and
		    e
		    (let ((e2 (car e))
			  (eprev nil)
			  res2)
		      (setq res nil)

		      ;; Try to solve for vv the list of equations e2.
		      (while (and e2
				  (setq res2 (or (and (eq (car e2) eprev)
						      res2)
						 (math-solve-for (car e2) 0 vv
								 solve-full))))
			(setq eprev (car e2)
			      res (cons (if (eq solve-full 'all)
					    (cdr res2)
					  (list res2))
					res)
			      e2 (cdr e2)))
		      (if e2
			  (setq res nil)

			;; Found a solution.  Now try other variables.
			(setq res (nreverse res)
			      res (math-solve-system-rec
				   (mapcar
				    'math-solve-system-subst
				    (delq (car e)
					  (copy-sequence eqn-list)))
				   (delq (car v) (copy-sequence var-list))
				   (let ((math-solve-simplifying nil)
					 (s (mapcar
					     (function
					      (lambda (x)
						(cons
						 (car x)
						 (math-solve-system-subst
						  (cdr x)))))
					     solns)))
				     (if elim
					 s
				       (cons (cons vv (apply 'append res))
					     s)))))
			(not res))))
		 (setq e (cdr e)))
	       (not res)))
	  (setq v (cdr v)))
	res)

    ;; Eliminated all variables, so now put solution into the proper format.
    (setq solns (sort solns
		      (function
		       (lambda (x y)
			 (not (memq (car x) (memq (car y) solve-vars)))))))
    (if (eq solve-full 'all)
	(math-transpose
	 (math-normalize
	  (cons 'vec
		(if solns
		    (mapcar (function (lambda (x) (cons 'vec (cdr x)))) solns)
		  (mapcar (function (lambda (x) (cons 'vec x))) eqn-list)))))
      (math-normalize
       (cons 'vec
	     (if solns
		 (mapcar (function (lambda (x) (cons 'calcFunc-eq x))) solns)
	       (mapcar 'car eqn-list)))))))

(defun math-solve-system-subst (x)    ; uses "res" and "v"
  (let ((accum nil)
	(res2 res))
    (while x
      (setq accum (nconc accum
			 (mapcar (function
				  (lambda (r)
				    (if math-solve-simplifying
					(math-simplify
					 (math-expr-subst (car x) vv r))
				      (math-expr-subst (car x) vv r))))
				 (car res2)))
	    x (cdr x)
	    res2 (cdr res2)))
    accum))


(defun math-get-from-counter (name)
  (let ((ctr (assq name calc-command-flags)))
    (if ctr
	(setcdr ctr (1+ (cdr ctr)))
      (setq ctr (cons name 1)
	    calc-command-flags (cons ctr calc-command-flags)))
    (cdr ctr)))

(defun math-solve-get-sign (val)
  (setq val (math-simplify val))
  (if (and (eq (car-safe val) '*)
	   (Math-numberp (nth 1 val)))
      (list '* (nth 1 val) (math-solve-get-sign (nth 2 val)))
    (and (eq (car-safe val) 'calcFunc-sqrt)
	 (eq (car-safe (nth 1 val)) '^)
	 (setq val (math-normalize (list '^
					 (nth 1 (nth 1 val))
					 (math-div (nth 2 (nth 1 val)) 2)))))
    (if solve-full
	(if (and (calc-var-value 'var-GenCount)
		 (Math-natnump var-GenCount)
		 (not (eq solve-full 'all)))
	    (prog1
		(math-mul (list 'calcFunc-as var-GenCount) val)
	      (setq var-GenCount (math-add var-GenCount 1))
	      (calc-refresh-evaltos 'var-GenCount))
	  (let* ((var (concat "s" (int-to-string (math-get-from-counter 'solve-sign))))
		 (var2 (list 'var (intern var) (intern (concat "var-" var)))))
	    (if (eq solve-full 'all)
		(setq math-solve-ranges (cons (list var2 1 -1)
					      math-solve-ranges)))
	    (math-mul var2 val)))
      (calc-record-why "*Choosing positive solution")
      val)))

(defun math-solve-get-int (val &optional range first)
  (if solve-full
      (if (and (calc-var-value 'var-GenCount)
	       (Math-natnump var-GenCount)
	       (not (eq solve-full 'all)))
	  (prog1
	      (math-mul val (list 'calcFunc-an var-GenCount))
	    (setq var-GenCount (math-add var-GenCount 1))
	    (calc-refresh-evaltos 'var-GenCount))
	(let* ((var (concat "n" (int-to-string
				 (math-get-from-counter 'solve-int))))
	       (var2 (list 'var (intern var) (intern (concat "var-" var)))))
	  (if (and range (eq solve-full 'all))
	      (setq math-solve-ranges (cons (cons var2
						  (cdr (calcFunc-index
							range (or first 0))))
					    math-solve-ranges)))
	  (math-mul val var2)))
    (calc-record-why "*Choosing 0 for arbitrary integer in solution")
    0))

(defun math-solve-sign (sign expr)
  (and sign
       (let ((s1 (math-possible-signs expr)))
	 (cond ((memq s1 '(4 6))
		sign)
	       ((memq s1 '(1 3))
		(- sign))))))

(defun math-looks-evenp (expr)
  (if (Math-integerp expr)
      (math-evenp expr)
    (if (memq (car expr) '(* /))
	(math-looks-evenp (nth 1 expr)))))

(defun math-solve-for (lhs rhs solve-var solve-full &optional sign)
  (if (math-expr-contains rhs solve-var)
      (math-solve-for (math-sub lhs rhs) 0 solve-var solve-full)
    (and (math-expr-contains lhs solve-var)
	 (math-with-extra-prec 1
	   (let* ((math-poly-base-variable solve-var)
		  (res (math-try-solve-for lhs rhs sign)))
	     (if (and (eq solve-full 'all)
		      (math-known-realp solve-var))
		 (let ((old-len (length res))
		       new-len)
		   (setq res (delq nil
				   (mapcar (function
					    (lambda (x)
					      (and (not (memq (car-safe x)
							      '(cplx polar)))
						   x)))
					   res))
			 new-len (length res))
		   (if (< new-len old-len)
		       (calc-record-why (if (= new-len 1)
					    "*All solutions were complex"
					  (format
					   "*Omitted %d complex solutions"
					   (- old-len new-len)))))))
	     res)))))

(defun math-solve-eqn (expr var full)
  (if (memq (car-safe expr) '(calcFunc-neq calcFunc-lt calcFunc-gt
					   calcFunc-leq calcFunc-geq))
      (let ((res (math-solve-for (cons '- (cdr expr))
				 0 var full
				 (if (eq (car expr) 'calcFunc-neq) nil 1))))
	(and res
	     (if (eq math-solve-sign 1)
		 (list (car expr) var res)
	       (if (eq math-solve-sign -1)
		   (list (car expr) res var)
		 (or (eq (car expr) 'calcFunc-neq)
		     (calc-record-why
		      "*Can't determine direction of inequality"))
		 (and (memq (car expr) '(calcFunc-neq calcFunc-lt calcFunc-gt))
		      (list 'calcFunc-neq var res))))))
    (let ((res (math-solve-for expr 0 var full)))
      (and res
	   (list 'calcFunc-eq var res)))))

(defun math-reject-solution (expr var func)
  (if (math-expr-contains expr var)
      (or (equal (car calc-next-why) '(* "Unable to find a symbolic solution"))
	  (calc-record-why "*Unable to find a solution")))
  (list func expr var))

(defun calcFunc-solve (expr var)
  (or (if (or (Math-vectorp expr) (Math-vectorp var))
	  (math-solve-system expr var nil)
	(math-solve-eqn expr var nil))
      (math-reject-solution expr var 'calcFunc-solve)))

(defun calcFunc-fsolve (expr var)
  (or (if (or (Math-vectorp expr) (Math-vectorp var))
	  (math-solve-system expr var t)
	(math-solve-eqn expr var t))
      (math-reject-solution expr var 'calcFunc-fsolve)))

(defun calcFunc-roots (expr var)
  (let ((math-solve-ranges nil))
    (or (if (or (Math-vectorp expr) (Math-vectorp var))
	    (math-solve-system expr var 'all)
	  (math-solve-for expr 0 var 'all))
      (math-reject-solution expr var 'calcFunc-roots))))

(defun calcFunc-finv (expr var)
  (let ((res (math-solve-for expr math-integ-var var nil)))
    (if res
	(math-normalize (math-expr-subst res math-integ-var var))
      (math-reject-solution expr var 'calcFunc-finv))))

(defun calcFunc-ffinv (expr var)
  (let ((res (math-solve-for expr math-integ-var var t)))
    (if res
	(math-normalize (math-expr-subst res math-integ-var var))
      (math-reject-solution expr var 'calcFunc-finv))))


(put 'calcFunc-inv 'math-inverse
     (function (lambda (x) (math-div 1 x))))
(put 'calcFunc-inv 'math-inverse-sign -1)

(put 'calcFunc-sqrt 'math-inverse
     (function (lambda (x) (math-sqr x))))

(put 'calcFunc-conj 'math-inverse
     (function (lambda (x) (list 'calcFunc-conj x))))

(put 'calcFunc-abs 'math-inverse
     (function (lambda (x) (math-solve-get-sign x))))

(put 'calcFunc-deg 'math-inverse
     (function (lambda (x) (list 'calcFunc-rad x))))
(put 'calcFunc-deg 'math-inverse-sign 1)

(put 'calcFunc-rad 'math-inverse
     (function (lambda (x) (list 'calcFunc-deg x))))
(put 'calcFunc-rad 'math-inverse-sign 1)

(put 'calcFunc-ln 'math-inverse
     (function (lambda (x) (list 'calcFunc-exp x))))
(put 'calcFunc-ln 'math-inverse-sign 1)

(put 'calcFunc-log10 'math-inverse
     (function (lambda (x) (list 'calcFunc-exp10 x))))
(put 'calcFunc-log10 'math-inverse-sign 1)

(put 'calcFunc-lnp1 'math-inverse
     (function (lambda (x) (list 'calcFunc-expm1 x))))
(put 'calcFunc-lnp1 'math-inverse-sign 1)

(put 'calcFunc-exp 'math-inverse
     (function (lambda (x) (math-add (math-normalize (list 'calcFunc-ln x))
				     (math-mul 2
					       (math-mul '(var pi var-pi)
							 (math-solve-get-int
							  '(var i var-i))))))))
(put 'calcFunc-exp 'math-inverse-sign 1)

(put 'calcFunc-expm1 'math-inverse
     (function (lambda (x) (math-add (math-normalize (list 'calcFunc-lnp1 x))
				     (math-mul 2
					       (math-mul '(var pi var-pi)
							 (math-solve-get-int
							  '(var i var-i))))))))
(put 'calcFunc-expm1 'math-inverse-sign 1)

(put 'calcFunc-sin 'math-inverse
     (function (lambda (x) (let ((n (math-solve-get-int 1)))
			     (math-add (math-mul (math-normalize
						  (list 'calcFunc-arcsin x))
						 (math-pow -1 n))
				       (math-mul (math-half-circle t)
						 n))))))

(put 'calcFunc-cos 'math-inverse
     (function (lambda (x) (math-add (math-solve-get-sign
				      (math-normalize
				       (list 'calcFunc-arccos x)))
				     (math-solve-get-int
				      (math-full-circle t))))))

(put 'calcFunc-tan 'math-inverse
     (function (lambda (x) (math-add (math-normalize (list 'calcFunc-arctan x))
				     (math-solve-get-int
				      (math-half-circle t))))))

(put 'calcFunc-arcsin 'math-inverse
     (function (lambda (x) (math-normalize (list 'calcFunc-sin x)))))

(put 'calcFunc-arccos 'math-inverse
     (function (lambda (x) (math-normalize (list 'calcFunc-cos x)))))

(put 'calcFunc-arctan 'math-inverse
     (function (lambda (x) (math-normalize (list 'calcFunc-tan x)))))

(put 'calcFunc-sinh 'math-inverse
     (function (lambda (x) (let ((n (math-solve-get-int 1)))
			     (math-add (math-mul (math-normalize
						  (list 'calcFunc-arcsinh x))
						 (math-pow -1 n))
				       (math-mul (math-half-circle t)
						 (math-mul
						  '(var i var-i)
						  n)))))))
(put 'calcFunc-sinh 'math-inverse-sign 1)

(put 'calcFunc-cosh 'math-inverse
     (function (lambda (x) (math-add (math-solve-get-sign
				      (math-normalize
				       (list 'calcFunc-arccosh x)))
				     (math-mul (math-full-circle t)
					       (math-solve-get-int
						'(var i var-i)))))))

(put 'calcFunc-tanh 'math-inverse
     (function (lambda (x) (math-add (math-normalize
				      (list 'calcFunc-arctanh x))
				     (math-mul (math-half-circle t)
					       (math-solve-get-int
						'(var i var-i)))))))
(put 'calcFunc-tanh 'math-inverse-sign 1)

(put 'calcFunc-arcsinh 'math-inverse
     (function (lambda (x) (math-normalize (list 'calcFunc-sinh x)))))
(put 'calcFunc-arcsinh 'math-inverse-sign 1)

(put 'calcFunc-arccosh 'math-inverse
     (function (lambda (x) (math-normalize (list 'calcFunc-cosh x)))))

(put 'calcFunc-arctanh 'math-inverse
     (function (lambda (x) (math-normalize (list 'calcFunc-tanh x)))))
(put 'calcFunc-arctanh 'math-inverse-sign 1)



(defun calcFunc-taylor (expr var num)
  (let ((x0 0) (v var))
    (if (memq (car-safe var) '(+ - calcFunc-eq))
	(setq x0 (if (eq (car var) '+) (math-neg (nth 2 var)) (nth 2 var))
	      v (nth 1 var)))
    (or (and (eq (car-safe v) 'var)
	     (math-expr-contains expr v)
	     (natnump num)
	     (let ((accum (math-expr-subst expr v x0))
		   (var2 (if (eq (car var) 'calcFunc-eq)
			     (cons '- (cdr var))
			   var))
		   (n 0)
		   (nfac 1)
		   (fprime expr))
	       (while (and (<= (setq n (1+ n)) num)
			   (setq fprime (calcFunc-deriv fprime v nil t)))
		 (setq fprime (math-simplify fprime)
		       nfac (math-mul nfac n)
		       accum (math-add accum
				       (math-div (math-mul (math-pow var2 n)
							   (math-expr-subst
							    fprime v x0))
						 nfac))))
	       (and fprime
		    (math-normalize accum))))
	(list 'calcFunc-taylor expr var num))))

;;; arch-tag: f2932ec8-dd63-418b-a542-11a644b9d4c4
;;; calcalg2.el ends here