home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / solve.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  32.7 KB  |  965 lines
  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. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module solve)
  13.  
  14. (load-macsyma-macros ratmac strmac)
  15.  
  16. (DECLARE-TOP (GENPREFIX V_)
  17.      (SPECIAL VAR-LIST EXPSUMSPLIT $DISPFLAG $NOLABELS CHECKFACTORS *G
  18.           $ALGEBRAIC EQUATIONS ;List of E-labels
  19.           *POWER *VARB *FLG $DERIVSUBST $NUMER $FLOAT
  20.           $%EMODE WFLAG GENVAR GENPAIRS VARLIST BROKEN-NOT-FREEOF
  21.           $FACTORFLAG
  22.           MULT ;Some crock which tracks multiplicities.
  23.           *ROOTS ;alternating list of solutions and multiplicities
  24.           *FAILURES ;alternating list of equations and multiplicities
  25.           *MYVAR $LISTCONSTVARS
  26.           *HAS*VAR *VAR $DONTFACTOR $LINENUM $LINECHAR
  27.           LINELABLE $KEEPFLOAT $RATFAC
  28.           ERRRJFFLAG ;A substitute for condition binding.
  29.           LSOLVEFLAG XM* XN* MUL* SOLVEXP)
  30.      (ARRAY* (NOTYPE XA* 2))
  31.      (FIXNUM THISN $LINENUM))
  32.  
  33. (DEFMVAR $BREAKUP T
  34.      "Causes solutions to cubic and quartic equations to be expressed in
  35.      terms of common subexpressions.")
  36.  
  37. (DEFMVAR $MULTIPLICITIES '$NOT_SET_YET
  38.      "Set to a list of the multiplicities of the individual solutions
  39.      returned by SOLVE, REALROOTS, or ALLROOTS.")
  40.  
  41. (DEFMVAR $LINSOLVEWARN T
  42.      "Needs to be documented.")
  43.  
  44. (DEFMVAR $SOLVE_INCONSISTENT_ERROR T
  45.      "If T gives an MAXIMA-ERROR if SOLVE meets up with inconsistent linear
  46.      equations.  If NIL, returns ((MLIST SIMP)) in this case.")
  47.  
  48. (DEFMVAR $PROGRAMMODE T
  49.      "Causes SOLVE to return its answers explicitly as elements
  50.      in a list rather than printing E-labels.")
  51.  
  52. (DEFMVAR $SOLVEDECOMPOSES T
  53.      "Causes SOLVE to use POLYDECOMP in attempting to solve polynomials.")
  54.  
  55. (DEFMVAR $SOLVEEXPLICIT NIL
  56.      "Causes SOLVE to return implicit solutions i.e. of the form F(x)=0.")
  57.  
  58. (DEFMVAR $SOLVEFACTORS T
  59.      "If T, then SOLVE will try to factor the expression.  The FALSE
  60.      setting may be desired in zl-SOME cases where factoring is not
  61.      necessary.")
  62.  
  63. (DEFMVAR $SOLVENULLWARN T
  64.      "Causes the user will be warned if SOLVE is called with either a
  65.      null equation list or a null variable list.  For example,
  66.      SOLVE([],[]); would print two warning messages and return [].")
  67.  
  68. (DEFMVAR $SOLVETRIGWARN T
  69.      "Causes SOLVE to print a warning message when it is uses
  70.      inverse trigonometric functions to solve an equation,
  71.      thereby losing solutions.")
  72.  
  73. (DEFMVAR $SOLVERADCAN NIL
  74.      "SOLVE will use RADCAN which will make SOLVE slower but will allow
  75.      certain problems containing exponentials and logs to be solved.")
  76.  
  77. ;; Utility macros
  78.  
  79. ;; In MacLisp, this turns into SUBRCALL if we are compiling, FUNCALL if
  80. ;; interpreted.  In LMLisp and other random systems, just turn into FUNCALL.
  81.  
  82. #+MacLisp
  83. (DEFMACRO SUBR-FUNCALL (FUNCTION . ARGS)
  84.       (COND ((STATUS FEATURE COMPLR) `(SUBRCALL NIL ,FUNCTION . ,ARGS))
  85.         (T `(FUNCALL ,FUNCTION . ,ARGS))))
  86.  
  87. #-MacLisp
  88. (DEFMACRO SUBR-FUNCALL (FUNCTION . ARGS) `(FUNCALL ,FUNCTION . ,ARGS))
  89.  
  90. ;; This macro returns the number of trivial equations.  It counts up the
  91. ;; number of zeros in a list.
  92.  
  93. (DEFMACRO NZLIST (lLIST)
  94.   `(DO ((L ,lLIST (CDR L))
  95.     (ZCOUNT 0))
  96.        ((NULL L) ZCOUNT)
  97.      (IF (AND (INTEGERP (CAR L)) (ZEROP (CAR L)))
  98.      (INCREMENT ZCOUNT))))
  99.  
  100. ;; This is only called on a variable.
  101.  
  102. (DEFMACRO ALLROOT (EXP)
  103.   `(SETQ *FAILURES (LIST* (MAKE-MEQUAL-SIMP ,EXP ,EXP) 1 *FAILURES)))
  104.  
  105.  
  106. ;; Finds variables, changes equations into expressions without MEQUAL.
  107. ;; Checks for consistency between the number of unknowns and equations.
  108. ;; Calls SOLVEX for simultaneous equations and SSOLVE for a single equation.
  109.  
  110. (DEFMFUN $SOLVE (*EQL &OPTIONAL (VARL NIL VARL-P))
  111.   (SETQ $MULTIPLICITIES (MAKE-MLIST))
  112.   (PROG (EQL               ;Expressions and variables being solved
  113.      $KEEPFLOAT $RATFAC       ;In case the user has set these
  114.      *ROOTS *FAILURES       ;*roots gets solutions, *failures "roots of"
  115.      BROKEN-NOT-FREEOF)       ;Has something to do with spliting up roots
  116.    (SETQ EQL
  117.      (COND ((ATOM *EQL) (NCONS *EQL))
  118.            ((EQ (G-REP-OPERATOR *EQL) 'MLIST)
  119.         (MAPCAR 'MEQHK (MAPCAR 'MEVAL (CDR *EQL))))
  120.            ((MEMQ (G-REP-OPERATOR *EQL)
  121.               '(MNOTEQUAL MGREATERP MLESSP MGEQP MLEQP))
  122.         (MERROR "Cannot solve inequalities. -SOLVE"))
  123.            (T (NCONS (MEQHK *EQL)))))
  124.  
  125.    (COND ((NULL VARL-P)       ;If the variable list wasn't supplied
  126.       (SETQ VARL       ;we have to supply it ourselves.
  127.         (LET (($LISTCONSTVARS NIL))
  128.              (CDR ($LISTOFVARS
  129.                ;If some trivial then use original equations
  130.                ;(primarily for case of X=X etc.)
  131.                (COND ((ZEROP (NZLIST EQL)) *EQL)
  132.                  (T EQL))))))         ;Usually throw trivia out!
  133.       (IF VARL (SETQ VARL (REMC VARL))))          ;Remove all constants
  134.      (T (SETQ VARL
  135.           (COND (($LISTP VARL) (MAPCAR #'MEVAL (CDR VARL)))
  136.             (T (LIST VARL))))))
  137.  
  138.       (IF (AND (NULL VARL) $SOLVENULLWARN)
  139.       (MTELL "~&Got a null variable list, continuing - SOLVE~%"))
  140.       (IF (AND (NULL EQL) $SOLVENULLWARN)
  141.       (MTELL "~&Got a null equation list, continuing - SOLVE~%"))
  142.       (IF (ORMAPC #'MNUMP VARL)
  143.       (MERROR "A number was found where a variable was expected -SOLVE"))
  144.  
  145.       (COND ((EQUAL EQL '(0)) (RETURN '$ALL))
  146.         ((OR (NULL VARL) (NULL EQL)) (RETURN (MAKE-MLIST-SIMP)))
  147.         ((AND (NULL (CDR VARL)) (NULL (CDR EQL)))
  148.          (RETURN (SSOLVE (CAR EQL) (CAR VARL))))
  149.         ((OR VARL-P (= (LENGTH VARL) (LENGTH EQL)))
  150.          (SETQ EQL (SOLVEX EQL VARL (NOT $PROGRAMMODE) T))
  151.          (RETURN (COND ((AND (CDR EQL) (NOT ($LISTP (CADR EQL))))
  152.                 (MAKE-MLIST EQL))
  153.                (T EQL)))))
  154.     (LET ((U (MAKE-MLIST-L VARL))
  155.       (E (COND (($LISTP *EQL) *EQL)
  156.            (T (MAKE-MLIST *EQL)))))
  157.       ;; MFORMAT doesn't have ~:[~] yet, so I just change this to
  158.       ;; make one of two possible calls to MERROR. Smaller codesize
  159.       ;; then what was here before anyway.
  160.       (IF (> (LENGTH VARL) (LENGTH EQL))
  161.       (MERROR "More unknowns than equations -SOLVE~
  162.           ~%Unknowns given :  ~%~M~
  163.           ~%Equations given:  ~%~M"
  164.           U E)
  165.       (MERROR "More equations than unknowns -SOLVE~
  166.           ~%Unknowns given :  ~%~M~
  167.           ~%Equations given:  ~%~M"
  168.           U E)))))
  169.  
  170.  
  171. ;; Removes anything from its list arg which solve considers not to be a
  172. ;; variable, i.e.  constants, functions or subscripted variables without
  173. ;; numeric args.
  174.  
  175. (DEFUN REMC (LST)
  176.   (DO ((L LST (CDR L)) (FL) (VL)) ((NULL L) VL)
  177.     (COND ((ATOM (SETQ FL (CAR L)))
  178.         (OR (MAXIMA-CONSTANTP FL) (SETQ VL (CONS FL VL))))
  179.           ((ANDMAPC #'$CONSTANTP (CDR FL)) (SETQ VL (CONS FL VL))))))
  180.  
  181. ;; List of multiplicities.  Why is this special?
  182.  
  183. (DECLARE-TOP (SPECIAL MULTI)) 
  184.  
  185. ;; Solve a single equation for a single unknown.
  186. ;; Obtains roots via solve and prints them.
  187.  
  188. (DEFUN SSOLVE (EXP *VAR &AUX EQUATIONS MULTI)
  189.   (LET (($SOLVETRIGWARN $SOLVETRIGWARN))
  190.     (COND ((NULL *VAR) '$ALL)
  191.       (T (SOLVE EXP *VAR 1)
  192.          (COND ((NOT (OR *ROOTS *FAILURES)) (MAKE-MLIST))
  193.            ($PROGRAMMODE
  194.             (PROG1 (MAKE-MLIST-L
  195.                 (NREVERSE
  196.                  (MAP2C #'(LAMBDA (EQN MULT) (PUSH MULT MULTI) EQN)
  197.                     (IF $SOLVEEXPLICIT
  198.                     *ROOTS
  199.                     (NCONC *ROOTS *FAILURES)))))
  200.                (SETQ $MULTIPLICITIES
  201.                  (MAKE-MLIST-L (NREVERSE MULTI)))))
  202.            (T (WHEN (AND *FAILURES (NOT $SOLVEEXPLICIT))
  203.                 (IF $DISPFLAG (MTELL "The roots of:~%"))
  204.                 (SOLVE2 *FAILURES))
  205.               (WHEN *ROOTS
  206.                 (IF $DISPFLAG (MTELL "Solution:~%"))
  207.                 (SOLVE2 *ROOTS))
  208.               (MAKE-MLIST-L EQUATIONS)))))))
  209.  
  210. ;; Solve takes three arguments, the expression to solve for zero, the variable
  211. ;; to solve for, and what multiplicity this solution is assumed to have (from
  212. ;; higher-level Solve's).  Solve returns NIL.  Isn't that useful?  The lists
  213. ;; *roots and *failures are special variables to which Solve prepends solutions
  214. ;; and their multiplicities in that order: *roots contains explicit solutions
  215. ;; of the form <var>=<function of independent variables>, and *failures
  216. ;; contains equations which if solved would yield additional solutions.
  217.  
  218. ;; Factors expression and reduces exponents by their gcd (via solventhp)
  219.  
  220. (DEFMFUN SOLVE (*EXP *VAR MULT
  221.            &AUX (GENVAR NIL)
  222.             ($DERIVSUBST NIL)
  223.             (EXP (FLOAT2RAT (MRATCHECK *EXP)))
  224.             (*MYVAR *VAR)
  225.             ($SAVEFACTORS T))
  226.  (PROG (FACTORS *HAS*VAR GENPAIRS $DONTFACTOR TEMP SYMBOL *G CHECKFACTORS 
  227.     VARLIST EXPSUMSPLIT)
  228.      (LET (($RATFAC T)) (SETQ EXP (RATDISREP (RATF EXP))))
  229.         ; Cancel out any simple 
  230.         ; (non-algebraic) common factors in numerator and 
  231.         ; denominator without altering the structure of the 
  232.         ; expression too much.
  233.         ; Also, RJFPROB in TEST;SOLVE TEST is now solved.
  234.         ; - JPG
  235.    A (COND ((ATOM EXP)
  236.         (COND ((EQ EXP *VAR)
  237.            (SOLVE3 0 MULT))
  238.           ((EQUAL EXP 0) (ALLROOT *VAR))
  239.           (T NIL)))
  240.        (T (SETQ EXP (MEQHK EXP))
  241.           (COND ((EQUAL EXP '(0))
  242.              (RETURN (ALLROOT *VAR)))
  243.             ((FREE EXP *VAR)
  244.              (RETURN NIL)))
  245.           (COND ((NOT (ATOM *VAR))
  246.              (SETQ SYMBOL (GENSYM))
  247.              (SETQ EXP (MAXIMA-SUBSTITUTE SYMBOL *VAR EXP))
  248.              (SETQ TEMP *VAR)
  249.              (SETQ *VAR SYMBOL)
  250.              (SETQ *MYVAR *VAR)))  ;keep *MYVAR up-to-date
  251.           
  252.           (COND ($SOLVERADCAN (SETQ EXP (RADCAN1 EXP))
  253.                   (IF (ATOM EXP) (GO A))))
  254.           
  255.           (COND ((EASY-CASES EXP *VAR)
  256.              (COND (SYMBOL (SETQ *ROOTS (SUBST TEMP *VAR *ROOTS))
  257.                    (SETQ *FAILURES (SUBST TEMP *VAR *FAILURES))))
  258.              (ROOTSORT *ROOTS)
  259.              (ROOTSORT *FAILURES)
  260.              (RETURN NIL)))
  261.           
  262.           (COND ((SETQ FACTORS (FIRST-ORDER-P EXP *VAR))
  263.              (SOLVE3 (RATDISREP
  264.                    (RATF (MAKE-MTIMES -1 (DIV* (CDR FACTORS)
  265.                                (CAR FACTORS)))))
  266.                  MULT))
  267.             
  268.             (T (SETQ VARLIST (LIST *VAR))
  269.                (FNEWVAR EXP)
  270.                (SETQ VARLIST (VARSORT VARLIST))
  271.                (LET ((VARTEMP)
  272.                  (RATNUMER (MRAT-NUMER (RATREP* EXP)))
  273.                  (NUMER-VARLIST VARLIST)
  274.                  (SUBST-LIST (TRIG-SUBST-P VARLIST)))
  275.               (SETQ VARLIST (NCONS *VAR))
  276.              (COND (SUBST-LIST
  277.                  (SETQ EXP (TRIG-SUBST EXP SUBST-LIST))
  278.                  (FNEWVAR EXP)
  279.                  (SETQ VARLIST (VARSORT VARLIST))
  280.                  (SETQ EXP (MRAT-NUMER (RATREP* EXP)))
  281.                  (SETQ VARTEMP VARLIST))
  282.                    (T (SETQ VARTEMP NUMER-VARLIST)
  283.                   (SETQ EXP RATNUMER)))
  284.              (SETQ VARLIST VARTEMP))
  285.                
  286.                (COND ((ATOM EXP) (GO A))
  287.                  ((SPECASEP EXP) (SOLVE1A EXP MULT))
  288.                  ((AND (NOT (PCOEFP EXP))
  289.                    (CDDR EXP)
  290.                    (NOT (EQUAL 1 (SETQ *G
  291.                                (SOLVENTHP (CDDDR EXP) (CADR EXP))))))
  292.                   (SOLVENTH EXP *G))
  293.                  (T (MAP2C #'SOLVE1A
  294.                        (COND ($SOLVEFACTORS (PFACTOR EXP))
  295.                          (T (LIST EXP 1))))))))))
  296.  
  297.       (COND (SYMBOL (SETQ *ROOTS (SUBST TEMP *VAR *ROOTS))
  298.             (SETQ *FAILURES (SUBST TEMP *VAR *FAILURES))))
  299.       (ROOTSORT *ROOTS)
  300.       (ROOTSORT *FAILURES)
  301.       (RETURN NIL)))
  302.  
  303. (DEFUN FLOAT2RAT (EXP)
  304.  (COND ((FLOATP EXP) (SETQ EXP (PREP1 EXP)) (MAKE-RAT-SIMP (CAR EXP) (CDR EXP)))
  305.        ((OR (ATOM EXP) (SPECREPP EXP)) EXP)
  306.        (T (RECUR-APPLY #'FLOAT2RAT EXP))))
  307.  
  308. ;;; The following takes care of cases where the expression is already in 
  309. ;;; factored form. This can introduce spurious roots if one of the factors
  310. ;;; is an expression that can be undefined or infinity for certain values of
  311. ;;; the variable in question. But soon this will be no worry because I will
  312. ;;; add a list of  "possible bad roots" to what $SOLVE returns.
  313. ;;; Solve is not fully recursive when it due to globals, $MULTIPLICIES
  314. ;;; may be screwed here. (Solve should be made recursive)
  315.  
  316. (DEFUN EASY-CASES (*EXP *VAR)
  317.        (COND ((EQ (CAAR *EXP) 'MTIMES)
  318.           (DO ((TERMS (CDR *EXP) (CDR TERMS)))
  319.           ((NULL TERMS))
  320.           (SOLVE (CAR TERMS) *VAR 1))
  321.           'MTIMES)
  322.          ((EQ (CAAR *EXP) 'MEXP)
  323.           (COND ((AND (INTEGERP  (CADDR *EXP))
  324.               (PLUSP (CADDR *EXP)))
  325.              (SOLVE (CADR *EXP) *VAR (CADDR *EXP))
  326.              'MEXPRAT)))))
  327.  
  328. ;;; Predicate to test for presence of troublesome trig functions to be
  329. ;;; canonicalized.  A  table of when to make substitutions should
  330. ;;; be used here. 
  331. ;;;  trig kind                     => SIN | COS | TAN ...   subst to make
  332. ;;; number around in expression ->     1     1     0         ......
  333. ;;; what you want to be able to do for example is to see if SIN and COS^2 
  334. ;;; are around and then make a reasonable substitution.
  335.  
  336. (DEFUN TRIG-SUBST-P (VLIST)
  337.   (AND (NOT (TRIG-NOT-SUBST-P VLIST))
  338.        (DO ((VAR (CAR VLIST) (CAR VLIST))
  339.         (VLIST (CDR VLIST) (CDR VLIST))
  340.         (SUBST-LIST))
  341.        ((NULL VAR) SUBST-LIST)
  342.      (COND ((AND (NOT (ATOM VAR))
  343.              (TRIG-CANNON (G-REP-OPERATOR VAR))
  344.              (NOT (FREE VAR *VAR)))
  345.         (PUSH VAR SUBST-LIST))))))
  346.  
  347. ;; Predicate to see when obviously not to substitute for trigs.
  348. ;; A hack in the direction of expression properties-table driven
  349. ;; substition. The "measure" of the expression is the total number
  350. ;; of different kinds of trig functions in the expression.
  351.  
  352. (DEFUN TRIG-NOT-SUBST-P (VLIST)
  353.   (LET ((TRIGS '(%SIN %COS %TAN %COT %CSC %SEC)))
  354.     (< (MEASURE #'SIGN-GJC (OPERATOR-FREQUENCY-TABLE VLIST TRIGS) TRIGS)
  355.        2)))
  356.  
  357. ;; To get the total "value" of things in a table, this case an assoc list.
  358. ;; (MEASURE FUNCTION ASSOCIATION-LIST SET) where FUNCTION is a function mapping
  359. ;; the range of the ASSOCIATION-LIST viewed as a function on the SET, to the
  360. ;; integers.
  361.  
  362. (DEFUN MEASURE (F ALIST SET &AUX (SUM 0))
  363.   (DOLIST (ELEMENT SET)
  364.     (INCREMENT SUM (FUNCALL F (CDR (ASSQ ELEMENT ALIST)))))
  365.   SUM)
  366.  
  367. ;; (defun MEASURE (F AL S)
  368. ;;        (do ((j 0 (f1+ j))
  369. ;;           (sum 0))
  370. ;;       ((= j (length S))  sum)
  371. ;;       (setq sum (f+ sum (funcall F (cdr (assoc (nth j S) al)))))))
  372.  
  373. ;; Named for uniqueness only
  374.  
  375. (DEFUN SIGN-GJC (X)
  376.   (COND ((OR (NULL X) (= X 0)) 0)
  377.     ((< 0 X) 1)
  378.     (T -1)))
  379.  
  380. ;; A function that can EXTEND a function
  381. ;; over two association lists. Note that I have been using association lists
  382. ;; as mere functions (that is, as sets of ordered pairs).
  383. ;; (EXTEND '+ L1 L2 S) could also be to take the union of two multi-sets in the
  384. ;; sample space S. (what the '&%%#?& has this got to do with SOLVE?) 
  385.  
  386. (DEFUN EXTEND (F L1 L2 S)
  387.   (DO ((J 0 (f1+ J))
  388.        (VALUE NIL))
  389.       ((= J (LENGTH S)) VALUE)
  390.     (SETQ VALUE (CONS (CONS (NTH J S)
  391.                 (FUNCALL F (CDR (zl-ASSOC (NTH J S) L1))
  392.                      (CDR (zl-ASSOC (NTH J S) L2))))
  393.               VALUE))))
  394.  
  395. ;; For the case where the value of assoc is NIL, we will need a special "+"
  396.  
  397. (DEFUN +MSET (A B) (f+ (OR A 0) (OR B 0)))
  398.  
  399. ;; To recursively looks through a list
  400. ;; structure (the VLIST) for members of the SET appearing in the MACSYMA 
  401. ;; functional position (caar list). Returning an assoc. list of appearence
  402. ;; frequencies. Notice the use of EXTEND.
  403.  
  404. (DEFUN OPERATOR-FREQUENCY-TABLE (VLIST SET)
  405.   (DO ((J 0 (f1+ J))
  406.        (IT)
  407.        (ASSL (DO ((K 0 (f1+ K))
  408.           (MADE NIL))
  409.          ((= K (LENGTH SET)) MADE)
  410.            (SETQ MADE (CONS (CONS (NTH K SET) 0)
  411.                 MADE)))))
  412.       ((= J (LENGTH VLIST)) ASSL)
  413.     (SETQ IT (NTH J VLIST))
  414.     (COND ((ATOM IT))
  415.       (T (SETQ ASSL (EXTEND #'+MSET (CONS (CONS (CAAR IT) 1) NIL)
  416.                 ASSL SET))
  417.          (SETQ ASSL (EXTEND #'+MSET ASSL
  418.                 (OPERATOR-FREQUENCY-TABLE (CDR IT) SET)
  419.                 SET))))))
  420.  
  421. (DEFUN TRIG-SUBST (EXP SUB-LIST)
  422.   (DO ((EXP EXP)
  423.        (SUB-LIST (CDR SUB-LIST) (CDR SUB-LIST))
  424.        (VAR (CAR SUB-LIST) (CAR SUB-LIST)))
  425.       ((NULL VAR) EXP)
  426.     (SETQ EXP
  427.       (MAXIMA-SUBSTITUTE (FUNCALL (TRIG-CANNON (G-REP-OPERATOR VAR))
  428.                    (MAKE-MLIST-L (G-REP-OPERANDS VAR)))
  429.               VAR EXP))))
  430.  
  431. ;; Here are the canonical trig substitutions.
  432.  
  433. (DEFUN-prop (%SEC TRIG-CANNON) (X)
  434.   (INV* (MAKE-G-REP '%COS (G-REP-FIRST-OPERAND X))))
  435.  
  436. (DEFUN-prop (%CSC TRIG-CANNON) (X)
  437.   (INV* (MAKE-G-REP '%SIN (G-REP-FIRST-OPERAND X))))
  438.  
  439. (DEFUN-prop (%TAN TRIG-CANNON) (X)
  440.   (DIV* (MAKE-G-REP '%SIN (G-REP-FIRST-OPERAND X))
  441.     (MAKE-G-REP '%COS (G-REP-FIRST-OPERAND X))))
  442.  
  443. (DEFUN-prop (%COT TRIG-CANNON) (X)
  444.   (DIV* (MAKE-G-REP '%COS (G-REP-FIRST-OPERAND X))
  445.     (MAKE-G-REP '%SIN (G-REP-FIRST-OPERAND X))))
  446.  
  447. (DEFUN-prop (%SECH TRIG-CANNON) (X)
  448.   (INV* (MAKE-G-REP '%COSH (G-REP-FIRST-OPERAND X))))
  449.  
  450. (DEFUN-prop (%CSCH TRIG-CANNON) (X)
  451.   (INV* (MAKE-G-REP '%SINH (G-REP-FIRST-OPERAND X))))
  452.  
  453. (DEFUN-prop (%TANH TRIG-CANNON) (X)
  454.   (DIV* (MAKE-G-REP '%SINH (G-REP-FIRST-OPERAND X))
  455.     (MAKE-G-REP '%COSH (G-REP-FIRST-OPERAND X))))
  456.  
  457. (DEFUN-prop (%COTH TRIG-CANNON) (X)
  458.   (DIV* (MAKE-G-REP '%COSH (G-REP-FIRST-OPERAND X))
  459.     (MAKE-G-REP '%SINH (G-REP-FIRST-OPERAND X))))
  460.  
  461. ;; Predicate to replace ISLINEAR....Returns NIL if not of for A*X+B, A and B
  462. ;; freeof X, else returns (A . B)
  463.  
  464. (DEFUN FIRST-ORDER-P (EXP VAR &AUX TEMP)
  465.   ;; Expand the expression at one level, i.e. distribute products
  466.   ;; over sums, but leave exponentiations alone.
  467.   ;; (X+1)^2*(X+Y) --> X*(X+1)^2 + Y*(X+1)^2
  468.   (SETQ EXP (EXPAND1 EXP 1 1))
  469.   (COND ((ATOM EXP) NIL)
  470.     (T (CASE (G-REP-OPERATOR EXP)
  471.           (MTIMES
  472.             (COND ((SETQ TEMP (LINEAR-TERM-P EXP VAR))
  473.                (MAKE-LINEQ TEMP 0))
  474.               (T NIL)))
  475.           (MPLUS
  476.             (DO ((ARG  (CAR (G-REP-OPERANDS EXP)) (CAR REST))
  477.              (REST (CDR (G-REP-OPERANDS EXP)) (CDR REST))
  478.              (LINEAR-TERM-LIST)
  479.              (CONSTANT-TERM-LIST)
  480.              (TEMP))
  481.             ((NULL ARG)
  482.              (IF LINEAR-TERM-LIST
  483.                  (MAKE-LINEQ (MAKE-MPLUS-L LINEAR-TERM-LIST)
  484.                      (IF CONSTANT-TERM-LIST
  485.                          (MAKE-MPLUS-L CONSTANT-TERM-LIST)
  486.                          0))))
  487.               (COND ((SETQ TEMP (LINEAR-TERM-P ARG VAR))
  488.                  (PUSH TEMP LINEAR-TERM-LIST))
  489.                 ((BROKEN-FREEOF VAR ARG)
  490.                  (PUSH ARG CONSTANT-TERM-LIST))
  491.                 (T (RETURN NIL)))))
  492.           (T NIL)))))
  493.  
  494. ;; Function to test if a term from an expanded expression is a linear term
  495. ;; check and see that exactly one item in the product is the main var and
  496. ;; all others are free of the main var.  Returns NIL or a G-REP expression.
  497.  
  498. (DEFUN LINEAR-TERM-P (EXP VAR)
  499.   (COND ((ATOM EXP)
  500.      (COND ((EQ EXP VAR) 1)
  501.            (T NIL)))
  502.     (T (CASE (G-REP-OPERATOR EXP)
  503.           (MTIMES
  504.             (DO ((FACTOR (CAR (G-REP-OPERANDS EXP))  ;individual factors
  505.                  (CAR REST))
  506.              (REST (CDR (G-REP-OPERANDS EXP))    ;factors yet to be done
  507.                    (CDR REST))
  508.              (MAIN-VAR-P)  ;nt -> main-var seen at top level
  509.              (LIST-OF-FACTORS))  ;accumulate our factors
  510.             ((NULL FACTOR)  ;for all factors
  511.              (AND MAIN-VAR-P
  512.                   ;no-main-var at top level -=> not linear
  513.                   (MAKE-MTIMES-L LIST-OF-FACTORS)))
  514.               (COND ((EQ FACTOR VAR)  ;if it's our main var
  515.                  ;note it...it has to be there to be a linear term
  516.                  (SETQ MAIN-VAR-P T))
  517.                 ((BROKEN-FREEOF VAR FACTOR)  ;if 
  518.                  (PUSH FACTOR LIST-OF-FACTORS))
  519.                 (T (RETURN NIL)))))
  520.           (T NIL)))))
  521.  
  522.  
  523. ;;; DISPATCHING FUNCTION ON DEGREE OF EXPRESSION
  524. ;;; This is a crock of shit, it should be data driven and be able to
  525. ;;; dispatch to all manner of special cases that are in a table.
  526. ;;; EXP here is a polynomial in MRAT form.  All of this well-structured,
  527. ;;; intelligently-designed code works by side effect.  SOLVECUBIC
  528. ;;; takes something that looks like (G0003 3 4 1 1 0 10) as an argument
  529. ;;; and returns something like ((MEQUAL) $X ((MTIMES) ...)).  You figure
  530. ;;; out where the $X comes from.
  531.  
  532. ;;; It comes from GENVARS/VARLIST, of course.  Isn't this wonderful rational
  533. ;;; function package irrational?  If you don't know about GENVARS and
  534. ;;; VARLIST, you'd better bite the bullet and learn...everything depends
  535. ;;; on them.  The canonical example of mis-use of special variables!
  536. ;;; --RWK
  537.  
  538. (DEFUN SOLVE1A (EXP MULT) 
  539.   (LET ((*MYVAR *MYVAR)
  540.     (*G NIL)) 
  541.     (COND ((ATOM EXP) NIL)
  542.       ((NOT (MEMALIKE (SETQ *MYVAR (PDIS (LIST (CAR EXP) 1 1)))
  543.               *HAS*VAR))
  544.        NIL)
  545.       ((EQUAL (CADR EXP) 1) (SOLVELIN EXP))
  546.       ((SPECASEP EXP) (SOLVESPEC EXP T))
  547.       ((EQUAL (CADR EXP) 2) (SOLVEQUAD EXP))
  548.       ((NOT (EQUAL 1 (SETQ *G (SOLVENTHP (CDDDR EXP) (CADR EXP)))))
  549.        (SOLVENTH EXP *G))
  550.       ((EQUAL (CADR EXP) 3) (SOLVECUBIC EXP))
  551.       ((EQUAL (CADR EXP) 4) (SOLVEQUARTIC EXP))
  552.       (T (LET ((TT (SOLVE-BY-DECOMPOSITION EXP *MYVAR)))
  553.            (SETQ *FAILURES (APPEND (SOLUTION-LOSSES TT) *FAILURES))
  554.            (SETQ *ROOTS    (APPEND (SOLUTION-WINS TT) *ROOTS)))))))
  555.  
  556. (DEFUN SOLVE-SIMPLIST (LIST-OF-THINGS)
  557.   (G-REP-OPERANDS (SIMPLIFYA (MAKE-MLIST-L LIST-OF-THINGS) NIL)))
  558.  
  559. ;; The Solve-by-decomposition program returns the cons of (ROOTS . FAILURES).
  560. ;; It returns a "Solution" object, that is, a CONS with the CAR being the
  561. ;; failures and the CDR being the successes.
  562. ;; It takes a POLY as an argument and returns a SOLUTION.
  563.  
  564. (DEFUN SOLVE-BY-DECOMPOSITION (POLY *$VAR)
  565.  (LET ((DECOMP))
  566.   (COND ((OR (NOT $SOLVEDECOMPOSES)
  567.          (= (LENGTH (SETQ DECOMP (POLYDECOMP POLY (POLY-VAR POLY)))) 1))
  568.      (MAKE-SOLUTION NIL `(,(MAKE-MEQUAL 0 (PDIS POLY)) 1)))
  569.     (T (DECOMP-TRACE (MAKE-MEQUAL 0 (RDIS (CAR DECOMP)))
  570.              DECOMP
  571.              (POLY-VAR POLY) *$VAR 1)))))
  572.  
  573. ;; DECOMP-TRACE is the recursive function which maps itself down the
  574. ;; intermediate solutions until the end is reached.  If it encounters
  575. ;; non-solvable equations it stops.  It returns a SOLUTION object, that is, a
  576. ;; CONS with the CAR being the failures and the CDR being the successes.
  577.  
  578. (DEFUN DECOMP-TRACE (EQN DECOMP VAR *$VAR MULT &AUX SOL CHAIN-SOL WINS LOSSES)
  579.   (SETQ SOL (IF DECOMP
  580.         (RE-SOLVE EQN *$VAR MULT)
  581.         (MAKE-SOLUTION `(,EQN 1) NIL)))
  582.   (COND ((SOLUTION-LOSSES SOL) SOL)
  583.     ;; End test
  584.     ((NULL DECOMP) SOL)
  585.     (T (DO ((L (SOLUTION-WINS SOL) (CDDR L)))
  586.            ((NULL L))
  587.            (SETQ CHAIN-SOL
  588.              (DECOMP-CHAIN (CAR L) (CDR DECOMP) VAR *$VAR (CADR L)))
  589.            (SETQ WINS (NCONC WINS
  590.                  (COPY-TOP-LEVEL (SOLUTION-WINS CHAIN-SOL))))
  591.            (SETQ LOSSES (NCONC LOSSES
  592.                    (COPY-TOP-LEVEL (SOLUTION-LOSSES CHAIN-SOL)))))
  593.        (MAKE-SOLUTION WINS LOSSES))))
  594.  
  595. ;; Decomp-chain is the function which formats the mess for the recursive call.
  596. ;; It returns a "Solution" object, that is, a CONS with the CAR being the
  597. ;; failures and the CDR being the successes.
  598.  
  599. (DEFUN DECOMP-CHAIN (RSOL DECOMP VAR *$VAR MULT)
  600.   (LET ((SOL (SIMPLIFY (MAKE-MEQUAL (RDIS (IF DECOMP (CAR DECOMP)
  601.                  ;; Include the var itself in the decomposition
  602.                     (MAKE-MRAT-BODY (MAKE-MRAT-POLY VAR '(1 1)) 1)))
  603.                     (MEQUAL-RHS RSOL)))))
  604.     (DECOMP-TRACE SOL DECOMP VAR *$VAR MULT)))
  605.  
  606. ;; RE-SOLVE calls SOLVE recursively, returning a SOLUTION object.
  607. ;; Will not decompose or factor.
  608.  
  609. (DEFUN RE-SOLVE (EQN VAR MULT)
  610.   (LET ((*ROOTS NIL)
  611.     (*FAILURES NIL)
  612.     ;; We've already decomposed and factored
  613.     ($SOLVEDECOMPOSES)
  614.     ($SOLVEFACTORS))
  615.     (SOLVE EQN VAR MULT)
  616.     (MAKE-SOLUTION *ROOTS *FAILURES)))
  617.  
  618. ;; SOLVENTH programs test to see if the variable of interest appears 
  619. ;; to some power in all terms.  If so, a new variable is substituted for it
  620. ;; and the simpler expression solved with the multiplicity
  621. ;; adjusted accordingly.
  622. ;; SOLVENTHP returns gcd of exponents.
  623.  
  624. (DEFUN SOLVENTHP (L GCD) 
  625.   (COND ((NULL L) GCD)
  626.     ((EQUAL GCD 1) 1)
  627.     (T (SOLVENTHP (CDDR L)
  628.               (GCD (CAR L) GCD)))))
  629.  
  630. ;; Reduces exponents by their gcd.
  631.  
  632. (DEFUN SOLVENTH (EXP *G) 
  633.   (LET ((*VARB (PDIS (MAKE-MRAT-POLY (POLY-VAR EXP) '(1 1))))
  634.     (EXP   (MAKE-MRAT-POLY (POLY-VAR EXP) (SOLVENTH1 (POLY-TERMS EXP)))))
  635.     (LET* ((RTS (RE-SOLVE-FULL (PDIS EXP) *VARB))
  636.        (FAILS (SOLUTION-LOSSES RTS))
  637.        (WINS (SOLUTION-WINS RTS))
  638.        (*POWER (MAKE-MEXPT *VARB *G)))
  639.       (MAP2C #'(LAMBDA (W Z)
  640.          (COND ((ATOM *VARB)
  641.             (SOLVE (MAKE-MEQUAL *POWER (MEQUAL-RHS W)) *VARB Z))
  642.                (T (LET ((RTS (RE-SOLVE-FULL
  643.                        (MAKE-MEQUAL *POWER (MEQUAL-RHS W))
  644.                        *VARB)))
  645.                 (MAP2C #'(LAMBDA (ROOT MULT)
  646.                        (SOLVE (MAKE-MEQUAL (MEQUAL-RHS ROOT) 0)
  647.                           *MYVAR MULT))
  648.                    (SOLUTION-WINS RTS))))))
  649.          WINS)
  650.       (MAP2C #'(LAMBDA (W Z)
  651.          (PUSH Z *FAILURES)
  652.          (PUSH (SOLVENTH3 W *POWER *VARB) *FAILURES))
  653.          FAILS)
  654.       *ROOTS)))
  655.  
  656. (DEFUN SOLVENTH3 (W *POWER *VARB &AUX VARLIST GENVAR *FLG W1 W2)
  657.   (COND ((BROKEN-FREEOF *VARB W) W)
  658.     (T (SETQ W1 (RATF (CADR W)))
  659.        (SETQ W2 (RATF (CADDR W)))
  660.        (SETQ VARLIST
  661.          (MAPCAR #'(LAMBDA (H) 
  662.                  (COND (*FLG H)
  663.                    ((ALIKE1 H *VARB)
  664.                     (SETQ *FLG T)
  665.                     *POWER)
  666.                    (T H)))
  667.              VARLIST))
  668.        (LIST (CAR W) (RDIS (CDR W1)) (RDIS (CDR W2))))))
  669.  
  670. (DECLARE-TOP (MUZZLED T))
  671. (DEFUN SOLVENTH1 (L) 
  672.   (COND ((NULL L) NIL)
  673.     (T (CONS (QUOTIENT (CAR L) *G)
  674.          (CONS (CADR L) (SOLVENTH1 (CDDR L)))))))
  675. (DECLARE-TOP (MUZZLED NIL))
  676.  
  677. ;; Will decompose or factor
  678.  
  679. (DEFUN RE-SOLVE-FULL (X VAR &AUX *ROOTS *FAILURES)
  680.   (SOLVE X VAR MULT)
  681.   (MAKE-SOLUTION *ROOTS *FAILURES))
  682.  
  683. ;; Sees if expression is of the form A*F(X)^N+B.
  684.  
  685. (DEFUN SPECASEP (E)
  686.        (AND (MEMALIKE (PDIS (LIST (CAR E) 1 1)) *HAS*VAR)
  687.         (OR (ATOM (CADDR E))
  688.         (NOT (MEMALIKE (PDIS (LIST (CAADDR E) 1 1))
  689.                  *HAS*VAR)))
  690.         (OR (NULL (CDDDR E)) (EQUAL (CADDDR E) 0))))
  691.  
  692. ;; Solves the special case A*F(X)^N+B.
  693.  
  694. (DECLARE-TOP (MUZZLED T))
  695. (DEFUN SOLVESPEC (EXP $%EMODE) 
  696.   (PROG (A B C) 
  697.     (SETQ A (PDIS (CADDR EXP)))
  698.     (SETQ C (PDIS (LIST (CAR EXP) 1 1)))
  699.     (COND ((NULL (CDDDR EXP))
  700.            (RETURN (SOLVE C *VAR (TIMES (CADR EXP) MULT)))))
  701.     (SETQ B (PDIS (PMINUS (CADDDR (CDR EXP)))))
  702.     (RETURN (SOLVESPEC1 C
  703.                 (SIMPNRT (DIV* B A) (CADR EXP))
  704.                 (MAKE-RAT 1 (CADR EXP))
  705.                 (CADR EXP)))))
  706. (DECLARE-TOP (MUZZLED NIL))
  707.  
  708. (DEFUN SOLVESPEC1 (VAR ROOT N THISN) 
  709.  (DO ((THISN THISN (f1- THISN))) ((ZEROP THISN))
  710.      (SOLVE (ADD* VAR (MUL* -1 ROOT
  711.                 (POWER* '$%E (MUL* 2 '$%PI '$%I THISN N))))
  712.         *VAR MULT)))
  713.  
  714.  
  715. ;; ADISPLINE displays a line like DISPLINE, and in addition, notes that it is
  716. ;; not free of *VAR if it isn't.
  717.  
  718. (DEFUN ADISPLINE (LINE)
  719.   ;; This may be redundant, but nice if ADISPLINE gets used where not needed.
  720.   (COND ((AND $BREAKUP (NOT $PROGRAMMODE))
  721.      (LET ((LINELABEL (DISPLINE LINE)))
  722.        (COND ((BROKEN-FREEOF *VAR LINE))
  723.          (T (SETQ BROKEN-NOT-FREEOF
  724.               (CONS LINELABEL BROKEN-NOT-FREEOF))))
  725.        LINELABEL))
  726.     (T (DISPLINE LINE))))
  727.  
  728. ;; Predicate to check if an expression which may be broken up
  729. ;; is freeof
  730.  
  731. (SETQ BROKEN-NOT-FREEOF NIL)
  732.  
  733. ;; For consistency, use backwards args.
  734.  
  735. (DEFUN BROKEN-FREEOF (VAR EXP)
  736.   (COND ($BREAKUP
  737.       (DO ((B-N-FO VAR (CAR B-N-FO-L))
  738.            (B-N-FO-L BROKEN-NOT-FREEOF (CDR B-N-FO-L)))
  739.           ((NULL B-N-FO) T)
  740.         (AND (NOT (ARGSFREEOF B-N-FO EXP))
  741.          (RETURN NIL))))
  742.     (T (ARGSFREEOF VAR EXP))))
  743.  
  744. ;; Adds solutions to roots list.
  745. ;; Solves for inverse of functions (via USOLVE)
  746.  
  747. (DEFUN SOLVE3 (EXP MULT) 
  748.        (SETQ EXP (SIMPLIFY EXP))
  749.        (COND ((NOT (BROKEN-FREEOF *VAR EXP))
  750.           (PUSH MULT *FAILURES)
  751.           (PUSH (MAKE-MEQUAL-SIMP (SIMPLIFY *MYVAR) EXP) *FAILURES))
  752.          (T (COND ((EQ *MYVAR *VAR)
  753.                (PUSH MULT *ROOTS)
  754.                (PUSH (MAKE-MEQUAL-SIMP *VAR EXP) *ROOTS))
  755.               ((ATOM *MYVAR)
  756.                (PUSH MULT *FAILURES)
  757.                (PUSH (MAKE-MEQUAL-SIMP *MYVAR EXP) *FAILURES))
  758.               (T (USOLVE EXP (G-REP-OPERATOR *MYVAR)))))))
  759.  
  760.  
  761. ;; Solve a linear equation.  Argument is a polynomial in pseudo-cre form.
  762. ;; This function is called for side-effect only.
  763.  
  764. (DEFUN SOLVELIN (EXP) 
  765.        (COND ((EQUAL 0. (PTERM (CDR EXP) 0.))
  766.           (SOLVE1A (CADDR EXP) MULT)))
  767.        (SOLVE3 (RDIS (RATREDUCE (PMINUS (PTERM (CDR EXP) 0.))
  768.                 (CADDR EXP)))
  769.            MULT))
  770.  
  771. ;; Solve a quadratic equation.  Argument is a polynomial in pseudo-cre form.
  772. ;; This function is called for side-effect only.
  773. ;; The code for handling the case where the discriminant = 0 seems to never
  774. ;; be run.  Presumably, the expression is factored higher up.
  775.  
  776. (DECLARE-TOP (MUZZLED T))
  777. (DEFUN SOLVEQUAD (EXP &AUX DISCRIM A B C)
  778.    (SETQ A (CADDR EXP))
  779.    (SETQ B (PTERM (CDR EXP) 1.))
  780.    (SETQ C (PTERM (CDR EXP) 0.))
  781.    (SETQ DISCRIM (SIMPLIFY (PDIS (PPLUS (PEXPT B 2.)
  782.                        (PMINUS (PTIMES 4. (PTIMES A C)))))))
  783.    (SETQ B (PDIS (PMINUS B)))
  784.    (SETQ A (PDIS (PTIMES 2. A)))
  785.    ;; At this point, everything is back in general representation.
  786.    (let ((varlist nil)) ;;2/6/2002 RJF
  787.      (COND ((EQUAL 0. DISCRIM)
  788.      (SOLVE3 (FULLRATSIMP `((MQUOTIENT) ,B ,A))
  789.          (TIMES 2. MULT)))
  790.     (T (SETQ DISCRIM (SIMPNRT DISCRIM 2.))
  791.        (SOLVE3 (FULLRATSIMP `((MQUOTIENT) ((MPLUS) ,B ,DISCRIM) ,A))
  792.            MULT)
  793.        (SOLVE3 (FULLRATSIMP `((MQUOTIENT) ((MPLUS) ,B ((MMINUS) ,DISCRIM)) ,A))
  794.            MULT)))))
  795. (DECLARE-TOP (MUZZLED NIL))
  796.  
  797. ;; Reorders V so that members which contain the variable of
  798. ;; interest come first.
  799.  
  800. (DEFUN VARSORT (V)
  801.   (LET ((*U NIL)
  802.     (*V (COPY-TOP-LEVEL V)))
  803.      (MAPC #'(LAMBDA (Z) 
  804.            (COND ((BROKEN-FREEOF *VAR Z)
  805.               (SETQ *U (CONS Z *U))
  806.               (SETQ *V (zl-DELETE Z *V 1)))))
  807.        V)
  808.      (SETQ $DONTFACTOR *U)
  809.      (SETQ *HAS*VAR *V)
  810.      (APPEND *U *V)))
  811.  
  812. ;; Solves for variable when it occurs within a function by taking the inverse.
  813. ;; When this code is fixed, the `((mplus) ,x ,y) forms should be rewritten as
  814. ;; (MAKE-MPLUS X Y).  I didn't do this because the code was buggy and it should
  815. ;; be fixed first.  - cwh
  816. ;; You mean you didn't do it because you were buggy.  Hope you're fixed soon!
  817. ;; --RWK
  818.  
  819. (DEFUN USOLVE (EXP OP) 
  820.    (PROG (INVERSE) 
  821.      (SETQ INVERSE
  822.            (COND
  823.         ((EQ OP 'MEXPT)
  824.          (COND ((BROKEN-FREEOF *VAR
  825.                 (CADR *MYVAR))
  826.             (COND ((EQUAL EXP 0)
  827.                    (GO FAIL)))
  828.             `((mplus) ((mminus) ,(caddr *myvar))
  829.                   ,(div* `((%log) ,exp)
  830.                      `((%log) ,(cadr *myvar)))))
  831.                ((BROKEN-FREEOF *VAR
  832.                 (CADDR *MYVAR))
  833.             (COND ((EQUAL EXP 0)
  834.                    (COND ((MNEGP (CADDR *MYVAR))
  835.                       (GO FAIL))
  836.                      (T (CADR *MYVAR))))
  837.                   ;; There is a bug right here.
  838.                   ;; SOLVE(SQRT(U)+1) should return U=1
  839.                   ;; This code is entered with EXP = -1, OP = MEXPT
  840.                   ;; *VAR = U, and *MYVAR = ((MEXPT) U ((RAT) 1 2))
  841.                   ;; BULLSHIT -- RWK.  That is precisely the bug
  842.                   ;; this code was added to fix!
  843.                   ((and (not (eq (ask-integer (caddr *myvar)
  844.                               '$INTEGER)
  845.                          '$yes))
  846.                     (free exp '$%i)
  847.                     (eq ($asksign exp) '$neg))
  848.                    (go fail))                    
  849.                   (T `((mplus) ,(cadr *myvar)
  850.                        ((mminus)
  851.                         ((mexpt) ,exp
  852.                              ,(div* 1 (caddr *myvar))))))))
  853.                (T (GO FAIL))))
  854.         ((SETQ INVERSE (GET OP '$INVERSE))
  855.          (WHEN (AND $SOLVETRIGWARN
  856.                 (MEMQ OP '(%SIN %COS %TAN %SEC
  857.                        %CSC %COT %COSH %SECH)))
  858.            (MTELL "~&SOLVE is using arc-trig functions to get ~
  859.                 a solution.~%Some solutions will be lost.~%")
  860.            (SETQ $SOLVETRIGWARN NIL))
  861.          `((MPLUS) ((MMINUS) ,(CADR *MYVAR))
  862.                ((,INVERSE) ,EXP)))
  863.         ((EQ OP '%LOG)
  864.          `((MPLUS) ((MMINUS) ,(CADR *MYVAR))
  865.                ((MEXPT) $%E ,EXP)))
  866.         (T (GO FAIL))))
  867.      (RETURN (SOLVE (SIMPLIFY INVERSE) *VAR MULT))
  868.     FAIL (RETURN (SETQ *FAILURES
  869.                (CONS (SIMPLIFY `((MEQUAL) ,*MYVAR ,EXP))
  870.                  (CONS MULT *FAILURES))))))
  871.  
  872. ;; Predicate for determining if an expression is messy enough to 
  873. ;; generate a new linelabel for it.
  874. ;; Expression must be in general form.
  875.  
  876. (DEFUN COMPLICATED (EXP)
  877.        (AND $BREAKUP
  878.         (NOT $PROGRAMMODE)
  879.         (NOT (FREE EXP 'MPLUS))))
  880.  
  881. (DECLARE-TOP (MUZZLED T))
  882. (DEFUN ROOTSORT (L) 
  883.        (PROG (A FM FM1) 
  884.     G1   (COND ((NULL L) (RETURN NIL)))
  885.          (SETQ A (CAR (SETQ FM L)))
  886.          (SETQ FM1 (CDR FM))
  887.     LOOP (COND ((NULL (CDDR FM)) (SETQ L (CDDR L)) (GO G1))
  888.            ((ALIKE1 (CADDR FM) A)
  889.             (RPLACA FM1 (PLUS (CAR FM1) (CADDDR FM)))
  890.             (RPLACD (CDR FM) (CDDDDR FM))
  891.             (GO LOOP)))
  892.          (SETQ FM (CDDR FM))
  893.          (GO LOOP)))
  894. (DECLARE-TOP (MUZZLED NIL))
  895.  
  896. ;; Stuff moving in from MAT to get it out of core.
  897.  
  898. (DEFMFUN $LINSOLVE (EQL VARL)
  899.   (LET (($RATFAC))
  900.     (SETQ EQL (COND (($LISTP EQL) (CDR EQL))
  901.             (T (NCONS EQL))))
  902.     (SETQ VARL (COND (($LISTP VARL) (REMRED (CDR VARL)))
  903.              (T (NCONS VARL))))
  904.     (DO ((VARL VARL (CDR VARL))) ((NULL VARL))
  905.     (COND ((MNUMP (CAR VARL))
  906.            (MERROR "Unacceptable variable to SOLVE: ~M"
  907.                (CAR VARL)
  908.                ))))
  909.     (COND ((NULL VARL) (MAKE-MLIST-SIMP))
  910.       (T (SOLVEX (MAPCAR 'MEQHK EQL) VARL (NOT $PROGRAMMODE) NIL)))
  911.     ))
  912.  
  913. ;; REMRED removes any repetition that may be in the variables list
  914. ;; The NREVERSE is significant here for some reason?
  915.  
  916. (DEFUN REMRED (L) (IF L (NREVERSE (UNION1 L NIL))))
  917.  
  918. (DEFUN SOLVEX (EQL VARL IND FLAG &AUX ($ALGEBRAIC $ALGEBRAIC))
  919.   (declare (special xa*))
  920.   (PROG (*VARL ANS VARLIST GENVAR LSOLVEFLAG XM* XN* MUL* SOLVEXP)
  921.     (SETQ *VARL VARL)
  922.     (SETQ SOLVEXP FLAG)
  923.     (SETQ LSOLVEFLAG T)
  924.     (SETQ EQL
  925.           (MAPCAR #'(LAMBDA (X) ($RATDISREP ($RATNUMER X)))
  926.               EQL))
  927.     (COND ((ATOM (LET ((ERRRJFFLAG  T))
  928.                (CATCH 'RATERR (FORMX FLAG 'XA* EQL VARL))))
  929.            ;; This flag is T if called from SOLVE
  930.            ;; and NIL if called from LINSOLVE.
  931.            (COND (FLAG (RETURN ($ALGSYS (MAKE-MLIST-L EQL)
  932.                         (MAKE-MLIST-L VARL))))
  933.              (T (MERROR "LINSOLVE ran into a nonlinear equation.")))))
  934.     (SETQ ANS (TFGELI 'XA* XN* XM*))
  935.     (IF (AND $LINSOLVEWARN (CAR ANS))
  936.         (MTELL "~&Dependent equations eliminated:  ~A~%" (CAR ANS)))
  937.     (IF (CADR ANS)
  938.         (IF $SOLVE_INCONSISTENT_ERROR
  939.         (MERROR "Inconsistent equations:  ~A" (CADR ANS))
  940.         (RETURN '((MLIST SIMP)))))
  941.     (DO ((J 0 (f1+ J)))
  942.         ((> J XM*))
  943.       ;;I put this in the value cell--wfs 
  944.        ; (STORE ( XA* 0 J) NIL))
  945.       (STORE (ARRAYCALL T XA* 0 J) NIL))
  946.     (PTORAT 'XA* XN* XM*)
  947.     (SETQ VARL
  948.           (XRUTOUT 'XA* XN* XM* 
  949.                (MAPCAR #'(LAMBDA (X) (ITH VARL X))
  950.                    (CADDR ANS))
  951.                IND))
  952.     (*REARRAY 'XA*)
  953.     (IF $PROGRAMMODE
  954.         (SETQ VARL (MAKE-MLIST-L (LINSORT (CDR VARL) *VARL))))
  955.     (RETURN VARL)))
  956.  
  957. ;; (LINSORT '(((MEQUAL) A2 FOO) ((MEQUAL) A3 BAR)) '(A3 A2))
  958. ;; returns (((MEQUAL) A3 BAR) ((MEQUAL) A2 FOO)) .
  959.  
  960. (DEFUN LINSORT (MEQ-LIST VAR-LIST)
  961.        (MAPCAR #'(LAMBDA (X) (CONS (CAAR MEQ-LIST) X))
  962.            (SORTCAR (MAPCAR #'CDR MEQ-LIST)
  963.             #'(LAMBDA (X Y)
  964.                 (zl-MEMBER Y (zl-MEMBER X VAR-LIST))))))
  965.