1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancments. ;;;;;
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 $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
25 *has
*var
*var $dontfactor
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
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.")
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))
78 ; (if (and (integerp (car l)) (zerop (car l)))
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)
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
)))))
116 ;; If the variable list wasn't supplied we have to supply it
117 ;; ourselves. Also remove constants like $%pi from the list.
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
129 (cond (($listp varl
) (remove-duplicates
130 (mapcar #'meval
(cdr 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.
145 ;; Trivially true equations for any set of variables.
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.
160 (= (length varl
) (length eql
)))
161 (setq eql
(solvex eql varl
(not $programmode
) t
))
163 (cond ((and (cdr eql
)
164 (not ($listp
(cadr 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
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
))
179 (intl:gettext
"solve: more unknowns than equations.~
180 ~%Unknowns given : ~%~M~
181 ~%Equations given: ~%~M")
184 (intl:gettext
"solve: more equations than unknowns.~
185 ~%Unknowns given : ~%~M~
186 ~%Equations given: ~%~M")
190 ;; Removes anything from its list arg which solve considers not to be a
191 ;; variable, i.e. constants, functions or subscripted variables without
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))
210 (make-mlist-l (nreverse (map2c #'(lambda (eqn mult
) (push mult multi
) eqn
)
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:~%")))
219 (when $dispflag
(mtell (intl:gettext
"solve: solution:~%")))
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
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.
249 ((equal exp
0) (allroot *var
))
251 (t (setq exp
(meqhk exp
))
252 (cond ((equal exp
'(0))
253 (return (allroot *var
)))
256 (cond ((not (atom *var
))
257 (setq symbol
(gensym))
258 (setq exp
(maxima-substitute symbol
*var exp
))
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
))))
273 (cond ((setq factors
(first-order-p exp
*var
))
275 (ratf (make-mtimes -
1 (div* (cdr factors
)
279 (t (setq varlist
(list *var
))
281 (setq varlist
(varsort varlist
))
283 (ratnumer (mrat-numer (ratrep* exp
)))
284 (numer-varlist varlist
)
285 (subst-list (trig-subst-p varlist
)))
286 (setq varlist
(ncons *var
))
288 (setq exp
(trig-subst exp subst-list
))
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
))
301 (not (equal 1 (setq *g
(solventhp (cdddr exp
) (cadr exp
))))))
303 (t (cond ($solvefactors
304 (map2c (lambda (x y
) (solve1a x
(m* mult y
)))
306 (t (solve1a exp mult
)))))))))
308 (cond (symbol (setq *roots
(subst temp
*var
*roots
))
309 (setq *failures
(subst temp
*var
*failures
))))
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
)))
331 (solve (car terms
) *var mult
))
334 ((eq (caar *exp
) 'mabs
) ;; abs(x) = 0 <=> x = 0
335 (solve (cadr *exp
) *var mult
)
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
344 ((and (integerp (caddr *exp
))
345 (plusp (caddr *exp
)))
346 (solve (cadr *exp
) *var
(m* mult
(caddr *exp
)))
349 ;;; Predicate to test for presence of troublesome trig functions to be
350 ;;; canonicalized. A table of when to make substitutions should
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
))
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
)
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
383 (defun measure (f alist set
&aux
(sum 0))
384 (dolist (element set
)
385 (incf sum
(funcall f
(cdr (assoc element alist
:test
#'eq
)))))
388 ;; Named for uniqueness only
391 (cond ((or (null x
) (= x
0)) 0)
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
)
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
))))
410 ;; For the case where the value of assoc is NIL, we will need a special "+"
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
)
423 (assl (do ((k 0 (1+ k
))
425 ((= k
(length set
)) made
)
426 (setq made
(cons (cons (nth k set
) 0)
428 ((= j
(length vlist
)) assl
)
429 (setq it
(nth j vlist
))
431 (t (setq assl
(extend #'+mset
(cons (cons (caar it
) 1) nil
)
433 (setq assl
(extend #'+mset assl
434 (operator-frequency-table (cdr it
) set
)
437 (defun trig-subst (exp sub-list
)
439 (sub-list (cdr sub-list
) (cdr sub-list
))
440 (var (car sub-list
) (car sub-list
)))
443 (maxima-substitute (funcall (trig-cannon (g-rep-operator var
))
444 (make-mlist-l (g-rep-operands var
)))
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
)
488 (cond ((setq temp
(linear-term-p exp var
))
492 (do ((arg (car (g-rep-operands exp
)) (car rest
))
493 (rest (cdr (g-rep-operands exp
)) (cdr rest
))
499 (make-lineq (make-mplus-l linear-term-list
)
500 (if constant-term-list
501 (make-mplus-l constant-term-list
)
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
))
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
)
516 (cond ((eq exp var
) 1)
518 (t (case (g-rep-operator exp
)
520 (do ((factor (car (g-rep-operands exp
)) ;individual factors
522 (rest (cdr (g-rep-operands exp
)) ;factors yet to be done
524 (main-var-p) ;nt -> main-var seen at top level
525 (list-of-factors)) ;accumulate our factors
526 ((null factor
) ;for all factors
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
533 ((broken-freeof var factor
) ;if
534 (push factor list-of-factors
))
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!
554 (defun solve1a (exp mult
)
555 (let ((*myvar
*myvar
)
557 (cond ((atom exp
) nil
)
558 ((not (memalike (setq *myvar
(simplify (pdis (list (car exp
) 1 1))))
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
)))))
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
)
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
)))
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
)
596 (re-solve eqn
*$var mult
)
597 (make-solution `(,eqn
1) nil
)))
598 (cond ((solution-losses sol
) sol
)
601 (t (do ((l (solution-wins sol
) (cddr l
)))
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
)
626 ;; We've already decomposed and factored
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
)
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
)
655 (solve (make-mequal *power
(mequal-rhs w
)) *varb z
))
656 (t (let ((rts (re-solve-full
657 (make-mequal *power
(mequal-rhs w
))
659 (map2c #'(lambda (root mult
)
660 (solve (make-mequal (mequal-rhs root
) 0)
662 (solution-wins rts
))))))
664 (map2c #'(lambda (w z
)
666 (push (solventh3 w
*power
*varb
) *failures
))
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
)))
675 (mapcar #'(lambda (h)
682 (list (car w
) (rdis (cdr w1
)) (rdis (cdr w2
))))))
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
)
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
)
700 (not (memalike (simplify (pdis (list (caaddr e
) 1 1)))
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
)
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
))
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
))))
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
))))
735 (t (displine line
))))
737 ;; Predicate to check if an expression which may be broken up
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
)
746 (do ((b-n-fo var
(car b-n-fo-l
))
747 (b-n-fo-l broken-not-freeof
(cdr b-n-fo-l
)))
749 (and (not (argsfreeof b-n-fo exp
))
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
)
763 (push (make-mequal-simp *var exp
) *roots
))
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))
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
)
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
))
798 (t (setq discrim
(simpnrt discrim
2))
799 (solve3 (fullratsimp `((mquotient) ((mplus) ,b
,discrim
) ,a
))
801 (solve3 (fullratsimp `((mquotient) ((mplus) ,b
((mminus) ,discrim
)) ,a
))
804 ;; Reorders V so that members which contain the variable of
805 ;; interest come first.
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
))
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!
826 ;; Solve <exp> = <*myvar> for <*var>, where <*myvar>=<op>(...)
827 (defun usolve (exp op
)
832 (cond ((broken-freeof *var
836 `((mplus) ((mminus) ,(caddr *myvar
))
837 ,(div* `((%log
) ,exp
)
838 `((%log
) ,(cadr *myvar
)))))
842 (cond ((mnegp (caddr *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
)
855 (eq ($asksign exp
) '$neg
))
857 (t `((mplus) ,(cadr *myvar
)
860 ,(div* 1 (caddr *myvar
))))))))
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
))
870 `((mplus) ((mminus) ,(cadr *myvar
))
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)
885 (not (free exp
'mplus
))))
889 g1
(cond ((null l
) (return nil
)))
890 (setq a
(car (setq fm l
)))
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
))
900 (defmfun $linsolve
(eql varl
)
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
)
906 (do ((varl varl
(cdr varl
)))
908 (when (mnump (car varl
))
909 (merror (intl:gettext
"solve: variable must not be a number; found: ~M") (car varl
))))
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
*)
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
)))
929 (return '((mlist simp
))))
932 ;;I put this in the value cell--wfs
933 (setf (aref xa
* 0 j
) nil
))
934 (ptorat 'xa
* xn
* xm
*)
936 (xrutout 'xa
* xn
* xm
*
937 (mapcar #'(lambda (x) (ith varl x
))
941 (setq varl
(make-mlist-l (linsort (cdr varl
) *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
)))