comparison lisp/calc/calcalg2.el @ 41047:73f364fd8aaa

Style cleanup; don't put closing parens on their own line, add "foo.el ends here" to each file, and update copyright date.

41046:14b73d89514a

;; Calculator for GNU Emacs, part II [calc-alg-2.el]
;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc.
;; Written by Dave Gillespie, [e-mail].

;; This file is part of GNU Emacs.

;; GNU Emacs is distributed in the hope that it will be useful,

	   (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 "$"))

     (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 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)

		    (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 "$"))

		  (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)

						  (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))))

	(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

							  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

							 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))

			     (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;

    (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 (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))

	     (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)

	;; 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))

	((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 (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

			    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)

						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)))

		  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))

	   (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))

	      (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)

	(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)))

			(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")

	     (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) '+)

	       (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-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)

	  (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))

      (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)))))
)

(setq math-tabulate-initial nil)
(setq math-tabulate-function nil)

(defun math-scan-for-limits (x)

		(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)))))
)


(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))
)
(setq math-disable-sums nil)

(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))

		  (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)

						       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))

		   (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))))
)

;; 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)))

		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))
)
(setq math-sum-int-pow-cache (list '(0 1)))

(defun math-to-exponentials (expr)
  (and (consp expr)
       (= (length expr) 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)))))
)


(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))
)
(setq math-disable-prods nil)

(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)))

	    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))))
)




;;; Attempt to reduce lhs = rhs to solve-var = rhs', where solve-var appears

	   t2)
	  ((setq t1 (math-expand-formula lhs))
	   (math-try-solve-for t1 rhs sign))
	  (t
	   (calc-record-why "*No inverse known" lhs)
	   nil)))
)

(setq math-solve-ranges nil)

(defun math-try-solve-prod ()
  (cond ((eq (car lhs) '*)

				  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)

	(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)

  (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)))

		    (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

    (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)

	(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))

						     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))

			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)
)

(defun math-poly-all-roots (var p &optional math-factoring)
  (catch 'ouch
    (let* ((math-symbolic-solve calc-symbolic-mode)
	   (roots nil)

	      (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)))))
)
(setq math-symbolic-solve nil)

(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)

		(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)

				 (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)

							   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))))))))))
)
(setq math-int-coefs nil)

;;; The following routine is from Numerical Recipes, section 9.5.
(defun math-poly-laguerre-root (p x polish)
  (let* ((calc-prefer-frac nil)

		 (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 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)

		       (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

		  (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

					 (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)))

	    (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)

						  (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)

		       (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))

		      "*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)

							   (math-expr-subst
							    fprime v x0))
						 nfac))))
	       (and fprime
		    (math-normalize accum))))
	(list 'calcFunc-taylor expr var num)))
)

41047:73f364fd8aaa

;; Calculator for GNU Emacs, part II [calc-alg-2.el]
;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.
;; Written by Dave Gillespie, [e-mail].

;; This file is part of GNU Emacs.

;; GNU Emacs is distributed in the hope that it will be useful,

	   (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 "$"))

     (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 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)

		    (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 "$"))

		  (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)

						  (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))))

	(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

							  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

							 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))

			     (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;

    (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 (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))

	     (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)

	;; 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))

	((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 (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

			    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)

						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)))

		  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))

	   (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))

	      (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)

	(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)))

			(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")

	     (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) '+)

	       (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-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)

	  (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))

      (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))))))

(setq math-tabulate-initial nil)
(setq math-tabulate-function nil)

(defun math-scan-for-limits (x)

		(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))))))


(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)))
(setq math-disable-sums nil)

(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))

		  (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)

						       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))

		   (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)))))

;; 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)))

		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)))
(setq math-sum-int-pow-cache (list '(0 1)))

(defun math-to-exponentials (expr)
  (and (consp expr)
       (= (length expr) 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))))))


(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)))
(setq math-disable-prods nil)

(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)))

	    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)))))




;;; Attempt to reduce lhs = rhs to solve-var = rhs', where solve-var appears

	   t2)
	  ((setq t1 (math-expand-formula lhs))
	   (math-try-solve-for t1 rhs sign))
	  (t
	   (calc-record-why "*No inverse known" lhs)
	   nil))))

(setq math-solve-ranges nil)

(defun math-try-solve-prod ()
  (cond ((eq (car lhs) '*)

				  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)

	(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)

  (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)))

		    (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

    (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)

	(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))

						     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))

			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))

(defun math-poly-all-roots (var p &optional math-factoring)
  (catch 'ouch
    (let* ((math-symbolic-solve calc-symbolic-mode)
	   (roots nil)

	      (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))))))
(setq math-symbolic-solve nil)

(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)

		(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)

				 (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)

							   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)))))))))))
(setq math-int-coefs nil)

;;; The following routine is from Numerical Recipes, section 9.5.
(defun math-poly-laguerre-root (p x polish)
  (let* ((calc-prefer-frac nil)

		 (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 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)

		       (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

		  (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

					 (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)))

	    (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)

						  (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)

		       (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))

		      "*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)

							   (math-expr-subst
							    fprime v x0))
						 nfac))))
	       (and fprime
		    (math-normalize accum))))
	(list 'calcFunc-taylor expr var num))))

;;; calcalg2.el ends here