aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lisp/calc/calcalg3.el339
1 files changed, 195 insertions, 144 deletions
diff --git a/lisp/calc/calcalg3.el b/lisp/calc/calcalg3.el
index 56ee7f69878..fd204e23b2b 100644
--- a/lisp/calc/calcalg3.el
+++ b/lisp/calc/calcalg3.el
@@ -3,8 +3,7 @@
3;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc. 3;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.
4 4
5;; Author: David Gillespie <daveg@synaptics.com> 5;; Author: David Gillespie <daveg@synaptics.com>
6;; Maintainers: D. Goel <deego@gnufans.org> 6;; Maintainer: Jay Belanger <belanger@truman.edu>
7;; Colin Walters <walters@debian.org>
8 7
9;; This file is part of GNU Emacs. 8;; This file is part of GNU Emacs.
10 9
@@ -99,8 +98,15 @@
99 (calc-enter-result 1 "poli" (list 'calcFunc-polint data 98 (calc-enter-result 1 "poli" (list 'calcFunc-polint data
100 (calc-top 1))))))) 99 (calc-top 1)))))))
101 100
101;; The variables calc-curve-nvars, calc-curve-varnames, calc-curve-model and calc-curve-coefnames are local to calc-curve-fit, but are
102;; used by calc-get-fit-variables which is called by calc-curve-fit.
103(defvar calc-curve-nvars)
104(defvar calc-curve-varnames)
105(defvar calc-curve-model)
106(defvar calc-curve-coefnames)
102 107
103(defun calc-curve-fit (arg &optional model coefnames varnames) 108(defun calc-curve-fit (arg &optional calc-curve-model
109 calc-curve-coefnames calc-curve-varnames)
104 (interactive "P") 110 (interactive "P")
105 (calc-slow-wrapper 111 (calc-slow-wrapper
106 (setq calc-aborted-prefix nil) 112 (setq calc-aborted-prefix nil)
@@ -108,7 +114,7 @@
108 (if (calc-is-hyperbolic) 'calcFunc-efit 114 (if (calc-is-hyperbolic) 'calcFunc-efit
109 'calcFunc-fit))) 115 'calcFunc-fit)))
110 key (which 0) 116 key (which 0)
111 n nvars temp data 117 n calc-curve-nvars temp data
112 (homog nil) 118 (homog nil)
113 (msgs '( "(Press ? for help)" 119 (msgs '( "(Press ? for help)"
114 "1 = linear or multilinear" 120 "1 = linear or multilinear"
@@ -120,7 +126,7 @@
120 "g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)" 126 "g = (a/b sqrt(2 pi)) exp(-0.5*((x-c)/b)^2)"
121 "h prefix = homogeneous model (no constant term)" 127 "h prefix = homogeneous model (no constant term)"
122 "' = alg entry, $ = stack, u = Model1, U = Model2"))) 128 "' = alg entry, $ = stack, u = Model1, U = Model2")))
123 (while (not model) 129 (while (not calc-curve-model)
124 (message "Fit to model: %s:%s" 130 (message "Fit to model: %s:%s"
125 (nth which msgs) 131 (nth which msgs)
126 (if homog " h" "")) 132 (if homog " h" ""))
@@ -150,44 +156,50 @@
150 (t (error "Bad prefix argument"))) 156 (t (error "Bad prefix argument")))
151 (or (math-matrixp data) (not (cdr (cdr data))) 157 (or (math-matrixp data) (not (cdr (cdr data)))
152 (error "Data matrix is not a matrix!")) 158 (error "Data matrix is not a matrix!"))
153 (setq nvars (- (length data) 2) 159 (setq calc-curve-nvars (- (length data) 2)
154 coefnames nil 160 calc-curve-coefnames nil
155 varnames nil) 161 calc-curve-varnames nil)
156 nil)) 162 nil))
157 ((= key ?1) ; linear or multilinear 163 ((= key ?1) ; linear or multilinear
158 (calc-get-fit-variables nvars (1+ nvars) (and homog 0)) 164 (calc-get-fit-variables calc-curve-nvars
159 (setq model (math-mul coefnames 165 (1+ calc-curve-nvars) (and homog 0))
160 (cons 'vec (cons 1 (cdr varnames)))))) 166 (setq calc-curve-model (math-mul calc-curve-coefnames
167 (cons 'vec (cons 1 (cdr calc-curve-varnames))))))
161 ((and (>= key ?2) (<= key ?9)) ; polynomial 168 ((and (>= key ?2) (<= key ?9)) ; polynomial
162 (calc-get-fit-variables 1 (- key ?0 -1) (and homog 0)) 169 (calc-get-fit-variables 1 (- key ?0 -1) (and homog 0))
163 (setq model (math-build-polynomial-expr (cdr coefnames) 170 (setq calc-curve-model
164 (nth 1 varnames)))) 171 (math-build-polynomial-expr (cdr calc-curve-coefnames)
172 (nth 1 calc-curve-varnames))))
165 ((= key ?i) ; exact polynomial 173 ((= key ?i) ; exact polynomial
166 (calc-get-fit-variables 1 (1- (length (nth 1 data))) 174 (calc-get-fit-variables 1 (1- (length (nth 1 data)))
167 (and homog 0)) 175 (and homog 0))
168 (setq model (math-build-polynomial-expr (cdr coefnames) 176 (setq calc-curve-model
169 (nth 1 varnames)))) 177 (math-build-polynomial-expr (cdr calc-curve-coefnames)
178 (nth 1 calc-curve-varnames))))
170 ((= key ?p) ; power law 179 ((= key ?p) ; power law
171 (calc-get-fit-variables nvars (1+ nvars) (and homog 1)) 180 (calc-get-fit-variables calc-curve-nvars
172 (setq model (math-mul (nth 1 coefnames) 181 (1+ calc-curve-nvars) (and homog 1))
182 (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
173 (calcFunc-reduce 183 (calcFunc-reduce
174 '(var mul var-mul) 184 '(var mul var-mul)
175 (calcFunc-map 185 (calcFunc-map
176 '(var pow var-pow) 186 '(var pow var-pow)
177 varnames 187 calc-curve-varnames
178 (cons 'vec (cdr (cdr coefnames)))))))) 188 (cons 'vec (cdr (cdr calc-curve-coefnames))))))))
179 ((= key ?^) ; exponential law 189 ((= key ?^) ; exponential law
180 (calc-get-fit-variables nvars (1+ nvars) (and homog 1)) 190 (calc-get-fit-variables calc-curve-nvars
181 (setq model (math-mul (nth 1 coefnames) 191 (1+ calc-curve-nvars) (and homog 1))
192 (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
182 (calcFunc-reduce 193 (calcFunc-reduce
183 '(var mul var-mul) 194 '(var mul var-mul)
184 (calcFunc-map 195 (calcFunc-map
185 '(var pow var-pow) 196 '(var pow var-pow)
186 (cons 'vec (cdr (cdr coefnames))) 197 (cons 'vec (cdr (cdr calc-curve-coefnames)))
187 varnames))))) 198 calc-curve-varnames)))))
188 ((memq key '(?e ?E)) 199 ((memq key '(?e ?E))
189 (calc-get-fit-variables nvars (1+ nvars) (and homog 1)) 200 (calc-get-fit-variables calc-curve-nvars
190 (setq model (math-mul (nth 1 coefnames) 201 (1+ calc-curve-nvars) (and homog 1))
202 (setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
191 (calcFunc-reduce 203 (calcFunc-reduce
192 '(var mul var-mul) 204 '(var mul var-mul)
193 (calcFunc-map 205 (calcFunc-map
@@ -198,45 +210,50 @@
198 (^ 10 (var a var-a)))) 210 (^ 10 (var a var-a))))
199 (calcFunc-map 211 (calcFunc-map
200 '(var mul var-mul) 212 '(var mul var-mul)
201 (cons 'vec (cdr (cdr coefnames))) 213 (cons 'vec (cdr (cdr calc-curve-coefnames)))
202 varnames)))))) 214 calc-curve-varnames))))))
203 ((memq key '(?x ?X)) 215 ((memq key '(?x ?X))
204 (calc-get-fit-variables nvars (1+ nvars) (and homog 0)) 216 (calc-get-fit-variables calc-curve-nvars
205 (setq model (math-mul coefnames 217 (1+ calc-curve-nvars) (and homog 0))
206 (cons 'vec (cons 1 (cdr varnames))))) 218 (setq calc-curve-model (math-mul calc-curve-coefnames
207 (setq model (if (eq key ?x) 219 (cons 'vec (cons 1 (cdr calc-curve-varnames)))))
208 (list 'calcFunc-exp model) 220 (setq calc-curve-model (if (eq key ?x)
209 (list '^ 10 model)))) 221 (list 'calcFunc-exp calc-curve-model)
222 (list '^ 10 calc-curve-model))))
210 ((memq key '(?l ?L)) 223 ((memq key '(?l ?L))
211 (calc-get-fit-variables nvars (1+ nvars) (and homog 0)) 224 (calc-get-fit-variables calc-curve-nvars
212 (setq model (math-mul coefnames 225 (1+ calc-curve-nvars) (and homog 0))
226 (setq calc-curve-model (math-mul calc-curve-coefnames
213 (cons 'vec 227 (cons 'vec
214 (cons 1 (cdr (calcFunc-map 228 (cons 1 (cdr (calcFunc-map
215 (if (eq key ?l) 229 (if (eq key ?l)
216 '(var ln var-ln) 230 '(var ln var-ln)
217 '(var log10 231 '(var log10
218 var-log10)) 232 var-log10))
219 varnames))))))) 233 calc-curve-varnames)))))))
220 ((= key ?q) 234 ((= key ?q)
221 (calc-get-fit-variables nvars (1+ (* 2 nvars)) (and homog 0)) 235 (calc-get-fit-variables calc-curve-nvars
222 (let ((c coefnames) 236 (1+ (* 2 calc-curve-nvars)) (and homog 0))
223 (v varnames)) 237 (let ((c calc-curve-coefnames)
224 (setq model (nth 1 c)) 238 (v calc-curve-varnames))
239 (setq calc-curve-model (nth 1 c))
225 (while (setq v (cdr v) c (cdr (cdr c))) 240 (while (setq v (cdr v) c (cdr (cdr c)))
226 (setq model (math-add 241 (setq calc-curve-model (math-add
227 model 242 calc-curve-model
228 (list '* 243 (list '*
229 (car c) 244 (car c)
230 (list '^ 245 (list '^
231 (list '- (car v) (nth 1 c)) 246 (list '- (car v) (nth 1 c))
232 2))))))) 247 2)))))))
233 ((= key ?g) 248 ((= key ?g)
234 (setq model (math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)") 249 (setq calc-curve-model
235 varnames '(vec (var XFit var-XFit)) 250 (math-read-expr "(AFit / BFit sqrt(2 pi)) exp(-0.5 * ((XFit - CFit) / BFit)^2)")
236 coefnames '(vec (var AFit var-AFit) 251 calc-curve-varnames '(vec (var XFit var-XFit))
252 calc-curve-coefnames '(vec (var AFit var-AFit)
237 (var BFit var-BFit) 253 (var BFit var-BFit)
238 (var CFit var-CFit))) 254 (var CFit var-CFit)))
239 (calc-get-fit-variables 1 (1- (length coefnames)) (and homog 1))) 255 (calc-get-fit-variables 1 (1- (length calc-curve-coefnames))
256 (and homog 1)))
240 ((memq key '(?\$ ?\' ?u ?U)) 257 ((memq key '(?\$ ?\' ?u ?U))
241 (let* ((defvars nil) 258 (let* ((defvars nil)
242 (record-entry nil)) 259 (record-entry nil))
@@ -244,74 +261,78 @@
244 (let* ((calc-dollar-values calc-arg-values) 261 (let* ((calc-dollar-values calc-arg-values)
245 (calc-dollar-used 0) 262 (calc-dollar-used 0)
246 (calc-hashes-used 0)) 263 (calc-hashes-used 0))
247 (setq model (calc-do-alg-entry "" "Model formula: ")) 264 (setq calc-curve-model (calc-do-alg-entry "" "Model formula: "))
248 (if (/= (length model) 1) 265 (if (/= (length calc-curve-model) 1)
249 (error "Bad format")) 266 (error "Bad format"))
250 (setq model (car model) 267 (setq calc-curve-model (car calc-curve-model)
251 record-entry t) 268 record-entry t)
252 (if (> calc-dollar-used 0) 269 (if (> calc-dollar-used 0)
253 (setq coefnames 270 (setq calc-curve-coefnames
254 (cons 'vec 271 (cons 'vec
255 (nthcdr (- (length calc-arg-values) 272 (nthcdr (- (length calc-arg-values)
256 calc-dollar-used) 273 calc-dollar-used)
257 (reverse calc-arg-values)))) 274 (reverse calc-arg-values))))
258 (if (> calc-hashes-used 0) 275 (if (> calc-hashes-used 0)
259 (setq coefnames 276 (setq calc-curve-coefnames
260 (cons 'vec (calc-invent-args 277 (cons 'vec (calc-invent-args
261 calc-hashes-used)))))) 278 calc-hashes-used))))))
262 (progn 279 (progn
263 (setq model (cond ((eq key ?u) 280 (setq calc-curve-model (cond ((eq key ?u)
264 (calc-var-value 'var-Model1)) 281 (calc-var-value 'var-Model1))
265 ((eq key ?U) 282 ((eq key ?U)
266 (calc-var-value 'var-Model2)) 283 (calc-var-value 'var-Model2))
267 (t (calc-top 1)))) 284 (t (calc-top 1))))
268 (or model (error "User model not yet defined")) 285 (or calc-curve-model (error "User model not yet defined"))
269 (if (math-vectorp model) 286 (if (math-vectorp calc-curve-model)
270 (if (and (memq (length model) '(3 4)) 287 (if (and (memq (length calc-curve-model) '(3 4))
271 (not (math-objvecp (nth 1 model))) 288 (not (math-objvecp (nth 1 calc-curve-model)))
272 (math-vectorp (nth 2 model)) 289 (math-vectorp (nth 2 calc-curve-model))
273 (or (null (nth 3 model)) 290 (or (null (nth 3 calc-curve-model))
274 (math-vectorp (nth 3 model)))) 291 (math-vectorp (nth 3 calc-curve-model))))
275 (setq varnames (nth 2 model) 292 (setq calc-curve-varnames (nth 2 calc-curve-model)
276 coefnames (or (nth 3 model) 293 calc-curve-coefnames
277 (cons 'vec 294 (or (nth 3 calc-curve-model)
278 (math-all-vars-but 295 (cons 'vec
279 model varnames))) 296 (math-all-vars-but
280 model (nth 1 model)) 297 calc-curve-model calc-curve-varnames)))
298 calc-curve-model (nth 1 calc-curve-model))
281 (error "Incorrect model specifier"))))) 299 (error "Incorrect model specifier")))))
282 (or varnames 300 (or calc-curve-varnames
283 (let ((with-y (eq (car-safe model) 'calcFunc-eq))) 301 (let ((with-y (eq (car-safe calc-curve-model) 'calcFunc-eq)))
284 (if coefnames 302 (if calc-curve-coefnames
285 (calc-get-fit-variables (if with-y (1+ nvars) nvars) 303 (calc-get-fit-variables
286 (1- (length coefnames)) 304 (if with-y (1+ calc-curve-nvars) calc-curve-nvars)
287 (math-all-vars-but 305 (1- (length calc-curve-coefnames))
288 model coefnames) 306 (math-all-vars-but
289 nil with-y) 307 calc-curve-model calc-curve-coefnames)
290 (let* ((coefs (math-all-vars-but model nil)) 308 nil with-y)
309 (let* ((coefs (math-all-vars-but calc-curve-model nil))
291 (vars nil) 310 (vars nil)
292 (n (- (length coefs) nvars (if with-y 2 1))) 311 (n (- (length coefs) calc-curve-nvars (if with-y 2 1)))
293 p) 312 p)
294 (if (< n 0) 313 (if (< n 0)
295 (error "Not enough variables in model")) 314 (error "Not enough variables in model"))
296 (setq p (nthcdr n coefs)) 315 (setq p (nthcdr n coefs))
297 (setq vars (cdr p)) 316 (setq vars (cdr p))
298 (setcdr p nil) 317 (setcdr p nil)
299 (calc-get-fit-variables (if with-y (1+ nvars) nvars) 318 (calc-get-fit-variables
300 (length coefs) 319 (if with-y (1+ calc-curve-nvars) calc-curve-nvars)
301 vars coefs with-y))))) 320 (length coefs)
321 vars coefs with-y)))))
302 (if record-entry 322 (if record-entry
303 (calc-record (list 'vec model varnames coefnames) 323 (calc-record (list 'vec calc-curve-model
324 calc-curve-varnames calc-curve-coefnames)
304 "modl")))) 325 "modl"))))
305 (t (beep)))) 326 (t (beep))))
306 (let ((calc-fit-to-trail t)) 327 (let ((calc-fit-to-trail t))
307 (calc-enter-result n (substring (symbol-name func) 9) 328 (calc-enter-result n (substring (symbol-name func) 9)
308 (list func model 329 (list func calc-curve-model
309 (if (= (length varnames) 2) 330 (if (= (length calc-curve-varnames) 2)
310 (nth 1 varnames) 331 (nth 1 calc-curve-varnames)
311 varnames) 332 calc-curve-varnames)
312 (if (= (length coefnames) 2) 333 (if (= (length calc-curve-coefnames) 2)
313 (nth 1 coefnames) 334 (nth 1 calc-curve-coefnames)
314 coefnames) 335 calc-curve-coefnames)
315 data)) 336 data))
316 (if (consp calc-fit-to-trail) 337 (if (consp calc-fit-to-trail)
317 (calc-record (calc-normalize calc-fit-to-trail) "parm")))))) 338 (calc-record (calc-normalize calc-fit-to-trail) "parm"))))))
@@ -340,7 +361,7 @@
340 (calc-invent-variables num but t base)))) 361 (calc-invent-variables num but t base))))
341 362
342(defun calc-get-fit-variables (nv nc &optional defv defc with-y homog) 363(defun calc-get-fit-variables (nv nc &optional defv defc with-y homog)
343 (or (= nv (if with-y (1+ nvars) nvars)) 364 (or (= nv (if with-y (1+ calc-curve-nvars) calc-curve-nvars))
344 (error "Wrong number of data vectors for this type of model")) 365 (error "Wrong number of data vectors for this type of model"))
345 (if (integerp defv) 366 (if (integerp defv)
346 (setq homog defv 367 (setq homog defv
@@ -388,12 +409,12 @@
388 (error "Expected %d parameter variable%s" nc (if (= nc 1) "" "s"))) 409 (error "Expected %d parameter variable%s" nc (if (= nc 1) "" "s")))
389 (if homog 410 (if homog
390 (setq coefs (cons 'vec (cons homog (cdr coefs))))) 411 (setq coefs (cons 'vec (cons homog (cdr coefs)))))
391 (if varnames 412 (if calc-curve-varnames
392 (setq model (math-multi-subst model (cdr varnames) (cdr vars)))) 413 (setq calc-curve-model (math-multi-subst calc-curve-model (cdr calc-curve-varnames) (cdr vars))))
393 (if coefnames 414 (if calc-curve-coefnames
394 (setq model (math-multi-subst model (cdr coefnames) (cdr coefs)))) 415 (setq calc-curve-model (math-multi-subst calc-curve-model (cdr calc-curve-coefnames) (cdr coefs))))
395 (setq varnames vars 416 (setq calc-curve-varnames vars
396 coefnames coefs))) 417 calc-curve-coefnames coefs)))
397 418
398 419
399 420
@@ -401,6 +422,9 @@
401;;; The following algorithms are from Numerical Recipes chapter 9. 422;;; The following algorithms are from Numerical Recipes chapter 9.
402 423
403;;; "rtnewt" with safety kludges 424;;; "rtnewt" with safety kludges
425
426(defvar var-DUMMY)
427
404(defun math-newton-root (expr deriv guess orig-guess limit) 428(defun math-newton-root (expr deriv guess orig-guess limit)
405 (math-working "newton" guess) 429 (math-working "newton" guess)
406 (let* ((var-DUMMY guess) 430 (let* ((var-DUMMY guess)
@@ -494,14 +518,20 @@
494 low vlow high vhigh)))))) 518 low vlow high vhigh))))))
495 519
496;;; Search for a root in an interval with no overt zero crossing. 520;;; Search for a root in an interval with no overt zero crossing.
521
522;; The variable math-root-widen is local to math-find-root, but
523;; is used by math-search-root, which is called (directly and
524;; indirectly) by math-find-root.
525(defvar math-root-widen)
526
497(defun math-search-root (expr deriv low vlow high vhigh) 527(defun math-search-root (expr deriv low vlow high vhigh)
498 (let (found) 528 (let (found)
499 (if root-widen 529 (if math-root-widen
500 (let ((iters 0) 530 (let ((iters 0)
501 (iterlim (if (eq root-widen 'point) 531 (iterlim (if (eq math-root-widen 'point)
502 (+ calc-internal-prec 10) 532 (+ calc-internal-prec 10)
503 20)) 533 20))
504 (factor (if (eq root-widen 'point) 534 (factor (if (eq math-root-widen 'point)
505 '(float 9 0) 535 '(float 9 0)
506 '(float 16 -1))) 536 '(float 16 -1)))
507 (prev nil) vprev waslow 537 (prev nil) vprev waslow
@@ -600,6 +630,9 @@
600 (list 'vec mid vmid))) 630 (list 'vec mid vmid)))
601 631
602;;; "mnewt" 632;;; "mnewt"
633
634(defvar math-root-vars [(var DUMMY var-DUMMY)])
635
603(defun math-newton-multi (expr jacob n guess orig-guess limit) 636(defun math-newton-multi (expr jacob n guess orig-guess limit)
604 (let ((m -1) 637 (let ((m -1)
605 (p guess) 638 (p guess)
@@ -624,9 +657,8 @@
624 (math-reject-arg nil "*Newton's method failed to converge")) 657 (math-reject-arg nil "*Newton's method failed to converge"))
625 (list 'vec next expr-val)))) 658 (list 'vec next expr-val))))
626 659
627(defvar math-root-vars [(var DUMMY var-DUMMY)])
628 660
629(defun math-find-root (expr var guess root-widen) 661(defun math-find-root (expr var guess math-root-widen)
630 (if (eq (car-safe expr) 'vec) 662 (if (eq (car-safe expr) 'vec)
631 (let ((n (1- (length expr))) 663 (let ((n (1- (length expr)))
632 (calc-symbolic-mode nil) 664 (calc-symbolic-mode nil)
@@ -710,7 +742,7 @@
710 var-DUMMY guess 742 var-DUMMY guess
711 vlow (math-evaluate-expr expr) 743 vlow (math-evaluate-expr expr)
712 vhigh vlow 744 vhigh vlow
713 root-widen 'point) 745 math-root-widen 'point)
714 (if (eq (car guess) 'intv) 746 (if (eq (car guess) 'intv)
715 (progn 747 (progn
716 (or (math-constp guess) (math-reject-arg guess 'constp)) 748 (or (math-constp guess) (math-reject-arg guess 'constp))
@@ -752,6 +784,8 @@
752 784
753;;; The following algorithms come from Numerical Recipes, chapter 10. 785;;; The following algorithms come from Numerical Recipes, chapter 10.
754 786
787(defvar math-min-vars [(var DUMMY var-DUMMY)])
788
755(defun math-min-eval (expr a) 789(defun math-min-eval (expr a)
756 (if (Math-vectorp a) 790 (if (Math-vectorp a)
757 (let ((m -1)) 791 (let ((m -1))
@@ -894,7 +928,7 @@
894 (tol (list 'float 1 (- -1 prec))) 928 (tol (list 'float 1 (- -1 prec)))
895 (zeps (list 'float 1 (- -5 prec))) 929 (zeps (list 'float 1 (- -5 prec)))
896 (e '(float 0 0)) 930 (e '(float 0 0))
897 u vu xm tol1 tol2 etemp p q r xv xw) 931 d u vu xm tol1 tol2 etemp p q r xv xw)
898 (while (progn 932 (while (progn
899 (setq xm (math-mul-float '(float 5 -1) 933 (setq xm (math-mul-float '(float 5 -1)
900 (math-add-float a b)) 934 (math-add-float a b))
@@ -1056,8 +1090,6 @@
1056 (list (math-add line-p xi) xi (nth 2 res)))) 1090 (list (math-add line-p xi) xi (nth 2 res))))
1057 1091
1058 1092
1059(defvar math-min-vars [(var DUMMY var-DUMMY)])
1060
1061(defun math-find-minimum (expr var guess min-widen) 1093(defun math-find-minimum (expr var guess min-widen)
1062 (let* ((calc-symbolic-mode nil) 1094 (let* ((calc-symbolic-mode nil)
1063 (n 0) 1095 (n 0)
@@ -1072,7 +1104,7 @@
1072 (math-dimension-error)) 1104 (math-dimension-error))
1073 (while (setq var (cdr var) guess (cdr guess)) 1105 (while (setq var (cdr var) guess (cdr guess))
1074 (or (eq (car-safe (car var)) 'var) 1106 (or (eq (car-safe (car var)) 'var)
1075 (math-reject-arg (car vg) "*Expected a variable")) 1107 (math-reject-arg (car var) "*Expected a variable"))
1076 (or (math-expr-contains expr (car var)) 1108 (or (math-expr-contains expr (car var))
1077 (math-reject-arg (car var) 1109 (math-reject-arg (car var)
1078 "*Formula does not contain specified variable")) 1110 "*Formula does not contain specified variable"))
@@ -1314,6 +1346,12 @@
1314 1346
1315 1347
1316;;; Open Romberg method; "qromo" in section 4.4. 1348;;; Open Romberg method; "qromo" in section 4.4.
1349
1350;; The variable math-ninteg-temp is local to math-ninteg-romberg,
1351;; but is used by math-ninteg-midpoint, which is used by
1352;; math-ninteg-romberg.
1353(defvar math-ninteg-temp)
1354
1317(defun math-ninteg-romberg (func expr lo hi mode) 1355(defun math-ninteg-romberg (func expr lo hi mode)
1318 (let ((curh '(float 1 0)) 1356 (let ((curh '(float 1 0))
1319 (h nil) 1357 (h nil)
@@ -1321,7 +1359,7 @@
1321 (j 0) 1359 (j 0)
1322 (ss nil) 1360 (ss nil)
1323 (prec calc-internal-prec) 1361 (prec calc-internal-prec)
1324 (integ-temp nil)) 1362 (math-ninteg-temp nil))
1325 (math-with-extra-prec 2 1363 (math-with-extra-prec 2
1326 ;; Limit on "j" loop must be 14 or less to keep "it" from overflowing. 1364 ;; Limit on "j" loop must be 14 or less to keep "it" from overflowing.
1327 (or (while (and (null ss) (<= (setq j (1+ j)) 8)) 1365 (or (while (and (null ss) (<= (setq j (1+ j)) 8))
@@ -1332,8 +1370,7 @@
1332 (if (math-lessp (math-abs (nth 1 res)) 1370 (if (math-lessp (math-abs (nth 1 res))
1333 (calcFunc-scf (math-abs (car res)) 1371 (calcFunc-scf (math-abs (car res))
1334 (- prec))) 1372 (- prec)))
1335 (setq math-ninteg-convergence j 1373 (setq ss (car res)))))
1336 ss (car res)))))
1337 (if (>= j 5) 1374 (if (>= j 5)
1338 (setq s (cdr s) 1375 (setq s (cdr s)
1339 h (cdr h))) 1376 h (cdr h)))
@@ -1354,15 +1391,15 @@
1354 res)) 1391 res))
1355 1392
1356 1393
1357(defun math-ninteg-midpoint (expr lo hi mode) ; uses "integ-temp" 1394(defun math-ninteg-midpoint (expr lo hi mode) ; uses "math-ninteg-temp"
1358 (if (eq mode 'inf) 1395 (if (eq mode 'inf)
1359 (let ((math-infinite-mode t) temp) 1396 (let ((math-infinite-mode t) temp)
1360 (setq temp (math-div 1 lo) 1397 (setq temp (math-div 1 lo)
1361 lo (math-div 1 hi) 1398 lo (math-div 1 hi)
1362 hi temp))) 1399 hi temp)))
1363 (if integ-temp 1400 (if math-ninteg-temp
1364 (let* ((it3 (* 3 (car integ-temp))) 1401 (let* ((it3 (* 3 (car math-ninteg-temp)))
1365 (math-working-step-2 (* 2 (car integ-temp))) 1402 (math-working-step-2 (* 2 (car math-ninteg-temp)))
1366 (math-working-step 0) 1403 (math-working-step 0)
1367 (range (math-sub hi lo)) 1404 (range (math-sub hi lo))
1368 (del (math-div range (math-float it3))) 1405 (del (math-div range (math-float it3)))
@@ -1371,7 +1408,7 @@
1371 (x (math-add lo (math-mul '(float 5 -1) del))) 1408 (x (math-add lo (math-mul '(float 5 -1) del)))
1372 (sum '(float 0 0)) 1409 (sum '(float 0 0))
1373 (j 0) temp) 1410 (j 0) temp)
1374 (while (<= (setq j (1+ j)) (car integ-temp)) 1411 (while (<= (setq j (1+ j)) (car math-ninteg-temp))
1375 (setq math-working-step (1+ math-working-step) 1412 (setq math-working-step (1+ math-working-step)
1376 temp (math-ninteg-evaluate expr x mode) 1413 temp (math-ninteg-evaluate expr x mode)
1377 math-working-step (1+ math-working-step) 1414 math-working-step (1+ math-working-step)
@@ -1379,17 +1416,17 @@
1379 expr (math-add x del2) 1416 expr (math-add x del2)
1380 mode))) 1417 mode)))
1381 x (math-add x del3))) 1418 x (math-add x del3)))
1382 (setq integ-temp (list it3 1419 (setq math-ninteg-temp (list it3
1383 (math-add (math-div (nth 1 integ-temp) 1420 (math-add (math-div (nth 1 math-ninteg-temp)
1384 '(float 3 0)) 1421 '(float 3 0))
1385 (math-mul sum del))))) 1422 (math-mul sum del)))))
1386 (setq integ-temp (list 1 (math-mul 1423 (setq math-ninteg-temp (list 1 (math-mul
1387 (math-sub hi lo) 1424 (math-sub hi lo)
1388 (math-ninteg-evaluate 1425 (math-ninteg-evaluate
1389 expr 1426 expr
1390 (math-mul (math-add lo hi) '(float 5 -1)) 1427 (math-mul (math-add lo hi) '(float 5 -1))
1391 mode))))) 1428 mode)))))
1392 (nth 1 integ-temp)) 1429 (nth 1 math-ninteg-temp))
1393 1430
1394 1431
1395 1432
@@ -1427,13 +1464,21 @@
1427 (math-with-extra-prec 2 1464 (math-with-extra-prec 2
1428 (math-general-fit expr vars coefs data 'full)))) 1465 (math-general-fit expr vars coefs data 'full))))
1429 1466
1467;; The variables math-fit-first-var, math-fit-first-coef and
1468;; math-fit-new-coefs are local to math-general-fit, but are used by
1469;; calcFunc-fitvar, calcFunc-fitparam and calcFunc-fitdummy
1470;; (respectively), which are used by math-general-fit.
1471(defvar math-fit-first-var)
1472(defvar math-fit-first-coef)
1473(defvar math-fit-new-coefs)
1474
1430(defun math-general-fit (expr vars coefs data mode) 1475(defun math-general-fit (expr vars coefs data mode)
1431 (let ((calc-simplify-mode nil) 1476 (let ((calc-simplify-mode nil)
1432 (math-dummy-counter math-dummy-counter) 1477 (math-dummy-counter math-dummy-counter)
1433 (math-in-fit 1) 1478 (math-in-fit 1)
1434 (extended (eq mode 'full)) 1479 (extended (eq mode 'full))
1435 (first-coef math-dummy-counter) 1480 (math-fit-first-coef math-dummy-counter)
1436 first-var 1481 math-fit-first-var
1437 (plain-expr expr) 1482 (plain-expr expr)
1438 orig-expr 1483 orig-expr
1439 have-sdevs need-chisq chisq 1484 have-sdevs need-chisq chisq
@@ -1441,7 +1486,7 @@
1441 (y-filter nil) 1486 (y-filter nil)
1442 y-dummy 1487 y-dummy
1443 (coef-filters nil) 1488 (coef-filters nil)
1444 new-coefs 1489 math-fit-new-coefs
1445 (xy-values nil) 1490 (xy-values nil)
1446 (weights nil) 1491 (weights nil)
1447 (var-YVAL nil) (var-YVALX nil) 1492 (var-YVAL nil) (var-YVALX nil)
@@ -1496,8 +1541,8 @@
1496 (setq dummy (math-dummy-variable) 1541 (setq dummy (math-dummy-variable)
1497 expr (math-expr-subst expr (car p) 1542 expr (math-expr-subst expr (car p)
1498 (list 'calcFunc-fitparam 1543 (list 'calcFunc-fitparam
1499 (- math-dummy-counter first-coef))))) 1544 (- math-dummy-counter math-fit-first-coef)))))
1500 (setq first-var math-dummy-counter 1545 (setq math-fit-first-var math-dummy-counter
1501 p vars) 1546 p vars)
1502 (while (setq p (cdr p)) 1547 (while (setq p (cdr p))
1503 (or (eq (car-safe (car p)) 'var) 1548 (or (eq (car-safe (car p)) 'var)
@@ -1505,8 +1550,8 @@
1505 (setq dummy (math-dummy-variable) 1550 (setq dummy (math-dummy-variable)
1506 expr (math-expr-subst expr (car p) 1551 expr (math-expr-subst expr (car p)
1507 (list 'calcFunc-fitvar 1552 (list 'calcFunc-fitvar
1508 (- math-dummy-counter first-var))))) 1553 (- math-dummy-counter math-fit-first-var)))))
1509 (if (< math-dummy-counter (+ first-var v)) 1554 (if (< math-dummy-counter (+ math-fit-first-var v))
1510 (setq dummy (math-dummy-variable))) ; dependent variable may be unnamed 1555 (setq dummy (math-dummy-variable))) ; dependent variable may be unnamed
1511 (setq y-dummy dummy 1556 (setq y-dummy dummy
1512 orig-expr expr) 1557 orig-expr expr)
@@ -1565,7 +1610,7 @@
1565 (setq sigmasqr (math-add (math-sqr (nth 2 xval)) 1610 (setq sigmasqr (math-add (math-sqr (nth 2 xval))
1566 (or sigmasqr 0)) 1611 (or sigmasqr 0))
1567 xval (nth 1 xval)))) 1612 xval (nth 1 xval))))
1568 (set (nth 2 (aref math-dummy-vars (+ first-var j))) xval) 1613 (set (nth 2 (aref math-dummy-vars (+ math-fit-first-var j))) xval)
1569 (setq j (1+ j))) 1614 (setq j (1+ j)))
1570 1615
1571 ;; Compute Y value for this data point. 1616 ;; Compute Y value for this data point.
@@ -1656,8 +1701,8 @@
1656 xy-values (cdr xy-values))))) 1701 xy-values (cdr xy-values)))))
1657 1702
1658 ;; Convert coefficients back into original terms. 1703 ;; Convert coefficients back into original terms.
1659 (setq new-coefs (copy-sequence beta)) 1704 (setq math-fit-new-coefs (copy-sequence beta))
1660 (let* ((bp new-coefs) 1705 (let* ((bp math-fit-new-coefs)
1661 (cp covar) 1706 (cp covar)
1662 (sigdat 1) 1707 (sigdat 1)
1663 (math-in-fit 3) 1708 (math-in-fit 3)
@@ -1673,9 +1718,9 @@
1673 (math-sqrt (math-mul (nth (setq j (1+ j)) 1718 (math-sqrt (math-mul (nth (setq j (1+ j))
1674 (car (setq cp (cdr cp)))) 1719 (car (setq cp (cdr cp))))
1675 sigdat)))))) 1720 sigdat))))))
1676 (setq new-coefs (math-evaluate-expr coef-filters)) 1721 (setq math-fit-new-coefs (math-evaluate-expr coef-filters))
1677 (if calc-fit-to-trail 1722 (if calc-fit-to-trail
1678 (let ((bp new-coefs) 1723 (let ((bp math-fit-new-coefs)
1679 (cp coefs) 1724 (cp coefs)
1680 (vec nil)) 1725 (vec nil))
1681 (while (setq bp (cdr bp) cp (cdr cp)) 1726 (while (setq bp (cdr bp) cp (cdr cp))
@@ -1695,7 +1740,7 @@
1695 (setq vec (cons (list 'calcFunc-fitparam n) vec) 1740 (setq vec (cons (list 'calcFunc-fitparam n) vec)
1696 n (1- n))) 1741 n (1- n)))
1697 vec) 1742 vec)
1698 (append (cdr new-coefs) (cdr vars)))) 1743 (append (cdr math-fit-new-coefs) (cdr vars))))
1699 1744
1700 ;; Package the result. 1745 ;; Package the result.
1701 (math-normalize 1746 (math-normalize
@@ -1719,20 +1764,20 @@
1719(defun calcFunc-fitvar (x) 1764(defun calcFunc-fitvar (x)
1720 (if (>= math-in-fit 2) 1765 (if (>= math-in-fit 2)
1721 (progn 1766 (progn
1722 (setq x (aref math-dummy-vars (+ first-var x -1))) 1767 (setq x (aref math-dummy-vars (+ math-fit-first-var x -1)))
1723 (or (calc-var-value (nth 2 x)) x)) 1768 (or (calc-var-value (nth 2 x)) x))
1724 (math-reject-arg x))) 1769 (math-reject-arg x)))
1725 1770
1726(defun calcFunc-fitparam (x) 1771(defun calcFunc-fitparam (x)
1727 (if (>= math-in-fit 2) 1772 (if (>= math-in-fit 2)
1728 (progn 1773 (progn
1729 (setq x (aref math-dummy-vars (+ first-coef x -1))) 1774 (setq x (aref math-dummy-vars (+ math-fit-first-coef x -1)))
1730 (or (calc-var-value (nth 2 x)) x)) 1775 (or (calc-var-value (nth 2 x)) x))
1731 (math-reject-arg x))) 1776 (math-reject-arg x)))
1732 1777
1733(defun calcFunc-fitdummy (x) 1778(defun calcFunc-fitdummy (x)
1734 (if (= math-in-fit 3) 1779 (if (= math-in-fit 3)
1735 (nth x new-coefs) 1780 (nth x math-fit-new-coefs)
1736 (math-reject-arg x))) 1781 (math-reject-arg x)))
1737 1782
1738(defun calcFunc-hasfitvars (expr) 1783(defun calcFunc-hasfitvars (expr)
@@ -1759,19 +1804,25 @@
1759 (sort (mapcar 'car vars) 1804 (sort (mapcar 'car vars)
1760 (function (lambda (x y) (string< (nth 1 x) (nth 1 y))))))) 1805 (function (lambda (x y) (string< (nth 1 x) (nth 1 y)))))))
1761 1806
1807;; The variables math-all-vars-vars (the vars for math-all-vars) and
1808;; math-all-vars-found are local to math-all-vars-in, but are used by
1809;; math-all-vars-rec which is called by math-all-vars-in.
1810(defvar math-all-vars-vars)
1811(defvar math-all-vars-found)
1812
1762(defun math-all-vars-in (expr) 1813(defun math-all-vars-in (expr)
1763 (let ((vars nil) 1814 (let ((math-all-vars-vars nil)
1764 found) 1815 math-all-vars-found)
1765 (math-all-vars-rec expr) 1816 (math-all-vars-rec expr)
1766 vars)) 1817 math-all-vars-vars))
1767 1818
1768(defun math-all-vars-rec (expr) 1819(defun math-all-vars-rec (expr)
1769 (if (Math-primp expr) 1820 (if (Math-primp expr)
1770 (if (eq (car-safe expr) 'var) 1821 (if (eq (car-safe expr) 'var)
1771 (or (math-const-var expr) 1822 (or (math-const-var expr)
1772 (if (setq found (assoc expr vars)) 1823 (if (setq math-all-vars-found (assoc expr math-all-vars-vars))
1773 (setcdr found (1+ (cdr found))) 1824 (setcdr math-all-vars-found (1+ (cdr math-all-vars-found)))
1774 (setq vars (cons (cons expr 1) vars))))) 1825 (setq math-all-vars-vars (cons (cons expr 1) math-all-vars-vars)))))
1775 (while (setq expr (cdr expr)) 1826 (while (setq expr (cdr expr))
1776 (math-all-vars-rec (car expr))))) 1827 (math-all-vars-rec (car expr)))))
1777 1828
generated by cgit v1.2.1 (git 2.18.0) at 2026-01-24 07:14:37 +0000