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

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