Add note that the lapack package needs to loaded to get the functions.
[maxima.git] / src / solve.lisp
blobe8a3e193fe9751f2e8cc5006f28f9e686e217c13
1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
4 ;;; ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
11 (in-package :maxima)
13 (macsyma-module solve)
15 (load-macsyma-macros ratmac strmac)
17 (declare-top (special expsumsplit *g
18 equations ;List of E-labels
19 *power *varb *flg
20 broken-not-freeof
21 mult ;Some crock which tracks multiplicities.
22 *roots ;alternating list of solutions and multiplicities
23 *failures ;alternating list of equations and multiplicities
24 *myvar
25 *has*var *var
26 xm* xn* mul*))
28 (defmvar $linsolvewarn t
29 "Needs to be documented.")
31 (defmvar $solvedecomposes t
32 "Causes `solve' to use `polydecomp' in attempting to solve polynomials.")
34 (defmvar $solveexplicit nil
35 "Causes `solve' to return implicit solutions i.e. of the form F(x)=0.")
37 (defmvar $solvenullwarn t
38 "Causes the user will be warned if SOLVE is called with either a
39 null equation list or a null variable list. For example,
40 SOLVE([],[]); would print two warning messages and return [].")
42 ;; Utility macros
44 ;; This macro returns the number of trivial equations. It counts up the
45 ;; number of zeros in a list.
47 ;(defmacro nzlist (llist)
48 ; `(do ((l ,llist (cdr l))
49 ; (zcount 0))
50 ; ((null l) zcount)
51 ; (if (and (integerp (car l)) (zerop (car l)))
52 ; (incf zcount))))
54 ;; This is only called on a variable.
56 (defmacro allroot (exp)
57 `(setq *failures (list* (make-mequal-simp ,exp ,exp) 1 *failures)))
59 ;; Finds variables, changes equations into expressions without MEQUAL.
60 ;; Checks for consistency between the number of unknowns and equations.
61 ;; Calls SOLVEX for simultaneous equations and SSOLVE for a single equation.
63 (defmfun $solve (*eql &optional (varl nil varl-p))
64 (setq $multiplicities (make-mlist))
65 (prog (eql ; Equations to solve
66 $keepfloat $ratfac ; In case user has set these
67 *roots ; *roots gets solutions,
68 *failures ; *failures "roots of"
69 broken-not-freeof) ;Has something to do with splitting up roots
71 ;; Create the equation list (this is a lisp list, not 'MLIST)
72 (setq eql
73 (cond
74 ;; If an atom, cons it.
75 ((atom *eql) (ncons *eql))
76 ;; If we have a list of equations, move everything over
77 ;; to one side, so x=5 -> x-5=0.
78 ((eq (g-rep-operator *eql) 'mlist)
79 (mapcar 'meqhk (mapcar 'meval (cdr *eql))))
80 ;; We can't solve inequalities
81 ((member (g-rep-operator *eql)
82 '(mnotequal mgreaterp mlessp mgeqp mleqp) :test #'eq)
83 (merror (intl:gettext "solve: cannot solve inequalities.")))
84 ;; Finally, assume we have just one equation, and put it
85 ;; on one side again.
86 (t (ncons (meqhk *eql)))))
88 (cond
89 ;; If the variable list wasn't supplied we have to supply it
90 ;; ourselves. Also remove constants like $%pi from the list.
91 ((null varl-p)
92 (setq varl
93 (let (($listconstvars nil))
94 (cdr ($listofvars *eql))))
95 (if varl (setq varl (remc varl)))) ; Remove all constants
97 ;; If we have got a variable list then if it's a list apply
98 ;; meval to each entry and then weed out duplicates. Else, just
99 ;; cons it.
101 (setq varl
102 (cond (($listp varl) (remove-duplicates
103 (mapcar #'meval (cdr varl))))
104 (t (list varl))))))
106 ;; Some sanity checks and warning messages.
107 (when (and (null varl) $solvenullwarn)
108 (mtell (intl:gettext "~&solve: variable list is empty, continuing anyway.~%")))
110 (when (and (null eql) $solvenullwarn)
111 (mtell (intl:gettext "~&solve: equation list is empty, continuing anyway.~%")))
113 (when (some #'mnump varl)
114 (merror (intl:gettext "solve: all variables must not be numbers.")))
116 ;; Deal with special cases.
117 (cond
118 ;; Trivially true equations for any set of variables.
119 ((equal eql '(0))
120 (return '$all))
122 ;; Trivially false equations: return []
123 ((or (null varl) (null eql))
124 (return (make-mlist-simp)))
126 ;; One equation in one variable: SSOLVE
127 ((and (null (cdr varl)) (null (cdr eql)))
128 (return (ssolve (car eql) (car varl))))
130 ;; We were given a variable list, or there are same # of eqns
131 ;; as unknowns: SOLVEX.
132 ((or varl-p
133 (= (length varl) (length eql)))
134 (setq eql (solvex eql varl (not $programmode) t))
135 (return
136 (cond ((and (cdr eql)
137 (not ($listp (cadr eql))))
138 (make-mlist eql))
139 (t eql)))))
141 ;; We don't know what to do, so complain. The let sets u to varl
142 ;; but as an MLIST list and e to the original eqns coerced to a
143 ;; list.
144 (let ((u (make-mlist-l varl))
145 (e (cond (($listp *eql) *eql)
146 (t (make-mlist *eql)))))
147 ;; MFORMAT doesn't have ~:[~] yet, so I just change this to
148 ;; make one of two possible calls to MERROR. Smaller codesize
149 ;; then what was here before anyway.
150 (if (> (length varl) (length eql))
151 (merror
152 (intl:gettext "solve: more unknowns than equations.~
153 ~%Unknowns given : ~%~M~
154 ~%Equations given: ~%~M")
155 u e)
156 (merror
157 (intl:gettext "solve: more equations than unknowns.~
158 ~%Unknowns given : ~%~M~
159 ~%Equations given: ~%~M")
160 u e)))))
163 ;; Removes anything from its list arg which solve considers not to be a
164 ;; variable, i.e. constants, functions or subscripted variables without
165 ;; numeric args.
167 (defun remc (lst)
168 (do ((l lst (cdr l)) (fl) (vl)) ((null l) vl)
169 (cond ((atom (setq fl (car l)))
170 (unless (maxima-constantp fl) (push fl vl)))
171 ((every #'$constantp (cdr fl)) (push fl vl)))))
173 ;; Solve a single equation for a single unknown.
174 ;; Obtains roots via solve and prints them.
176 (defun ssolve (exp *var)
177 (let (($solvetrigwarn $solvetrigwarn)
178 equations multi)
179 (cond ((null *var) '$all)
180 (t (solve exp *var 1)
181 (cond ((not (or *roots *failures)) (make-mlist))
182 ($programmode
183 (prog1
184 (make-mlist-l (nreverse (map2c #'(lambda (eqn mult) (push mult multi) eqn)
185 (if $solveexplicit
186 *roots
187 (nconc *roots *failures)))))
188 (setq $multiplicities (make-mlist-l (nreverse multi)))))
190 (let (soln)
191 (when (and *failures (not $solveexplicit))
192 (when $dispflag (mtell (intl:gettext "solve: the roots of:~%")))
193 (multiple-value-setq (soln equations)
194 (solve2 *failures equations)))
195 (when *roots
196 (when $dispflag (mtell (intl:gettext "solve: solution:~%")))
197 (multiple-value-setq (soln equations)
198 (solve2 *roots equations)))
199 (make-mlist-l equations))))))))
201 ;; Solve takes three arguments, the expression to solve for zero, the variable
202 ;; to solve for, and what multiplicity this solution is assumed to have (from
203 ;; higher-level Solve's). Solve returns NIL. Isn't that useful? The lists
204 ;; *roots and *failures are special variables to which Solve prepends solutions
205 ;; and their multiplicities in that order: *roots contains explicit solutions
206 ;; of the form <var>=<function of independent variables>, and *failures
207 ;; contains equations which if solved would yield additional solutions.
209 ;; Factors expression and reduces exponents by their gcd (via solventhp)
211 (defun solve (*exp *var mult &aux (genvar nil) ($derivsubst nil)
212 (exp (float2rat (mratcheck *exp)))
213 (*myvar *var) ($savefactors t))
214 (prog (factors *has*var genpairs $dontfactor temp symbol *g *checkfactors*
215 varlist expsumsplit)
216 (let (($ratfac t))
217 (setq exp (ratdisrep (ratf exp))))
218 ;; Cancel out any simple
219 ;; (non-algebraic) common factors in numerator and
220 ;; denominator without altering the structure of the
221 ;; expression too much.
222 ;; Also, RJFPROB in TEST;SOLVE TEST is now solved.
223 ;; - JPG
224 a (cond ((atom exp)
225 (cond ((eq exp *var)
226 (solve3 0 mult))
227 ((equal exp 0) (allroot *var))
228 (t nil)))
229 (t (setq exp (meqhk exp))
230 (cond ((equal exp '(0))
231 (return (allroot *var)))
232 ((free exp *var)
233 (return nil)))
234 (cond ((not (atom *var))
235 (setq symbol (gensym))
236 (setq exp (maxima-substitute symbol *var exp))
237 (setq temp *var)
238 (setq *var symbol)
239 (setq *myvar *var))) ;keep *MYVAR up-to-date
241 (cond ($solveradcan (setq exp (radcan1 exp *var))
242 (if (atom exp) (go a))))
244 (cond ((easy-cases exp *var mult)
245 (cond (symbol (setq *roots (subst temp *var *roots))
246 (setq *failures (subst temp *var *failures))))
247 (rootsort *roots)
248 (rootsort *failures)
249 (return nil)))
251 (cond ((setq factors (first-order-p exp *var))
252 (solve3 (ratdisrep
253 (ratf (make-mtimes -1 (div* (cdr factors)
254 (car factors)))))
255 mult))
257 (t (setq varlist (list *var))
258 (fnewvar exp)
259 (setq varlist (varsort varlist))
260 (let ((vartemp)
261 (ratnumer (mrat-numer (ratrep* exp)))
262 (numer-varlist varlist)
263 (subst-list (trig-subst-p varlist)))
264 (setq varlist (ncons *var))
265 (cond (subst-list
266 (setq exp (trig-subst exp subst-list))
267 (fnewvar exp)
268 (setq varlist (varsort varlist))
269 (setq exp (mrat-numer (ratrep* exp)))
270 (setq vartemp varlist))
271 (t (setq vartemp numer-varlist)
272 (setq exp ratnumer)))
273 (setq varlist vartemp))
275 (cond ((atom exp) (go a))
276 ((of-form-A*F<X>^N+B exp) (solve1a exp mult))
277 ((and (not (pcoefp exp))
278 (cddr exp)
279 (not (equal 1 (setq *g (solventhp (cdddr exp) (cadr exp))))))
280 (solventh exp *g))
281 (t (cond ($solvefactors
282 (map2c (lambda (x y) (solve1a x (m* mult y)))
283 (pfactor exp)))
284 (t (solve1a exp mult)))))))))
286 (cond (symbol (setq *roots (subst temp *var *roots))
287 (setq *failures (subst temp *var *failures))))
288 (rootsort *roots)
289 (rootsort *failures)
290 (return nil)))
292 (defun float2rat (exp)
293 (cond ((floatp exp) (setq exp (prep1 exp)) (make-rat-simp (car exp) (cdr exp)))
294 ((or (atom exp) (specrepp exp)) exp)
295 (t (recur-apply #'float2rat exp))))
297 ;;; The following takes care of cases where the expression is already in
298 ;;; factored form. This can introduce spurious roots if one of the factors
299 ;;; is an expression that can be undefined or infinity for certain values of
300 ;;; the variable in question. But soon this will be no worry because I will
301 ;;; add a list of "possible bad roots" to what $SOLVE returns.
302 ;;; Passes multiplicity to recursive calls to solve.
304 (defun easy-cases (*exp *var mult)
305 (cond ((or (atom *exp) (atom (car *exp))) nil)
306 ((eq (caar *exp) 'mtimes)
307 (do ((terms (cdr *exp) (cdr terms)))
308 ((null terms))
309 (solve (car terms) *var mult))
310 'mtimes)
312 ((eq (caar *exp) 'mabs) ;; abs(x) = 0 <=> x = 0
313 (solve (cadr *exp) *var mult)
314 'mabs)
316 ((eq (caar *exp) 'mexpt)
317 (cond ((and (freeof *var (cadr *exp))
318 (not (zerop1 (cadr *exp))))
319 ;; no solutions: c^x is never zero
320 'mexpt)
322 ((and (integerp (caddr *exp))
323 (plusp (caddr *exp)))
324 (solve (cadr *exp) *var (m* mult (caddr *exp)))
325 'mexprat)))))
327 ;;; Predicate to test for presence of troublesome trig functions to be
328 ;;; canonicalized. A table of when to make substitutions should
329 ;;; be used here.
330 ;;; trig kind => SIN | COS | TAN ... subst to make
331 ;;; number around in expression -> 1 1 0 ......
332 ;;; what you want to be able to do for example is to see if SIN and COS^2
333 ;;; are around and then make a reasonable substitution.
335 (defun trig-subst-p (vlist)
336 (and (not (trig-not-subst-p vlist))
337 (do ((var (car vlist) (car vlist))
338 (vlist (cdr vlist) (cdr vlist))
339 (subst-list))
340 ((null var) subst-list)
341 (cond ((and (not (atom var))
342 (trig-cannon (g-rep-operator var))
343 (not (free var *var)))
344 (push var subst-list))))))
346 ;; Predicate to see when obviously not to substitute for trigs.
347 ;; A hack in the direction of expression properties-table driven
348 ;; substitution. The "measure" of the expression is the total number
349 ;; of different kinds of trig functions in the expression.
351 (defun trig-not-subst-p (vlist)
352 (let ((trigs '(%sin %cos %tan %cot %csc %sec)))
353 (< (measure #'sign-gjc (operator-frequency-table vlist trigs) trigs)
354 2)))
356 ;; To get the total "value" of things in a table, this case an assoc list.
357 ;; (MEASURE FUNCTION ASSOCIATION-LIST SET) where FUNCTION is a function mapping
358 ;; the range of the ASSOCIATION-LIST viewed as a function on the SET, to the
359 ;; integers.
361 (defun measure (f alist set &aux (sum 0))
362 (dolist (element set)
363 (incf sum (funcall f (cdr (assoc element alist :test #'eq)))))
364 sum)
366 ;; Named for uniqueness only
368 (defun sign-gjc (x)
369 (cond ((or (null x) (= x 0)) 0)
370 ((< 0 x) 1)
371 (t -1)))
373 ;; A function that can EXTEND a function
374 ;; over two association lists. Note that I have been using association lists
375 ;; as mere functions (that is, as sets of ordered pairs).
376 ;; (EXTEND '+ L1 L2 S) could also be to take the union of two multi-sets in the
377 ;; sample space S. (what the '&%%#?& has this got to do with SOLVE?)
379 (defun extend (f l1 l2 s)
380 (do ((j 0 (1+ j))
381 (value nil))
382 ((= j (length s)) value)
383 (setq value (cons (cons (nth j s)
384 (funcall f (cdr (assoc (nth j s) l1 :test #'equal))
385 (cdr (assoc (nth j s) l2 :test #'equal))))
386 value))))
388 ;; For the case where the value of assoc is NIL, we will need a special "+"
390 (defun +mset (a b)
391 (+ (or a 0) (or b 0)))
393 ;; To recursively looks through a list
394 ;; structure (the VLIST) for members of the SET appearing in the MACSYMA
395 ;; functional position (caar list). Returning an assoc. list of appearance
396 ;; frequencies. Notice the use of EXTEND.
398 (defun operator-frequency-table (vlist set)
399 (do ((j 0 (1+ j))
400 (it)
401 (assl (do ((k 0 (1+ k))
402 (made nil))
403 ((= k (length set)) made)
404 (setq made (cons (cons (nth k set) 0)
405 made)))))
406 ((= j (length vlist)) assl)
407 (setq it (nth j vlist))
408 (cond ((atom it))
409 (t (setq assl (extend #'+mset (cons (cons (caar it) 1) nil)
410 assl set))
411 (setq assl (extend #'+mset assl
412 (operator-frequency-table (cdr it) set)
413 set))))))
415 (defun trig-subst (exp sub-list)
416 (do ((exp exp)
417 (sub-list (cdr sub-list) (cdr sub-list))
418 (var (car sub-list) (car sub-list)))
419 ((null var) exp)
420 (setq exp
421 (maxima-substitute (funcall (trig-cannon (g-rep-operator var))
422 (make-mlist-l (g-rep-operands var)))
423 var exp))))
425 ;; Here are the canonical trig substitutions.
427 (defun-prop (%sec trig-cannon) (x)
428 (inv* (make-g-rep '%cos (g-rep-first-operand x))))
430 (defun-prop (%csc trig-cannon) (x)
431 (inv* (make-g-rep '%sin (g-rep-first-operand x))))
433 (defun-prop (%tan trig-cannon) (x)
434 (div* (make-g-rep '%sin (g-rep-first-operand x))
435 (make-g-rep '%cos (g-rep-first-operand x))))
437 (defun-prop (%cot trig-cannon) (x)
438 (div* (make-g-rep '%cos (g-rep-first-operand x))
439 (make-g-rep '%sin (g-rep-first-operand x))))
441 (defun-prop (%sech trig-cannon) (x)
442 (inv* (make-g-rep '%cosh (g-rep-first-operand x))))
444 (defun-prop (%csch trig-cannon) (x)
445 (inv* (make-g-rep '%sinh (g-rep-first-operand x))))
447 (defun-prop (%tanh trig-cannon) (x)
448 (div* (make-g-rep '%sinh (g-rep-first-operand x))
449 (make-g-rep '%cosh (g-rep-first-operand x))))
451 (defun-prop (%coth trig-cannon) (x)
452 (div* (make-g-rep '%cosh (g-rep-first-operand x))
453 (make-g-rep '%sinh (g-rep-first-operand x))))
455 ;; Predicate to replace ISLINEAR....Returns NIL if not of for A*X+B, A and B
456 ;; freeof X, else returns (A . B)
458 (defun first-order-p (exp var &aux temp)
459 ;; Expand the expression at one level, i.e. distribute products
460 ;; over sums, but leave exponentiations alone.
461 ;; (X+1)^2*(X+Y) --> X*(X+1)^2 + Y*(X+1)^2
462 (setq exp (expand1 exp 1 1))
463 (cond ((atom exp) nil)
464 (t (case (g-rep-operator exp)
465 (mtimes
466 (cond ((setq temp (linear-term-p exp var))
467 (make-lineq temp 0))
468 (t nil)))
469 (mplus
470 (do ((arg (car (g-rep-operands exp)) (car rest))
471 (rest (cdr (g-rep-operands exp)) (cdr rest))
472 (linear-term-list)
473 (constant-term-list)
474 (temp))
475 ((null arg)
476 (if linear-term-list
477 (make-lineq (make-mplus-l linear-term-list)
478 (if constant-term-list
479 (make-mplus-l constant-term-list)
480 0))))
481 (cond ((setq temp (linear-term-p arg var))
482 (push temp linear-term-list))
483 ((broken-freeof var arg)
484 (push arg constant-term-list))
485 (t (return nil)))))
486 (t nil)))))
488 ;; Function to test if a term from an expanded expression is a linear term
489 ;; check and see that exactly one item in the product is the main var and
490 ;; all others are free of the main var. Returns NIL or a G-REP expression.
492 (defun linear-term-p (exp var)
493 (cond ((atom exp)
494 (cond ((eq exp var) 1)
495 (t nil)))
496 (t (case (g-rep-operator exp)
497 (mtimes
498 (do ((factor (car (g-rep-operands exp)) ;individual factors
499 (car rest))
500 (rest (cdr (g-rep-operands exp)) ;factors yet to be done
501 (cdr rest))
502 (main-var-p) ;nt -> main-var seen at top level
503 (list-of-factors)) ;accumulate our factors
504 ((null factor) ;for all factors
505 (and main-var-p
506 ;no-main-var at top level -=> not linear
507 (make-mtimes-l list-of-factors)))
508 (cond ((eq factor var) ;if it's our main var
509 ;note it...it has to be there to be a linear term
510 (setq main-var-p t))
511 ((broken-freeof var factor) ;if
512 (push factor list-of-factors))
513 (t (return nil)))))
514 (t nil)))))
517 ;;; DISPATCHING FUNCTION ON DEGREE OF EXPRESSION
518 ;;; This is a crock of shit, it should be data driven and be able to
519 ;;; dispatch to all manner of special cases that are in a table.
520 ;;; EXP here is a polynomial in MRAT form. All of this well-structured,
521 ;;; intelligently-designed code works by side effect. SOLVECUBIC
522 ;;; takes something that looks like (G0003 3 4 1 1 0 10) as an argument
523 ;;; and returns something like ((MEQUAL) $X ((MTIMES) ...)). You figure
524 ;;; out where the $X comes from.
526 ;;; It comes from GENVARS/VARLIST, of course. Isn't this wonderful rational
527 ;;; function package irrational? If you don't know about GENVARS and
528 ;;; VARLIST, you'd better bite the bullet and learn...everything depends
529 ;;; on them. The canonical example of mis-use of special variables!
530 ;;; --RWK
532 (defun solve1a (exp mult)
533 (let ((*myvar *myvar)
534 (*g nil))
535 (cond ((atom exp) nil)
536 ((not (memalike (setq *myvar (simplify (pdis (list (car exp) 1 1))))
537 *has*var))
538 nil)
539 ((equal (cadr exp) 1) (solvelin exp))
540 ((of-form-A*F<X>^N+B exp) (solve-A*F<X>^N+B exp t))
541 ((equal (cadr exp) 2) (solvequad exp))
542 ((not (equal 1 (setq *g (solventhp (cdddr exp) (cadr exp)))))
543 (solventh exp *g))
544 ((equal (cadr exp) 3) (solvecubic exp))
545 ((equal (cadr exp) 4) (solvequartic exp))
546 (t (let ((tt (solve-by-decomposition exp *myvar)))
547 (setq *failures (append (solution-losses tt) *failures))
548 (setq *roots (append (solution-wins tt) *roots)))))))
550 (defun solve-simplist (list-of-things)
551 (g-rep-operands (simplifya (make-mlist-l list-of-things) nil)))
553 ;; The Solve-by-decomposition program returns the cons of (ROOTS . FAILURES).
554 ;; It returns a "Solution" object, that is, a CONS with the CAR being the
555 ;; failures and the CDR being the successes.
556 ;; It takes a POLY as an argument and returns a SOLUTION.
558 (defun solve-by-decomposition (poly *$var)
559 (let ((decomp))
560 (cond ((or (not $solvedecomposes)
561 (= (length (setq decomp (polydecomp poly (poly-var poly)))) 1))
562 (make-solution nil `(,(make-mequal 0 (pdis poly)) 1)))
563 (t (decomp-trace (make-mequal 0 (rdis (car decomp)))
564 decomp
565 (poly-var poly) *$var 1)))))
567 ;; DECOMP-TRACE is the recursive function which maps itself down the
568 ;; intermediate solutions until the end is reached. If it encounters
569 ;; non-solvable equations it stops. It returns a SOLUTION object, that is, a
570 ;; CONS with the CAR being the failures and the CDR being the successes.
572 (defun decomp-trace (eqn decomp var *$var mult &aux sol chain-sol wins losses)
573 (setq sol (if decomp
574 (re-solve eqn *$var mult)
575 (make-solution `(,eqn 1) nil)))
576 (cond ((solution-losses sol) sol)
577 ;; End test
578 ((null decomp) sol)
579 (t (do ((l (solution-wins sol) (cddr l)))
580 ((null l))
581 (setq chain-sol
582 (decomp-chain (car l) (cdr decomp) var *$var (cadr l)))
583 (setq wins (nconc wins (copy-list (solution-wins chain-sol))))
584 (setq losses (nconc losses (copy-list (solution-losses chain-sol)))))
585 (make-solution wins losses))))
587 ;; Decomp-chain is the function which formats the mess for the recursive call.
588 ;; It returns a "Solution" object, that is, a CONS with the CAR being the
589 ;; failures and the CDR being the successes.
591 (defun decomp-chain (rsol decomp var *$var mult)
592 (let ((sol (simplify (make-mequal (rdis (if decomp (car decomp)
593 ;; Include the var itself in the decomposition
594 (make-mrat-body (make-mrat-poly var '(1 1)) 1)))
595 (mequal-rhs rsol)))))
596 (decomp-trace sol decomp var *$var mult)))
598 ;; RE-SOLVE calls SOLVE recursively, returning a SOLUTION object.
599 ;; Will not decompose or factor.
601 (defun re-solve (eqn var mult)
602 (let ((*roots nil)
603 (*failures nil)
604 ;; We've already decomposed and factored
605 ($solvedecomposes)
606 ($solvefactors))
607 (solve eqn var mult)
608 (make-solution *roots *failures)))
610 ;; SOLVENTH programs test to see if the variable of interest appears
611 ;; to some power in all terms. If so, a new variable is substituted for it
612 ;; and the simpler expression solved with the multiplicity
613 ;; adjusted accordingly.
614 ;; SOLVENTHP returns gcd of exponents.
616 (defun solventhp (l gcd)
617 (cond ((null l) gcd)
618 ((equal gcd 1) 1)
619 (t (solventhp (cddr l)
620 (gcd (car l) gcd)))))
622 ;; Reduces exponents by their gcd.
624 (defun solventh (exp *g)
625 (let ((*varb (pdis (make-mrat-poly (poly-var exp) '(1 1))))
626 (exp (make-mrat-poly (poly-var exp) (solventh1 (poly-terms exp)))))
627 (let* ((rts (re-solve-full (pdis exp) *varb))
628 (fails (solution-losses rts))
629 (wins (solution-wins rts))
630 (*power (make-mexpt *varb *g)))
631 (map2c #'(lambda (w z)
632 (cond ((atom *varb)
633 (solve (make-mequal *power (mequal-rhs w)) *varb z))
634 (t (let ((rts (re-solve-full
635 (make-mequal *power (mequal-rhs w))
636 *varb)))
637 (map2c #'(lambda (root mult)
638 (solve (make-mequal (mequal-rhs root) 0)
639 *myvar mult))
640 (solution-wins rts))))))
641 wins)
642 (map2c #'(lambda (w z)
643 (push z *failures)
644 (push (solventh3 w *power *varb) *failures))
645 fails)
646 *roots)))
648 (defun solventh3 (w *power *varb &aux varlist genvar *flg w1 w2)
649 (cond ((broken-freeof *varb w) w)
650 (t (setq w1 (ratf (cadr w)))
651 (setq w2 (ratf (caddr w)))
652 (setq varlist
653 (mapcar #'(lambda (h)
654 (cond (*flg h)
655 ((alike1 h *varb)
656 (setq *flg t)
657 *power)
658 (t h)))
659 varlist))
660 (list (car w) (rdis (cdr w1)) (rdis (cdr w2))))))
662 (defun solventh1 (l)
663 (cond ((null l) nil)
664 (t (cons (quotient (car l) *g)
665 (cons (cadr l) (solventh1 (cddr l)))))))
667 ;; Will decompose or factor
669 (defun re-solve-full (x var &aux *roots *failures)
670 (solve x var mult)
671 (make-solution *roots *failures))
673 ;; Sees if expression is of the form A*F<X>^N+B.
675 (defun of-form-A*F<X>^N+B (e)
676 (and (memalike (simplify (pdis (list (car e) 1 1))) *has*var)
677 (or (atom (caddr e))
678 (not (memalike (simplify (pdis (list (caaddr e) 1 1)))
679 *has*var)))
680 (or (null (cdddr e)) (equal (cadddr e) 0))))
682 ;; Solves the special case A*F<X>^N+B.
684 (defun solve-A*F<X>^N+B (exp $%emode)
685 (prog (a b c)
686 (setq a (pdis (caddr exp)))
687 (setq c (pdis (list (car exp) 1 1)))
688 (cond ((null (cdddr exp))
689 (return (solve c *var (* (cadr exp) mult)))))
690 (setq b (pdis (pminus (cadddr (cdr exp)))))
691 (return (solve-A*F<X>^N+B1 c
692 (simpnrt (div* b a) (cadr exp))
693 (make-rat 1 (cadr exp))
694 (cadr exp)))))
696 (defun solve-A*F<X>^N+B1 (var root n thisn)
697 (do ((thisn thisn (1- thisn))) ((zerop thisn))
698 (solve (add* var (mul* -1 root (power* '$%e (mul* 2 '$%pi '$%i thisn n))))
699 *var mult)))
702 ;; ADISPLINE displays a line like DISPLINE, and in addition, notes that it is
703 ;; not free of *VAR if it isn't.
705 (defun adispline (line)
706 ;; This may be redundant, but nice if ADISPLINE gets used where not needed.
707 (cond ((and $breakup (not $programmode))
708 (let ((linelabel (displine line)))
709 (cond ((broken-freeof *var line))
710 (t (setq broken-not-freeof
711 (cons linelabel broken-not-freeof))))
712 linelabel))
713 (t (displine line))))
715 ;; Predicate to check if an expression which may be broken up
716 ;; is freeof
718 (setq broken-not-freeof nil)
720 ;; For consistency, use backwards args.
721 ;; == (freeof var exp) but works even if solution is broken up ($breakup=t)
722 (defun broken-freeof (var exp)
723 (cond ($breakup
724 (do ((b-n-fo var (car b-n-fo-l))
725 (b-n-fo-l broken-not-freeof (cdr b-n-fo-l)))
726 ((null b-n-fo) t)
727 (and (not (argsfreeof b-n-fo exp))
728 (return nil))))
729 (t (argsfreeof var exp))))
731 ;; Adds solutions to roots list.
732 ;; Solves for inverse of functions (via USOLVE)
734 (defun solve3 (exp mult)
735 (setq exp (simplify exp))
736 (cond ((not (broken-freeof *var exp))
737 (push mult *failures)
738 (push (make-mequal-simp (simplify *myvar) exp) *failures))
739 (t (cond ((eq *myvar *var)
740 (push mult *roots)
741 (push (make-mequal-simp *var exp) *roots))
742 ((atom *myvar)
743 (push mult *failures)
744 (push (make-mequal-simp *myvar exp) *failures))
745 (t (usolve exp (g-rep-operator *myvar)))))))
748 ;; Solve a linear equation. Argument is a polynomial in pseudo-cre form.
749 ;; This function is called for side-effect only.
751 (defun solvelin (exp)
752 (cond ((equal 0 (ptterm (cdr exp) 0))
753 (solve1a (caddr exp) mult)))
754 (solve3 (rdis (ratreduce (pminus (ptterm (cdr exp) 0))
755 (caddr exp)))
756 mult))
758 ;; Solve a quadratic equation. Argument is a polynomial in pseudo-cre form.
759 ;; This function is called for side-effect only.
760 ;; The code for handling the case where the discriminant = 0 seems to never
761 ;; be run. Presumably, the expression is factored higher up.
763 (defun solvequad (exp &aux discrim a b c)
764 (setq a (caddr exp))
765 (setq b (ptterm (cdr exp) 1.))
766 (setq c (ptterm (cdr exp) 0.))
767 (setq discrim (simplify (pdis (pplus (pexpt b 2.)
768 (pminus (ptimes 4. (ptimes a c)))))))
769 (setq b (pdis (pminus b)))
770 (setq a (pdis (ptimes 2. a)))
771 ;; At this point, everything is back in general representation.
772 (let ((varlist nil)) ;;2/6/2002 RJF
773 (cond ((equal 0 discrim)
774 (solve3 (fullratsimp `((mquotient) ,b ,a))
775 (* 2 mult)))
776 (t (setq discrim (simpnrt discrim 2))
777 (solve3 (fullratsimp `((mquotient) ((mplus) ,b ,discrim) ,a))
778 mult)
779 (solve3 (fullratsimp `((mquotient) ((mplus) ,b ((mminus) ,discrim)) ,a))
780 mult)))))
782 ;; Reorders V so that members which contain the variable of
783 ;; interest come first.
785 (defun varsort (v)
786 (let ((*u nil)
787 (*v (copy-list v)))
788 (mapc #'(lambda (z)
789 (cond ((broken-freeof *var z)
790 (setq *u (cons z *u))
791 (setq *v (delete z *v :count 1 :test #'equal)))))
793 (setq $dontfactor *u)
794 (setq *has*var (mapcar #'resimplify *v))
795 (append *u *v)))
797 ;; Solves for variable when it occurs within a function by taking the inverse.
798 ;; When this code is fixed, the `((mplus) ,x ,y) forms should be rewritten as
799 ;; (MAKE-MPLUS X Y). I didn't do this because the code was buggy and it should
800 ;; be fixed first. - cwh
801 ;; You mean you didn't do it because you were buggy. Hope you're fixed soon!
802 ;; --RWK
804 ;; Solve <exp> = <*myvar> for <*var>, where <*myvar>=<op>(...)
805 (defun usolve (exp op)
806 (prog (inverse)
807 (setq inverse
808 (cond
809 ((eq op 'mexpt)
810 (cond ((broken-freeof *var
811 (cadr *myvar))
812 (cond ((equal exp 0)
813 (go fail)))
814 `((mplus) ((mminus) ,(caddr *myvar))
815 ,(div* `((%log) ,exp)
816 `((%log) ,(cadr *myvar)))))
817 ((broken-freeof *var
818 (caddr *myvar))
819 (cond ((equal exp 0)
820 (cond ((mnegp (caddr *myvar))
821 (go fail))
822 (t (cadr *myvar))))
823 ;; There is a bug right here.
824 ;; SOLVE(SQRT(U)+1) should return U=1
825 ;; This code is entered with EXP = -1, OP = MEXPT
826 ;; *VAR = U, and *MYVAR = ((MEXPT) U ((RAT) 1 2))
827 ;; BULLSHIT -- RWK. That is precisely the bug
828 ;; this code was added to fix!
829 ((and (not (eq (ask-integer (caddr *myvar)
830 '$integer)
831 '$yes))
832 (free exp '$%i)
833 (eq ($asksign exp) '$neg))
834 (go fail))
835 (t `((mplus) ,(cadr *myvar)
836 ((mminus)
837 ((mexpt) ,exp
838 ,(div* 1 (caddr *myvar))))))))
839 (t (go fail))))
840 ((setq inverse (get op '$inverse))
841 (when (and $solvetrigwarn
842 (member op '(%sin %cos %tan %sec %csc %cot %cosh %sech) :test #'eq))
843 (mtell (intl:gettext "~&solve: using arc-trig functions to get a solution.~%Some solutions will be lost.~%"))
844 (setq $solvetrigwarn nil))
845 `((mplus) ((mminus) ,(cadr *myvar))
846 ((,inverse) ,exp)))
847 ((eq op '%log)
848 `((mplus) ((mminus) ,(cadr *myvar))
849 ((mexpt) $%e ,exp)))
850 (t (go fail))))
851 (return (solve (simplify inverse) *var mult))
852 fail (return (setq *failures
853 (cons (simplify `((mequal) ,*myvar ,exp))
854 (cons mult *failures))))))
856 ;; Predicate for determining if an expression is messy enough to
857 ;; generate a new linelabel for it.
858 ;; Expression must be in general form.
860 (defun complicated (exp)
861 (and $breakup
862 (not $programmode)
863 (not (free exp 'mplus))))
865 (defun rootsort (l)
866 (prog (a fm fm1)
867 g1 (cond ((null l) (return nil)))
868 (setq a (car (setq fm l)))
869 (setq fm1 (cdr fm))
870 loop (cond ((null (cddr fm)) (setq l (cddr l)) (go g1))
871 ((alike1 (caddr fm) a)
872 (rplaca fm1 (+ (car fm1) (cadddr fm)))
873 (rplacd (cdr fm) (cddddr fm))
874 (go loop)))
875 (setq fm (cddr fm))
876 (go loop)))
878 (defmfun $linsolve (eql varl)
879 (let (($ratfac))
880 (setq eql (if ($listp eql) (cdr eql) (ncons eql)))
881 (setq varl (if ($listp varl)
882 (delete-duplicates (cdr varl) :test #'equal :from-end t)
883 (ncons varl)))
884 (do ((varl varl (cdr varl)))
885 ((null varl))
886 (when (mnump (car varl))
887 (merror (intl:gettext "solve: variable must not be a number; found: ~M") (car varl))))
888 (if (null varl)
889 (make-mlist-simp)
890 (solvex (mapcar 'meqhk eql) varl (not $programmode) nil))))
892 (defun solvex (eql varl ind flag &aux ($algebraic $algebraic))
893 (declare (special xa*))
894 (prog (*varl ans varlist genvar xm* xn* mul*)
895 (setq *varl varl)
896 (setq eql (mapcar #'(lambda (x) ($ratdisrep ($ratnumer x))) eql))
897 (cond ((atom (ignore-rat-err (formx flag 'xa* eql varl)))
898 ;; This flag is T if called from SOLVE
899 ;; and NIL if called from LINSOLVE.
900 (cond (flag (return ($algsys (make-mlist-l eql)
901 (make-mlist-l varl))))
902 (t (merror (intl:gettext "linsolve: cannot solve a nonlinear equation."))))))
903 (setq ans (tfgeli 'xa* xn* xm*))
904 (if (and $linsolvewarn (car ans))
905 (mtell (intl:gettext "~&solve: dependent equations eliminated: ~A~%") (car ans)))
906 (if (cadr ans)
907 (return '((mlist simp))))
908 (do ((j 0 (1+ j)))
909 ((> j xm*))
910 ;;I put this in the value cell--wfs
911 (setf (aref xa* 0 j) nil))
912 (ptorat 'xa* xn* xm*)
913 (setq varl
914 (xrutout 'xa* xn* xm*
915 (mapcar #'(lambda (x) (ith varl x))
916 (caddr ans))
917 ind))
918 (if $programmode
919 (setq varl (make-mlist-l (linsort (cdr varl) *varl))))
920 (return varl)))
922 ;; (LINSORT '(((MEQUAL) A2 FOO) ((MEQUAL) A3 BAR)) '(A3 A2))
923 ;; returns (((MEQUAL) A3 BAR) ((MEQUAL) A2 FOO)) .
925 (defun linsort (meq-list var-list)
926 (mapcar #'(lambda (x) (cons (caar meq-list) x))
927 (sort (mapcar #'cdr meq-list)
928 #'(lambda (x y) (member y (member x var-list :test #'equal) :test #'equal)) :key #'car)))