lisp/calc/calcalg3.el

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