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