lisp/calc/calcalg3.el

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