code.delx.au - gnu-emacs/blob - lisp/calc/calcalg3.el

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