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