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 1835