cgit.sxemacs.org Git - packages/blob - xemacs-packages/calc/calc-alg-3.el

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