lisp/calc/calcalg3.el

   1 ;;; calcalg3.el --- more algebraic functions for Calc
   2
   3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.
   4
   5 ;; Author: David Gillespie <daveg@synaptics.com>
   6 ;; Maintainers: D. Goel <deego@gnufans.org>
   7 ;;              Colin Walters <walters@debian.org>
   8
   9 ;; This file is part of GNU Emacs.
  10
  11 ;; GNU Emacs is distributed in the hope that it will be useful,
  12 ;; but WITHOUT ANY WARRANTY.  No author or distributor
  13 ;; accepts responsibility to anyone for the consequences of using it
  14 ;; or for whether it serves any particular purpose or works at all,
  15 ;; unless he says so in writing.  Refer to the GNU Emacs General Public
  16 ;; License for full details.
  17
  18 ;; Everyone is granted permission to copy, modify and redistribute
  19 ;; GNU Emacs, but only under the conditions described in the
  20 ;; GNU Emacs General Public License.   A copy of this license is
  21 ;; supposed to have been given to you along with GNU Emacs so you
  22 ;; can know your rights and responsibilities.  It should be in a
  23 ;; file named COPYING.  Among other things, the copyright notice
  24 ;; and this notice must be preserved on all copies.
  25
  26 ;;; Commentary:
  27
  28 ;;; Code:
  29
  30 ;; This file is autoloaded from calc-ext.el.
  31 (require 'calc-ext)
  32
  33 (require 'calc-macs)
  34
  35 (defun calc-Need-calc-alg-3 () nil)
  36
  37
  38 (defun calc-find-root (var)
  39   (interactive "sVariable(s) to solve for: ")
  40   (calc-slow-wrapper
  41    (let ((func (if (calc-is-hyperbolic) 'calcFunc-wroot 'calcFunc-root)))
  42      (if (or (equal var "") (equal var "$"))
  43          (calc-enter-result 2 "root" (list func
  44                                            (calc-top-n 3)
  45                                            (calc-top-n 1)
  46                                            (calc-top-n 2)))
  47        (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
  48                            (not (string-match "\\[" var)))
  49                       (math-read-expr (concat "[" var "]"))
  50                     (math-read-expr var))))
  51          (if (eq (car-safe var) 'error)
  52              (error "Bad format in expression: %s" (nth 1 var)))
  53          (calc-enter-result 1 "root" (list func
  54                                            (calc-top-n 2)
  55                                            var
  56                                            (calc-top-n 1))))))))
  57
  58 (defun calc-find-minimum (var)
  59   (interactive "sVariable(s) to minimize over: ")
  60   (calc-slow-wrapper
  61    (let ((func (if (calc-is-inverse)
  62                    (if (calc-is-hyperbolic)
  63                        'calcFunc-wmaximize 'calcFunc-maximize)
  64                  (if (calc-is-hyperbolic)
  65                      'calcFunc-wminimize 'calcFunc-minimize)))
  66          (tag (if (calc-is-inverse) "max" "min")))
  67      (if (or (equal var "") (equal var "$"))
  68          (calc-enter-result 2 tag (list func
  69                                         (calc-top-n 3)
  70                                         (calc-top-n 1)
  71                                         (calc-top-n 2)))
  72        (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
  73                            (not (string-match "\\[" var)))
  74                       (math-read-expr (concat "[" var "]"))
  75                     (math-read-expr var))))
  76          (if (eq (car-safe var) 'error)
  77              (error "Bad format in expression: %s" (nth 1 var)))
  78          (calc-enter-result 1 tag (list func
  79                                         (calc-top-n 2)
  80                                         var
  81                                         (calc-top-n 1))))))))
  82
  83 (defun calc-find-maximum (var)
  84   (interactive "sVariable to maximize over: ")
  85   (calc-invert-func)
  86   (calc-find-minimum var))
  87
  88
  89 (defun calc-poly-interp (arg)
  90   (interactive "P")
  91   (calc-slow-wrapper
  92    (let ((data (calc-top 2)))
  93      (if (or (consp arg) (eq arg 0) (eq arg 2))
  94          (setq data (cons 'vec (calc-top-list 2 2)))
  95        (or (null arg)
  96            (error "Bad prefix argument")))
  97      (if (calc-is-hyperbolic)
  98          (calc-enter-result 1 "rati" (list 'calcFunc-ratint data (calc-top 1)))
  99        (calc-enter-result 1 "poli" (list 'calcFunc-polint data
 100                                          (calc-top 1)))))))
 101
 102
 103 (defun calc-curve-fit (arg &optional model coefnames varnames)
 104   (interactive "P")
 105   (calc-slow-wrapper
 106    (setq calc-aborted-prefix nil)
 107    (let ((func (if (calc-is-inverse) 'calcFunc-xfit
 108                  (if (calc-is-hyperbolic) 'calcFunc-efit
 109                    'calcFunc-fit)))
 110          key (which 0)
 111          n nvars temp data
 112          (homog nil)
 113          (msgs '( "(Press ? for help)"
 114                   "1 = linear or multilinear"
 115                   "2-9 = polynomial fits; i = interpolating polynomial"
 116                   "p = a x^b, ^ = a b^x"
 117                   "e = a exp(b x), x = exp(a + b x), l = a + b ln(x)"
 118                   "E = a 10^(b x), X = 10^(a + b x), L = a + b log10(x)"
 119                   "q = a + b (x-c)^2"
 120                   "g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)"
 121                   "h prefix = homogeneous model (no constant term)"
 122                   "' = alg entry, $ = stack, u = Model1, U = Model2")))
 123      (while (not model)
 124        (message "Fit to model: %s:%s"
 125                 (nth which msgs)
 126                 (if homog " h" ""))
 127        (setq key (read-char))
 128        (cond ((= key ?\C-g)
 129               (keyboard-quit))
 130              ((= key ??)
 131               (setq which (% (1+ which) (length msgs))))
 132              ((memq key '(?h ?H))
 133               (setq homog (not homog)))
 134              ((progn
 135                 (if (eq key ?\$)
 136                     (setq n 1)
 137                   (setq n 0))
 138                 (cond ((null arg)
 139                        (setq n (1+ n)
 140                              data (calc-top n)))
 141                       ((or (consp arg) (eq arg 0))
 142                        (setq n (+ n 2)
 143                              data (calc-top n)
 144                              data (if (math-matrixp data)
 145                                       (append data (list (calc-top (1- n))))
 146                                     (list 'vec data (calc-top (1- n))))))
 147                       ((> (setq arg (prefix-numeric-value arg)) 0)
 148                        (setq data (cons 'vec (calc-top-list arg (1+ n)))
 149                              n (+ n arg)))
 150                       (t (error "Bad prefix argument")))
 151                 (or (math-matrixp data) (not (cdr (cdr data)))
 152                     (error "Data matrix is not a matrix!"))
 153                 (setq nvars (- (length data) 2)
 154                       coefnames nil
 155                       varnames nil)
 156                 nil))
 157              ((= key ?1)  ; linear or multilinear
 158               (calc-get-fit-variables nvars (1+ nvars) (and homog 0))
 159               (setq model (math-mul coefnames
 160                                     (cons 'vec (cons 1 (cdr varnames))))))
 161              ((and (>= key ?2) (<= key ?9))   ; polynomial
 162               (calc-get-fit-variables 1 (- key ?0 -1) (and homog 0))
 163               (setq model (math-build-polynomial-expr (cdr coefnames)
 164                                                       (nth 1 varnames))))
 165              ((= key ?i)  ; exact polynomial
 166               (calc-get-fit-variables 1 (1- (length (nth 1 data)))
 167                                       (and homog 0))
 168               (setq model (math-build-polynomial-expr (cdr coefnames)
 169                                                       (nth 1 varnames))))
 170              ((= key ?p)  ; power law
 171               (calc-get-fit-variables nvars (1+ nvars) (and homog 1))
 172               (setq model (math-mul (nth 1 coefnames)
 173                                     (calcFunc-reduce
 174                                      '(var mul var-mul)
 175                                      (calcFunc-map
 176                                       '(var pow var-pow)
 177                                       varnames
 178                                       (cons 'vec (cdr (cdr coefnames))))))))
 179              ((= key ?^)  ; exponential law
 180               (calc-get-fit-variables nvars (1+ nvars) (and homog 1))
 181               (setq model (math-mul (nth 1 coefnames)
 182                                     (calcFunc-reduce
 183                                      '(var mul var-mul)
 184                                      (calcFunc-map
 185                                       '(var pow var-pow)
 186                                       (cons 'vec (cdr (cdr coefnames)))
 187                                       varnames)))))
 188              ((memq key '(?e ?E))
 189               (calc-get-fit-variables nvars (1+ nvars) (and homog 1))
 190               (setq model (math-mul (nth 1 coefnames)
 191                                     (calcFunc-reduce
 192                                      '(var mul var-mul)
 193                                      (calcFunc-map
 194                                       (if (eq key ?e)
 195                                           '(var exp var-exp)
 196                                         '(calcFunc-lambda
 197                                           (var a var-a)
 198                                           (^ 10 (var a var-a))))
 199                                       (calcFunc-map
 200                                        '(var mul var-mul)
 201                                        (cons 'vec (cdr (cdr coefnames)))
 202                                        varnames))))))
 203              ((memq key '(?x ?X))
 204               (calc-get-fit-variables nvars (1+ nvars) (and homog 0))
 205               (setq model (math-mul coefnames
 206                                     (cons 'vec (cons 1 (cdr varnames)))))
 207               (setq model (if (eq key ?x)
 208                               (list 'calcFunc-exp model)
 209                             (list '^ 10 model))))
 210              ((memq key '(?l ?L))
 211               (calc-get-fit-variables nvars (1+ nvars) (and homog 0))
 212               (setq model (math-mul coefnames
 213                                     (cons 'vec
 214                                           (cons 1 (cdr (calcFunc-map
 215                                                         (if (eq key ?l)
 216                                                             '(var ln var-ln)
 217                                                           '(var log10
 218                                                                 var-log10))
 219                                                         varnames)))))))
 220              ((= key ?q)
 221               (calc-get-fit-variables nvars (1+ (* 2 nvars)) (and homog 0))
 222               (let ((c coefnames)
 223                     (v varnames))
 224                 (setq model (nth 1 c))
 225                 (while (setq v (cdr v) c (cdr (cdr c)))
 226                   (setq model (math-add
 227                                model
 228                                (list '*
 229                                      (car c)
 230                                      (list '^
 231                                            (list '- (car v) (nth 1 c))
 232                                            2)))))))
 233              ((= key ?g)
 234               (setq model (math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)")
 235                     varnames '(vec (var XFit var-XFit))
 236                     coefnames '(vec (var AFit var-AFit)
 237                                     (var BFit var-BFit)
 238                                     (var CFit var-CFit)))
 239               (calc-get-fit-variables 1 (1- (length coefnames)) (and homog 1)))
 240              ((memq key '(?\$ ?\' ?u ?U))
 241               (let* ((defvars nil)
 242                      (record-entry nil))
 243                 (if (eq key ?\')
 244                     (let* ((calc-dollar-values calc-arg-values)
 245                            (calc-dollar-used 0)
 246                            (calc-hashes-used 0))
 247                       (setq model (calc-do-alg-entry "" "Model formula: "))
 248                       (if (/= (length model) 1)
 249                           (error "Bad format"))
 250                       (setq model (car model)
 251                             record-entry t)
 252                       (if (> calc-dollar-used 0)
 253                           (setq coefnames
 254                                 (cons 'vec
 255                                       (nthcdr (- (length calc-arg-values)
 256                                                  calc-dollar-used)
 257                                               (reverse calc-arg-values))))
 258                         (if (> calc-hashes-used 0)
 259                             (setq coefnames
 260                                   (cons 'vec (calc-invent-args
 261                                               calc-hashes-used))))))
 262                   (progn
 263                     (setq model (cond ((eq key ?u)
 264                                        (calc-var-value 'var-Model1))
 265                                       ((eq key ?U)
 266                                        (calc-var-value 'var-Model2))
 267                                       (t (calc-top 1))))
 268                     (or model (error "User model not yet defined"))
 269                     (if (math-vectorp model)
 270                         (if (and (memq (length model) '(3 4))
 271                                  (not (math-objvecp (nth 1 model)))
 272                                  (math-vectorp (nth 2 model))
 273                                  (or (null (nth 3 model))
 274                                      (math-vectorp (nth 3 model))))
 275                             (setq varnames (nth 2 model)
 276                                   coefnames (or (nth 3 model)
 277                                                 (cons 'vec
 278                                                       (math-all-vars-but
 279                                                        model varnames)))
 280                                   model (nth 1 model))
 281                           (error "Incorrect model specifier")))))
 282                 (or varnames
 283                     (let ((with-y (eq (car-safe model) 'calcFunc-eq)))
 284                       (if coefnames
 285                           (calc-get-fit-variables (if with-y (1+ nvars) nvars)
 286                                                   (1- (length coefnames))
 287                                                   (math-all-vars-but
 288                                                    model coefnames)
 289                                                   nil with-y)
 290                         (let* ((coefs (math-all-vars-but model nil))
 291                                (vars nil)
 292                                (n (- (length coefs) nvars (if with-y 2 1)))
 293                                p)
 294                           (if (< n 0)
 295                               (error "Not enough variables in model"))
 296                           (setq p (nthcdr n coefs))
 297                           (setq vars (cdr p))
 298                           (setcdr p nil)
 299                           (calc-get-fit-variables (if with-y (1+ nvars) nvars)
 300                                                   (length coefs)
 301                                                   vars coefs with-y)))))
 302                 (if record-entry
 303                     (calc-record (list 'vec model varnames coefnames)
 304                                  "modl"))))
 305              (t (beep))))
 306      (let ((calc-fit-to-trail t))
 307        (calc-enter-result n (substring (symbol-name func) 9)
 308                           (list func model
 309                                 (if (= (length varnames) 2)
 310                                     (nth 1 varnames)
 311                                   varnames)
 312                                 (if (= (length coefnames) 2)
 313                                     (nth 1 coefnames)
 314                                   coefnames)
 315                                 data))
 316        (if (consp calc-fit-to-trail)
 317            (calc-record (calc-normalize calc-fit-to-trail) "parm"))))))
 318
 319 (defun calc-invent-independent-variables (n &optional but)
 320   (calc-invent-variables n but '(x y z t) "x"))
 321
 322 (defun calc-invent-parameter-variables (n &optional but)
 323   (calc-invent-variables n but '(a b c d) "a"))
 324
 325 (defun calc-invent-variables (num but names base)
 326   (let ((vars nil)
 327         (n num) (nn 0)
 328         var)
 329     (while (and (> n 0) names)
 330       (setq var (math-build-var-name (if (consp names)
 331                                          (car names)
 332                                        (concat base (int-to-string
 333                                                      (setq nn (1+ nn)))))))
 334       (or (math-expr-contains (cons 'vec but) var)
 335           (setq vars (cons var vars)
 336                 n (1- n)))
 337       (or (symbolp names) (setq names (cdr names))))
 338     (if (= n 0)
 339         (nreverse vars)
 340       (calc-invent-variables num but t base))))
 341
 342 (defun calc-get-fit-variables (nv nc &optional defv defc with-y homog)
 343   (or (= nv (if with-y (1+ nvars) nvars))
 344       (error "Wrong number of data vectors for this type of model"))
 345   (if (integerp defv)
 346       (setq homog defv
 347             defv nil))
 348   (if homog
 349       (setq nc (1- nc)))
 350   (or defv
 351       (setq defv (calc-invent-independent-variables nv)))
 352   (or defc
 353       (setq defc (calc-invent-parameter-variables nc defv)))
 354   (let ((vars (read-string (format "Fitting variables: (default %s; %s) "
 355                                    (mapconcat 'symbol-name
 356                                               (mapcar (function (lambda (v)
 357                                                                   (nth 1 v)))
 358                                                       defv)
 359                                               ",")
 360                                    (mapconcat 'symbol-name
 361                                               (mapcar (function (lambda (v)
 362                                                                   (nth 1 v)))
 363                                                       defc)
 364                                               ","))))
 365         (coefs nil))
 366     (setq vars (if (string-match "\\[" vars)
 367                    (math-read-expr vars)
 368                  (math-read-expr (concat "[" vars "]"))))
 369     (if (eq (car-safe vars) 'error)
 370         (error "Bad format in expression: %s" (nth 2 vars)))
 371     (or (math-vectorp vars)
 372         (error "Expected a variable or vector of variables"))
 373     (if (equal vars '(vec))
 374         (setq vars (cons 'vec defv)
 375               coefs (cons 'vec defc))
 376       (if (math-vectorp (nth 1 vars))
 377           (if (and (= (length vars) 3)
 378                    (math-vectorp (nth 2 vars)))
 379               (setq coefs (nth 2 vars)
 380                     vars (nth 1 vars))
 381             (error
 382              "Expected independent variables vector, then parameters vector"))
 383         (setq coefs (cons 'vec defc))))
 384     (or (= nv (1- (length vars)))
 385         (and (not with-y) (= (1+ nv) (1- (length vars))))
 386         (error "Expected %d independent variable%s" nv (if (= nv 1) "" "s")))
 387     (or (= nc (1- (length coefs)))
 388         (error "Expected %d parameter variable%s" nc (if (= nc 1) "" "s")))
 389     (if homog
 390         (setq coefs (cons 'vec (cons homog (cdr coefs)))))
 391     (if varnames
 392         (setq model (math-multi-subst model (cdr varnames) (cdr vars))))
 393     (if coefnames
 394         (setq model (math-multi-subst model (cdr coefnames) (cdr coefs))))
 395     (setq varnames vars
 396           coefnames coefs)))
 397
 398
 399
 400
 401 ;;; The following algorithms are from Numerical Recipes chapter 9.
 402
 403 ;;; "rtnewt" with safety kludges
 404 (defun math-newton-root (expr deriv guess orig-guess limit)
 405   (math-working "newton" guess)
 406   (let* ((var-DUMMY guess)
 407          next dval)
 408     (setq next (math-evaluate-expr expr)
 409           dval (math-evaluate-expr deriv))
 410     (if (and (Math-numberp next)
 411              (Math-numberp dval)
 412              (not (Math-zerop dval)))
 413         (progn
 414           (setq next (math-sub guess (math-div next dval)))
 415           (if (math-nearly-equal guess (setq next (math-float next)))
 416               (progn
 417                 (setq var-DUMMY next)
 418                 (list 'vec next (math-evaluate-expr expr)))
 419             (if (Math-lessp (math-abs-approx (math-sub next orig-guess))
 420                             limit)
 421                 (math-newton-root expr deriv next orig-guess limit)
 422               (math-reject-arg next "*Newton's method failed to converge"))))
 423       (math-reject-arg next "*Newton's method encountered a singularity"))))
 424
 425 ;;; Inspired by "rtsafe"
 426 (defun math-newton-search-root (expr deriv guess vguess ostep oostep
 427                                      low vlow high vhigh)
 428   (let ((var-DUMMY guess)
 429         (better t)
 430         pos step next vnext)
 431     (if guess
 432         (math-working "newton" (list 'intv 0 low high))
 433       (math-working "bisect" (list 'intv 0 low high))
 434       (setq ostep (math-mul-float (math-sub-float high low)
 435                                   '(float 5 -1))
 436             guess (math-add-float low ostep)
 437             var-DUMMY guess
 438             vguess (math-evaluate-expr expr))
 439       (or (Math-realp vguess)
 440           (progn
 441             (setq ostep (math-mul-float ostep '(float 6 -1))
 442                   guess (math-add-float low ostep)
 443                   var-DUMMY guess
 444                   vguess (math-evaluate-expr expr))
 445             (or (math-realp vguess)
 446                 (progn
 447                   (setq ostep (math-mul-float ostep '(float 123456 -5))
 448                         guess (math-add-float low ostep)
 449                         var-DUMMY guess
 450                         vguess nil))))))
 451     (or vguess
 452         (setq vguess (math-evaluate-expr expr)))
 453     (or (Math-realp vguess)
 454         (math-reject-arg guess "*Newton's method encountered a singularity"))
 455     (setq vguess (math-float vguess))
 456     (if (eq (Math-negp vlow) (setq pos (Math-posp vguess)))
 457         (setq high guess
 458               vhigh vguess)
 459       (if (eq (Math-negp vhigh) pos)
 460           (setq low guess
 461                 vlow vguess)
 462         (setq better nil)))
 463     (if (or (Math-zerop vguess)
 464             (math-nearly-equal low high))
 465         (list 'vec guess vguess)
 466       (setq step (math-evaluate-expr deriv))
 467       (if (and (Math-realp step)
 468                (not (Math-zerop step))
 469                (setq step (math-div-float vguess (math-float step))
 470                      next (math-sub-float guess step))
 471                (not (math-lessp-float high next))
 472                (not (math-lessp-float next low)))
 473           (progn
 474             (setq var-DUMMY next
 475                   vnext (math-evaluate-expr expr))
 476             (if (or (Math-zerop vnext)
 477                     (math-nearly-equal next guess))
 478                 (list 'vec next vnext)
 479               (if (and better
 480                        (math-lessp-float (math-abs (or oostep
 481                                                        (math-sub-float
 482                                                         high low)))
 483                                          (math-abs
 484                                           (math-mul-float '(float 2 0)
 485                                                           step))))
 486                   (math-newton-search-root expr deriv nil nil nil ostep
 487                                            low vlow high vhigh)
 488                 (math-newton-search-root expr deriv next vnext step ostep
 489                                          low vlow high vhigh))))
 490         (if (or (and (Math-posp vlow) (Math-posp vhigh))
 491                 (and (Math-negp vlow) (Math-negp vhigh)))
 492             (math-search-root expr deriv low vlow high vhigh)
 493           (math-newton-search-root expr deriv nil nil nil ostep
 494                                    low vlow high vhigh))))))
 495
 496 ;;; Search for a root in an interval with no overt zero crossing.
 497 (defun math-search-root (expr deriv low vlow high vhigh)
 498   (let (found)
 499     (if root-widen
 500         (let ((iters 0)
 501               (iterlim (if (eq root-widen 'point)
 502                            (+ calc-internal-prec 10)
 503                          20))
 504               (factor (if (eq root-widen 'point)
 505                           '(float 9 0)
 506                         '(float 16 -1)))
 507               (prev nil) vprev waslow
 508               diff)
 509           (while (or (and (math-posp vlow) (math-posp vhigh))
 510                      (and (math-negp vlow) (math-negp vhigh)))
 511             (math-working "widen" (list 'intv 0 low high))
 512             (if (> (setq iters (1+ iters)) iterlim)
 513                 (math-reject-arg (list 'intv 0 low high)
 514                                  "*Unable to bracket root"))
 515             (if (= iters calc-internal-prec)
 516                 (setq factor '(float 16 -1)))
 517             (setq diff (math-mul-float (math-sub-float high low) factor))
 518             (if (Math-zerop diff)
 519                 (setq high (calcFunc-incr high 10))
 520               (if (math-lessp-float (math-abs vlow) (math-abs vhigh))
 521                   (setq waslow t
 522                         prev low
 523                         low (math-sub low diff)
 524                         var-DUMMY low
 525                         vprev vlow
 526                         vlow (math-evaluate-expr expr))
 527                 (setq waslow nil
 528                       prev high
 529                       high (math-add high diff)
 530                       var-DUMMY high
 531                       vprev vhigh
 532                       vhigh (math-evaluate-expr expr)))))
 533           (if prev
 534               (if waslow
 535                   (setq high prev vhigh vprev)
 536                 (setq low prev vlow vprev)))
 537           (setq found t))
 538       (or (Math-realp vlow)
 539           (math-reject-arg vlow 'realp))
 540       (or (Math-realp vhigh)
 541           (math-reject-arg vhigh 'realp))
 542       (let ((xvals (list low high))
 543             (yvals (list vlow vhigh))
 544             (pos (Math-posp vlow))
 545             (levels 0)
 546             (step (math-sub-float high low))
 547             xp yp var-DUMMY)
 548         (while (and (<= (setq levels (1+ levels)) 5)
 549                     (not found))
 550           (setq xp xvals
 551                 yp yvals
 552                 step (math-mul-float step '(float 497 -3)))
 553           (while (and (cdr xp) (not found))
 554             (if (Math-realp (car yp))
 555                 (setq low (car xp)
 556                       vlow (car yp)))
 557             (setq high (math-add-float (car xp) step)
 558                   var-DUMMY high
 559                   vhigh (math-evaluate-expr expr))
 560             (math-working "search" high)
 561             (if (and (Math-realp vhigh)
 562                      (eq (math-negp vhigh) pos))
 563                 (setq found t)
 564               (setcdr xp (cons high (cdr xp)))
 565               (setcdr yp (cons vhigh (cdr yp)))
 566               (setq xp (cdr (cdr xp))
 567                     yp (cdr (cdr yp))))))))
 568     (if found
 569         (if (Math-zerop vhigh)
 570             (list 'vec high vhigh)
 571           (if (Math-zerop vlow)
 572               (list 'vec low vlow)
 573             (if deriv
 574                 (math-newton-search-root expr deriv nil nil nil nil
 575                                          low vlow high vhigh)
 576               (math-bisect-root expr low vlow high vhigh))))
 577       (math-reject-arg (list 'intv 3 low high)
 578                        "*Unable to find a sign change in this interval"))))
 579
 580 ;;; "rtbis"  (but we should be using Brent's method)
 581 (defun math-bisect-root (expr low vlow high vhigh)
 582   (let ((step (math-sub-float high low))
 583         (pos (Math-posp vhigh))
 584         var-DUMMY
 585         mid vmid)
 586     (while (not (or (math-nearly-equal low
 587                                        (setq step (math-mul-float
 588                                                    step '(float 5 -1))
 589                                              mid (math-add-float low step)))
 590                     (progn
 591                       (setq var-DUMMY mid
 592                             vmid (math-evaluate-expr expr))
 593                       (Math-zerop vmid))))
 594       (math-working "bisect" mid)
 595       (if (eq (Math-posp vmid) pos)
 596           (setq high mid
 597                 vhigh vmid)
 598         (setq low mid
 599               vlow vmid)))
 600     (list 'vec mid vmid)))
 601
 602 ;;; "mnewt"
 603 (defun math-newton-multi (expr jacob n guess orig-guess limit)
 604   (let ((m -1)
 605         (p guess)
 606         p2 expr-val jacob-val next)
 607     (while (< (setq p (cdr p) m (1+ m)) n)
 608       (set (nth 2 (aref math-root-vars m)) (car p)))
 609     (setq expr-val (math-evaluate-expr expr)
 610           jacob-val (math-evaluate-expr jacob))
 611     (unless (and (math-constp expr-val)
 612                  (math-constp jacob-val))
 613       (math-reject-arg guess "*Newton's method encountered a singularity"))
 614     (setq next (math-add guess (math-div (math-float (math-neg expr-val))
 615                                          (math-float jacob-val)))
 616           p guess p2 next)
 617     (math-working "newton" next)
 618     (while (and (setq p (cdr p) p2 (cdr p2))
 619                 (math-nearly-equal (car p) (car p2))))
 620     (if p
 621         (if (Math-lessp (math-abs-approx (math-sub next orig-guess))
 622                         limit)
 623             (math-newton-multi expr jacob n next orig-guess limit)
 624           (math-reject-arg nil "*Newton's method failed to converge"))
 625       (list 'vec next expr-val))))
 626
 627 (defvar math-root-vars [(var DUMMY var-DUMMY)])
 628
 629 (defun math-find-root (expr var guess root-widen)
 630   (if (eq (car-safe expr) 'vec)
 631       (let ((n (1- (length expr)))
 632             (calc-symbolic-mode nil)
 633             (var-DUMMY nil)
 634             (jacob (list 'vec))
 635             p p2 m row)
 636         (unless (eq (car-safe var) 'vec)
 637           (math-reject-arg var 'vectorp))
 638         (unless (= (length var) (1+ n))
 639           (math-dimension-error))
 640         (setq expr (copy-sequence expr))
 641         (while (>= n (length math-root-vars))
 642           (let ((symb (intern (concat "math-root-v"
 643                                       (int-to-string
 644                                        (length math-root-vars))))))
 645             (setq math-root-vars (vconcat math-root-vars
 646                                           (vector (list 'var symb symb))))))
 647         (setq m -1)
 648         (while (< (setq m (1+ m)) n)
 649           (set (nth 2 (aref math-root-vars m)) nil))
 650         (setq m -1 p var)
 651         (while (setq m (1+ m) p (cdr p))
 652           (or (eq (car-safe (car p)) 'var)
 653               (math-reject-arg var "*Expected a variable"))
 654           (setq p2 expr)
 655           (while (setq p2 (cdr p2))
 656             (setcar p2 (math-expr-subst (car p2) (car p)
 657                                         (aref math-root-vars m)))))
 658         (unless (eq (car-safe guess) 'vec)
 659           (math-reject-arg guess 'vectorp))
 660         (unless (= (length guess) (1+ n))
 661           (math-dimension-error))
 662         (setq guess (copy-sequence guess)
 663               p guess)
 664         (while (setq p (cdr p))
 665           (or (Math-numberp (car guess))
 666               (math-reject-arg guess 'numberp))
 667           (setcar p (math-float (car p))))
 668         (setq p expr)
 669         (while (setq p (cdr p))
 670           (if (assq (car-safe (car p)) calc-tweak-eqn-table)
 671               (setcar p (math-sub (nth 1 (car p)) (nth 2 (car p)))))
 672           (setcar p (math-evaluate-expr (car p)))
 673           (setq row (list 'vec)
 674                 m -1)
 675           (while (< (setq m (1+ m)) n)
 676             (nconc row (list (math-evaluate-expr
 677                               (or (calcFunc-deriv (car p)
 678                                                   (aref math-root-vars m)
 679                                                   nil t)
 680                                   (math-reject-arg
 681                                    expr
 682                                    "*Formulas must be differentiable"))))))
 683           (nconc jacob (list row)))
 684         (setq m (math-abs-approx guess))
 685         (math-newton-multi expr jacob n guess guess
 686                            (if (math-zerop m) '(float 1 3) (math-mul m 10))))
 687     (unless (eq (car-safe var) 'var)
 688       (math-reject-arg var "*Expected a variable"))
 689     (unless (math-expr-contains expr var)
 690       (math-reject-arg expr "*Formula does not contain specified variable"))
 691     (if (assq (car expr) calc-tweak-eqn-table)
 692         (setq expr (math-sub (nth 1 expr) (nth 2 expr))))
 693     (math-with-extra-prec 2
 694       (setq expr (math-expr-subst expr var '(var DUMMY var-DUMMY)))
 695       (let* ((calc-symbolic-mode nil)
 696              (var-DUMMY nil)
 697              (expr (math-evaluate-expr expr))
 698              (deriv (calcFunc-deriv expr '(var DUMMY var-DUMMY) nil t))
 699              low high vlow vhigh)
 700         (and deriv (setq deriv (math-evaluate-expr deriv)))
 701         (setq guess (math-float guess))
 702         (if (and (math-numberp guess)
 703                  deriv)
 704             (math-newton-root expr deriv guess guess
 705                               (if (math-zerop guess) '(float 1 6)
 706                                 (math-mul (math-abs-approx guess) 100)))
 707           (if (Math-realp guess)
 708               (setq low guess
 709                     high guess
 710                     var-DUMMY guess
 711                     vlow (math-evaluate-expr expr)
 712                     vhigh vlow
 713                     root-widen 'point)
 714             (if (eq (car guess) 'intv)
 715                 (progn
 716                   (or (math-constp guess) (math-reject-arg guess 'constp))
 717                   (setq low (nth 2 guess)
 718                         high (nth 3 guess))
 719                   (if (memq (nth 1 guess) '(0 1))
 720                       (setq low (calcFunc-incr low 1 high)))
 721                   (if (memq (nth 1 guess) '(0 2))
 722                       (setq high (calcFunc-incr high -1 low)))
 723                   (setq var-DUMMY low
 724                         vlow (math-evaluate-expr expr)
 725                         var-DUMMY high
 726                         vhigh (math-evaluate-expr expr)))
 727               (if (math-complexp guess)
 728                   (math-reject-arg "*Complex root finder must have derivative")
 729                 (math-reject-arg guess 'realp))))
 730           (if (Math-zerop vlow)
 731               (list 'vec low vlow)
 732             (if (Math-zerop vhigh)
 733                 (list 'vec high vhigh)
 734               (if (and deriv (Math-numberp vlow) (Math-numberp vhigh))
 735                   (math-newton-search-root expr deriv nil nil nil nil
 736                                            low vlow high vhigh)
 737                 (if (or (and (Math-posp vlow) (Math-posp vhigh))
 738                         (and (Math-negp vlow) (Math-negp vhigh))
 739                         (not (Math-numberp vlow))
 740                         (not (Math-numberp vhigh)))
 741                     (math-search-root expr deriv low vlow high vhigh)
 742                   (math-bisect-root expr low vlow high vhigh))))))))))
 743
 744 (defun calcFunc-root (expr var guess)
 745   (math-find-root expr var guess nil))
 746
 747 (defun calcFunc-wroot (expr var guess)
 748   (math-find-root expr var guess t))
 749
 750
 751
 752
 753 ;;; The following algorithms come from Numerical Recipes, chapter 10.
 754
 755 (defun math-min-eval (expr a)
 756   (if (Math-vectorp a)
 757       (let ((m -1))
 758         (while (setq m (1+ m) a (cdr a))
 759           (set (nth 2 (aref math-min-vars m)) (car a))))
 760     (setq var-DUMMY a))
 761   (setq a (math-evaluate-expr expr))
 762   (if (Math-ratp a)
 763       (math-float a)
 764     (if (eq (car a) 'float)
 765         a
 766       (math-reject-arg a 'realp))))
 767
 768 (defvar math-min-or-max "minimum")
 769
 770 ;;; A bracket for a minimum is a < b < c where f(b) < f(a) and f(b) < f(c).
 771
 772 ;;; "mnbrak"
 773 (defun math-widen-min (expr a b)
 774   (let ((done nil)
 775         (iters 30)
 776         incr c va vb vc u vu r q ulim bc ba qr)
 777     (or b (setq b (math-mul a '(float 101 -2))))
 778     (setq va (math-min-eval expr a)
 779           vb (math-min-eval expr b))
 780     (if (math-lessp-float va vb)
 781         (setq u a a b b u
 782               vu va va vb vb vu))
 783     (setq c (math-add-float b (math-mul-float '(float 161803 -5)
 784                                               (math-sub-float b a)))
 785           vc (math-min-eval expr c))
 786     (while (and (not done) (math-lessp-float vc vb))
 787       (math-working "widen" (list 'intv 0 a c))
 788       (if (= (setq iters (1- iters)) 0)
 789           (math-reject-arg nil (format "*Unable to find a %s near the interval"
 790                                        math-min-or-max)))
 791       (setq bc (math-sub-float b c)
 792             ba (math-sub-float b a)
 793             r (math-mul-float ba (math-sub-float vb vc))
 794             q (math-mul-float bc (math-sub-float vb va))
 795             qr (math-sub-float q r))
 796       (if (math-lessp-float (math-abs qr) '(float 1 -20))
 797           (setq qr (if (math-negp qr) '(float -1 -20) '(float 1 -20))))
 798       (setq u (math-sub-float
 799                b
 800                (math-div-float (math-sub-float (math-mul-float bc q)
 801                                                (math-mul-float ba r))
 802                                (math-mul-float '(float 2 0) qr)))
 803             ulim (math-add-float b (math-mul-float '(float -1 2) bc))
 804             incr (math-negp bc))
 805       (if (if incr (math-lessp-float b u) (math-lessp-float u b))
 806           (if (if incr (math-lessp-float u c) (math-lessp-float c u))
 807               (if (math-lessp-float (setq vu (math-min-eval expr u)) vc)
 808                   (setq a b  va vb
 809                         b u  vb vu
 810                         done t)
 811                 (if (math-lessp-float vb vu)
 812                     (setq c u  vc vu
 813                           done t)
 814                   (setq u (math-add-float c (math-mul-float '(float -161803 -5)
 815                                                             bc))
 816                         vu (math-min-eval expr u))))
 817             (if (if incr (math-lessp-float u ulim) (math-lessp-float ulim u))
 818                 (if (math-lessp-float (setq vu (math-min-eval expr u)) vc)
 819                     (setq b c  vb vc
 820                           c u  vc vu
 821                           u (math-add-float c (math-mul-float
 822                                                '(float -161803 -5)
 823                                                (math-sub-float b c)))
 824                           vu (math-min-eval expr u)))
 825               (setq u ulim
 826                     vu (math-min-eval expr u))))
 827         (setq u (math-add-float c (math-mul-float '(float -161803 -5)
 828                                                   bc))
 829               vu (math-min-eval expr u)))
 830       (setq a b  va vb
 831             b c  vb vc
 832             c u  vc vu))
 833     (if (math-lessp-float a c)
 834         (list a va b vb c vc)
 835       (list c vc b vb a va))))
 836
 837 (defun math-narrow-min (expr a c intv)
 838   (let ((xvals (list a c))
 839         (yvals (list (math-min-eval expr a)
 840                      (math-min-eval expr c)))
 841         (levels 0)
 842         (step (math-sub-float c a))
 843         (found nil)
 844         xp yp b)
 845     (while (and (<= (setq levels (1+ levels)) 5)
 846                 (not found))
 847       (setq xp xvals
 848             yp yvals
 849             step (math-mul-float step '(float 497 -3)))
 850       (while (and (cdr xp) (not found))
 851         (setq b (math-add-float (car xp) step))
 852         (math-working "search" b)
 853         (setcdr xp (cons b (cdr xp)))
 854         (setcdr yp (cons (math-min-eval expr b) (cdr yp)))
 855         (if (and (math-lessp-float (nth 1 yp) (car yp))
 856                  (math-lessp-float (nth 1 yp) (nth 2 yp)))
 857             (setq found t)
 858           (setq xp (cdr xp)
 859                 yp (cdr yp))
 860           (if (and (cdr (cdr yp))
 861                    (math-lessp-float (nth 1 yp) (car yp))
 862                    (math-lessp-float (nth 1 yp) (nth 2 yp)))
 863               (setq found t)
 864             (setq xp (cdr xp)
 865                   yp (cdr yp))))))
 866     (if found
 867         (list (car xp) (car yp)
 868               (nth 1 xp) (nth 1 yp)
 869               (nth 2 xp) (nth 2 yp))
 870       (or (if (math-lessp-float (car yvals) (nth 1 yvals))
 871               (and (memq (nth 1 intv) '(2 3))
 872                    (let ((min (car yvals)))
 873                      (while (and (setq yvals (cdr yvals))
 874                                  (math-lessp-float min (car yvals))))
 875                      (and (not yvals)
 876                           (list (nth 2 intv) min))))
 877             (and (memq (nth 1 intv) '(1 3))
 878                  (setq yvals (nreverse yvals))
 879                  (let ((min (car yvals)))
 880                    (while (and (setq yvals (cdr yvals))
 881                                (math-lessp-float min (car yvals))))
 882                    (and (not yvals)
 883                         (list (nth 3 intv) min)))))
 884           (math-reject-arg nil (format "*Unable to find a %s in the interval"
 885                                        math-min-or-max))))))
 886
 887 ;;; "brent"
 888 (defun math-brent-min (expr prec a va x vx b vb)
 889   (let ((iters (+ 20 (* 5 prec)))
 890         (w x)
 891         (vw vx)
 892         (v x)
 893         (vv vx)
 894         (tol (list 'float 1 (- -1 prec)))
 895         (zeps (list 'float 1 (- -5 prec)))
 896         (e '(float 0 0))
 897         u vu xm tol1 tol2 etemp p q r xv xw)
 898     (while (progn
 899              (setq xm (math-mul-float '(float 5 -1)
 900                                       (math-add-float a b))
 901                    tol1 (math-add-float
 902                          zeps
 903                          (math-mul-float tol (math-abs x)))
 904                    tol2 (math-mul-float tol1 '(float 2 0)))
 905              (math-lessp-float (math-sub-float tol2
 906                                                (math-mul-float
 907                                                 '(float 5 -1)
 908                                                 (math-sub-float b a)))
 909                                (math-abs (math-sub-float x xm))))
 910       (if (= (setq iters (1- iters)) 0)
 911           (math-reject-arg nil (format "*Unable to converge on a %s"
 912                                        math-min-or-max)))
 913       (math-working "brent" x)
 914       (if (math-lessp-float (math-abs e) tol1)
 915           (setq e (if (math-lessp-float x xm)
 916                       (math-sub-float b x)
 917                     (math-sub-float a x))
 918                 d (math-mul-float '(float 381966 -6) e))
 919         (setq xw (math-sub-float x w)
 920               r (math-mul-float xw (math-sub-float vx vv))
 921               xv (math-sub-float x v)
 922               q (math-mul-float xv (math-sub-float vx vw))
 923               p (math-sub-float (math-mul-float xv q)
 924                                 (math-mul-float xw r))
 925               q (math-mul-float '(float 2 0) (math-sub-float q r)))
 926         (if (math-posp q)
 927             (setq p (math-neg-float p))
 928           (setq q (math-neg-float q)))
 929         (setq etemp e
 930               e d)
 931         (if (and (math-lessp-float (math-abs p)
 932                                    (math-abs (math-mul-float
 933                                               '(float 5 -1)
 934                                               (math-mul-float q etemp))))
 935                  (math-lessp-float (math-mul-float
 936                                     q (math-sub-float a x)) p)
 937                  (math-lessp-float p (math-mul-float
 938                                       q (math-sub-float b x))))
 939             (progn
 940               (setq d (math-div-float p q)
 941                     u (math-add-float x d))
 942               (if (or (math-lessp-float (math-sub-float u a) tol2)
 943                       (math-lessp-float (math-sub-float b u) tol2))
 944                   (setq d (if (math-lessp-float xm x)
 945                               (math-neg-float tol1)
 946                             tol1))))
 947           (setq e (if (math-lessp-float x xm)
 948                       (math-sub-float b x)
 949                     (math-sub-float a x))
 950                 d (math-mul-float '(float 381966 -6) e))))
 951       (setq u (math-add-float x
 952                               (if (math-lessp-float (math-abs d) tol1)
 953                                   (if (math-negp d)
 954                                       (math-neg-float tol1)
 955                                     tol1)
 956                                 d))
 957             vu (math-min-eval expr u))
 958       (if (math-lessp-float vx vu)
 959           (progn
 960             (if (math-lessp-float u x)
 961                 (setq a u)
 962               (setq b u))
 963             (if (or (equal w x)
 964                     (not (math-lessp-float vw vu)))
 965                 (setq v w  vv vw
 966                       w u  vw vu)
 967               (if (or (equal v x)
 968                       (equal v w)
 969                       (not (math-lessp-float vv vu)))
 970                   (setq v u  vv vu))))
 971         (if (math-lessp-float u x)
 972             (setq b x)
 973           (setq a x))
 974         (setq v w  vv vw
 975               w x  vw vx
 976               x u  vx vu)))
 977     (list 'vec x vx)))
 978
 979 ;;; "powell"
 980 (defun math-powell-min (expr n guesses prec)
 981   (let* ((f1dim (math-line-min-func expr n))
 982          (xi (calcFunc-idn 1 n))
 983          (p (cons 'vec (mapcar 'car guesses)))
 984          (pt p)
 985          (ftol (list 'float 1 (- prec)))
 986          (fret (math-min-eval expr p))
 987          fp ptt fptt xit i ibig del diff res)
 988     (while (progn
 989              (setq fp fret
 990                    ibig 0
 991                    del '(float 0 0)
 992                    i 0)
 993              (while (<= (setq i (1+ i)) n)
 994                (setq fptt fret
 995                      res (math-line-min f1dim p
 996                                         (math-mat-col xi i)
 997                                         n prec)
 998                      p (let ((calc-internal-prec prec))
 999                          (math-normalize (car res)))
1000                      fret (nth 2 res)
1001                      diff (math-abs (math-sub-float fptt fret)))
1002                (if (math-lessp-float del diff)
1003                    (setq del diff
1004                          ibig i)))
1005              (math-lessp-float
1006               (math-mul-float ftol
1007                               (math-add-float (math-abs fp)
1008                                               (math-abs fret)))
1009               (math-mul-float '(float 2 0)
1010                               (math-abs (math-sub-float fp
1011                                                         fret)))))
1012       (setq ptt (math-sub (math-mul '(float 2 0) p) pt)
1013             xit (math-sub p pt)
1014             pt p
1015             fptt (math-min-eval expr ptt))
1016       (if (and (math-lessp-float fptt fp)
1017                (math-lessp-float
1018                 (math-mul-float
1019                  (math-mul-float '(float 2 0)
1020                                  (math-add-float
1021                                   (math-sub-float fp
1022                                                   (math-mul-float '(float 2 0)
1023                                                                   fret))
1024                                   fptt))
1025                  (math-sqr-float (math-sub-float
1026                                   (math-sub-float fp fret) del)))
1027                 (math-mul-float del
1028                                 (math-sqr-float (math-sub-float fp fptt)))))
1029           (progn
1030             (setq res (math-line-min f1dim p xit n prec)
1031                   p (car res)
1032                   fret (nth 2 res)
1033                   i 0)
1034             (while (<= (setq i (1+ i)) n)
1035               (setcar (nthcdr ibig (nth i xi))
1036                       (nth i (nth 1 res)))))))
1037     (list 'vec p fret)))
1038
1039 (defun math-line-min-func (expr n)
1040   (let ((m -1))
1041     (while (< (setq m (1+ m)) n)
1042       (set (nth 2 (aref math-min-vars m))
1043            (list '+
1044                  (list '*
1045                        '(var DUMMY var-DUMMY)
1046                        (list 'calcFunc-mrow '(var line-xi line-xi) (1+ m)))
1047                  (list 'calcFunc-mrow '(var line-p line-p) (1+ m)))))
1048     (math-evaluate-expr expr)))
1049
1050 (defun math-line-min (f1dim line-p line-xi n prec)
1051   (let* ((var-DUMMY nil)
1052          (expr (math-evaluate-expr f1dim))
1053          (params (math-widen-min expr '(float 0 0) '(float 1 0)))
1054          (res (apply 'math-brent-min expr prec params))
1055          (xi (math-mul (nth 1 res) line-xi)))
1056     (list (math-add line-p xi) xi (nth 2 res))))
1057
1058
1059 (defvar math-min-vars [(var DUMMY var-DUMMY)])
1060
1061 (defun math-find-minimum (expr var guess min-widen)
1062   (let* ((calc-symbolic-mode nil)
1063          (n 0)
1064          (var-DUMMY nil)
1065          (isvec (math-vectorp var))
1066          g guesses)
1067     (or (math-vectorp var)
1068         (setq var (list 'vec var)))
1069     (or (math-vectorp guess)
1070         (setq guess (list 'vec guess)))
1071     (or (= (length var) (length guess))
1072         (math-dimension-error))
1073     (while (setq var (cdr var) guess (cdr guess))
1074       (or (eq (car-safe (car var)) 'var)
1075           (math-reject-arg (car vg) "*Expected a variable"))
1076       (or (math-expr-contains expr (car var))
1077           (math-reject-arg (car var)
1078                            "*Formula does not contain specified variable"))
1079       (while (>= (1+ n) (length math-min-vars))
1080         (let ((symb (intern (concat "math-min-v"
1081                                     (int-to-string
1082                                      (length math-min-vars))))))
1083           (setq math-min-vars (vconcat math-min-vars
1084                                        (vector (list 'var symb symb))))))
1085       (set (nth 2 (aref math-min-vars n)) nil)
1086       (set (nth 2 (aref math-min-vars (1+ n))) nil)
1087       (if (math-complexp (car guess))
1088           (setq expr (math-expr-subst expr
1089                                       (car var)
1090                                       (list '+ (aref math-min-vars n)
1091                                             (list '*
1092                                                   (aref math-min-vars (1+ n))
1093                                                   '(cplx 0 1))))
1094                 guesses (let ((g (math-float (math-complex (car guess)))))
1095                           (cons (list (nth 2 g) nil nil)
1096                                 (cons (list (nth 1 g) nil nil t)
1097                                       guesses)))
1098                 n (+ n 2))
1099         (setq expr (math-expr-subst expr
1100                                     (car var)
1101                                     (aref math-min-vars n))
1102               guesses (cons (if (math-realp (car guess))
1103                                 (list (math-float (car guess)) nil nil)
1104                               (if (and (eq (car-safe (car guess)) 'intv)
1105                                        (math-constp (car guess)))
1106                                   (list (math-mul
1107                                          (math-add (nth 2 (car guess))
1108                                                    (nth 3 (car guess)))
1109                                          '(float 5 -1))
1110                                         (math-float (nth 2 (car guess)))
1111                                         (math-float (nth 3 (car guess)))
1112                                         (car guess))
1113                                 (math-reject-arg (car guess) 'realp)))
1114                             guesses)
1115               n (1+ n))))
1116     (setq guesses (nreverse guesses)
1117           expr (math-evaluate-expr expr))
1118     (if (= n 1)
1119         (let* ((params (if (nth 1 (car guesses))
1120                            (if min-widen
1121                                (math-widen-min expr
1122                                                (nth 1 (car guesses))
1123                                                (nth 2 (car guesses)))
1124                              (math-narrow-min expr
1125                                               (nth 1 (car guesses))
1126                                               (nth 2 (car guesses))
1127                                               (nth 3 (car guesses))))
1128                          (math-widen-min expr
1129                                          (car (car guesses))
1130                                          nil)))
1131                (prec calc-internal-prec)
1132                (res (if (cdr (cdr params))
1133                         (math-with-extra-prec (+ calc-internal-prec 2)
1134                           (apply 'math-brent-min expr prec params))
1135                       (cons 'vec params))))
1136           (if isvec
1137               (list 'vec (list 'vec (nth 1 res)) (nth 2 res))
1138             res))
1139       (let* ((prec calc-internal-prec)
1140              (res (math-with-extra-prec (+ calc-internal-prec 2)
1141                     (math-powell-min expr n guesses prec)))
1142              (p (nth 1 res))
1143              (vec (list 'vec)))
1144         (while (setq p (cdr p))
1145           (if (nth 3 (car guesses))
1146               (progn
1147                 (nconc vec (list (math-normalize
1148                                   (list 'cplx (car p) (nth 1 p)))))
1149                 (setq p (cdr p)
1150                       guesses (cdr guesses)))
1151             (nconc vec (list (car p))))
1152           (setq guesses (cdr guesses)))
1153         (if isvec
1154             (list 'vec vec (nth 2 res))
1155           (list 'vec (nth 1 vec) (nth 2 res)))))))
1156
1157 (defun calcFunc-minimize (expr var guess)
1158   (let ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
1159         (math-min-or-max "minimum"))
1160     (math-find-minimum (math-normalize expr)
1161                        (math-normalize var)
1162                        (math-normalize guess) nil)))
1163
1164 (defun calcFunc-wminimize (expr var guess)
1165   (let ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
1166         (math-min-or-max "minimum"))
1167     (math-find-minimum (math-normalize expr)
1168                        (math-normalize var)
1169                        (math-normalize guess) t)))
1170
1171 (defun calcFunc-maximize (expr var guess)
1172   (let* ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
1173          (math-min-or-max "maximum")
1174          (res (math-find-minimum (math-normalize (math-neg expr))
1175                                  (math-normalize var)
1176                                  (math-normalize guess) nil)))
1177     (list 'vec (nth 1 res) (math-neg (nth 2 res)))))
1178
1179 (defun calcFunc-wmaximize (expr var guess)
1180   (let* ((calc-internal-prec (max (/ calc-internal-prec 2) 3))
1181          (math-min-or-max "maximum")
1182          (res (math-find-minimum (math-normalize (math-neg expr))
1183                                  (math-normalize var)
1184                                  (math-normalize guess) t)))
1185     (list 'vec (nth 1 res) (math-neg (nth 2 res)))))
1186
1187
1188
1189
1190 ;;; The following algorithms come from Numerical Recipes, chapter 3.
1191
1192 (defun calcFunc-polint (data x)
1193   (or (math-matrixp data) (math-reject-arg data 'matrixp))
1194   (or (= (length data) 3)
1195       (math-reject-arg data "*Wrong number of data rows"))
1196   (or (> (length (nth 1 data)) 2)
1197       (math-reject-arg data "*Too few data points"))
1198   (if (and (math-vectorp x) (or (math-constp x) math-expand-formulas))
1199       (cons 'vec (mapcar (function (lambda (x) (calcFunc-polint data x)))
1200                          (cdr x)))
1201     (or (math-objectp x) math-expand-formulas (math-reject-arg x 'objectp))
1202     (math-with-extra-prec 2
1203       (cons 'vec (math-poly-interp (cdr (nth 1 data)) (cdr (nth 2 data)) x
1204                                    nil)))))
1205 (put 'calcFunc-polint 'math-expandable t)
1206
1207
1208 (defun calcFunc-ratint (data x)
1209   (or (math-matrixp data) (math-reject-arg data 'matrixp))
1210   (or (= (length data) 3)
1211       (math-reject-arg data "*Wrong number of data rows"))
1212   (or (> (length (nth 1 data)) 2)
1213       (math-reject-arg data "*Too few data points"))
1214   (if (and (math-vectorp x) (or (math-constp x) math-expand-formulas))
1215       (cons 'vec (mapcar (function (lambda (x) (calcFunc-ratint data x)))
1216                          (cdr x)))
1217     (or (math-objectp x) math-expand-formulas (math-reject-arg x 'objectp))
1218     (math-with-extra-prec 2
1219       (cons 'vec (math-poly-interp (cdr (nth 1 data)) (cdr (nth 2 data)) x
1220                                    (cdr (cdr (cdr (nth 1 data)))))))))
1221 (put 'calcFunc-ratint 'math-expandable t)
1222
1223
1224 (defun math-poly-interp (xa ya x ratp)
1225   (let ((n (length xa))
1226         (dif nil)
1227         (ns nil)
1228         (xax nil)
1229         (c (copy-sequence ya))
1230         (d (copy-sequence ya))
1231         (i 0)
1232         (m 0)
1233         y dy (xp xa) xpm cp dp temp)
1234     (while (<= (setq i (1+ i)) n)
1235       (setq xax (cons (math-sub (car xp) x) xax)
1236             xp (cdr xp)
1237             temp (math-abs (car xax)))
1238       (if (or (null dif) (math-lessp temp dif))
1239           (setq dif temp
1240                 ns i)))
1241     (setq xax (nreverse xax)
1242           ns (1- ns)
1243           y (nth ns ya))
1244     (if (math-zerop dif)
1245         (list y 0)
1246       (while (< (setq m (1+ m)) n)
1247         (setq i 0
1248               xp xax
1249               xpm (nthcdr m xax)
1250               cp c
1251               dp d)
1252         (while (<= (setq i (1+ i)) (- n m))
1253           (if ratp
1254               (let ((t2 (math-div (math-mul (car xp) (car dp)) (car xpm))))
1255                 (setq temp (math-div (math-sub (nth 1 cp) (car dp))
1256                                      (math-sub t2 (nth 1 cp))))
1257                 (setcar dp (math-mul (nth 1 cp) temp))
1258                 (setcar cp (math-mul t2 temp)))
1259             (if (math-equal (car xp) (car xpm))
1260                 (math-reject-arg (cons 'vec xa) "*Duplicate X values"))
1261             (setq temp (math-div (math-sub (nth 1 cp) (car dp))
1262                                  (math-sub (car xp) (car xpm))))
1263             (setcar dp (math-mul (car xpm) temp))
1264             (setcar cp (math-mul (car xp) temp)))
1265           (setq cp (cdr cp)
1266                 dp (cdr dp)
1267                 xp (cdr xp)
1268                 xpm (cdr xpm)))
1269         (if (< (+ ns ns) (- n m))
1270             (setq dy (nth ns c))
1271           (setq ns (1- ns)
1272                 dy (nth ns d)))
1273         (setq y (math-add y dy)))
1274       (list y dy))))
1275
1276
1277
1278 ;;; The following algorithms come from Numerical Recipes, chapter 4.
1279
1280 (defun calcFunc-ninteg (expr var lo hi)
1281   (setq lo (math-evaluate-expr lo)
1282         hi (math-evaluate-expr hi))
1283   (or (math-numberp lo) (math-infinitep lo) (math-reject-arg lo 'numberp))
1284   (or (math-numberp hi) (math-infinitep hi) (math-reject-arg hi 'numberp))
1285   (if (math-lessp hi lo)
1286       (math-neg (calcFunc-ninteg expr var hi lo))
1287     (setq expr (math-expr-subst expr var '(var DUMMY var-DUMMY)))
1288     (let ((var-DUMMY nil)
1289           (calc-symbolic-mode nil)
1290           (calc-prefer-frac nil)
1291           (sum 0))
1292       (setq expr (math-evaluate-expr expr))
1293       (if (equal lo '(neg (var inf var-inf)))
1294           (let ((thi (if (math-lessp hi '(float -2 0))
1295                          hi '(float -2 0))))
1296             (setq sum (math-ninteg-romberg
1297                        'math-ninteg-midpoint expr
1298                          (math-float lo) (math-float thi) 'inf)
1299                   lo thi)))
1300       (if (equal hi '(var inf var-inf))
1301           (let ((tlo (if (math-lessp '(float 2 0) lo)
1302                          lo '(float 2 0))))
1303             (setq sum (math-add sum
1304                                 (math-ninteg-romberg
1305                                  'math-ninteg-midpoint expr
1306                                  (math-float tlo) (math-float hi) 'inf))
1307                   hi tlo)))
1308       (or (math-equal lo hi)
1309           (setq sum (math-add sum
1310                               (math-ninteg-romberg
1311                                'math-ninteg-midpoint expr
1312                                (math-float lo) (math-float hi) nil))))
1313       sum)))
1314
1315
1316 ;;; Open Romberg method; "qromo" in section 4.4.
1317 (defun math-ninteg-romberg (func expr lo hi mode)
1318   (let ((curh '(float 1 0))
1319         (h nil)
1320         (s nil)
1321         (j 0)
1322         (ss nil)
1323         (prec calc-internal-prec)
1324         (integ-temp nil))
1325     (math-with-extra-prec 2
1326       ;; Limit on "j" loop must be 14 or less to keep "it" from overflowing.
1327       (or (while (and (null ss) (<= (setq j (1+ j)) 8))
1328             (setq s (nconc s (list (funcall func expr lo hi mode)))
1329                   h (nconc h (list curh)))
1330             (if (>= j 3)
1331                 (let ((res (math-poly-interp h s '(float 0 0) nil)))
1332                   (if (math-lessp (math-abs (nth 1 res))
1333                                   (calcFunc-scf (math-abs (car res))
1334                                                 (- prec)))
1335                       (setq math-ninteg-convergence j
1336                             ss (car res)))))
1337             (if (>= j 5)
1338                 (setq s (cdr s)
1339                       h (cdr h)))
1340             (setq curh (math-div-float curh '(float 9 0))))
1341           ss
1342           (math-reject-arg nil (format "*Integral failed to converge"))))))
1343
1344
1345 (defun math-ninteg-evaluate (expr x mode)
1346   (if (eq mode 'inf)
1347       (setq x (math-div '(float 1 0) x)))
1348   (let* ((var-DUMMY x)
1349          (res (math-evaluate-expr expr)))
1350     (or (Math-numberp res)
1351         (math-reject-arg res "*Integrand does not evaluate to a number"))
1352     (if (eq mode 'inf)
1353         (setq res (math-mul res (math-sqr x))))
1354     res))
1355
1356
1357 (defun math-ninteg-midpoint (expr lo hi mode)    ; uses "integ-temp"
1358   (if (eq mode 'inf)
1359       (let ((math-infinite-mode t) temp)
1360         (setq temp (math-div 1 lo)
1361               lo (math-div 1 hi)
1362               hi temp)))
1363   (if integ-temp
1364       (let* ((it3 (* 3 (car integ-temp)))
1365              (math-working-step-2 (* 2 (car integ-temp)))
1366              (math-working-step 0)
1367              (range (math-sub hi lo))
1368              (del (math-div range (math-float it3)))
1369              (del2 (math-add del del))
1370              (del3 (math-add del del2))
1371              (x (math-add lo (math-mul '(float 5 -1) del)))
1372              (sum '(float 0 0))
1373              (j 0) temp)
1374         (while (<= (setq j (1+ j)) (car integ-temp))
1375           (setq math-working-step (1+ math-working-step)
1376                 temp (math-ninteg-evaluate expr x mode)
1377                 math-working-step (1+ math-working-step)
1378                 sum (math-add sum (math-add temp (math-ninteg-evaluate
1379                                                   expr (math-add x del2)
1380                                                   mode)))
1381                 x (math-add x del3)))
1382         (setq integ-temp (list it3
1383                                (math-add (math-div (nth 1 integ-temp)
1384                                                    '(float 3 0))
1385                                          (math-mul sum del)))))
1386     (setq integ-temp (list 1 (math-mul
1387                               (math-sub hi lo)
1388                               (math-ninteg-evaluate
1389                                expr
1390                                (math-mul (math-add lo hi) '(float 5 -1))
1391                                mode)))))
1392   (nth 1 integ-temp))
1393
1394
1395
1396
1397
1398 ;;; The following algorithms come from Numerical Recipes, chapter 14.
1399
1400 (defvar math-dummy-vars [(var DUMMY var-DUMMY)])
1401 (defvar math-dummy-counter 0)
1402 (defun math-dummy-variable ()
1403   (if (= math-dummy-counter (length math-dummy-vars))
1404       (let ((symb (intern (format "math-dummy-%d" math-dummy-counter))))
1405         (setq math-dummy-vars (vconcat math-dummy-vars
1406                                        (vector (list 'var symb symb))))))
1407   (set (nth 2 (aref math-dummy-vars math-dummy-counter)) nil)
1408   (prog1
1409       (aref math-dummy-vars math-dummy-counter)
1410     (setq math-dummy-counter (1+ math-dummy-counter))))
1411
1412 (defvar math-in-fit 0)
1413 (defvar calc-fit-to-trail nil)
1414
1415 (defun calcFunc-fit (expr vars &optional coefs data)
1416   (let ((math-in-fit 10))
1417     (math-with-extra-prec 2
1418       (math-general-fit expr vars coefs data nil))))
1419
1420 (defun calcFunc-efit (expr vars &optional coefs data)
1421   (let ((math-in-fit 10))
1422     (math-with-extra-prec 2
1423       (math-general-fit expr vars coefs data 'sdev))))
1424
1425 (defun calcFunc-xfit (expr vars &optional coefs data)
1426   (let ((math-in-fit 10))
1427     (math-with-extra-prec 2
1428       (math-general-fit expr vars coefs data 'full))))
1429
1430 (defun math-general-fit (expr vars coefs data mode)
1431   (let ((calc-simplify-mode nil)
1432         (math-dummy-counter math-dummy-counter)
1433         (math-in-fit 1)
1434         (extended (eq mode 'full))
1435         (first-coef math-dummy-counter)
1436         first-var
1437         (plain-expr expr)
1438         orig-expr
1439         have-sdevs need-chisq chisq
1440         (x-funcs nil)
1441         (y-filter nil)
1442         y-dummy
1443         (coef-filters nil)
1444         new-coefs
1445         (xy-values nil)
1446         (weights nil)
1447         (var-YVAL nil) (var-YVALX nil)
1448         covar beta
1449         n nn m mm v dummy p)
1450
1451     ;; Validate and parse arguments.
1452     (or data
1453         (if coefs
1454             (setq data coefs
1455                   coefs nil)
1456           (if (math-vectorp expr)
1457               (if (memq (length expr) '(3 4))
1458                   (setq data vars
1459                         vars (nth 2 expr)
1460                         coefs (nth 3 expr)
1461                         expr (nth 1 expr))
1462                 (math-dimension-error))
1463             (setq data vars
1464                   vars nil
1465                   coefs nil))))
1466     (or (math-matrixp data) (math-reject-arg data 'matrixp))
1467     (setq v (1- (length data))
1468           n (1- (length (nth 1 data))))
1469     (or (math-vectorp vars) (null vars)
1470         (setq vars (list 'vec vars)))
1471     (or (math-vectorp coefs) (null coefs)
1472         (setq coefs (list 'vec coefs)))
1473     (or coefs
1474         (setq coefs (cons 'vec (math-all-vars-but expr vars))))
1475     (or vars
1476         (if (<= (1- (length coefs)) v)
1477             (math-reject-arg coefs "*Not enough variables in model")
1478           (setq coefs (copy-sequence coefs))
1479           (let ((p (nthcdr (- (length coefs) v
1480                               (if (eq (car-safe expr) 'calcFunc-eq) 1 0))
1481                            coefs)))
1482             (setq vars (cons 'vec (cdr p)))
1483             (setcdr p nil))))
1484     (or (= (1- (length vars)) v)
1485         (= (length vars) v)
1486         (math-reject-arg vars "*Number of variables does not match data"))
1487     (setq m (1- (length coefs)))
1488     (if (< m 1)
1489         (math-reject-arg coefs "*Need at least one parameter"))
1490
1491     ;; Rewrite expr in terms of fitparam and fitvar, make into an equation.
1492     (setq p coefs)
1493     (while (setq p (cdr p))
1494       (or (eq (car-safe (car p)) 'var)
1495           (math-reject-arg (car p) "*Expected a variable"))
1496       (setq dummy (math-dummy-variable)
1497             expr (math-expr-subst expr (car p)
1498                                   (list 'calcFunc-fitparam
1499                                         (- math-dummy-counter first-coef)))))
1500     (setq first-var math-dummy-counter
1501           p vars)
1502     (while (setq p (cdr p))
1503       (or (eq (car-safe (car p)) 'var)
1504           (math-reject-arg (car p) "*Expected a variable"))
1505       (setq dummy (math-dummy-variable)
1506             expr (math-expr-subst expr (car p)
1507                                   (list 'calcFunc-fitvar
1508                                         (- math-dummy-counter first-var)))))
1509     (if (< math-dummy-counter (+ first-var v))
1510         (setq dummy (math-dummy-variable))) ; dependent variable may be unnamed
1511     (setq y-dummy dummy
1512           orig-expr expr)
1513     (or (eq (car-safe expr) 'calcFunc-eq)
1514         (setq expr (list 'calcFunc-eq (list 'calcFunc-fitvar v) expr)))
1515
1516     (let ((calc-symbolic-mode nil))
1517
1518       ;; Apply rewrites to put expr into a linear-like form.
1519       (setq expr (math-evaluate-expr expr)
1520             expr (math-rewrite (list 'calcFunc-fitmodel expr)
1521                                '(var FitRules var-FitRules))
1522             math-in-fit 2
1523             expr (math-evaluate-expr expr))
1524       (or (and (eq (car-safe expr) 'calcFunc-fitsystem)
1525                (= (length expr) 4)
1526                (math-vectorp (nth 2 expr))
1527                (math-vectorp (nth 3 expr))
1528                (> (length (nth 2 expr)) 1)
1529                (= (length (nth 3 expr)) (1+ m)))
1530           (math-reject-arg plain-expr "*Model expression is too complex"))
1531       (setq y-filter (nth 1 expr)
1532             x-funcs (vconcat (cdr (nth 2 expr)))
1533             coef-filters (nth 3 expr)
1534             mm (length x-funcs))
1535       (if (equal y-filter y-dummy)
1536           (setq y-filter nil))
1537
1538       ;; Build the (square) system of linear equations to be solved.
1539       (setq beta (cons 'vec (make-list mm 0))
1540             covar (cons 'vec (mapcar 'copy-sequence (make-list mm beta))))
1541       (let* ((ptrs (vconcat (cdr data)))
1542              (isigsq 1)
1543              (xvals (make-vector mm 0))
1544              (i 0)
1545              j k xval yval sigmasqr wt covj covjk covk betaj lud)
1546         (while (<= (setq i (1+ i)) n)
1547
1548           ;; Assign various independent variables for this data point.
1549           (setq j 0
1550                 sigmasqr nil)
1551           (while (< j v)
1552             (aset ptrs j (cdr (aref ptrs j)))
1553             (setq xval (car (aref ptrs j)))
1554             (if (= j (1- v))
1555                 (if sigmasqr
1556                     (progn
1557                       (if (eq (car-safe xval) 'sdev)
1558                           (setq sigmasqr (math-add (math-sqr (nth 2 xval))
1559                                                    sigmasqr)
1560                                 xval (nth 1 xval)))
1561                       (if y-filter
1562                           (setq xval (math-make-sdev xval
1563                                                      (math-sqrt sigmasqr))))))
1564               (if (eq (car-safe xval) 'sdev)
1565                   (setq sigmasqr (math-add (math-sqr (nth 2 xval))
1566                                            (or sigmasqr 0))
1567                         xval (nth 1 xval))))
1568             (set (nth 2 (aref math-dummy-vars (+ first-var j))) xval)
1569             (setq j (1+ j)))
1570
1571           ;; Compute Y value for this data point.
1572           (if y-filter
1573               (setq yval (math-evaluate-expr y-filter))
1574             (setq yval (symbol-value (nth 2 y-dummy))))
1575           (if (eq (car-safe yval) 'sdev)
1576               (setq sigmasqr (math-sqr (nth 2 yval))
1577                     yval (nth 1 yval)))
1578           (if (= i 1)
1579               (setq have-sdevs sigmasqr
1580                     need-chisq (or extended
1581                                    (and (eq mode 'sdev) (not have-sdevs)))))
1582           (if have-sdevs
1583               (if sigmasqr
1584                   (progn
1585                     (setq isigsq (math-div 1 sigmasqr))
1586                     (if need-chisq
1587                         (setq weights (cons isigsq weights))))
1588                 (math-reject-arg yval "*Mixed error forms and plain numbers"))
1589             (if sigmasqr
1590                 (math-reject-arg yval "*Mixed error forms and plain numbers")))
1591
1592           ;; Compute X values for this data point and update covar and beta.
1593           (if (eq (car-safe xval) 'sdev)
1594               (set (nth 2 y-dummy) (nth 1 xval)))
1595           (setq j 0
1596                 covj covar
1597                 betaj beta)
1598           (while (< j mm)
1599             (setq wt (math-evaluate-expr (aref x-funcs j)))
1600             (aset xvals j wt)
1601             (setq wt (math-mul wt isigsq)
1602                   betaj (cdr betaj)
1603                   covjk (car (setq covj (cdr covj)))
1604                   k 0)
1605             (while (<= k j)
1606               (setq covjk (cdr covjk))
1607               (setcar covjk (math-add (car covjk)
1608                                       (math-mul wt (aref xvals k))))
1609               (setq k (1+ k)))
1610             (setcar betaj (math-add (car betaj) (math-mul wt yval)))
1611             (setq j (1+ j)))
1612           (if need-chisq
1613               (setq xy-values (cons (append xvals (list yval)) xy-values))))
1614
1615         ;; Fill in symmetric half of covar matrix.
1616         (setq j 0
1617               covj covar)
1618         (while (< j (1- mm))
1619           (setq k j
1620                 j (1+ j)
1621                 covjk (nthcdr j (car (setq covj (cdr covj))))
1622                 covk (nthcdr j covar))
1623           (while (< (setq k (1+ k)) mm)
1624             (setq covjk (cdr covjk)
1625                   covk (cdr covk))
1626             (setcar covjk (nth j (car covk))))))
1627
1628       ;; Solve the linear system.
1629       (if mode
1630           (progn
1631             (setq covar (math-matrix-inv-raw covar))
1632             (if covar
1633                 (setq beta (math-mul covar beta))
1634               (if (math-zerop (math-abs beta))
1635                   (setq covar (calcFunc-diag 0 (1- (length beta))))
1636                 (math-reject-arg orig-expr "*Singular matrix")))
1637             (or (math-vectorp covar)
1638                 (setq covar (list 'vec (list 'vec covar)))))
1639         (setq beta (math-div beta covar)))
1640
1641       ;; Compute chi-square statistic if necessary.
1642       (if need-chisq
1643           (let (bp xp sum)
1644             (setq chisq 0)
1645             (while xy-values
1646               (setq bp beta
1647                     xp (car xy-values)
1648                     sum 0)
1649               (while (setq bp (cdr bp))
1650                 (setq sum (math-add sum (math-mul (car bp) (car xp)))
1651                       xp (cdr xp)))
1652               (setq sum (math-sqr (math-sub (car xp) sum)))
1653               (if weights (setq sum (math-mul sum (car weights))))
1654               (setq chisq (math-add chisq sum)
1655                     weights (cdr weights)
1656                     xy-values (cdr xy-values)))))
1657
1658       ;; Convert coefficients back into original terms.
1659       (setq new-coefs (copy-sequence beta))
1660       (let* ((bp new-coefs)
1661              (cp covar)
1662              (sigdat 1)
1663              (math-in-fit 3)
1664              (j 0))
1665         (and mode (not have-sdevs)
1666              (setq sigdat (if (<= n mm)
1667                               0
1668                             (math-div chisq (- n mm)))))
1669         (if mode
1670             (while (setq bp (cdr bp))
1671               (setcar bp (math-make-sdev
1672                           (car bp)
1673                           (math-sqrt (math-mul (nth (setq j (1+ j))
1674                                                     (car (setq cp (cdr cp))))
1675                                                sigdat))))))
1676         (setq new-coefs (math-evaluate-expr coef-filters))
1677         (if calc-fit-to-trail
1678             (let ((bp new-coefs)
1679                   (cp coefs)
1680                   (vec nil))
1681               (while (setq bp (cdr bp) cp (cdr cp))
1682                 (setq vec (cons (list 'calcFunc-eq (car cp) (car bp)) vec)))
1683               (setq calc-fit-to-trail (cons 'vec (nreverse vec)))))))
1684
1685     ;; Substitute best-fit coefficients back into original formula.
1686     (setq expr (math-multi-subst
1687                 orig-expr
1688                 (let ((n v)
1689                       (vec nil))
1690                   (while (>= n 1)
1691                     (setq vec (cons (list 'calcFunc-fitvar n) vec)
1692                           n (1- n)))
1693                   (setq n m)
1694                   (while (>= n 1)
1695                     (setq vec (cons (list 'calcFunc-fitparam n) vec)
1696                           n (1- n)))
1697                   vec)
1698                 (append (cdr new-coefs) (cdr vars))))
1699
1700     ;; Package the result.
1701     (math-normalize
1702      (if extended
1703          (list 'vec expr beta covar
1704                (let ((p coef-filters)
1705                      (n 0))
1706                  (while (and (setq n (1+ n) p (cdr p))
1707                              (eq (car-safe (car p)) 'calcFunc-fitdummy)
1708                              (eq (nth 1 (car p)) n)))
1709                  (if p
1710                      coef-filters
1711                    (list 'vec)))
1712                chisq
1713                (if (and have-sdevs (> n mm))
1714                    (list 'calcFunc-utpc chisq (- n mm))
1715                  '(var nan var-nan)))
1716        expr))))
1717
1718
1719 (defun calcFunc-fitvar (x)
1720   (if (>= math-in-fit 2)
1721       (progn
1722         (setq x (aref math-dummy-vars (+ first-var x -1)))
1723         (or (calc-var-value (nth 2 x)) x))
1724     (math-reject-arg x)))
1725
1726 (defun calcFunc-fitparam (x)
1727   (if (>= math-in-fit 2)
1728       (progn
1729         (setq x (aref math-dummy-vars (+ first-coef x -1)))
1730         (or (calc-var-value (nth 2 x)) x))
1731     (math-reject-arg x)))
1732
1733 (defun calcFunc-fitdummy (x)
1734   (if (= math-in-fit 3)
1735       (nth x new-coefs)
1736     (math-reject-arg x)))
1737
1738 (defun calcFunc-hasfitvars (expr)
1739   (if (Math-primp expr)
1740       0
1741     (if (eq (car expr) 'calcFunc-fitvar)
1742         (nth 1 expr)
1743       (apply 'max (mapcar 'calcFunc-hasfitvars (cdr expr))))))
1744
1745 (defun calcFunc-hasfitparams (expr)
1746   (if (Math-primp expr)
1747       0
1748     (if (eq (car expr) 'calcFunc-fitparam)
1749         (nth 1 expr)
1750       (apply 'max (mapcar 'calcFunc-hasfitparams (cdr expr))))))
1751
1752
1753 (defun math-all-vars-but (expr but)
1754   (let* ((vars (math-all-vars-in expr))
1755          (p but))
1756     (while p
1757       (setq vars (delq (assoc (car-safe p) vars) vars)
1758             p (cdr p)))
1759     (sort (mapcar 'car vars)
1760           (function (lambda (x y) (string< (nth 1 x) (nth 1 y)))))))
1761
1762 (defun math-all-vars-in (expr)
1763   (let ((vars nil)
1764         found)
1765     (math-all-vars-rec expr)
1766     vars))
1767
1768 (defun math-all-vars-rec (expr)
1769   (if (Math-primp expr)
1770       (if (eq (car-safe expr) 'var)
1771           (or (math-const-var expr)
1772               (if (setq found (assoc expr vars))
1773                   (setcdr found (1+ (cdr found)))
1774                 (setq vars (cons (cons expr 1) vars)))))
1775     (while (setq expr (cdr expr))
1776       (math-all-vars-rec (car expr)))))
1777
1778 ;;; arch-tag: ff9f2920-8111-48b5-b3fa-b0682c3e44a6
1779 ;;; calcalg3.el ends here