Newsgroups: comp.sources.misc
From: daveg@synaptics.com (David Gillespie)
Subject:  v24i054:  gnucalc - GNU Emacs Calculator, v2.00, Part06/56
Message-ID: <1991Oct29.042514.7256@sparky.imd.sterling.com>
X-Md4-Signature: 9f30cde948d4d7fb5331e55ee96320de
Date: Tue, 29 Oct 1991 04:25:14 GMT
Approved: kent@sparky.imd.sterling.com

Submitted-by: daveg@synaptics.com (David Gillespie)
Posting-number: Volume 24, Issue 54
Archive-name: gnucalc/part06
Environment: Emacs
Supersedes: gmcalc: Volume 13, Issue 27-45

---- Cut Here and unpack ----
#!/bin/sh
# this is Part.06 (part 6 of a multipart archive)
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc-alg-2.el continued
#
if test ! -r _shar_seq_.tmp; then
	echo 'Please unpack part 1 first!'
	exit 1
fi
(read Scheck
 if test "$Scheck" != 6; then
	echo Please unpack part "$Scheck" next!
	exit 1
 else
	exit 0
 fi
) < _shar_seq_.tmp || exit 1
if test ! -f _shar_wnt_.tmp; then
	echo 'x - still skipping calc-alg-2.el'
else
echo 'x - continuing file calc-alg-2.el'
sed 's/^X//' << 'SHAR_EOF' >> 'calc-alg-2.el' &&
X		 (mapcar (function (lambda (x) (cons 'calcFunc-eq x))) solns)
X	       (mapcar 'car eqn-list))))))
)
X
(defun math-solve-system-subst (x)    ; uses "res" and "v"
X  (let ((accum nil)
X	(res2 res))
X    (while x
X      (setq accum (nconc accum
X			 (mapcar (function
X				  (lambda (r)
X				    (if math-solve-simplifying
X					(math-simplify
X					 (math-expr-subst (car x) vv r))
X				      (math-expr-subst (car x) vv r))))
X				 (car res2)))
X	    x (cdr x)
X	    res2 (cdr res2)))
X    accum)
)
X
X
(defun math-get-from-counter (name)
X  (let ((ctr (assq name calc-command-flags)))
X    (if ctr
X	(setcdr ctr (1+ (cdr ctr)))
X      (setq ctr (cons name 1)
X	    calc-command-flags (cons ctr calc-command-flags)))
X    (cdr ctr))
)
X
(defun math-solve-get-sign (val)
X  (setq val (math-simplify val))
X  (if (and (eq (car-safe val) '*)
X	   (Math-numberp (nth 1 val)))
X      (list '* (nth 1 val) (math-solve-get-sign (nth 2 val)))
X    (and (eq (car-safe val) 'calcFunc-sqrt)
X	 (eq (car-safe (nth 1 val)) '^)
X	 (setq val (math-normalize (list '^
X					 (nth 1 (nth 1 val))
X					 (math-div (nth 2 (nth 1 val)) 2)))))
X    (if solve-full
X	(if (and (calc-var-value 'var-GenCount)
X		 (Math-natnump var-GenCount)
X		 (not (eq solve-full 'all)))
X	    (prog1
X		(math-mul (list 'calcFunc-as var-GenCount) val)
X	      (setq var-GenCount (math-add var-GenCount 1))
X	      (calc-refresh-evaltos 'var-GenCount))
X	  (let* ((var (concat "s" (math-get-from-counter 'solve-sign)))
X		 (var2 (list 'var (intern var) (intern (concat "var-" var)))))
X	    (if (eq solve-full 'all)
X		(setq math-solve-ranges (cons (list var2 1 -1)
X					      math-solve-ranges)))
X	    (math-mul var2 val)))
X      (calc-record-why "*Choosing positive solution")
X      val))
)
X
(defun math-solve-get-int (val &optional range first)
X  (if solve-full
X      (if (and (calc-var-value 'var-GenCount)
X	       (Math-natnump var-GenCount)
X	       (not (eq solve-full 'all)))
X	  (prog1
X	      (math-mul val (list 'calcFunc-an var-GenCount))
X	    (setq var-GenCount (math-add var-GenCount 1))
X	    (calc-refresh-evaltos 'var-GenCount))
X	(let* ((var (concat "n" (math-get-from-counter 'solve-int)))
X	       (var2 (list 'var (intern var) (intern (concat "var-" var)))))
X	  (if (and range (eq solve-full 'all))
X	      (setq math-solve-ranges (cons (cons var2
X						  (cdr (calcFunc-index
X							range (or first 0))))
X					    math-solve-ranges)))
X	  (math-mul val var2)))
X    (calc-record-why "*Choosing 0 for arbitrary integer in solution")
X    0)
)
X
(defun math-solve-sign (sign expr)
X  (and sign
X       (let ((s1 (math-possible-signs expr)))
X	 (cond ((memq s1 '(4 6))
X		sign)
X	       ((memq s1 '(1 3))
X		(- sign)))))
)
X
(defun math-looks-evenp (expr)
X  (if (Math-integerp expr)
X      (math-evenp expr)
X    (if (memq (car expr) '(* /))
X	(math-looks-evenp (nth 1 expr))))
)
X
(defun math-solve-for (lhs rhs solve-var solve-full &optional sign)
X  (if (math-expr-contains rhs solve-var)
X      (math-solve-for (math-sub lhs rhs) 0 solve-var solve-full)
X    (and (math-expr-contains lhs solve-var)
X	 (math-with-extra-prec 1
X	   (let* ((math-poly-base-variable solve-var)
X		  (res (math-try-solve-for lhs rhs sign)))
X	     (if (and (eq solve-full 'all)
X		      (math-known-realp solve-var))
X		 (let ((old-len (length res))
X		       new-len)
X		   (setq res (delq nil
X				   (mapcar (function
X					    (lambda (x)
X					      (and (not (memq (car-safe x)
X							      '(cplx polar)))
X						   x)))
X					   res))
X			 new-len (length res))
X		   (if (< new-len old-len)
X		       (calc-record-why (if (= new-len 1)
X					    "*All solutions were complex"
X					  (format
X					   "*Omitted %d complex solutions"
X					   (- old-len new-len)))))))
X	     res))))
)
X
(defun math-solve-eqn (expr var full)
X  (if (memq (car-safe expr) '(calcFunc-neq calcFunc-lt calcFunc-gt
X					   calcFunc-leq calcFunc-geq))
X      (let ((res (math-solve-for (cons '- (cdr expr))
X				 0 var full
X				 (if (eq (car expr) 'calcFunc-neq) nil 1))))
X	(and res
X	     (if (eq math-solve-sign 1)
X		 (list (car expr) var res)
X	       (if (eq math-solve-sign -1)
X		   (list (car expr) res var)
X		 (or (eq (car expr) 'calcFunc-neq)
X		     (calc-record-why
X		      "*Can't determine direction of inequality"))
X		 (and (memq (car expr) '(calcFunc-neq calcFunc-lt calcFunc-gt))
X		      (list 'calcFunc-neq var res))))))
X    (let ((res (math-solve-for expr 0 var full)))
X      (and res
X	   (list 'calcFunc-eq var res))))
)
X
(defun math-reject-solution (expr var func)
X  (if (math-expr-contains expr var)
X      (or (equal (car calc-next-why) '(* "Unable to find a symbolic solution"))
X	  (calc-record-why "*Unable to find a solution")))
X  (list func expr var)
)
X
(defun calcFunc-solve (expr var)
X  (or (if (or (Math-vectorp expr) (Math-vectorp var))
X	  (math-solve-system expr var nil)
X	(math-solve-eqn expr var nil))
X      (math-reject-solution expr var 'calcFunc-solve))
)
X
(defun calcFunc-fsolve (expr var)
X  (or (if (or (Math-vectorp expr) (Math-vectorp var))
X	  (math-solve-system expr var t)
X	(math-solve-eqn expr var t))
X      (math-reject-solution expr var 'calcFunc-fsolve))
)
X
(defun calcFunc-roots (expr var)
X  (let ((math-solve-ranges nil))
X    (or (if (or (Math-vectorp expr) (Math-vectorp var))
X	    (math-solve-system expr var 'all)
X	  (math-solve-for expr 0 var 'all))
X      (math-reject-solution expr var 'calcFunc-roots)))
)
X
(defun calcFunc-finv (expr var)
X  (let ((res (math-solve-for expr math-integ-var var nil)))
X    (if res
X	(math-normalize (math-expr-subst res math-integ-var var))
X      (math-reject-solution expr var 'calcFunc-finv)))
)
X
(defun calcFunc-ffinv (expr var)
X  (let ((res (math-solve-for expr math-integ-var var t)))
X    (if res
X	(math-normalize (math-expr-subst res math-integ-var var))
X      (math-reject-solution expr var 'calcFunc-finv)))
)
X
X
(put 'calcFunc-inv 'math-inverse
X     (function (lambda (x) (math-div 1 x))))
(put 'calcFunc-inv 'math-inverse-sign -1)
X
(put 'calcFunc-sqrt 'math-inverse
X     (function (lambda (x) (math-sqr x))))
X
(put 'calcFunc-conj 'math-inverse
X     (function (lambda (x) (list 'calcFunc-conj x))))
X
(put 'calcFunc-abs 'math-inverse
X     (function (lambda (x) (math-solve-get-sign x))))
X
(put 'calcFunc-deg 'math-inverse
X     (function (lambda (x) (list 'calcFunc-rad x))))
(put 'calcFunc-deg 'math-inverse-sign 1)
X
(put 'calcFunc-rad 'math-inverse
X     (function (lambda (x) (list 'calcFunc-deg x))))
(put 'calcFunc-rad 'math-inverse-sign 1)
X
(put 'calcFunc-ln 'math-inverse
X     (function (lambda (x) (list 'calcFunc-exp x))))
(put 'calcFunc-ln 'math-inverse-sign 1)
X
(put 'calcFunc-log10 'math-inverse
X     (function (lambda (x) (list 'calcFunc-exp10 x))))
(put 'calcFunc-log10 'math-inverse-sign 1)
X
(put 'calcFunc-lnp1 'math-inverse
X     (function (lambda (x) (list 'calcFunc-expm1 x))))
(put 'calcFunc-lnp1 'math-inverse-sign 1)
X
(put 'calcFunc-exp 'math-inverse
X     (function (lambda (x) (math-add (math-normalize (list 'calcFunc-ln x))
X				     (math-mul 2
X					       (math-mul '(var pi var-pi)
X							 (math-solve-get-int
X							  '(var i var-i))))))))
(put 'calcFunc-exp 'math-inverse-sign 1)
X
(put 'calcFunc-expm1 'math-inverse
X     (function (lambda (x) (math-add (math-normalize (list 'calcFunc-lnp1 x))
X				     (math-mul 2
X					       (math-mul '(var pi var-pi)
X							 (math-solve-get-int
X							  '(var i var-i))))))))
(put 'calcFunc-expm1 'math-inverse-sign 1)
X
(put 'calcFunc-sin 'math-inverse
X     (function (lambda (x) (let ((n (math-solve-get-int 1)))
X			     (math-add (math-mul (math-normalize
X						  (list 'calcFunc-arcsin x))
X						 (math-pow -1 n))
X				       (math-mul (math-half-circle t)
X						 n))))))
X
(put 'calcFunc-cos 'math-inverse
X     (function (lambda (x) (math-add (math-solve-get-sign
X				      (math-normalize
X				       (list 'calcFunc-arccos x)))
X				     (math-solve-get-int
X				      (math-full-circle t))))))
X
(put 'calcFunc-tan 'math-inverse
X     (function (lambda (x) (math-add (math-normalize (list 'calcFunc-arctan x))
X				     (math-solve-get-int
X				      (math-half-circle t))))))
X
(put 'calcFunc-arcsin 'math-inverse
X     (function (lambda (x) (math-normalize (list 'calcFunc-sin x)))))
X
(put 'calcFunc-arccos 'math-inverse
X     (function (lambda (x) (math-normalize (list 'calcFunc-cos x)))))
X
(put 'calcFunc-arctan 'math-inverse
X     (function (lambda (x) (math-normalize (list 'calcFunc-tan x)))))
X
(put 'calcFunc-sinh 'math-inverse
X     (function (lambda (x) (let ((n (math-solve-get-int 1)))
X			     (math-add (math-mul (math-normalize
X						  (list 'calcFunc-arcsinh x))
X						 (math-pow -1 n))
X				       (math-mul (math-half-circle t)
X						 (math-mul
X						  '(var i var-i)
X						  n)))))))
(put 'calcFunc-sinh 'math-inverse-sign 1)
X
(put 'calcFunc-cosh 'math-inverse
X     (function (lambda (x) (math-add (math-solve-get-sign
X				      (math-normalize
X				       (list 'calcFunc-arccosh x)))
X				     (math-mul (math-full-circle t)
X					       (math-solve-get-int
X						'(var i var-i)))))))
X
(put 'calcFunc-tanh 'math-inverse
X     (function (lambda (x) (math-add (math-normalize
X				      (list 'calcFunc-arctanh x))
X				     (math-mul (math-half-circle t)
X					       (math-solve-get-int
X						'(var i var-i)))))))
(put 'calcFunc-tanh 'math-inverse-sign 1)
X
(put 'calcFunc-arcsinh 'math-inverse
X     (function (lambda (x) (math-normalize (list 'calcFunc-sinh x)))))
(put 'calcFunc-arcsinh 'math-inverse-sign 1)
X
(put 'calcFunc-arccosh 'math-inverse
X     (function (lambda (x) (math-normalize (list 'calcFunc-cosh x)))))
X
(put 'calcFunc-arctanh 'math-inverse
X     (function (lambda (x) (math-normalize (list 'calcFunc-tanh x)))))
(put 'calcFunc-arctanh 'math-inverse-sign 1)
X
X
X
(defun calcFunc-taylor (expr var num)
X  (let ((x0 0) (v var))
X    (if (memq (car-safe var) '(+ - calcFunc-eq))
X	(setq x0 (if (eq (car var) '+) (math-neg (nth 2 var)) (nth 2 var))
X	      v (nth 1 var)))
X    (or (and (eq (car-safe v) 'var)
X	     (math-expr-contains expr v)
X	     (natnump num)
X	     (let ((accum (math-expr-subst expr v x0))
X		   (var2 (if (eq (car var) 'calcFunc-eq)
X			     (cons '- (cdr var))
X			   var))
X		   (n 0)
X		   (nfac 1)
X		   (fprime expr))
X	       (while (and (<= (setq n (1+ n)) num)
X			   (setq fprime (calcFunc-deriv fprime v nil t)))
X		 (setq fprime (math-simplify fprime)
X		       nfac (math-mul nfac n)
X		       accum (math-add accum
X				       (math-div (math-mul (math-pow var2 n)
X							   (math-expr-subst
X							    fprime v x0))
X						 nfac))))
X	       (and fprime
X		    (math-normalize accum))))
X	(list 'calcFunc-taylor expr var num)))
)
X
X
X
X
SHAR_EOF
echo 'File calc-alg-2.el is complete' &&
chmod 0644 calc-alg-2.el ||
echo 'restore of calc-alg-2.el failed'
Wc_c="`wc -c < 'calc-alg-2.el'`"
test 108099 -eq "$Wc_c" ||
	echo 'calc-alg-2.el: original size 108099, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= calc-alg-3.el ==============
if test -f 'calc-alg-3.el' -a X"$1" != X"-c"; then
	echo 'x - skipping calc-alg-3.el (File already exists)'
	rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting calc-alg-3.el (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'calc-alg-3.el' &&
;; Calculator for GNU Emacs, part II [calc-alg-3.el]
;; Copyright (C) 1990, 1991 Free Software Foundation, Inc.
;; Written by Dave Gillespie, daveg@csvax.cs.caltech.edu.
X
;; This file is part of GNU Emacs.
X
;; 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.
X
;; 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.
X
X
X
;; This file is autoloaded from calc-ext.el.
(require 'calc-ext)
X
(require 'calc-macs)
X
(defun calc-Need-calc-alg-3 () nil)
X
X
(defun calc-find-root (var)
X  (interactive "sVariable(s) to solve for: ")
X  (calc-slow-wrapper
X   (let ((func (if (calc-is-hyperbolic) 'calcFunc-wroot 'calcFunc-root)))
X     (if (or (equal var "") (equal var "$"))
X	 (calc-enter-result 2 "root" (list func
X					   (calc-top-n 3)
X					   (calc-top-n 1)
X					   (calc-top-n 2)))
X       (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
X			   (not (string-match "\\[" var)))
X		      (math-read-expr (concat "[" var "]"))
X		    (math-read-expr var))))
X	 (if (eq (car-safe var) 'error)
X	     (error "Bad format in expression: %s" (nth 1 var)))
X	 (calc-enter-result 1 "root" (list func
X					   (calc-top-n 2)
X					   var
X					   (calc-top-n 1)))))))
)
X
(defun calc-find-minimum (var)
X  (interactive "sVariable(s) to minimize over: ")
X  (calc-slow-wrapper
X   (let ((func (if (calc-is-inverse)
X		   (if (calc-is-hyperbolic)
X		       'calcFunc-wmaximize 'calcFunc-maximize)
X		 (if (calc-is-hyperbolic)
X		     'calcFunc-wminimize 'calcFunc-minimize)))
X	 (tag (if (calc-is-inverse) "max" "min")))
X     (if (or (equal var "") (equal var "$"))
X	 (calc-enter-result 2 tag (list func
X					(calc-top-n 3)
X					(calc-top-n 1)
X					(calc-top-n 2)))
X       (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
X			   (not (string-match "\\[" var)))
X		      (math-read-expr (concat "[" var "]"))
X		    (math-read-expr var))))
X	 (if (eq (car-safe var) 'error)
X	     (error "Bad format in expression: %s" (nth 1 var)))
X	 (calc-enter-result 1 tag (list func
X					(calc-top-n 2)
X					var
X					(calc-top-n 1)))))))
)
X
(defun calc-find-maximum (var)
X  (interactive "sVariable to maximize over: ")
X  (calc-invert-func)
X  (calc-find-minimum var)
)
X
X
(defun calc-poly-interp (arg)
X  (interactive "P")
X  (calc-slow-wrapper
X   (let ((data (calc-top 2)))
X     (if (or (consp arg) (eq arg 0) (eq arg 2))
X	 (setq data (cons 'vec (calc-top-list 2 2)))
X       (or (null arg)
X	   (error "Bad prefix argument")))
X     (if (calc-is-hyperbolic)
X	 (calc-enter-result 1 "rati" (list 'calcFunc-ratint data (calc-top 1)))
X       (calc-enter-result 1 "poli" (list 'calcFunc-polint data
X					 (calc-top 1))))))
)
X
X
(defun calc-curve-fit (arg &optional model coefnames varnames)
X  (interactive "P")
X  (calc-slow-wrapper
X   (setq calc-aborted-prefix nil)
X   (let ((func (if (calc-is-inverse) 'calcFunc-xfit
X		 (if (calc-is-hyperbolic) 'calcFunc-efit
X		   'calcFunc-fit)))
X	 key (which 0)
X	 n nvars temp data
X	 (homog nil)
X	 (msgs '( "(Press ? for help)"
X		  "1 = linear or multilinear"
X		  "2-9 = polynomial fits; i = interpolating polynomial"
X		  "p = a x^b, ^ = a b^x"
X		  "e = a exp(b x), x = exp(a + b x), l = a + b ln(x)"
X		  "E = a 10^(b x), X = 10^(a + b x), L = a + b log10(x)"
X		  "q = a + b (x-c)^2"
X		  "g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)"
X		  "h prefix = homogeneous model (no constant term)"
X		  "' = alg entry, $ = stack, u = Model1, U = Model2")))
X     (while (not model)
X       (message "Fit to model: %s:%s"
X		(nth which msgs)
X		(if homog " h" ""))
X       (setq key (read-char))
X       (cond ((= key ?\C-g)
X	      (keyboard-quit))
X	     ((= key ??)
X	      (setq which (% (1+ which) (length msgs))))
X	     ((memq key '(?h ?H))
X	      (setq homog (not homog)))
X	     ((progn
X		(if (eq key ?\$)
X		    (setq n 1)
X		  (setq n 0))
X		(cond ((null arg)
X		       (setq n (1+ n)
X			     data (calc-top n)))
X		      ((or (consp arg) (eq arg 0))
X		       (setq n (+ n 2)
X			     data (calc-top n)
X			     data (if (math-matrixp data)
X				      (append data (list (calc-top (1- n))))
X				    (list 'vec data (calc-top (1- n))))))
X		      ((> (setq arg (prefix-numeric-value arg)) 0)
X		       (setq data (cons 'vec (calc-top-list arg (1+ n)))
X			     n (+ n arg)))
X		      (t (error "Bad prefix argument")))
X		(or (math-matrixp data) (not (cdr (cdr data)))
X		    (error "Data matrix is not a matrix!"))
X		(setq nvars (- (length data) 2)
X		      coefnames nil
X		      varnames nil)
X		nil))
X	     ((= key ?1)  ; linear or multilinear
X	      (calc-get-fit-variables nvars (1+ nvars) (and homog 0))
X	      (setq model (math-mul coefnames
X				    (cons 'vec (cons 1 (cdr varnames))))))
X	     ((and (>= key ?2) (<= key ?9))   ; polynomial
X	      (calc-get-fit-variables 1 (- key ?0 -1) (and homog 0))
X	      (setq model (math-build-polynomial-expr (cdr coefnames)
X						      (nth 1 varnames))))
X	     ((= key ?i)  ; exact polynomial
X	      (calc-get-fit-variables 1 (1- (length (nth 1 data)))
X				      (and homog 0))
X	      (setq model (math-build-polynomial-expr (cdr coefnames)
X						      (nth 1 varnames))))
X	     ((= key ?p)  ; power law
X	      (calc-get-fit-variables nvars (1+ nvars) (and homog 1))
X	      (setq model (math-mul (nth 1 coefnames)
X				    (calcFunc-reduce
X				     '(var mul var-mul)
X				     (calcFunc-map
X				      '(var pow var-pow)
X				      varnames
X				      (cons 'vec (cdr (cdr coefnames))))))))
X	     ((= key ?^)  ; exponential law
X	      (calc-get-fit-variables nvars (1+ nvars) (and homog 1))
X	      (setq model (math-mul (nth 1 coefnames)
X				    (calcFunc-reduce
X				     '(var mul var-mul)
X				     (calcFunc-map
X				      '(var pow var-pow)
X				      (cons 'vec (cdr (cdr coefnames)))
X				      varnames)))))
X	     ((memq key '(?e ?E))
X	      (calc-get-fit-variables nvars (1+ nvars) (and homog 1))
X	      (setq model (math-mul (nth 1 coefnames)
X				    (calcFunc-reduce
X				     '(var mul var-mul)
X				     (calcFunc-map
X				      (if (eq key ?e)
X					  '(var exp var-exp)
X					'(calcFunc-lambda
X					  (var a var-a)
X					  (^ 10 (var a var-a))))
X				      (calcFunc-map
X				       '(var mul var-mul)
X				       (cons 'vec (cdr (cdr coefnames)))
X				       varnames))))))
X	     ((memq key '(?x ?X))
X	      (calc-get-fit-variables nvars (1+ nvars) (and homog 0))
X	      (setq model (math-mul coefnames
X				    (cons 'vec (cons 1 (cdr varnames)))))
X	      (setq model (if (eq key ?x)
X			      (list 'calcFunc-exp model)
X			    (list '^ 10 model))))
X	     ((memq key '(?l ?L))
X	      (calc-get-fit-variables nvars (1+ nvars) (and homog 0))
X	      (setq model (math-mul coefnames
X				    (cons 'vec
X					  (cons 1 (cdr (calcFunc-map
X							(if (eq key ?l)
X							    '(var ln var-ln)
X							  '(var log10
X								var-log10))
X							varnames)))))))
X	     ((= key ?q)
X	      (calc-get-fit-variables nvars (1+ (* 2 nvars)) (and homog 0))
X	      (let ((c coefnames)
X		    (v varnames))
X		(setq model (nth 1 c))
X		(while (setq v (cdr v) c (cdr (cdr c)))
X		  (setq model (math-add
X			       model
X			       (list '*
X				     (car c)
X				     (list '^
X					   (list '- (car v) (nth 1 c))
X					   2)))))))
X	     ((= key ?g)
X	      (setq model (math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)")
X		    varnames '(vec (var XFit var-XFit))
X		    coefnames '(vec (var AFit var-AFit)
X				    (var BFit var-BFit)
X				    (var CFit var-CFit)))
X	      (calc-get-fit-variables 1 (1- (length coefnames)) (and homog 1)))
X	     ((memq key '(?\$ ?\' ?u ?U))
X	      (let* ((defvars nil)
X		     (record-entry nil))
X		(if (eq key ?\')
X		    (let* ((calc-dollar-values calc-arg-values)
X			   (calc-dollar-used 0)
X			   (calc-hashes-used 0))
X		      (setq model (calc-do-alg-entry "" "Model formula: "))
X		      (if (/= (length model) 1)
X			  (error "Bad format"))
X		      (setq model (car model)
X			    record-entry t)
X		      (if (> calc-dollar-used 0)
X			  (setq coefnames
X				(cons 'vec
X				      (nthcdr (- (length calc-arg-values)
X						 calc-dollar-used)
X					      (reverse calc-arg-values))))
X			(if (> calc-hashes-used 0)
X			    (setq coefnames
X				  (cons 'vec (calc-invent-args
X					      calc-hashes-used))))))
X		  (progn
X		    (setq model (cond ((eq key ?u)
X				       (calc-var-value 'var-Model1))
X				      ((eq key ?U)
X				       (calc-var-value 'var-Model2))
X				      (t (calc-top 1))))
X		    (or model (error "User model not yet defined"))
X		    (if (math-vectorp model)
X			(if (and (memq (length model) '(3 4))
X				 (not (math-objvecp (nth 1 model)))
X				 (math-vectorp (nth 2 model))
X				 (or (null (nth 3 model))
X				     (math-vectorp (nth 3 model))))
X			    (setq varnames (nth 2 model)
X				  coefnames (or (nth 3 model)
X						(cons 'vec
X						      (math-all-vars-but
X						       model varnames)))
X				  model (nth 1 model))
X			  (error "Incorrect model specifier")))))
X		(or varnames
X		    (let ((with-y (eq (car-safe model) 'calcFunc-eq)))
X		      (if coefnames
X			  (calc-get-fit-variables (if with-y (1+ nvars) nvars)
X						  (1- (length coefnames))
X						  (math-all-vars-but
X						   model coefnames)
X						  nil with-y)
X			(let* ((coefs (math-all-vars-but model nil))
X			       (vars nil)
X			       (n (- (length coefs) nvars (if with-y 2 1)))
X			       p)
X			  (if (< n 0)
X			      (error "Not enough variables in model"))
X			  (setq p (nthcdr n coefs))
X			  (setq vars (cdr p))
X			  (setcdr p nil)
X			  (calc-get-fit-variables (if with-y (1+ nvars) nvars)
X						  (length coefs)
X						  vars coefs with-y)))))
X		(if record-entry
X		    (calc-record (list 'vec model varnames coefnames)
X				 "modl"))))
X	     (t (beep))))
X     (let ((calc-fit-to-trail t))
X       (calc-enter-result n (substring (symbol-name func) 9)
X			  (list func model
X				(if (= (length varnames) 2)
X				    (nth 1 varnames)
X				  varnames)
X				(if (= (length coefnames) 2)
X				    (nth 1 coefnames)
X				  coefnames)
X				data))
X       (if (consp calc-fit-to-trail)
X	   (calc-record (calc-normalize calc-fit-to-trail) "parm")))))
)
X
(defun calc-invent-independent-variables (n &optional but)
X  (calc-invent-variables n but '(x y z t) "x")
)
X
(defun calc-invent-parameter-variables (n &optional but)
X  (calc-invent-variables n but '(a b c d) "a")
)
X
(defun calc-invent-variables (num but names base)
X  (let ((vars nil)
X	(n num) (nn 0)
X	var)
X    (while (and (> n 0) names)
X      (setq var (math-build-var-name (if (consp names)
X					 (car names)
X				       (concat base (setq nn (1+ nn))))))
X      (or (math-expr-contains (cons 'vec but) var)
X	  (setq vars (cons var vars)
X		n (1- n)))
X      (or (symbolp names) (setq names (cdr names))))
X    (if (= n 0)
X	(nreverse vars)
X      (calc-invent-variables num but t base)))
)
X
(defun calc-get-fit-variables (nv nc &optional defv defc with-y homog)
X  (or (= nv (if with-y (1+ nvars) nvars))
X      (error "Wrong number of data vectors for this type of model"))
X  (if (integerp defv)
X      (setq homog defv
X	    defv nil))
X  (if homog
X      (setq nc (1- nc)))
X  (or defv
X      (setq defv (calc-invent-independent-variables nv)))
X  (or defc
X      (setq defc (calc-invent-parameter-variables nc defv)))
X  (let ((vars (read-string (format "Fitting variables: (default %s; %s) "
X				   (mapconcat 'symbol-name
X					      (mapcar (function (lambda (v)
X								  (nth 1 v)))
X						      defv)
X					      ",")
X				   (mapconcat 'symbol-name
X					      (mapcar (function (lambda (v)
X								  (nth 1 v)))
X						      defc)
X					      ","))))
X	(coefs nil))
X    (setq vars (if (string-match "\\[" vars)
X		   (math-read-expr vars)
X		 (math-read-expr (concat "[" vars "]"))))
X    (if (eq (car-safe vars) 'error)
X	(error "Bad format in expression: %s" (nth 2 vars)))
X    (or (math-vectorp vars)
X	(error "Expected a variable or vector of variables"))
X    (if (equal vars '(vec))
X	(setq vars (cons 'vec defv)
X	      coefs (cons 'vec defc))
X      (if (math-vectorp (nth 1 vars))
X	  (if (and (= (length vars) 3)
X		   (math-vectorp (nth 2 vars)))
X	      (setq coefs (nth 2 vars)
X		    vars (nth 1 vars))
X	    (error
X	     "Expected independent variables vector, then parameters vector"))
X	(setq coefs (cons 'vec defc))))
X    (or (= nv (1- (length vars)))
X	(and (not with-y) (= (1+ nv) (1- (length vars))))
X	(error "Expected %d independent variable%s" nv (if (= nv 1) "" "s")))
X    (or (= nc (1- (length coefs)))
X	(error "Expected %d parameter variable%s" nc (if (= nc 1) "" "s")))
X    (if homog
X	(setq coefs (cons 'vec (cons homog (cdr coefs)))))
X    (if varnames
X	(setq model (math-multi-subst model (cdr varnames) (cdr vars))))
X    (if coefnames
X	(setq model (math-multi-subst model (cdr coefnames) (cdr coefs))))
X    (setq varnames vars
X	  coefnames coefs))
)
X
X
X
X
;;; The following algorithms are from Numerical Recipes chapter 9.
X
;;; "rtnewt" with safety kludges
(defun math-newton-root (expr deriv guess orig-guess limit)
X  (math-working "newton" guess)
X  (let* ((var-DUMMY guess)
X	 next dval)
X    (setq next (math-evaluate-expr expr)
X	  dval (math-evaluate-expr deriv))
X    (if (and (Math-numberp next)
X	     (Math-numberp dval)
X	     (not (Math-zerop dval)))
X	(progn
X	  (setq next (math-sub guess (math-div next dval)))
X	  (if (math-nearly-equal guess (setq next (math-float next)))
X	      (progn
X		(setq var-DUMMY next)
X		(list 'vec next (math-evaluate-expr expr)))
X	    (if (Math-lessp (math-abs-approx (math-sub next orig-guess))
X			    limit)
X		(math-newton-root expr deriv next orig-guess limit)
X	      (math-reject-arg next "*Newton's method failed to converge"))))
X      (math-reject-arg next "*Newton's method encountered a singularity")))
)
X
;;; Inspired by "rtsafe"
(defun math-newton-search-root (expr deriv guess vguess ostep oostep
X				     low vlow high vhigh)
X  (let ((var-DUMMY guess)
X	(better t)
X	pos step next vnext)
X    (if guess
X	(math-working "newton" (list 'intv 0 low high))
X      (math-working "bisect" (list 'intv 0 low high))
X      (setq ostep (math-mul-float (math-sub-float high low)
X				  '(float 5 -1))
X	    guess (math-add-float low ostep)
X	    var-DUMMY guess
X	    vguess (math-evaluate-expr expr))
X      (or (Math-realp vguess)
X	  (progn
X	    (setq ostep (math-mul-float ostep '(float 6 -1))
X		  guess (math-add-float low ostep)
X		  var-DUMMY guess
X		  vguess (math-evaluate-expr expr))
X	    (or (math-realp vguess)
X		(progn
X		  (setq ostep (math-mul-float ostep '(float 123456 -5))
X			guess (math-add-float low ostep)
X			var-DUMMY guess
X			vguess nil))))))
X    (or vguess
X	(setq vguess (math-evaluate-expr expr)))
X    (or (Math-realp vguess)
X	(math-reject-arg guess "*Newton's method encountered a singularity"))
X    (setq vguess (math-float vguess))
X    (if (eq (Math-negp vlow) (setq pos (Math-posp vguess)))
X	(setq high guess
X	      vhigh vguess)
X      (if (eq (Math-negp vhigh) pos)
X	  (setq low guess
X		vlow vguess)
X	(setq better nil)))
X    (if (or (Math-zerop vguess)
X	    (math-nearly-equal low high))
X	(list 'vec guess vguess)
X      (setq step (math-evaluate-expr deriv))
X      (if (and (Math-realp step)
X	       (not (Math-zerop step))
X	       (setq step (math-div-float vguess (math-float step))
X		     next (math-sub-float guess step))
X	       (not (math-lessp-float high next))
X	       (not (math-lessp-float next low)))
X	  (progn
X	    (setq var-DUMMY next
X		  vnext (math-evaluate-expr expr))
X	    (if (or (Math-zerop vnext)
X		    (math-nearly-equal next guess))
X		(list 'vec next vnext)
X	      (if (and better
X		       (math-lessp-float (math-abs (or oostep
X						       (math-sub-float
X							high low)))
X					 (math-abs
X					  (math-mul-float '(float 2 0)
X							  step))))
X		  (math-newton-search-root expr deriv nil nil nil ostep
X					   low vlow high vhigh)
X		(math-newton-search-root expr deriv next vnext step ostep
X					 low vlow high vhigh))))
X	(if (or (and (Math-posp vlow) (Math-posp vhigh))
X		(and (Math-negp vlow) (Math-negp vhigh)))
X	    (math-search-root expr deriv low vlow high vhigh)
X	  (math-newton-search-root expr deriv nil nil nil ostep
X				   low vlow high vhigh)))))
)
X
;;; Search for a root in an interval with no overt zero crossing.
(defun math-search-root (expr deriv low vlow high vhigh)
X  (let (found)
X    (if root-widen
X	(let ((iters 0)
X	      (iterlim (if (eq root-widen 'point)
X			   (+ calc-internal-prec 10)
X			 20))
X	      (factor (if (eq root-widen 'point)
X			  '(float 9 0)
X			'(float 16 -1)))
X	      (prev nil) vprev waslow
X	      diff)
X	  (while (or (and (math-posp vlow) (math-posp vhigh))
X		     (and (math-negp vlow) (math-negp vhigh)))
X	    (math-working "widen" (list 'intv 0 low high))
X	    (if (> (setq iters (1+ iters)) iterlim)
X		(math-reject-arg (list 'intv 0 low high)
X				 "*Unable to bracket root"))
X	    (if (= iters calc-internal-prec)
X		(setq factor '(float 16 -1)))
X	    (setq diff (math-mul-float (math-sub-float high low) factor))
X	    (if (Math-zerop diff)
X		(setq high (calcFunc-incr high 10))
X	      (if (math-lessp-float (math-abs vlow) (math-abs vhigh))
X		  (setq waslow t
X			prev low
X			low (math-sub low diff)
X			var-DUMMY low
X			vprev vlow
X			vlow (math-evaluate-expr expr))
X		(setq waslow nil
X		      prev high
X		      high (math-add high diff)
X		      var-DUMMY high
X		      vprev vhigh
X		      vhigh (math-evaluate-expr expr)))))
X	  (if prev
X	      (if waslow
X		  (setq high prev vhigh vprev)
X		(setq low prev vlow vprev)))
X	  (setq found t))
X      (or (Math-realp vlow)
X	  (math-reject-arg vlow 'realp))
X      (or (Math-realp vhigh)
X	  (math-reject-arg vhigh 'realp))
X      (let ((xvals (list low high))
X	    (yvals (list vlow vhigh))
X	    (pos (Math-posp vlow))
X	    (levels 0)
X	    (step (math-sub-float high low))
X	    xp yp var-DUMMY)
X	(while (and (<= (setq levels (1+ levels)) 5)
X		    (not found))
X	  (setq xp xvals
X		yp yvals
X		step (math-mul-float step '(float 497 -3)))
X	  (while (and (cdr xp) (not found))
X	    (if (Math-realp (car yp))
X		(setq low (car xp)
X		      vlow (car yp)))
X	    (setq high (math-add-float (car xp) step)
X		  var-DUMMY high
X		  vhigh (math-evaluate-expr expr))
X	    (math-working "search" high)
X	    (if (and (Math-realp vhigh)
X		     (eq (math-negp vhigh) pos))
X		(setq found t)
X	      (setcdr xp (cons high (cdr xp)))
X	      (setcdr yp (cons vhigh (cdr yp)))
X	      (setq xp (cdr (cdr xp))
X		    yp (cdr (cdr yp))))))))
X    (if found
X	(if (Math-zerop vhigh)
X	    (list 'vec high vhigh)
X	  (if (Math-zerop vlow)
X	      (list 'vec low vlow)
X	    (if deriv
X		(math-newton-search-root expr deriv nil nil nil nil
X					 low vlow high vhigh)
X	      (math-bisect-root expr low vlow high vhigh))))
X      (math-reject-arg (list 'intv 3 low high)
X		       "*Unable to find a sign change in this interval")))
)
X
;;; "rtbis"  (but we should be using Brent's method)
(defun math-bisect-root (expr low vlow high vhigh)
X  (let ((step (math-sub-float high low))
X	(pos (Math-posp vhigh))
X	var-DUMMY
X	mid vmid)
X    (while (not (or (math-nearly-equal low
X				       (setq step (math-mul-float
X						   step '(float 5 -1))
X					     mid (math-add-float low step)))
X		    (progn
X		      (setq var-DUMMY mid
X			    vmid (math-evaluate-expr expr))
X		      (Math-zerop vmid))))
X      (math-working "bisect" mid)
X      (if (eq (Math-posp vmid) pos)
X	  (setq high mid
X		vhigh vmid)
X	(setq low mid
X	      vlow vmid)))
X    (list 'vec mid vmid))
)
X
;;; "mnewt"
(defun math-newton-multi (expr jacob n guess orig-guess limit)
X  (let ((m -1)
X	(p guess)
X	p2 expr-val jacob-val next)
X    (while (< (setq p (cdr p) m (1+ m)) n)
X      (set (nth 2 (aref math-root-vars m)) (car p)))
X    (setq expr-val (math-evaluate-expr expr)
X	  jacob-val (math-evaluate-expr jacob))
X    (or (and (math-constp expr-val)
X	     (math-constp jacob-val))
X	(math-reject-arg guess "*Newton's method encountered a singularity"))
X    (setq next (math-add guess (math-div (math-float (math-neg expr-val))
X					 (math-float jacob-val)))
X	  p guess p2 next)
X    (math-working "newton" next)
X    (while (and (setq p (cdr p) p2 (cdr p2))
X		(math-nearly-equal (car p) (car p2))))
X    (if p
X	(if (Math-lessp (math-abs-approx (math-sub next orig-guess))
X			limit)
X	    (math-newton-multi expr jacob n next orig-guess limit)
X	  (math-reject-arg nil "*Newton's method failed to converge"))
X      (list 'vec next expr-val)))
)
X
(defvar math-root-vars [(var DUMMY var-DUMMY)])
X
(defun math-find-root (expr var guess root-widen)
X  (if (eq (car-safe expr) 'vec)
X      (let ((n (1- (length expr)))
X	    (calc-symbolic-mode nil)
X	    (var-DUMMY nil)
X	    (jacob (list 'vec))
X	    p p2 m row)
X	(or (eq (car-safe var) 'vec)
X	    (math-reject-arg var 'vectorp))
X	(or (= (length var) (1+ n))
X	    (math-dimension-error))
X	(setq expr (copy-sequence expr))
X	(while (>= n (length math-root-vars))
X	  (let ((symb (intern (concat "math-root-v"
X				      (int-to-string
X				       (length math-root-vars))))))
X	    (setq math-root-vars (vconcat math-root-vars
X					  (vector (list 'var symb symb))))))
X	(setq m -1)
X	(while (< (setq m (1+ m)) n)
X	  (set (nth 2 (aref math-root-vars m)) nil))
X	(setq m -1 p var)
X	(while (setq m (1+ m) p (cdr p))
X	  (or (eq (car-safe (car p)) 'var)
X	      (math-reject-arg var "*Expected a variable"))
X	  (setq p2 expr)
X	  (while (setq p2 (cdr p2))
X	    (setcar p2 (math-expr-subst (car p2) (car p)
X					(aref math-root-vars m)))))
X	(or (eq (car-safe guess) 'vec)
X	    (math-reject-arg guess 'vectorp))
X	(or (= (length guess) (1+ n))
X	    (math-dimension-error))
X	(setq guess (copy-sequence guess)
X	      p guess)
X	(while (setq p (cdr p))
X	  (or (Math-numberp (car guess))
X	      (math-reject-arg guess 'numberp))
X	  (setcar p (math-float (car p))))
X	(setq p expr)
X	(while (setq p (cdr p))
X	  (if (assq (car-safe (car p)) calc-tweak-eqn-table)
X	      (setcar p (math-sub (nth 1 (car p)) (nth 2 (car p)))))
X	  (setcar p (math-evaluate-expr (car p)))
X	  (setq row (list 'vec)
X		m -1)
X	  (while (< (setq m (1+ m)) n)
X	    (nconc row (list (math-evaluate-expr
X			      (or (calcFunc-deriv (car p)
X						  (aref math-root-vars m)
X						  nil t)
X				  (math-reject-arg
X				   expr
X				   "*Formulas must be differentiable"))))))
X	  (nconc jacob (list row)))
X	(setq m (math-abs-approx guess))
X	(math-newton-multi expr jacob n guess guess
X			   (if (math-zerop m) '(float 1 3) (math-mul m 10))))
X    (or (eq (car-safe var) 'var)
X	(math-reject-arg var "*Expected a variable"))
X    (or (math-expr-contains expr var)
X	(math-reject-arg expr "*Formula does not contain specified variable"))
X    (if (assq (car expr) calc-tweak-eqn-table)
X	(setq expr (math-sub (nth 1 expr) (nth 2 expr))))
X    (math-with-extra-prec 2
X      (setq expr (math-expr-subst expr var '(var DUMMY var-DUMMY)))
X      (let* ((calc-symbolic-mode nil)
X	     (var-DUMMY nil)
X	     (expr (math-evaluate-expr expr))
X	     (deriv (calcFunc-deriv expr '(var DUMMY var-DUMMY) nil t))
X	     low high vlow vhigh)
X	(and deriv (setq deriv (math-evaluate-expr deriv)))
X	(setq guess (math-float guess))
X	(if (and (math-numberp guess)
X		 deriv)
X	    (math-newton-root expr deriv guess guess
X			      (if (math-zerop guess) '(float 1 6)
X				(math-mul (math-abs-approx guess) 100)))
X	  (if (Math-realp guess)
X	      (setq low guess
X		    high guess
X		    var-DUMMY guess
X		    vlow (math-evaluate-expr expr)
X		    vhigh vlow
X		    root-widen 'point)
X	    (if (eq (car guess) 'intv)
X		(progn
X		  (or (math-constp guess) (math-reject-arg guess 'constp))
X		  (setq low (nth 2 guess)
X			high (nth 3 guess))
X		  (if (memq (nth 1 guess) '(0 1))
X		      (setq low (calcFunc-incr low 1 high)))
X		  (if (memq (nth 1 guess) '(0 2))
X		      (setq high (calcFunc-incr high -1 low)))
X		  (setq var-DUMMY low
X			vlow (math-evaluate-expr expr)
X			var-DUMMY high
X			vhigh (math-evaluate-expr expr)))
X	      (if (math-complexp guess)
X		  (math-reject-arg "*Complex root finder must have derivative")
X		(math-reject-arg guess 'realp))))
X	  (if (Math-zerop vlow)
X	      (list 'vec low vlow)
X	    (if (Math-zerop vhigh)
X		(list 'vec high vhigh)
X	      (if (and deriv (Math-numberp vlow) (Math-numberp vhigh))
X		  (math-newton-search-root expr deriv nil nil nil nil
X					   low vlow high vhigh)
X		(if (or (and (Math-posp vlow) (Math-posp vhigh))
X			(and (Math-negp vlow) (Math-negp vhigh))
X			(not (Math-numberp vlow))
X			(not (Math-numberp vhigh)))
X		    (math-search-root expr deriv low vlow high vhigh)
X		  (math-bisect-root expr low vlow high vhigh)))))))))
)
X
(defun calcFunc-root (expr var guess)
X  (math-find-root expr var guess nil)
)
X
(defun calcFunc-wroot (expr var guess)
X  (math-find-root expr var guess t)
)
X
X
X
X
;;; The following algorithms come from Numerical Recipes, chapter 10.
X
(defun math-min-eval (expr a)
X  (if (Math-vectorp a)
X      (let ((m -1))
X	(while (setq m (1+ m) a (cdr a))
X	  (set (nth 2 (aref math-min-vars m)) (car a))))
X    (setq var-DUMMY a))
X  (setq a (math-evaluate-expr expr))
X  (if (Math-ratp a)
X      (math-float a)
X    (if (eq (car a) 'float)
X	a
X      (math-reject-arg a 'realp)))
)
X
X
;;; A bracket for a minimum is a < b < c where f(b) < f(a) and f(b) < f(c).
X
;;; "mnbrak"
(defun math-widen-min (expr a b)
X  (let ((done nil)
X	(iters 30)
X	incr c va vb vc u vu r q ulim bc ba qr)
X    (or b (setq b (math-mul a '(float 101 -2))))
X    (setq va (math-min-eval expr a)
X	  vb (math-min-eval expr b))
X    (if (math-lessp-float va vb)
X	(setq u a a b b u
X	      vu va va vb vb vu))
X    (setq c (math-add-float b (math-mul-float '(float 161803 -5)
X					      (math-sub-float b a)))
X	  vc (math-min-eval expr c))
X    (while (and (not done) (math-lessp-float vc vb))
X      (math-working "widen" (list 'intv 0 a c))
X      (if (= (setq iters (1- iters)) 0)
X	  (math-reject-arg nil (format "*Unable to find a %s near the interval"
X				       math-min-or-max)))
X      (setq bc (math-sub-float b c)
X	    ba (math-sub-float b a)
X	    r (math-mul-float ba (math-sub-float vb vc))
X	    q (math-mul-float bc (math-sub-float vb va))
X	    qr (math-sub-float q r))
X      (if (math-lessp-float (math-abs qr) '(float 1 -20))
X	  (setq qr (if (math-negp qr) '(float -1 -20) '(float 1 -20))))
X      (setq u (math-sub-float
X	       b
X	       (math-div-float (math-sub-float (math-mul-float bc q)
X					       (math-mul-float ba r))
X			       (math-mul-float '(float 2 0) qr)))
X	    ulim (math-add-float b (math-mul-float '(float -1 2) bc))
X	    incr (math-negp bc))
X      (if (if incr (math-lessp-float b u) (math-lessp-float u b))
X	  (if (if incr (math-lessp-float u c) (math-lessp-float c u))
X	      (if (math-lessp-float (setq vu (math-min-eval expr u)) vc)
X		  (setq a b  va vb
X			b u  vb vu
X			done t)
X		(if (math-lessp-float vb vu)
X		    (setq c u  vc vu
X			  done t)
X		  (setq u (math-add-float c (math-mul-float '(float -161803 -5)
X							    bc))
X			vu (math-min-eval expr u))))
X	    (if (if incr (math-lessp-float u ulim) (math-lessp-float ulim u))
X		(if (math-lessp-float (setq vu (math-min-eval expr u)) vc)
X		    (setq b c  vb vc
X			  c u  vc vu
X			  u (math-add-float c (math-mul-float
X					       '(float -161803 -5)
X					       (math-sub-float b c)))
X			  vu (math-min-eval expr u)))
X	      (setq u ulim
X		    vu (math-min-eval expr u))))
X	(setq u (math-add-float c (math-mul-float '(float -161803 -5)
X						  bc))
X	      vu (math-min-eval expr u)))
X      (setq a b  va vb
X	    b c  vb vc
X	    c u  vc vu))
X    (if (math-lessp-float a c)
X	(list a va b vb c vc)
X      (list c vc b vb a va)))
)
X
(defun math-narrow-min (expr a c intv)
X  (let ((xvals (list a c))
X	(yvals (list (math-min-eval expr a)
X		     (math-min-eval expr c)))
X	(levels 0)
X	(step (math-sub-float c a))
X	(found nil)
X	xp yp b)
X    (while (and (<= (setq levels (1+ levels)) 5)
X		(not found))
X      (setq xp xvals
X	    yp yvals
X	    step (math-mul-float step '(float 497 -3)))
X      (while (and (cdr xp) (not found))
X	(setq b (math-add-float (car xp) step))
X	(math-working "search" b)
X	(setcdr xp (cons b (cdr xp)))
X	(setcdr yp (cons (math-min-eval expr b) (cdr yp)))
X	(if (and (math-lessp-float (nth 1 yp) (car yp))
X		 (math-lessp-float (nth 1 yp) (nth 2 yp)))
X	    (setq found t)
X	  (setq xp (cdr xp)
X		yp (cdr yp))
X	  (if (and (cdr (cdr yp))
X		   (math-lessp-float (nth 1 yp) (car yp))
X		   (math-lessp-float (nth 1 yp) (nth 2 yp)))
X	      (setq found t)
X	    (setq xp (cdr xp)
X		  yp (cdr yp))))))
X    (if found
X	(list (car xp) (car yp)
X	      (nth 1 xp) (nth 1 yp)
X	      (nth 2 xp) (nth 2 yp))
X      (or (if (math-lessp-float (car yvals) (nth 1 yvals))
X	      (and (memq (nth 1 intv) '(2 3))
X		   (let ((min (car yvals)))
X		     (while (and (setq yvals (cdr yvals))
X				 (math-lessp-float min (car yvals))))
X		     (and (not yvals)
X			  (list (nth 2 intv) min))))
X	    (and (memq (nth 1 intv) '(1 3))
X		 (setq yvals (nreverse yvals))
X		 (let ((min (car yvals)))
X		   (while (and (setq yvals (cdr yvals))
X			       (math-lessp-float min (car yvals))))
X		   (and (not yvals)
X			(list (nth 3 intv) min)))))
X	  (math-reject-arg nil (format "*Unable to find a %s in the interval"
X				       math-min-or-max)))))
)
X
;;; "brent"
(defun math-brent-min (expr prec a va x vx b vb)
X  (let ((iters (+ 20 (* 5 prec)))
X	(w x)
X	(vw vx)
X	(v x)
X	(vv vx)
X	(tol (list 'float 1 (- -1 prec)))
X	(zeps (list 'float 1 (- -5 prec)))
X	(e '(float 0 0))
X	u vu xm tol1 tol2 etemp p q r xv xw)
X    (while (progn
X	     (setq xm (math-mul-float '(float 5 -1)
X				      (math-add-float a b))
X		   tol1 (math-add-float
X			 zeps
X			 (math-mul-float tol (math-abs x)))
X		   tol2 (math-mul-float tol1 '(float 2 0)))
X	     (math-lessp-float (math-sub-float tol2
X					       (math-mul-float
X						'(float 5 -1)
X						(math-sub-float b a)))
X			       (math-abs (math-sub-float x xm))))
X      (if (= (setq iters (1- iters)) 0)
X	  (math-reject-arg nil (format "*Unable to converge on a %s"
X				       math-min-or-max)))
X      (math-working "brent" x)
X      (if (math-lessp-float (math-abs e) tol1)
X	  (setq e (if (math-lessp-float x xm)
X		      (math-sub-float b x)
X		    (math-sub-float a x))
X		d (math-mul-float '(float 381966 -6) e))
X	(setq xw (math-sub-float x w)
X	      r (math-mul-float xw (math-sub-float vx vv))
X	      xv (math-sub-float x v)
X	      q (math-mul-float xv (math-sub-float vx vw))
X	      p (math-sub-float (math-mul-float xv q)
X				(math-mul-float xw r))
X	      q (math-mul-float '(float 2 0) (math-sub-float q r)))
X	(if (math-posp q)
X	    (setq p (math-neg-float p))
X	  (setq q (math-neg-float q)))
X	(setq etemp e
X	      e d)
X	(if (and (math-lessp-float (math-abs p)
X				   (math-abs (math-mul-float
X					      '(float 5 -1)
X					      (math-mul-float q etemp))))
X		 (math-lessp-float (math-mul-float
X				    q (math-sub-float a x)) p)
X		 (math-lessp-float p (math-mul-float
X				      q (math-sub-float b x))))
X	    (progn
X	      (setq d (math-div-float p q)
X		    u (math-add-float x d))
X	      (if (or (math-lessp-float (math-sub-float u a) tol2)
X		      (math-lessp-float (math-sub-float b u) tol2))
X		  (setq d (if (math-lessp-float xm x)
X			      (math-neg-float tol1)
X			    tol1))))
X	  (setq e (if (math-lessp-float x xm)
X		      (math-sub-float b x)
X		    (math-sub-float a x))
X		d (math-mul-float '(float 381966 -6) e))))
X      (setq u (math-add-float x
X			      (if (math-lessp-float (math-abs d) tol1)
X				  (if (math-negp d)
X				      (math-neg-float tol1)
X				    tol1)
X				d))
X	    vu (math-min-eval expr u))
X      (if (math-lessp-float vx vu)
X	  (progn
X	    (if (math-lessp-float u x)
X		(setq a u)
X	      (setq b u))
X	    (if (or (equal w x)
X		    (not (math-lessp-float vw vu)))
X		(setq v w  vv vw
X		      w u  vw vu)
X	      (if (or (equal v x)
X		      (equal v w)
X		      (not (math-lessp-float vv vu)))
X		  (setq v u  vv vu))))
X	(if (math-lessp-float u x)
X	    (setq b x)
X	  (setq a x))
X	(setq v w  vv vw
X	      w x  vw vx
X	      x u  vx vu)))
X    (list 'vec x vx))
)
X
;;; "powell"
(defun math-powell-min (expr n guesses prec)
X  (let* ((f1dim (math-line-min-func expr n))
X	 (xi (calcFunc-idn 1 n))
X	 (p (cons 'vec (mapcar 'car guesses)))
X	 (pt p)
X	 (ftol (list 'float 1 (- prec)))
X	 (fret (math-min-eval expr p))
X	 fp ptt fptt xit i ibig del diff res)
X    (while (progn
X	     (setq fp fret
X		   ibig 0
X		   del '(float 0 0)
X		   i 0)
X	     (while (<= (setq i (1+ i)) n)
X	       (setq fptt fret
X		     res (math-line-min f1dim p
X					(math-mat-col xi i)
X					n prec)
X		     p (let ((calc-internal-prec prec))
X			 (math-normalize (car res)))
X		     fret (nth 2 res)
X		     diff (math-abs (math-sub-float fptt fret)))
X	       (if (math-lessp-float del diff)
X		   (setq del diff
X			 ibig i)))
X	     (math-lessp-float
X	      (math-mul-float ftol
X			      (math-add-float (math-abs fp)
X					      (math-abs fret)))
X	      (math-mul-float '(float 2 0)
X			      (math-abs (math-sub-float fp
X							fret)))))
X      (setq ptt (math-sub (math-mul '(float 2 0) p) pt)
X	    xit (math-sub p pt)
X	    pt p
X	    fptt (math-min-eval expr ptt))
X      (if (and (math-lessp-float fptt fp)
X	       (math-lessp-float
X		(math-mul-float
X		 (math-mul-float '(float 2 0)
X				 (math-add-float
X				  (math-sub-float fp
X						  (math-mul-float '(float 2 0)
X								  fret))
X				  fptt))
X		 (math-sqr-float (math-sub-float
X				  (math-sub-float fp fret) del)))
X		(math-mul-float del
X				(math-sqr-float (math-sub-float fp fptt)))))
X	  (progn
X	    (setq res (math-line-min f1dim p xit n prec)
X		  p (car res)
X		  fret (nth 2 res)
X		  i 0)
X	    (while (<= (setq i (1+ i)) n)
X	      (setcar (nthcdr ibig (nth i xi))
X		      (nth i (nth 1 res)))))))
X    (list 'vec p fret))
)
X
(defun math-line-min-func (expr n)
X  (let ((m -1))
X    (while (< (setq m (1+ m)) n)
X      (set (nth 2 (aref math-min-vars m))
X	   (list '+
X		 (list '*
X		       '(var DUMMY var-DUMMY)
X		       (list 'calcFunc-mrow '(var line-xi line-xi) (1+ m)))
X		 (list 'calcFunc-mrow '(var line-p line-p) (1+ m)))))
X    (math-evaluate-expr expr))
)
X
(defun math-line-min (f1dim line-p line-xi n prec)
X  (let* ((var-DUMMY nil)
X	 (expr (math-evaluate-expr f1dim))
X	 (params (math-widen-min expr '(float 0 0) '(float 1 0)))
X	 (res (apply 'math-brent-min expr prec params))
X	 (xi (math-mul (nth 1 res) line-xi)))
X    (list (math-add line-p xi) xi (nth 2 res)))
)
X
X
(defvar math-min-vars [(var DUMMY var-DUMMY)])
X
(defun math-find-minimum (expr var guess min-widen)
X  (let* ((calc-symbolic-mode nil)
X	 (n 0)
X	 (var-DUMMY nil)
X	 (isvec (math-vectorp var))
X	 g guesses)
X    (or (math-vectorp var)
X	(setq var (list 'vec var)))
X    (or (math-vectorp guess)
X	(setq guess (list 'vec guess)))
X    (or (= (length var) (length guess))
X	(math-dimension-error))
X    (while (setq var (cdr var) guess (cdr guess))
X      (or (eq (car-safe (car var)) 'var)
X	  (math-reject-arg (car vg) "*Expected a variable"))
X      (or (math-expr-contains expr (car var))
X	  (math-reject-arg (car var)
X			   "*Formula does not contain specified variable"))
X      (while (>= (1+ n) (length math-min-vars))
X	(let ((symb (intern (concat "math-min-v"
X				    (int-to-string
X				     (length math-min-vars))))))
X	  (setq math-min-vars (vconcat math-min-vars
X				       (vector (list 'var symb symb))))))
X      (set (nth 2 (aref math-min-vars n)) nil)
X      (set (nth 2 (aref math-min-vars (1+ n))) nil)
X      (if (math-complexp (car guess))
X	  (setq expr (math-expr-subst expr
X				      (car var)
X				      (list '+ (aref math-min-vars n)
X					    (list '*
X						  (aref math-min-vars (1+ n))
X						  '(cplx 0 1))))
X		guesses (let ((g (math-float (math-complex (car guess)))))
X			  (cons (list (nth 2 g) nil nil)
X				(cons (list (nth 1 g) nil nil t)
X				      guesses)))
X		n (+ n 2))
X	(setq expr (math-expr-subst expr
X				    (car var)
X				    (aref math-min-vars n))
X	      guesses (cons (if (math-realp (car guess))
X				(list (math-float (car guess)) nil nil)
X			      (if (and (eq (car-safe (car guess)) 'intv)
X				       (math-constp (car guess)))
X				  (list (math-mul
X					 (math-add (nth 2 (car guess))
X						   (nth 3 (car guess)))
X					 '(float 5 -1))
X					(math-float (nth 2 (car guess)))
X					(math-float (nth 3 (car guess)))
X					(car guess))
X				(math-reject-arg (car guess) 'realp)))
X			    guesses)
X	      n (1+ n))))
X    (setq guesses (nreverse guesses)
X	  expr (math-evaluate-expr expr))
X    (if (= n 1)
X	(let* ((params (if (nth 1 (car guesses))
X			   (if min-widen
X			       (math-widen-min expr
X					       (nth 1 (car guesses))
X					       (nth 2 (car guesses)))
X			     (math-narrow-min expr
X					      (nth 1 (car guesses))
X					      (nth 2 (car guesses))
X					      (nth 3 (car guesses))))
X			 (math-widen-min expr
X					 (car (car guesses))
X					 nil)))
X	       (prec calc-internal-prec)
X	       (res (if (cdr (cdr params))
X			(math-with-extra-prec (+ calc-internal-prec 2)
X			  (apply 'math-brent-min expr prec params))
X		      (cons 'vec params))))
X	  (if isvec
X	      (list 'vec (list 'vec (nth 1 res)) (nth 2 res))
X	    res))
X      (let* ((prec calc-internal-prec)
X	     (res (math-with-extra-prec (+ calc-internal-prec 2)
X		    (math-powell-min expr n guesses prec)))
X	     (p (nth 1 res))
X	     (vec (list 'vec)))
X	(while (setq p (cdr p))
X	  (if (nth 3 (car guesses))
X	      (progn
X		(nconc vec (list (math-normalize
X				  (list 'cplx (car p) (nth 1 p)))))
X		(setq p (cdr p)
X		      guesses (cdr guesses)))
X	    (nconc vec (list (car p))))
X	  (setq guesses (cdr guesses)))
X	(if isvec
X	    (list 'vec vec (nth 2 res))
X	  (list 'vec (nth 1 vec) (nth 2 res))))))
)
(setq math-min-or-max "minimum")
X
(defun calcFunc-minimize (expr var guess)
X  (let ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
X	(math-min-or-max "minimum"))
X    (math-find-minimum (math-normalize expr)
X		       (math-normalize var)
X		       (math-normalize guess) nil))
)
X
(defun calcFunc-wminimize (expr var guess)
X  (let ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
X	(math-min-or-max "minimum"))
X    (math-find-minimum (math-normalize expr)
X		       (math-normalize var)
X		       (math-normalize guess) t))
)
X
(defun calcFunc-maximize (expr var guess)
X  (let* ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
X	 (math-min-or-max "maximum")
X	 (res (math-find-minimum (math-normalize (math-neg expr))
X				 (math-normalize var)
X				 (math-normalize guess) nil)))
X    (list 'vec (nth 1 res) (math-neg (nth 2 res))))
)
X
(defun calcFunc-wmaximize (expr var guess)
X  (let* ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
X	 (math-min-or-max "maximum")
X	 (res (math-find-minimum (math-normalize (math-neg expr))
X				 (math-normalize var)
X				 (math-normalize guess) t)))
X    (list 'vec (nth 1 res) (math-neg (nth 2 res))))
)
X
X
X
X
;;; The following algorithms come from Numerical Recipes, chapter 3.
X
(defun calcFunc-polint (data x)
X  (or (math-matrixp data) (math-reject-arg data 'matrixp))
X  (or (= (length data) 3)
X      (math-reject-arg data "*Wrong number of data rows"))
X  (or (> (length (nth 1 data)) 2)
X      (math-reject-arg data "*Too few data points"))
X  (if (and (math-vectorp x) (or (math-constp x) math-expand-formulas))
X      (cons 'vec (mapcar (function (lambda (x) (calcFunc-polint data x)))
X			 (cdr x)))
X    (or (math-objectp x) math-expand-formulas (math-reject-arg x 'objectp))
X    (math-with-extra-prec 2
X      (cons 'vec (math-poly-interp (cdr (nth 1 data)) (cdr (nth 2 data)) x
X				   nil))))
)
(put 'calcFunc-polint 'math-expandable t)
X
X
(defun calcFunc-ratint (data x)
X  (or (math-matrixp data) (math-reject-arg data 'matrixp))
X  (or (= (length data) 3)
X      (math-reject-arg data "*Wrong number of data rows"))
X  (or (> (length (nth 1 data)) 2)
X      (math-reject-arg data "*Too few data points"))
X  (if (and (math-vectorp x) (or (math-constp x) math-expand-formulas))
X      (cons 'vec (mapcar (function (lambda (x) (calcFunc-ratint data x)))
X			 (cdr x)))
X    (or (math-objectp x) math-expand-formulas (math-reject-arg x 'objectp))
X    (math-with-extra-prec 2
X      (cons 'vec (math-poly-interp (cdr (nth 1 data)) (cdr (nth 2 data)) x
X				   (cdr (cdr (cdr (nth 1 data))))))))
)
(put 'calcFunc-ratint 'math-expandable t)
X
X
(defun math-poly-interp (xa ya x ratp)
X  (let ((n (length xa))
X	(dif nil)
X	(ns nil)
X	(xax nil)
X	(c (copy-sequence ya))
X	(d (copy-sequence ya))
X	(i 0)
X	(m 0)
X	y dy (xp xa) xpm cp dp temp)
X    (while (<= (setq i (1+ i)) n)
X      (setq xax (cons (math-sub (car xp) x) xax)
X	    xp (cdr xp)
X	    temp (math-abs (car xax)))
X      (if (or (null dif) (math-lessp temp dif))
X	  (setq dif temp
X		ns i)))
X    (setq xax (nreverse xax)
X	  ns (1- ns)
X	  y (nth ns ya))
X    (if (math-zerop dif)
X	(list y 0)
X      (while (< (setq m (1+ m)) n)
X	(setq i 0
X	      xp xax
X	      xpm (nthcdr m xax)
X	      cp c
X	      dp d)
X	(while (<= (setq i (1+ i)) (- n m))
X	  (if ratp
X	      (let ((t2 (math-div (math-mul (car xp) (car dp)) (car xpm))))
X		(setq temp (math-div (math-sub (nth 1 cp) (car dp))
X				     (math-sub t2 (nth 1 cp))))
X		(setcar dp (math-mul (nth 1 cp) temp))
X		(setcar cp (math-mul t2 temp)))
X	    (if (math-equal (car xp) (car xpm))
X		(math-reject-arg (cons 'vec xa) "*Duplicate X values"))
X	    (setq temp (math-div (math-sub (nth 1 cp) (car dp))
X				 (math-sub (car xp) (car xpm))))
X	    (setcar dp (math-mul (car xpm) temp))
X	    (setcar cp (math-mul (car xp) temp)))
X	  (setq cp (cdr cp)
X		dp (cdr dp)
X		xp (cdr xp)
X		xpm (cdr xpm)))
X	(if (< (+ ns ns) (- n m))
X	    (setq dy (nth ns c))
X	  (setq ns (1- ns)
X		dy (nth ns d)))
X	(setq y (math-add y dy)))
X      (list y dy)))
)
X
X
X
;;; The following algorithms come from Numerical Recipes, chapter 4.
X
(defun calcFunc-ninteg (expr var lo hi)
X  (setq lo (math-evaluate-expr lo)
X	hi (math-evaluate-expr hi))
X  (or (math-numberp lo) (math-infinitep lo) (math-reject-arg lo 'numberp))
X  (or (math-numberp hi) (math-infinitep hi) (math-reject-arg hi 'numberp))
X  (if (math-lessp hi lo)
X      (math-neg (calcFunc-ninteg expr var hi lo))
X    (setq expr (math-expr-subst expr var '(var DUMMY var-DUMMY)))
X    (let ((var-DUMMY nil)
X	  (calc-symbolic-mode nil)
X	  (calc-prefer-frac nil)
X	  (sum 0))
X      (setq expr (math-evaluate-expr expr))
X      (if (equal lo '(neg (var inf var-inf)))
X	  (let ((thi (if (math-lessp hi '(float -2 0))
X			 hi '(float -2 0))))
X	    (setq sum (math-ninteg-romberg
X		       'math-ninteg-midpoint expr
X			 (math-float lo) (math-float thi) 'inf)
X		  lo thi)))
X      (if (equal hi '(var inf var-inf))
X	  (let ((tlo (if (math-lessp '(float 2 0) lo)
X			 lo '(float 2 0))))
X	    (setq sum (math-add sum
X				(math-ninteg-romberg
X				 'math-ninteg-midpoint expr
SHAR_EOF
true || echo 'restore of calc-alg-3.el failed'
fi
echo 'End of  part 6'
echo 'File calc-alg-3.el is continued in part 7'
echo 7 > _shar_seq_.tmp
exit 0
exit 0 # Just in case...
-- 
Kent Landfield                   INTERNET: kent@sparky.IMD.Sterling.COM
Sterling Software, IMD           UUCP:     uunet!sparky!kent
Phone:    (402) 291-8300         FAX:      (402) 291-4362
Please send comp.sources.misc-related mail to kent@uunet.uu.net.