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