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