]> code.delx.au - gnu-emacs/blob - lisp/calc/calcalg3.el
delx.au / code / gnu-emacs / blob


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