calcalg3.el

Revision 4220, 59.5 kB (checked in by miyoshi, 9 months ago)
Sync up with Emacs22.2.

Line
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 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, or (at your option)
14	;; 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; see the file COPYING. If not, write to the
23	;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
24	;; Boston, MA 02110-1301, USA.
25
26	;;; Commentary:
27
28	;;; Code:
29
30	;; This file is autoloaded from calc-ext.el.
31
32	(require 'calc-ext)
33	(require 'calc-macs)
34
35	(defun calc-find-root (var)
36	(interactive "sVariable(s) to solve for: ")
37	(calc-slow-wrapper
38	(let ((func (if (calc-is-hyperbolic) 'calcFunc-wroot 'calcFunc-root)))
39	(if (or (equal var "") (equal var "$"))
40	(calc-enter-result 2 "root" (list func
41	(calc-top-n 3)
42	(calc-top-n 1)
43	(calc-top-n 2)))
44	(let ((var (if (and (string-match ",\\\|[^ ] +[^ ]" var)
45	(not (string-match "\\[" var)))
46	(math-read-expr (concat "[" var "]"))
47	(math-read-expr var))))
48	(if (eq (car-safe var) 'error)
49	(error "Bad format in expression: %s" (nth 1 var)))
50	(calc-enter-result 1 "root" (list func
51	(calc-top-n 2)
52	var
53	(calc-top-n 1))))))))
54
55	(defun calc-find-minimum (var)
56	(interactive "sVariable(s) to minimize over: ")
57	(calc-slow-wrapper
58	(let ((func (if (calc-is-inverse)
59	(if (calc-is-hyperbolic)
60	'calcFunc-wmaximize 'calcFunc-maximize)
61	(if (calc-is-hyperbolic)
62	'calcFunc-wminimize 'calcFunc-minimize)))
63	(tag (if (calc-is-inverse) "max" "min")))
64	(if (or (equal var "") (equal var "$"))
65	(calc-enter-result 2 tag (list func
66	(calc-top-n 3)
67	(calc-top-n 1)
68	(calc-top-n 2)))
69	(let ((var (if (and (string-match ",\\\|[^ ] +[^ ]" var)
70	(not (string-match "\\[" var)))
71	(math-read-expr (concat "[" var "]"))
72	(math-read-expr var))))
73	(if (eq (car-safe var) 'error)
74	(error "Bad format in expression: %s" (nth 1 var)))
75	(calc-enter-result 1 tag (list func
76	(calc-top-n 2)
77	var
78	(calc-top-n 1))))))))
79
80	(defun calc-find-maximum (var)
81	(interactive "sVariable to maximize over: ")
82	(calc-invert-func)
83	(calc-find-minimum var))
84
85
86	(defun calc-poly-interp (arg)
87	(interactive "P")
88	(calc-slow-wrapper
89	(let ((data (calc-top 2)))
90	(if (or (consp arg) (eq arg 0) (eq arg 2))
91	(setq data (cons 'vec (calc-top-list 2 2)))
92	(or (null arg)
93	(error "Bad prefix argument")))
94	(if (calc-is-hyperbolic)
95	(calc-enter-result 1 "rati" (list 'calcFunc-ratint data (calc-top 1)))
96	(calc-enter-result 1 "poli" (list 'calcFunc-polint data
97	(calc-top 1)))))))
98
99	;; The variables calc-curve-nvars, calc-curve-varnames, calc-curve-model and calc-curve-coefnames are local to calc-curve-fit, but are
100	;; used by calc-get-fit-variables which is called by calc-curve-fit.
101	(defvar calc-curve-nvars)
102	(defvar calc-curve-varnames)
103	(defvar calc-curve-model)
104	(defvar calc-curve-coefnames)
105
106	(defvar calc-curve-fit-history nil
107	"History for calc-curve-fit.")
108
109	(defun calc-curve-fit (arg &optional calc-curve-model
110	calc-curve-coefnames calc-curve-varnames)
111	(interactive "P")
112	(calc-slow-wrapper
113	(setq calc-aborted-prefix nil)
114	(let ((func (if (calc-is-inverse) 'calcFunc-xfit
115	(if (calc-is-hyperbolic) 'calcFunc-efit
116	'calcFunc-fit)))
117	key (which 0)
118	n calc-curve-nvars temp data
119	(homog nil)
120	(msgs '( "(Press ? for help)"
121	"1 = linear or multilinear"
122	"2-9 = polynomial fits; i = interpolating polynomial"
123	"p = a x^b, ^ = a b^x"
124	"e = a exp(b x), x = exp(a + b x), l = a + b ln(x)"
125	"E = a 10^(b x), X = 10^(a + b x), L = a + b log10(x)"
126	"q = a + b (x-c)^2"
127	*"g = (a/b sqrt(2 pi)) exp(-0.5((x-c)/b)^2)"**
128	"h prefix = homogeneous model (no constant term)"
129	"' = alg entry, $ = stack, u = Model1, U = Model2")))
130	(while (not calc-curve-model)
131	(message "Fit to model: %s:%s"
132	(nth which msgs)
133	(if homog " h" ""))
134	(setq key (read-char))
135	(cond ((= key ?\C-g)
136	(keyboard-quit))
137	((= key ??)
138	(setq which (% (1+ which) (length msgs))))
139	((memq key '(?h ?H))
140	(setq homog (not homog)))
141	((progn
142	(if (eq key ?\$)
143	(setq n 1)
144	(setq n 0))
145	(cond ((null arg)
146	(setq n (1+ n)
147	data (calc-top n)))
148	((or (consp arg) (eq arg 0))
149	(setq n (+ n 2)
150	data (calc-top n)
151	data (if (math-matrixp data)
152	(append data (list (calc-top (1- n))))
153	(list 'vec data (calc-top (1- n))))))
154	((> (setq arg (prefix-numeric-value arg)) 0)
155	(setq data (cons 'vec (calc-top-list arg (1+ n)))
156	n (+ n arg)))
157	(t (error "Bad prefix argument")))
158	(or (math-matrixp data) (not (cdr (cdr data)))
159	(error "Data matrix is not a matrix!"))
160	(setq calc-curve-nvars (- (length data) 2)
161	calc-curve-coefnames nil
162	calc-curve-varnames nil)
163	nil))
164	((= key ?1) ; linear or multilinear
165	(calc-get-fit-variables calc-curve-nvars
166	(1+ calc-curve-nvars) (and homog 0))
167	(setq calc-curve-model (math-mul calc-curve-coefnames
168	(cons 'vec (cons 1 (cdr calc-curve-varnames))))))
169	((and (>= key ?2) (<= key ?9)) ; polynomial
170	(calc-get-fit-variables 1 (- key ?0 -1) (and homog 0))
171	(setq calc-curve-model
172	(math-build-polynomial-expr (cdr calc-curve-coefnames)
173	(nth 1 calc-curve-varnames))))
174	((= key ?i) ; exact polynomial
175	(calc-get-fit-variables 1 (1- (length (nth 1 data)))
176	(and homog 0))
177	(setq calc-curve-model
178	(math-build-polynomial-expr (cdr calc-curve-coefnames)
179	(nth 1 calc-curve-varnames))))
180	((= key ?p) ; power law
181	(calc-get-fit-variables calc-curve-nvars
182	(1+ calc-curve-nvars) (and homog 1))
183	(setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
184	(calcFunc-reduce
185	'(var mul var-mul)
186	(calcFunc-map
187	'(var pow var-pow)
188	calc-curve-varnames
189	(cons 'vec (cdr (cdr calc-curve-coefnames))))))))
190	((= key ?^) ; exponential law
191	(calc-get-fit-variables calc-curve-nvars
192	(1+ calc-curve-nvars) (and homog 1))
193	(setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
194	(calcFunc-reduce
195	'(var mul var-mul)
196	(calcFunc-map
197	'(var pow var-pow)
198	(cons 'vec (cdr (cdr calc-curve-coefnames)))
199	calc-curve-varnames)))))
200	((memq key '(?e ?E))
201	(calc-get-fit-variables calc-curve-nvars
202	(1+ calc-curve-nvars) (and homog 1))
203	(setq calc-curve-model (math-mul (nth 1 calc-curve-coefnames)
204	(calcFunc-reduce
205	'(var mul var-mul)
206	(calcFunc-map
207	(if (eq key ?e)
208	'(var exp var-exp)
209	'(calcFunc-lambda
210	(var a var-a)
211	(^ 10 (var a var-a))))
212	(calcFunc-map
213	'(var mul var-mul)
214	(cons 'vec (cdr (cdr calc-curve-coefnames)))
215	calc-curve-varnames))))))
216	((memq key '(?x ?X))
217	(calc-get-fit-variables calc-curve-nvars
218	(1+ calc-curve-nvars) (and homog 0))
219	(setq calc-curve-model (math-mul calc-curve-coefnames
220	(cons 'vec (cons 1 (cdr calc-curve-varnames)))))
221	(setq calc-curve-model (if (eq key ?x)
222	(list 'calcFunc-exp calc-curve-model)
223	(list '^ 10 calc-curve-model))))
224	((memq key '(?l ?L))
225	(calc-get-fit-variables calc-curve-nvars
226	(1+ calc-curve-nvars) (and homog 0))
227	(setq calc-curve-model (math-mul calc-curve-coefnames
228	(cons 'vec
229	(cons 1 (cdr (calcFunc-map
230	(if (eq key ?l)
231	'(var ln var-ln)
232	'(var log10
233	var-log10))
234	calc-curve-varnames)))))))
235	((= key ?q)
236	(calc-get-fit-variables calc-curve-nvars
237	(1+ (* 2 calc-curve-nvars)) (and homog 0))
238	(let ((c calc-curve-coefnames)
239	(v calc-curve-varnames))
240	(setq calc-curve-model (nth 1 c))
241	(while (setq v (cdr v) c (cdr (cdr c)))
242	(setq calc-curve-model (math-add
243	calc-curve-model
244	(list '*
245	(car c)
246	(list '^
247	(list '- (car v) (nth 1 c))
248	2)))))))
249	((= key ?g)
250	(setq calc-curve-model
251	(math-read-expr *"(AFit / BFit sqrt(2 pi)) exp(-0.5 ((XFit - CFit) / BFit)^2)"**)
252	calc-curve-varnames '(vec (var XFit var-XFit))
253	calc-curve-coefnames '(vec (var AFit var-AFit)
254	(var BFit var-BFit)
255	(var CFit var-CFit)))
256	(calc-get-fit-variables 1 (1- (length calc-curve-coefnames))
257	(and homog 1)))
258	((memq key '(?\$ ?\' ?u ?U))
259	(let* ((defvars nil)
260	(record-entry nil))
261	(if (eq key ?\')
262	(let* ((calc-dollar-values calc-arg-values)
263	(calc-dollar-used 0)
264	(calc-hashes-used 0))
265	(setq calc-curve-model (calc-do-alg-entry "" "Model formula: "
266	nil 'calc-curve-fit-history))
267	(if (/= (length calc-curve-model) 1)
268	(error "Bad format"))
269	(setq calc-curve-model (car calc-curve-model)
270	record-entry t)
271	(if (> calc-dollar-used 0)
272	(setq calc-curve-coefnames
273	(cons 'vec
274	(nthcdr (- (length calc-arg-values)
275	calc-dollar-used)
276	(reverse calc-arg-values))))
277	(if (> calc-hashes-used 0)
278	(setq calc-curve-coefnames
279	(cons 'vec (calc-invent-args
280	calc-hashes-used))))))
281	(progn
282	(setq calc-curve-model (cond ((eq key ?u)
283	(calc-var-value 'var-Model1))
284	((eq key ?U)
285	(calc-var-value 'var-Model2))
286	(t (calc-top 1))))
287	(or calc-curve-model (error "User model not yet defined"))
288	(if (math-vectorp calc-curve-model)
289	(if (and (memq (length calc-curve-model) '(3 4))
290	(not (math-objvecp (nth 1 calc-curve-model)))
291	(math-vectorp (nth 2 calc-curve-model))
292	(or (null (nth 3 calc-curve-model))
293	(math-vectorp (nth 3 calc-curve-model))))
294	(setq calc-curve-varnames (nth 2 calc-curve-model)
295	calc-curve-coefnames
296	(or (nth 3 calc-curve-model)
297	(cons 'vec
298	(math-all-vars-but
299	calc-curve-model calc-curve-varnames)))
300	calc-curve-model (nth 1 calc-curve-model))
301	(error "Incorrect model specifier")))))
302	(or calc-curve-varnames
303	(let ((with-y (eq (car-safe calc-curve-model) 'calcFunc-eq)))
304	(if calc-curve-coefnames
305	(calc-get-fit-variables
306	(if with-y (1+ calc-curve-nvars) calc-curve-nvars)
307	(1- (length calc-curve-coefnames))
308	(math-all-vars-but
309	calc-curve-model calc-curve-coefnames)
310	nil with-y)
311	(let* ((coefs (math-all-vars-but calc-curve-model nil))
312	(vars nil)
313	(n (- (length coefs) calc-curve-nvars (if with-y 2 1)))
314	p)
315	(if (< n 0)
316	(error "Not enough variables in model"))
317	(setq p (nthcdr n coefs))
318	(setq vars (cdr p))
319	(setcdr p nil)
320	(calc-get-fit-variables
321	(if with-y (1+ calc-curve-nvars) calc-curve-nvars)
322	(length coefs)
323	vars coefs with-y)))))
324	(if record-entry
325	(calc-record (list 'vec calc-curve-model
326	calc-curve-varnames calc-curve-coefnames)
327	"modl"))))
328	(t (beep))))
329	(let ((calc-fit-to-trail t))
330	(calc-enter-result n (substring (symbol-name func) 9)
331	(list func calc-curve-model
332	(if (= (length calc-curve-varnames) 2)
333	(nth 1 calc-curve-varnames)
334	calc-curve-varnames)
335	(if (= (length calc-curve-coefnames) 2)
336	(nth 1 calc-curve-coefnames)
337	calc-curve-coefnames)
338	data))
339	(if (consp calc-fit-to-trail)
340	(calc-record (calc-normalize calc-fit-to-trail) "parm"))))))
341
342	(defun calc-invent-independent-variables (n &optional but)
343	(calc-invent-variables n but '(x y z t) "x"))
344
345	(defun calc-invent-parameter-variables (n &optional but)
346	(calc-invent-variables n but '(a b c d) "a"))
347
348	(defun calc-invent-variables (num but names base)
349	(let ((vars nil)
350	(n num) (nn 0)
351	var)
352	(while (and (> n 0) names)
353	(setq var (math-build-var-name (if (consp names)
354	(car names)
355	(concat base (int-to-string
356	(setq nn (1+ nn)))))))
357	(or (math-expr-contains (cons 'vec but) var)
358	(setq vars (cons var vars)
359	n (1- n)))
360	(or (symbolp names) (setq names (cdr names))))
361	(if (= n 0)
362	(nreverse vars)
363	(calc-invent-variables num but t base))))
364
365	(defun calc-get-fit-variables (nv nc &optional defv defc with-y homog)
366	(or (= nv (if with-y (1+ calc-curve-nvars) calc-curve-nvars))
367	(error "Wrong number of data vectors for this type of model"))
368	(if (integerp defv)
369	(setq homog defv
370	defv nil))
371	(if homog
372	(setq nc (1- nc)))
373	(or defv
374	(setq defv (calc-invent-independent-variables nv)))
375	(or defc
376	(setq defc (calc-invent-parameter-variables nc defv)))
377	(let ((vars (read-string (format "Fitting variables (default %s; %s): "
378	(mapconcat 'symbol-name
379	(mapcar (function (lambda (v)
380	(nth 1 v)))
381	defv)
382	",")
383	(mapconcat 'symbol-name
384	(mapcar (function (lambda (v)
385	(nth 1 v)))
386	defc)
387	","))))
388	(coefs nil))
389	(setq vars (if (string-match "\\[" vars)
390	(math-read-expr vars)
391	(math-read-expr (concat "[" vars "]"))))
392	(if (eq (car-safe vars) 'error)
393	(error "Bad format in expression: %s" (nth 2 vars)))
394	(or (math-vectorp vars)
395	(error "Expected a variable or vector of variables"))
396	(if (equal vars '(vec))
397	(setq vars (cons 'vec defv)
398	coefs (cons 'vec defc))
399	(if (math-vectorp (nth 1 vars))
400	(if (and (= (length vars) 3)
401	(math-vectorp (nth 2 vars)))
402	(setq coefs (nth 2 vars)
403	vars (nth 1 vars))
404	(error
405	"Expected independent variables vector, then parameters vector"))
406	(setq coefs (cons 'vec defc))))
407	(or (= nv (1- (length vars)))
408	(and (not with-y) (= (1+ nv) (1- (length vars))))
409	(error "Expected %d independent variable%s" nv (if (= nv 1) "" "s")))
410	(or (= nc (1- (length coefs)))
411	(error "Expected %d parameter variable%s" nc (if (= nc 1) "" "s")))
412	(if homog
413	(setq coefs (cons 'vec (cons homog (cdr coefs)))))
414	(if calc-curve-varnames
415	(setq calc-curve-model (math-multi-subst calc-curve-model (cdr calc-curve-varnames) (cdr vars))))
416	(if calc-curve-coefnames
417	(setq calc-curve-model (math-multi-subst calc-curve-model (cdr calc-curve-coefnames) (cdr coefs))))
418	(setq calc-curve-varnames vars
419	calc-curve-coefnames coefs)))
420
421
422
423
424	;;; The following algorithms are from Numerical Recipes chapter 9.
425
426	;;; "rtnewt" with safety kludges
427
428	(defvar var-DUMMY)
429
430	(defun math-newton-root (expr deriv guess orig-guess limit)
431	(math-working "newton" guess)
432	(let* ((var-DUMMY guess)
433	next dval)
434	(setq next (math-evaluate-expr expr)
435	dval (math-evaluate-expr deriv))
436	(if (and (Math-numberp next)
437	(Math-numberp dval)
438	(not (Math-zerop dval)))
439	(progn
440	(setq next (math-sub guess (math-div next dval)))
441	(if (math-nearly-equal guess (setq next (math-float next)))
442	(progn
443	(setq var-DUMMY next)
444	(list 'vec next (math-evaluate-expr expr)))
445	(if (Math-lessp (math-abs-approx (math-sub next orig-guess))
446	limit)
447	(math-newton-root expr deriv next orig-guess limit)
448	(math-reject-arg next *"Newton's method failed to converge"**))))
449	(math-reject-arg next *"Newton's method encountered a singularity"**))))
450
451	;;; Inspired by "rtsafe"
452	(defun math-newton-search-root (expr deriv guess vguess ostep oostep
453	low vlow high vhigh)
454	(let ((var-DUMMY guess)
455	(better t)
456	pos step next vnext)
457	(if guess
458	(math-working "newton" (list 'intv 0 low high))
459	(math-working "bisect" (list 'intv 0 low high))
460	(setq ostep (math-mul-float (math-sub-float high low)
461	'(float 5 -1))
462	guess (math-add-float low ostep)
463	var-DUMMY guess
464	vguess (math-evaluate-expr expr))
465	(or (Math-realp vguess)
466	(progn
467	(setq ostep (math-mul-float ostep '(float 6 -1))
468	guess (math-add-float low ostep)
469	var-DUMMY guess
470	vguess (math-evaluate-expr expr))
471	(or (math-realp vguess)
472	(progn
473	(setq ostep (math-mul-float ostep '(float 123456 -5))
474	guess (math-add-float low ostep)
475	var-DUMMY guess
476	vguess nil))))))
477	(or vguess
478	(setq vguess (math-evaluate-expr expr)))
479	(or (Math-realp vguess)
480	(math-reject-arg guess *"Newton's method encountered a singularity"**))
481	(setq vguess (math-float vguess))
482	(if (eq (Math-negp vlow) (setq pos (Math-posp vguess)))
483	(setq high guess
484	vhigh vguess)
485	(if (eq (Math-negp vhigh) pos)
486	(setq low guess
487	vlow vguess)
488	(setq better nil)))
489	(if (or (Math-zerop vguess)
490	(math-nearly-equal low high))
491	(list 'vec guess vguess)
492	(setq step (math-evaluate-expr deriv))
493	(if (and (Math-realp step)
494	(not (Math-zerop step))
495	(setq step (math-div-float vguess (math-float step))
496	next (math-sub-float guess step))
497	(not (math-lessp-float high next))
498	(not (math-lessp-float next low)))
499	(progn
500	(setq var-DUMMY next
501	vnext (math-evaluate-expr expr))
502	(if (or (Math-zerop vnext)
503	(math-nearly-equal next guess))
504	(list 'vec next vnext)
505	(if (and better
506	(math-lessp-float (math-abs (or oostep
507	(math-sub-float
508	high low)))
509	(math-abs
510	(math-mul-float '(float 2 0)
511	step))))
512	(math-newton-search-root expr deriv nil nil nil ostep
513	low vlow high vhigh)
514	(math-newton-search-root expr deriv next vnext step ostep
515	low vlow high vhigh))))
516	(if (or (and (Math-posp vlow) (Math-posp vhigh))
517	(and (Math-negp vlow) (Math-negp vhigh)))
518	(math-search-root expr deriv low vlow high vhigh)
519	(math-newton-search-root expr deriv nil nil nil ostep
520	low vlow high vhigh))))))
521
522	;;; Search for a root in an interval with no overt zero crossing.
523
524	;; The variable math-root-widen is local to math-find-root, but
525	;; is used by math-search-root, which is called (directly and
526	;; indirectly) by math-find-root.
527	(defvar math-root-widen)
528
529	(defun math-search-root (expr deriv low vlow high vhigh)
530	(let (found)
531	(if math-root-widen
532	(let ((iters 0)
533	(iterlim (if (eq math-root-widen 'point)
534	(+ calc-internal-prec 10)
535	20))
536	(factor (if (eq math-root-widen 'point)
537	'(float 9 0)
538	'(float 16 -1)))
539	(prev nil) vprev waslow
540	diff)
541	(while (or (and (math-posp vlow) (math-posp vhigh))
542	(and (math-negp vlow) (math-negp vhigh)))
543	(math-working "widen" (list 'intv 0 low high))
544	(if (> (setq iters (1+ iters)) iterlim)
545	(math-reject-arg (list 'intv 0 low high)
546	*"Unable to bracket root"**))
547	(if (= iters calc-internal-prec)
548	(setq factor '(float 16 -1)))
549	(setq diff (math-mul-float (math-sub-float high low) factor))
550	(if (Math-zerop diff)
551	(setq high (calcFunc-incr high 10))
552	(if (math-lessp-float (math-abs vlow) (math-abs vhigh))
553	(setq waslow t
554	prev low
555	low (math-sub low diff)
556	var-DUMMY low
557	vprev vlow
558	vlow (math-evaluate-expr expr))
559	(setq waslow nil
560	prev high
561	high (math-add high diff)
562	var-DUMMY high
563	vprev vhigh
564	vhigh (math-evaluate-expr expr)))))
565	(if prev
566	(if waslow
567	(setq high prev vhigh vprev)
568	(setq low prev vlow vprev)))
569	(setq found t))
570	(or (Math-realp vlow)
571	(math-reject-arg vlow 'realp))
572	(or (Math-realp vhigh)
573	(math-reject-arg vhigh 'realp))
574	(let ((xvals (list low high))
575	(yvals (list vlow vhigh))
576	(pos (Math-posp vlow))
577	(levels 0)
578	(step (math-sub-float high low))
579	xp yp var-DUMMY)
580	(while (and (<= (setq levels (1+ levels)) 5)
581	(not found))
582	(setq xp xvals
583	yp yvals
584	step (math-mul-float step '(float 497 -3)))
585	(while (and (cdr xp) (not found))
586	(if (Math-realp (car yp))
587	(setq low (car xp)
588	vlow (car yp)))
589	(setq high (math-add-float (car xp) step)
590	var-DUMMY high
591	vhigh (math-evaluate-expr expr))
592	(math-working "search" high)
593	(if (and (Math-realp vhigh)
594	(eq (math-negp vhigh) pos))
595	(setq found t)
596	(setcdr xp (cons high (cdr xp)))
597	(setcdr yp (cons vhigh (cdr yp)))
598	(setq xp (cdr (cdr xp))
599	yp (cdr (cdr yp))))))))
600	(if found
601	(if (Math-zerop vhigh)
602	(list 'vec high vhigh)
603	(if (Math-zerop vlow)
604	(list 'vec low vlow)
605	(if deriv
606	(math-newton-search-root expr deriv nil nil nil nil
607	low vlow high vhigh)
608	(math-bisect-root expr low vlow high vhigh))))
609	(math-reject-arg (list 'intv 3 low high)
610	*"Unable to find a sign change in this interval"**))))
611
612	;;; "rtbis" (but we should be using Brent's method)
613	(defun math-bisect-root (expr low vlow high vhigh)
614	(let ((step (math-sub-float high low))
615	(pos (Math-posp vhigh))
616	var-DUMMY
617	mid vmid)
618	(while (not (or (math-nearly-equal low
619	(setq step (math-mul-float
620	step '(float 5 -1))
621	mid (math-add-float low step)))
622	(progn
623	(setq var-DUMMY mid
624	vmid (math-evaluate-expr expr))
625	(Math-zerop vmid))))
626	(math-working "bisect" mid)
627	(if (eq (Math-posp vmid) pos)
628	(setq high mid
629	vhigh vmid)
630	(setq low mid
631	vlow vmid)))
632	(list 'vec mid vmid)))
633
634	;;; "mnewt"
635
636	(defvar math-root-vars [(var DUMMY var-DUMMY)])
637
638	(defun math-newton-multi (expr jacob n guess orig-guess limit)
639	(let ((m -1)
640	(p guess)
641	p2 expr-val jacob-val next)
642	(while (< (setq p (cdr p) m (1+ m)) n)
643	(set (nth 2 (aref math-root-vars m)) (car p)))
644	(setq expr-val (math-evaluate-expr expr)
645	jacob-val (math-evaluate-expr jacob))
646	(unless (and (math-constp expr-val)
647	(math-constp jacob-val))
648	(math-reject-arg guess *"Newton's method encountered a singularity"**))
649	(setq next (math-add guess (math-div (math-float (math-neg expr-val))
650	(math-float jacob-val)))
651	p guess p2 next)
652	(math-working "newton" next)
653	(while (and (setq p (cdr p) p2 (cdr p2))
654	(math-nearly-equal (car p) (car p2))))
655	(if p
656	(if (Math-lessp (math-abs-approx (math-sub next orig-guess))
657	limit)
658	(math-newton-multi expr jacob n next orig-guess limit)
659	(math-reject-arg nil *"Newton's method failed to converge"**))
660	(list 'vec next expr-val))))
661
662
663	(defun math-find-root (expr var guess math-root-widen)
664	(if (eq (car-safe expr) 'vec)
665	(let ((n (1- (length expr)))
666	(calc-symbolic-mode nil)
667	(var-DUMMY nil)
668	(jacob (list 'vec))
669	p p2 m row)
670	(unless (eq (car-safe var) 'vec)
671	(math-reject-arg var 'vectorp))
672	(unless (= (length var) (1+ n))
673	(math-dimension-error))
674	(setq expr (copy-sequence expr))
675	(while (>= n (length math-root-vars))
676	(let ((symb (intern (concat "math-root-v"
677	(int-to-string
678	(length math-root-vars))))))
679	(setq math-root-vars (vconcat math-root-vars
680	(vector (list 'var symb symb))))))
681	(setq m -1)
682	(while (< (setq m (1+ m)) n)
683	(set (nth 2 (aref math-root-vars m)) nil))
684	(setq m -1 p var)
685	(while (setq m (1+ m) p (cdr p))
686	(or (eq (car-safe (car p)) 'var)
687	(math-reject-arg var *"Expected a variable"**))
688	(setq p2 expr)
689	(while (setq p2 (cdr p2))
690	(setcar p2 (math-expr-subst (car p2) (car p)
691	(aref math-root-vars m)))))
692	(unless (eq (car-safe guess) 'vec)
693	(math-reject-arg guess 'vectorp))
694	(unless (= (length guess) (1+ n))
695	(math-dimension-error))
696	(setq guess (copy-sequence guess)
697	p guess)
698	(while (setq p (cdr p))
699	(or (Math-numberp (car guess))
700	(math-reject-arg guess 'numberp))
701	(setcar p (math-float (car p))))
702	(setq p expr)
703	(while (setq p (cdr p))
704	(if (assq (car-safe (car p)) calc-tweak-eqn-table)
705	(setcar p (math-sub (nth 1 (car p)) (nth 2 (car p)))))
706	(setcar p (math-evaluate-expr (car p)))
707	(setq row (list 'vec)
708	m -1)
709	(while (< (setq m (1+ m)) n)
710	(nconc row (list (math-evaluate-expr
711	(or (calcFunc-deriv (car p)
712	(aref math-root-vars m)
713	nil t)
714	(math-reject-arg
715	expr
716	*"Formulas must be differentiable"**))))))
717	(nconc jacob (list row)))
718	(setq m (math-abs-approx guess))
719	(math-newton-multi expr jacob n guess guess
720	(if (math-zerop m) '(float 1 3) (math-mul m 10))))
721	(unless (eq (car-safe var) 'var)
722	(math-reject-arg var *"Expected a variable"**))
723	(unless (math-expr-contains expr var)
724	(math-reject-arg expr *"Formula does not contain specified var**