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