Usenet 1994 January/usenetsourcesnewsgroupsinfomagicjanuary1994.iso/sources/misc/volume24/gnucalc/part06

Text File | 1991-10-28 | 55.2 KB | 1,737 lines

Newsgroups: comp.sources.misc
From: [e-mail] (David Gillespie)
Subject:  v24i054:  gnucalc - GNU Emacs Calculator, v2.00, Part06/56
Message-ID: <[e-mail]>
X-Md4-Signature: 9f30cde948d4d7fb5331e55ee96320de
Date: Tue, 29 Oct 1991 04:25:14 GMT
Approved: [e-mail]

Submitted-by: [e-mail] (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, [e-mail].
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: [e-mail]
Sterling Software, IMD           UUCP:     uunet!sparky!kent
Phone:    [phone]         FAX:      [phone]
Please send comp.sources.misc-related mail to [e-mail].