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