1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
13 (macsyma-module solve
)
15 (load-macsyma-macros ratmac strmac
)
17 (declare-top (special expsumsplit
*g
18 equations
;List of E-labels
21 mult
;Some crock which tracks multiplicities.
22 *roots
;alternating list of solutions and multiplicities
23 *failures
;alternating list of equations and multiplicities
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 [].")
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))
51 ; (if (and (integerp (car l)) (zerop (car l)))
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)
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
86 (t (ncons (meqhk *eql
)))))
89 ;; If the variable list wasn't supplied we have to supply it
90 ;; ourselves. Also remove constants like $%pi from the list.
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
102 (cond (($listp varl
) (remove-duplicates
103 (mapcar #'meval
(cdr 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.
118 ;; Trivially true equations for any set of variables.
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.
133 (= (length varl
) (length eql
)))
134 (setq eql
(solvex eql varl
(not $programmode
) t
))
136 (cond ((and (cdr eql
)
137 (not ($listp
(cadr 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
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
))
152 (intl:gettext
"solve: more unknowns than equations.~
153 ~%Unknowns given : ~%~M~
154 ~%Equations given: ~%~M")
157 (intl:gettext
"solve: more equations than unknowns.~
158 ~%Unknowns given : ~%~M~
159 ~%Equations given: ~%~M")
163 ;; Removes anything from its list arg which solve considers not to be a
164 ;; variable, i.e. constants, functions or subscripted variables without
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
)
179 (cond ((null *var
) '$all
)
180 (t (solve exp
*var
1)
181 (cond ((not (or *roots
*failures
)) (make-mlist))
184 (make-mlist-l (nreverse (map2c #'(lambda (eqn mult
) (push mult multi
) eqn
)
187 (nconc *roots
*failures
)))))
188 (setq $multiplicities
(make-mlist-l (nreverse multi
)))))
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
)))
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
*
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.
227 ((equal exp
0) (allroot *var
))
229 (t (setq exp
(meqhk exp
))
230 (cond ((equal exp
'(0))
231 (return (allroot *var
)))
234 (cond ((not (atom *var
))
235 (setq symbol
(gensym))
236 (setq exp
(maxima-substitute symbol
*var exp
))
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
))))
251 (cond ((setq factors
(first-order-p exp
*var
))
253 (ratf (make-mtimes -
1 (div* (cdr factors
)
257 (t (setq varlist
(list *var
))
259 (setq varlist
(varsort varlist
))
261 (ratnumer (mrat-numer (ratrep* exp
)))
262 (numer-varlist varlist
)
263 (subst-list (trig-subst-p varlist
)))
264 (setq varlist
(ncons *var
))
266 (setq exp
(trig-subst exp subst-list
))
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
))
279 (not (equal 1 (setq *g
(solventhp (cdddr exp
) (cadr exp
))))))
281 (t (cond ($solvefactors
282 (map2c (lambda (x y
) (solve1a x
(m* mult y
)))
284 (t (solve1a exp mult
)))))))))
286 (cond (symbol (setq *roots
(subst temp
*var
*roots
))
287 (setq *failures
(subst temp
*var
*failures
))))
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
)))
309 (solve (car terms
) *var mult
))
312 ((eq (caar *exp
) 'mabs
) ;; abs(x) = 0 <=> x = 0
313 (solve (cadr *exp
) *var mult
)
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
322 ((and (integerp (caddr *exp
))
323 (plusp (caddr *exp
)))
324 (solve (cadr *exp
) *var
(m* mult
(caddr *exp
)))
327 ;;; Predicate to test for presence of troublesome trig functions to be
328 ;;; canonicalized. A table of when to make substitutions should
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
))
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
)
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
361 (defun measure (f alist set
&aux
(sum 0))
362 (dolist (element set
)
363 (incf sum
(funcall f
(cdr (assoc element alist
:test
#'eq
)))))
366 ;; Named for uniqueness only
369 (cond ((or (null x
) (= x
0)) 0)
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
)
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
))))
388 ;; For the case where the value of assoc is NIL, we will need a special "+"
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
)
401 (assl (do ((k 0 (1+ k
))
403 ((= k
(length set
)) made
)
404 (setq made
(cons (cons (nth k set
) 0)
406 ((= j
(length vlist
)) assl
)
407 (setq it
(nth j vlist
))
409 (t (setq assl
(extend #'+mset
(cons (cons (caar it
) 1) nil
)
411 (setq assl
(extend #'+mset assl
412 (operator-frequency-table (cdr it
) set
)
415 (defun trig-subst (exp sub-list
)
417 (sub-list (cdr sub-list
) (cdr sub-list
))
418 (var (car sub-list
) (car sub-list
)))
421 (maxima-substitute (funcall (trig-cannon (g-rep-operator var
))
422 (make-mlist-l (g-rep-operands var
)))
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
)
466 (cond ((setq temp
(linear-term-p exp var
))
470 (do ((arg (car (g-rep-operands exp
)) (car rest
))
471 (rest (cdr (g-rep-operands exp
)) (cdr rest
))
477 (make-lineq (make-mplus-l linear-term-list
)
478 (if constant-term-list
479 (make-mplus-l constant-term-list
)
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
))
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
)
494 (cond ((eq exp var
) 1)
496 (t (case (g-rep-operator exp
)
498 (do ((factor (car (g-rep-operands exp
)) ;individual factors
500 (rest (cdr (g-rep-operands exp
)) ;factors yet to be done
502 (main-var-p) ;nt -> main-var seen at top level
503 (list-of-factors)) ;accumulate our factors
504 ((null factor
) ;for all factors
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
511 ((broken-freeof var factor
) ;if
512 (push factor list-of-factors
))
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!
532 (defun solve1a (exp mult
)
533 (let ((*myvar
*myvar
)
535 (cond ((atom exp
) nil
)
536 ((not (memalike (setq *myvar
(simplify (pdis (list (car exp
) 1 1))))
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
)))))
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
)
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
)))
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
)
574 (re-solve eqn
*$var mult
)
575 (make-solution `(,eqn
1) nil
)))
576 (cond ((solution-losses sol
) sol
)
579 (t (do ((l (solution-wins sol
) (cddr l
)))
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
)
604 ;; We've already decomposed and factored
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
)
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
)
633 (solve (make-mequal *power
(mequal-rhs w
)) *varb z
))
634 (t (let ((rts (re-solve-full
635 (make-mequal *power
(mequal-rhs w
))
637 (map2c #'(lambda (root mult
)
638 (solve (make-mequal (mequal-rhs root
) 0)
640 (solution-wins rts
))))))
642 (map2c #'(lambda (w z
)
644 (push (solventh3 w
*power
*varb
) *failures
))
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
)))
653 (mapcar #'(lambda (h)
660 (list (car w
) (rdis (cdr w1
)) (rdis (cdr w2
))))))
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
)
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
)
678 (not (memalike (simplify (pdis (list (caaddr e
) 1 1)))
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
)
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
))
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
))))
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
))))
713 (t (displine line
))))
715 ;; Predicate to check if an expression which may be broken up
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
)
724 (do ((b-n-fo var
(car b-n-fo-l
))
725 (b-n-fo-l broken-not-freeof
(cdr b-n-fo-l
)))
727 (and (not (argsfreeof b-n-fo exp
))
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
)
741 (push (make-mequal-simp *var exp
) *roots
))
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))
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
)
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
))
776 (t (setq discrim
(simpnrt discrim
2))
777 (solve3 (fullratsimp `((mquotient) ((mplus) ,b
,discrim
) ,a
))
779 (solve3 (fullratsimp `((mquotient) ((mplus) ,b
((mminus) ,discrim
)) ,a
))
782 ;; Reorders V so that members which contain the variable of
783 ;; interest come first.
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
))
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!
804 ;; Solve <exp> = <*myvar> for <*var>, where <*myvar>=<op>(...)
805 (defun usolve (exp op
)
810 (cond ((broken-freeof *var
814 `((mplus) ((mminus) ,(caddr *myvar
))
815 ,(div* `((%log
) ,exp
)
816 `((%log
) ,(cadr *myvar
)))))
820 (cond ((mnegp (caddr *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
)
833 (eq ($asksign exp
) '$neg
))
835 (t `((mplus) ,(cadr *myvar
)
838 ,(div* 1 (caddr *myvar
))))))))
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
))
848 `((mplus) ((mminus) ,(cadr *myvar
))
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)
863 (not (free exp
'mplus
))))
867 g1
(cond ((null l
) (return nil
)))
868 (setq a
(car (setq fm l
)))
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
))
878 (defmfun $linsolve
(eql varl
)
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
)
884 (do ((varl varl
(cdr varl
)))
886 (when (mnump (car varl
))
887 (merror (intl:gettext
"solve: variable must not be a number; found: ~M") (car varl
))))
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
*)
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
)))
907 (return '((mlist simp
))))
910 ;;I put this in the value cell--wfs
911 (setf (aref xa
* 0 j
) nil
))
912 (ptorat 'xa
* xn
* xm
*)
914 (xrutout 'xa
* xn
* xm
*
915 (mapcar #'(lambda (x) (ith varl x
))
919 (setq varl
(make-mlist-l (linsort (cdr varl
) *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
)))