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