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