reduce-historical/r30/solve.red

COMMENT SOLVE MODULE;

%******************* Global Declarations ***************************;

SYMBOLIC;

FLAG('(!*SOLVEWRITE), 'SHARE);

ARRAY !!CF(12), !!INTERVAL(10,2), !!EXACT(10);

GLOBAL '(!!HIPOW !!GCD !*SOLVESINGULAR SM!* MP!* !*ALLBRANCH
         !*SOLVEWRITE !!ARBINT !*SOLVEINTERVAL !!INTERVALARRAY);

!*SOLVESINGULAR := T;  % Solves consistent, singular eqns (0=0) if T;
!*ALLBRANCH     := T;  % Returns all branches of solutions if T;
!*SOLVEWRITE    := T;  % Prints solutions if T;
%!*SOLVEINTERVAL = NIL;% Attempts to isolate insoluble, real roots if T;

!!INTERVALARRAY := '!!INTERVAL;  % Value is the name of an array used to
				 %   pass args to RealRoot routines;
!!ARBINT    := 0;                % Index for arbitrary constants;

%  !!HIPOW : SOLVECOEFF returns highest power of its arg in this
%  !!GCD   : SOLVECOEFF returns GCD of powers of its arg in this
%  !!CF    : Array of coeffs from SOLVECOEFF
%
%  SM!*      : List of solutions
%  MP!*      : List of multiplicities;

ALGEBRAIC MATRIX SOLN, MULTIPLICITY;

%******************* Utility Functions *****************************;

SYMBOLIC PROCEDURE RPLACX U;
BEGIN SCALAR CARU;
  CARU := CAR U;
  RETURN RPLACD(RPLACA(U,CDR U),CARU)
END;

SYMBOLIC PROCEDURE UNIVARIATEP F;
  % F is a standard form.  Non-nil iff F is univariate or a constant;
DOMAINP F OR
(DOMAINP LC F AND (DOMAINP RED F OR
                   ((MVAR F = MVAR RED F) AND UNIVARIATEP RED F) ));

SYMBOLIC SMACRO PROCEDURE SUBTRSQ(U,V);
   ADDSQ(U, NEGSQ V);

SYMBOLIC SMACRO PROCEDURE VARLIS U;
   %U is an r-polynomial.
   %value is an ordered list of variables in U;
   VARLIS1(U,NIL);

SYMBOLIC SMACRO PROCEDURE LFCTR U;
   COMMENT RETURNS LEFTFACTOR OF A PAIR.  USED BY
      SUMFACTORS IN IEQN.RED;
   CAAR U;

SYMBOLIC OPERATOR LCMD;

SYMBOLIC PROCEDURE LCMD(C,D);
   COMMENT C and D are prefix rational numbers.  Returns
      integer least-common-multiple of their denominators;
   LCM(DENR SIMP!* C, DENR SIMP!* D);

SYMBOLIC PROCEDURE VARLIS1(U,V);
   IF DOMAINP U THEN V
    ELSE VARLIS1(CDR U,VARLIS1(CDAR U,ORDAS(CAAAR U,V)));

SYMBOLIC PROCEDURE ORDAS(A,L);
   IF NULL L THEN LIST A
    ELSE IF A=CAR L THEN L
    ELSE IF ORDP(A,CAR L) THEN A . L
    ELSE CAR L . ORDAS(A,CDR L);

SYMBOLIC PROCEDURE RATNUMP X;
   COMMENT Returns T iff any prefix expression X is a rational
      number;
   ATOM NUMR(X:=SIMP!* X) AND ATOM DENR X;

FLAG ('(RATNUMP), 'DIRECT);

SYMBOLIC PROCEDURE KARGLIS(KNAME, KLIS);
   COMMENT KNAME evaluates to an atom and KLIS to a list of
      kernels.  Returns the list of kernels named KNAME in KLIS;
   IF NULL KLIS THEN NIL
   ELSE UNION(KARG1(KNAME, CAR KLIS), KARGLIS(KNAME,CDR KLIS));

SYMBOLIC PROCEDURE KARG1(KNAME, KRN);
   COMMENT KNAME evaluates to an atom and KRN to a kernel.
      Returns a list of kernels named KNAME in KRN;
   IF ATOM KRN THEN NIL
   ELSE IF CAR KRN=KNAME THEN UNION(KARGLIS(KNAME,CDR KRN),
      LIST(KRN))
   ELSE KARGLIS(KNAME, CDR KRN);

SYMBOLIC PROCEDURE ALLKERN ELST;
   COMMENT Returns list of all top-level kernels in the list of
      standard forms ELST;
   IF NULL ELST THEN NIL
   ELSE UNION(VARLIS CAR ELST, ALLKERN CDR ELST);

SYMBOLIC OPERATOR FREEOFKERN;

SYMBOLIC PROCEDURE FREEOFKERN(X,U);
   COMMENT Returns T iff any expression U is free of kernel X;
   IF ATOM X THEN FREEOF(U,X)
   ELSE FREEOF(SUBST('!!DUM,X,U),'!!DUM);

FLAG('(FREEOFKERN),'DIRECT);

SYMBOLIC PROCEDURE TOPKERN(EX, X);
   BEGIN COMMENT Returns list of toplevel kernels in the
     standard form EX that contain the kernel X;
   SCALAR ALLK, WITHX;
   ALLK := VARLIS EX;
   WHILE  ALLK DO<<
      IF NOT FREEOFKERN(X,CAR ALLK) THEN WITHX:=CAR ALLK.WITHX;
      ALLK:=CDR ALLK>>;
   RETURN WITHX
   END;

SYMBOLIC PROCEDURE COEFLIS(EX);
% EX is a standard form.
% Returns a list of the coefficients of the main variable
%   in ex in the form ((expon.coeff) (expon.coeff) ... ),
%   where the expon's occur in increasing order, and entries
%   do not occur of zero coefficients;
   BEGIN
      SCALAR X, ANS, OLDKORD, VAR;
      X := EX;
      IF DOMAINP(X) THEN
         RETURN (0 . X);
      VAR := MVAR(EX);
      WHILE (NOT DOMAINP(X)) AND MVAR(X)=VAR DO <<
         ANS := (LDEG(X) . LC(X)) . ANS;
         X := RED(X) >>;
      IF X THEN
         ANS := (0 . X) . ANS;
      RETURN ANS
   END;

%******************* Temporary Factoring Routine *******************;

% The following square free factoring routine, based on the Reduce
%   function SQFRF, will eventually be replaced by the Norman-Moore
%   complete factorization technique.;

FLUID '(!*GCD);

SYMBOLIC PROCEDURE FACTLIS(EX, KLIST);
% EX is a standard form.
% KLIST is a list of kernels.
% Returns a list of square free factors containing the elements of
% KLIST in the form ((integer exponent . standard form factor) ...).;
% Factors constant with respect to KLIST are discarded;
BEGIN
   SCALAR  FIRST, ANS, OLDGCD, OLDKORD; INTEGER EXPON;
   OLDGCD := !*GCD;                     
   !*GCD  := T;               % Must be on for SQFRF;
   OLDKORD := SETKORDER(KLIST);         
   EX := REORDER(EX);                   
   WHILE (NOT DOMAINP(EX)) AND (MVAR(EX) MEMBER KLIST) DO <<
      FIRST := PP(EX);
      IF NOT DOMAINP(FIRST) THEN <<
         % Non-zero roots;
         EX  := QUOTF(EX, FIRST);
         FIRST := SQFRF(FIRST);
         FOR EACH X IN FIRST DO
            IF NOT DOMAINP X THEN
               ANS :=  RPLACX X . ANS >>
      ELSE <<
         % Zero root (possibly multiple);
         ANS := (LDEG(EX) . !*K2F(MVAR(EX))) . ANS;
         EX  := QUOTF(EX, !*P2F(LPOW(EX))) >> >>;
   % Restore the state of the world;
   SETKORDER(OLDKORD);
   !*GCD  := OLDGCD;
   RETURN ANS
END;

%******************* SOLVE Statement ******************************;

SYMBOLIC PROCEDURE SIMPSOLVE ARGLIST;
    BEGIN
        INTEGER NARGS;
        NARGS := LENGTH(ARGLIST);       
        RETURN !*F2Q IF NARGS=1 THEN SOLVE0(CAR ARGLIST,NIL)
		      ELSE IF NARGS=2
		       THEN SOLVE0(CAR ARGLIST, CADR ARGLIST)
		      ELSE SOLVE0(CAR ARGLIST,'LST . CDR ARGLIST)
    END;

PUT ('SOLVE,'SIMPFN,'SIMPSOLVE);

%******************* Fundamental SOLVE Procedures ******************;

SYMBOLIC PROCEDURE SOLVE0(ELST, XLST);

   BEGIN COMMENT ELST is any prefix expression, including the
      kernel named LST with any number of arguments.  XLST is
      a kernel, perhaps named LST with any number of arguments.
      Solves eqns in ELST for vars in XLST, putting solutions
      and multiplicities in SOLN and MULTIPLICITIES.
      Prints SOLN if !*SOLVEWRITE is non-nil.
      Returns number of rows in global matrix SOLN;
   SCALAR FLST, VARS, NONLIN;  INTEGER NEQN, I;
   %/ MAYBELOADMATR();
   ALGEBRAIC CLEAR SOLN, MULTIPLICITY;
   SM!* := MP!* := NIL;
   IF NOT ATOM ELST  AND CAR ELST='LST THEN ELST:=CDR ELST
   ELSE ELST:=LIST ELST;
   NEQN:=0;
   WHILE  ELST DO <<FLST:= NUMR SIMP!* CAR ELST . FLST;
      NEQN:=NEQN+1;  ELST:= CDR ELST >>;
% Note that ELST and XLST are reversed from the order entered;
   IF NULL XLST THEN <<VARS := ALLKERN FLST;
	 WRITE "UNKNOWNS:";
	 MAPCAR(REVERSE VARS, FUNCTION MATHPRINT);
	 TERPRI()>>
   ELSE<<IF ATOM XLST OR NOT(CAR XLST='LST)THEN XLST:=LIST(XLST)
         ELSE XLST:=CDR XLST;
         WHILE  XLST DO<<
            VARS:=MVAR !*A2F CAR XLST.VARS;
	    XLST:= CDR XLST>>>>;
   IF NOT(NEQN=LENGTH VARS) THEN REDERR
    "SOLVE CALLED WITH UNEQUAL NUMBER OF EXPRESSIONS AND UNKNOWNS";
   IF NEQN=1 THEN
      IF NULL (FLST:=CAR FLST) THEN
        IF !*SOLVESINGULAR THEN <<!!ARBINT:=!!ARBINT+1;
	   CONSSMMP(SIMP!* LIST('ARBCOMPLEX,!!ARBINT), 1) >>
        ELSE RETURN 0
      ELSE <<VARS:=CAR VARS;
         SOLVE1(FLST./1, VARS, 1) >>
   COMMENT More than one equation;
   ELSE <<
      SM!* := TP1(SOLVESYS(FLST, VARS));
      MP!* := LIST(LIST(MK!*SQ(!*F2Q(1)))) >>;
   SM!* := MAPC2(SM!*, FUNCTION MK!*SQ);
   PUT('MULTIPLICITY, 'MATRIX, 'MAT . MP!*);
   PUT('SOLN, 'MATRIX, 'MAT . SM!*);
   IF !*SOLVEWRITE THEN
      MATPRI(SM!*, 'SOLN);
   RETURN LENGTH SM!*
   END;

SYMBOLIC PROCEDURE CONSSMMP(S, M);
   BEGIN COMMENT S is a standard quotient and M is an integer.
      Conses (S) to global variable SM!* and (M) to global
      variable MP!*;
   SM!* := LIST(S) . SM!*;
   MP!* := LIST(MK!*SQ(M./1)) . MP!*
   END;

SYMBOLIC PROCEDURE SOLVEF(F, V);
% F is a standard form, V is a kernel.  Returns a list of
% pairs, each of which car is a standard quotient and cdr an
% integer.  If the integer is positive, the SQ is a zero of
% F with multiplicity equal to the integer.  Otherwise it is
% an insoluble factor, with multiplicity the absolute value of
% the integer;
BEGIN SCALAR OLDSOLVEWRITE, ANS;
   OLDSOLVEWRITE := !*SOLVEWRITE;
   !*SOLVEWRITE := NIL;
   SOLVE0(MK!*SQ(!*F2Q(F)), V);
   ANS := PAIR(MAPCAR(SM!*, FUNCTION LAMBDA(X); SIMP!*(CAR(X))),
	       MAPCAR(MP!*, FUNCTION CAR) );
   !*SOLVEWRITE := OLDSOLVEWRITE;
   RETURN ANS
END;

%******************* Procedures for solving a single eqn ***********;

SYMBOLIC PROCEDURE SOLVE1 (EX, X, MUL);
   BEGIN COMMENT Factors standard quotient EX with respect to
      toplevel occurrences of X and kernels containing variable
      X.  Factors containing more than one such kernel
      are appended to SM!*, with a negative multiplicity
      indicating unsolvability, and SOLVE2 is applied
      to the other factors.  Integer MUL is the multiplicity
      passed from any previous factorizations.  Returns NIL;
   SCALAR E1, X1, TKLIST;  INTEGER MU;
   EX := NUMR EX;
   TKLIST := TOPKERN(EX,X);
   IF NULL TKLIST THEN RETURN NIL;
   EX := FACTLIS(EX, TKLIST);
   WHILE EX DO <<
      E1 := CDAR(EX);
      X1 := TOPKERN(E1, X);
      MU := MUL*CAAR EX;
      IF  X1 THEN
         IF NULL CDR X1 THEN
           SOLVE2(E1,CAR X1,X,MU)
	 ELSE IF SMEMQ('SOL,
            (X1:=SIMP!* LIST('SOL,MK!*SQ(E1./1), X))) THEN
	       CONSSMMP(E1./1, -MU)
         ELSE
           SOLVE1(X1,X,MU);
      EX := CDR(EX) >>
   END;

SYMBOLIC PROCEDURE SOLVE2(E1, X1, X, MU);
  BEGIN COMMENT E1 is a standard form, MU is an
      integer, X1 and X are kernels. Uses roots of unity, known
      inverses, together with quadratic, cubic and quartic
      formulas, treating other cases as unsolvable. Returns NIL;
  SCALAR B, C, D, F;  INTEGER N;
  F:= ERRORSET(SOLVECOEFF(E1, X1),NIL,NIL);
  N:= !!GCD;

  COMMENT Test for single power of X1;
  IF ATOM(F) THEN CONSSMMP(E1./1, -MU)
  ELSE IF (F:=CAR F)=-1 THEN <<
    B:= LIST('EXPT, MK!*SQ QUOTSQ(NEGSQ SIMP!* GETELV(LIST('!!CF,0)),
      SIMP!* GETELV(LIST('!!CF,1))), MK!*SQ(1 ./!!GCD));
    FOR K := 0:N-1 DO <<
      SETELV(LIST('!!CF,1), SIMP!* LIST('TIMES,B,
        MKEXP LIST('QUOTIENT,LIST('TIMES,K,2,'PI),N)));

      COMMENT  x = b;
      IF X1=X THEN CONSSMMP(GETELV(LIST('!!CF, 1)), MU)

      COMMENT  LOG(x) = b;
      ELSE IF CAR X1 = 'LOG THEN SOLVE1           
         (SUBTRSQ(SIMP!* CADR X1,SIMP!* LIST('EXPT,'E,MK!*SQ
         GETELV(LIST('!!CF, 1)))),X,MU)

      ELSE IF CAR X1 = 'EXPT THEN

        COMMENT c**(...) = b;
	IF FREEOF(CADR X1,X) THEN <<
          IF !*ALLBRANCH THEN <<!!ARBINT:=!!ARBINT+1;
            C:=LIST('TIMES,2,'I,'PI,LIST('ARBINT,!!ARBINT)) >>
          ELSE C:=0;
          SOLVE1(SUBTRSQ(SIMP!* CADDR X1,QUOTSQ(ADDSQ(
	    SIMP!* LIST('LOG,MK!*SQ GETELV(LIST('!!CF, 1))),SIMP!* C),
	    SIMP!* LIST('LOG,CADR X1))),X,MU) >>

	ELSE IF FREEOF(CADDR X1,X) THEN

          COMMENT  (...)**(m/n) = b;
          IF RATNUMP CADDR X1 THEN SOLVE1(SUBTRSQ(
	    EXPTSQ(SIMP!* CADR X1,NUMR SIMP!* CADDR X1),
            SIMP!* LIST('EXPT,MK!*SQ GETELV(LIST('!!CF, 1)),MK!*SQ(DENR
            SIMP!* CADDR X1./1))),X,MU)

          COMMENT (...)**c = b;
          ELSE <<
            IF !*ALLBRANCH THEN <<!!ARBINT:=!!ARBINT+1;
              C:=MKEXP LIST('TIMES,LIST
                ('ARBREAL,!!ARBINT)) >>
            ELSE C:=1;
            SOLVE1(SUBTRSQ(SIMP!* CADR X1,MULTSQ(SIMP!*
	      LIST('EXPT,MK!*SQ GETELV(LIST('!!CF, 1)), MK!*SQ INVSQ
	      SIMP!* CADDR X1),SIMP!* C)), X, MU) >>

        COMMENT (...)**(...) = b : transcendental;
	ELSE CONSSMMP(SUBTRSQ(SIMP!* X1,GETELV(LIST('!!CF, 1))), -MU)

      COMMENT SIN(...) = b;
      ELSE IF CAR X1='SIN THEN<<
        IF !*ALLBRANCH THEN <<
          !!ARBINT:=!!ARBINT+1;
          F:=LIST('TIMES,2,'PI,LIST('ARBINT,!!ARBINT)) >>
        ELSE
          F:=0;
        C:=SIMP!* CADR X1;
        D:=LIST('ASIN,MK!*SQ GETELV(LIST('!!CF, 1)));
        SOLVE1(SUBTRSQ(C,SIMP!* LIST('PLUS,D,F)),X,MU);
        IF !*ALLBRANCH THEN SOLVE1(SUBTRSQ(C,SIMP!* LIST
          ('PLUS,'PI,MK!*SQ
          SUBTRSQ(SIMP!* F,SIMP!* D))), X, MU) >>

      COMMENT COS(...) = b;
      ELSE IF CAR X1='COS THEN<<
        IF !*ALLBRANCH THEN<<!!ARBINT:=!!ARBINT+1;
              C:=LIST('TIMES,2,'PI,LIST('ARBINT,!!ARBINT))>>
        ELSE C:=0;
        C:=SUBTRSQ(SIMP!* CADR X1,SIMP!* C);
        D:=SIMP!* LIST('ACOS,MK!*SQ GETELV(LIST('!!CF,1)));
        SOLVE1(SUBTRSQ(C,D), X, MU);
        IF !*ALLBRANCH THEN SOLVE1(ADDSQ(C,D), X, MU) >>   

      COMMENT Unknown inverse;
      ELSE IF NULL(F:=GET(CAR X1,'INVERSE))THEN
	CONSSMMP(SUBTRSQ(SIMP!* X1,GETELV(LIST('!!CF,1))), -MU)

      COMMENT Other, known inverse;
      ELSE SOLVE1(SUBTRSQ(SIMP!* CADR X1,SIMP!*
        LIST(F,MK!*SQ GETELV(LIST('!!CF,1)))), X, MU)>> >>      

  COMMENT Test for 2 powers of X1;
  ELSE IF F>=0 THEN <<
      D:= SIMP!* GETELV(LIST('!!CF,2));
      C := QUOTSQ(SIMP!* GETELV(LIST('!!CF,0)),D);
      D := QUOTSQ(SIMP!* GETELV(LIST('!!CF,1)),MULTSQ((2 ./1),D));
      C:=SIMP!* LIST('EXPT, MK!*SQ SUBTRSQ(EXPTSQ(D,2),C),
        MK!*SQ(1 ./2));
      D := ADDSQ(D, EXPTSQ(SIMP!* X1, N));
      SOLVE1(SUBTRSQ(D,C), X, MU);
      SOLVE1(ADDSQ(D,C), X, MU) >>
  ELSE SOLVE22(E1,X1,X,MU)
 END;

SYMBOLIC PROCEDURE SOLVE22(E1,X1,X,MU);
   BEGIN SCALAR B,C,D,F; INTEGER N;
    COMMENT Test for reciprocal equation, cubic, or quartic;
      F:=(!!HIPOW+1)/2;  D:=EXPTSQ(SIMP!* X1,N);
      IF (FOR J:=0:F DO IF NOT(GETELV(LIST('!!CF,J))
                               =GETELV(LIST('!!CF,!!HIPOW-J)) )
                        THEN RETURN T)
        THEN IF (FOR J:=0:F DO IF  NUMR ADDSQ(SIMP!*
          GETELV(LIST('!!CF,J)), SIMP!* GETELV(LIST('!!CF,!!HIPOW-J)))
             THEN RETURN T)
          THEN IF !!HIPOW=3 THEN SOLVECUBIC(D,X,MU,T)
            ELSE IF !!HIPOW=4 THEN SOLVEQUARTIC(D,X,MU)
              ELSE IF !*SOLVEINTERVAL AND UNIVARIATEP E1 THEN
                     SOLVEINTERVAL(E1,MU)
		ELSE CONSSMMP(E1./1, -MU)

        COMMENT Antisymmetric reciprocal equation;
        ELSE <<  C:=ADDSQ(D,(-1 ./1));
          SOLVE1(C, X, MU);
          E1:= QUOTSQ(E1./1, C);
          IF F+F = !!HIPOW THEN <<C:=ADDSQ(D,(1 ./1));
            SOLVE1(C, X, MU);
            E1:= QUOTSQ(E1, C) >>;
          SOLVE1(E1, X, MU) >>

      COMMENT Symmetric reciprocal equation;
      ELSE IF F+F=!!HIPOW+1 THEN <<
          C:=ADDSQ(D, 1 ./1);
          SOLVE1(C,X,MU);
          SOLVE1(QUOTSQ(E1./1, C), X, MU) >>
        ELSE <<
	  B:=SM!*;
          SETELV(LIST('!!CF, 0), SIMP!* 2);
          SETELV(LIST('!!CF, 1), SIMP!* '!!X);
          C:=ADDSQ(MULTSQ(SIMP!* GETELV(LIST('!!CF,F+1)),
			  GETELV(LIST('!!CF,1))),
			  SIMP!* GETELV(LIST('!!CF,F)));
          FOR J:=2:F DO <<
	    SETELV(LIST('!!CF, J),
		   SUBTRSQ(MULTSQ(GETELV(LIST('!!CF,1)),
				  GETELV(LIST('!!CF,J-1))),
			   GETELV(LIST('!!CF,J-2))));
            C:=ADDSQ(C,MULTSQ(GETELV(LIST('!!CF,J)),
                              SIMP!* GETELV(LIST('!!CF,F+J)) )) >>;
          SOLVE1(C,'!!X,MU);  C:=F:=NIL;
	  WHILE NOT(SM!*=B) DO <<
	    C:=CAR SM!* . C;
	    F:=CAR MP!* . F;
	    SM!*:=CDR SM!*;
	    MP!*:=CDR MP!* >>;
          WHILE  C DO <<
            SOLVE1(ADDSQ(1 ./1,MULTSQ(D,SUBTRSQ(D,CAAR C))),
               X, !*A2F CAAR F*MU);
	    C:=CDR C >>  >>
  END;

SYMBOLIC PROCEDURE MKEXP U;
   (LAMBDA X;
      LIST('PLUS,LIST('COS,X),LIST('TIMES,'I,LIST('SIN,X))))
   REVAL U;

SYMBOLIC PROCEDURE SOLVECOEFF(EX, VAR);
% EX is a standard form.
% VAR is a kernel.
% Puts the coefficients (as prefix standard quotients) of
%    VAR in EX into the elements of !!CF, with index equal
%    to the exponent divided by the GCD of all the
%    exponents.  This GCD is put into !!GCD, and the
%    highest power divided by the GCD is put into
%    !!HIPOW.
% Returns the lowest power if the highest is equal to 2;
%    -1 if the highest power is less than 2, or -1 if
%    the highest power is greater than 2.
% This bizarre behaviour stems from the rewriting of the
%    Reduce COEFF function originally used by SOLVE.
%    Hopefully this will be rewritten someday without
%    the kludginess.
% Note that !!CF (an array), !!GCD, and !!HIPOW are globals.;
BEGIN
   SCALAR CLIST, X, OLDKORD;
   OLDKORD := SETKORDER(LIST(VAR));
   CLIST := REORDER (EX);
   SETKORDER(OLDKORD);
   !!HIPOW := LDEG(CLIST);
   CLIST := COEFLIS(CLIST);
   !!GCD := 0;
   X := CLIST;
   WHILE X DO <<
      !!GCD := GCDN(CAAR(X), !!GCD);
      X := CDR(X) >>;
   X := CLIST;
   FOR I := 0:(CAR(DIMENSION('!!CF))-1) DO
      SETELV(LIST('!!CF, I), NIL);
   WHILE X DO <<
      SETELV(LIST('!!CF, CAAR(X)/!!GCD), MK!*SQ(CDAR(X) ./ 1));
      X := CDR(X) >>;
   !!HIPOW := !!HIPOW/!!GCD;
   IF !!HIPOW=2 THEN
      RETURN CAAR(CLIST)/!!GCD
   ELSE IF !!HIPOW<2 THEN
      RETURN -1
   ELSE
      RETURN -2
END;

SYMBOLIC PROCEDURE SOLVEINTERVAL(EX, MUL);
% EX is a standard form,  MUL is an integer.   Isolates
% insoluble, real roots  of EX  in rational  intervals,
% stuffing result in the form  INTERVL(Lowlim,Highlim)
% into SM!* with multiplicity MUL put into MP!*.;
BEGIN  INTEGER I;
  REALROOT(PREPF EX,PREPSQ !*K2Q MVAR EX,!!INTERVALARRAY,'!!EXACT);
  FOR I := 1:GETELV LIST('!!EXACT,0) DO
    CONSSMMP(SIMP!* GETELV LIST('!!EXACT,I), MUL);
  FOR I := 1:GETELV LIST(!!INTERVALARRAY,0,0) DO
    CONSSMMP(SIMP!* LIST('INTERVL,
                         GETELV LIST(!!INTERVALARRAY,I,1),
                         GETELV LIST(!!INTERVALARRAY,I,2) ),
             MUL)
END;

SYMBOLIC PROCEDURE REALROOT(U,V,W,X);
   REDERR("Real root finding not yet implemented");


%***************** Procedures for solving Cubic and Quartic eqns ***;

SYMBOLIC PROCEDURE SOLVECUBIC(X1, X, MU, CUBE3) ;
   BEGIN COMMENT Solves !!CF(3)*X1**3 + !!CF(2)*X1**2 ...
      X1 and X are
      kernels, M and MU are integers, CUBE3 is T or NIL.
      Returns NIL;
   SCALAR A,B,C,D;
   D:=SIMP!* GETELV(LIST('!!CF,3));
   C:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,2)),D);
   B:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,1)),D);
   A:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,0)),D);
   A:=MULTSQ(ADDSQ(MULTSQ((9 ./1),MULTSQ(C,B)), ADDSQ(MULTSQ
      ((-27 ./1),A),MULTSQ((-2 ./1),EXPTSQ(C,3)))),(1 ./54));
   B := MULTSQ((-1 ./9),ADDSQ(EXPTSQ(C,2),MULTSQ((-3 ./1),B)));
   D := SIMP!* LIST('EXPT, MK!*SQ ADDSQ(EXPTSQ(B,3),
      EXPTSQ(A,2)), MK!*SQ(1 ./2));
   D := SIMP!* LIST('EXPT, MK!*SQ ADDSQ(A,D),MK!*SQ(1 ./3));
   A := NEGSQ QUOTSQ(B,D);
   B := ADDSQ(D, A);
   C := ADDSQ(X1, MULTSQ(C,(1 ./3)));
   SOLVE1(SUBTRSQ(C,B), X, MU);
   IF CUBE3 THEN <<C := ADDSQ(MULTSQ(B,(1 ./2)), C);
      D := MULTSQ(SIMP!* LIST('EXPT,MK!*SQ(-3 ./4),MK!*SQ
         (1 ./2)), SUBTRSQ(D,A));
      SOLVE1(ADDSQ(C,D), X, MU);
      SOLVE1(SUBTRSQ(C,D), X, MU)>>
   END;

SYMBOLIC PROCEDURE SOLVEQUARTIC(X1,X,MU) ;
   BEGIN COMMENT Solves !!CF(4)*X1**4 + !!CF(3)*X1**3 + ....
      X1 is a standard quotient, X is a kernel, MU is an integer,
      CUBE3 is T or NIL.  Returns NIL;
   SCALAR A,B,C,D,F;
   F:=SIMP!* GETELV(LIST('!!CF,4));
   A:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,0)),F);
   B:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,1)),F);
   C:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,2)),F);
   D:=QUOTSQ(SIMP!* GETELV(LIST('!!CF,3)),F);
   F := ADDSQ(EXPTSQ(D,2), MULTSQ((-4 ./1),C));
   SETELV(LIST('!!CF, 0), MK!*SQ NEGSQ ADDSQ(EXPTSQ(B,2),MULTSQ(A,F)));
   SETELV(LIST('!!CF, 1), MK!*SQ ADDSQ(MULTSQ(B,D),MULTSQ((-4 ./1),A)));
   SETELV(LIST('!!CF, 2), MK!*SQ NEGSQ C);
   SETELV(LIST('!!CF, 3), 1);
   SOLVECUBIC(SIMP!* X, X, MU, NIL);
   B := CAAR SM!*;
   SM!* := CDR SM!*;
   MP!*:= CDR MP!*;
   A := SIMP!* LIST('EXPT, MK!*SQ ADDSQ(EXPTSQ(B,2),MULTSQ(A,
      (-4 ./1))), MK!*SQ(1 ./2));
   F := SIMP!* LIST('EXPT, MK!*SQ ADDSQ(F,MULTSQ(B,(4 ./1))),
      MK!*SQ(1 ./2));
   SOLVE1(ADDSQ(EXPTSQ(X1,2),MULTSQ((1 ./2),ADDSQ(MULTSQ(X1,ADDSQ
      (D,F)), ADDSQ(B,A)))), X, MU);
   SOLVE1(ADDSQ(EXPTSQ(X1,2),MULTSQ((1 ./2),ADDSQ(MULTSQ(X1,
      SUBTRSQ(D,F)), SUBTRSQ(B,A)))), X, MU);
   END;

%******************* Procedures for solving a system of eqns *******;

SYMBOLIC PROCEDURE SOLVESYS(EXLIST,VARLIST);
% EXLIST is a list of standard forms.
% VARLIST is a list of kernels.
% If EXLIST and VARLIST are of the same length and the
%   elements of VARLIST are linear in the elements of
%   exlist, and further the system of linear eqns so
%   defined is non-singular, then SOLVESYS returns a
%   list of standard quotients which are solutions of
%   the system, ordered as in VARLIST.
% Otherwise an error results.;
BEGIN
   SCALAR MTRX, RHS;    % Coeffs and right side of system;
   SCALAR ROW, OLDKORD;
   IF LENGTH(EXLIST) NEQ LENGTH(VARLIST) THEN
      REDERR "SOLVESYS given unequal number of eqns & unknowns";
   OLDKORD := SETKORDER(VARLIST);
   EXLIST := MAPCAR(EXLIST, 'REORDER);
   FOR EACH EX IN EXLIST DO <<
      ROW := NIL;
      FOR EACH VAR IN VARLIST DO<<
         IF NOT DOMAINP(EX) AND
            (MVAR(EX)=VAR AND LDEG(EX)>1
             OR (NOT FREEOFKERN(VAR, LC(EX)))
             OR (NOT FREEOFKERN(VAR, RED(EX))) ) THEN
               REDERR
       "SOLVE given system of non linear-fractional equations";
         IF NOT DOMAINP(EX) AND MVAR(EX)=VAR THEN <<
            ROW := !*F2Q(LC(EX)) . ROW;
            EX := RED(EX) >>
         ELSE
            ROW := !*F2Q(NIL) . ROW >>;
      RHS := LIST(!*F2Q(NEGF(EX))) . RHS;
      MTRX := ROW . MTRX >>;
   SETKORDER(OLDKORD);
   RETURN SOLVELNRSOLVE(MTRX, RHS)
END;

SYMBOLIC PROCEDURE SOLVELNRSOLVE(U,V);
% U is a matrix canonical form, V a compatible matrix form.
% Result is the solution,y, to the matrix equation U*y=V.
% If !*SOLVESINGULAR is non-nil, introduces arbitrary constants
% if necessary.  Returns an error if the system represented is
% inconsistent or if !*SOLVESINGULAR is nil and U is singular.;
   BEGIN INTEGER N, K; SCALAR X,!*S!*, PERM;
        X := !*EXP; !*EXP := T; N := LENGTH U; PERM := INDEXLIS(1, N);
        U := CAR NORMMAT AUGMENT(U,V);
        IF NOT !*SOLVESINGULAR THEN
           U := BAREISS U
        ELSE <<
           U := SOLVEBAREISS(U, PERM);
           IF U THEN
              U := INSERTARBCONSTS(CDR(U),
                                   CAR(U)+1,
                                   FUNCTION MAKEARBCOMPLEX) >>;
        !*S!* := BACKSUB(U,N);
        U := MAPC2(RHSIDE(CAR !*S!*,N),
                   FUNCTION (LAMBDA J; CANCEL(J . CDR !*S!*)));
        !*EXP := X;
        RETURN PERMUTE(U, PERM);
   END;

SYMBOLIC PROCEDURE SOLVEBAREISS(U, PERM);
  %The 2-step integer preserving elimination method of Bareiss
  %based on the implementation of Lipson;
  %This is based on the Bareiss function in the Reduce matrix package,
  %modified to reduce singular matrices.  If PERM is nil, behaves
  %as BAREISS, except a pair is returned for non-singular U, whose
  %cdr is the triangularized U.  The car is the rank of U, which in
  %this case is always LENGTH(U).
  %Otherwise PERM is a list of the integers 1,2...length(U).
  %As columns are interchanged, then so are the elements of PERM.
  %In this case a pair is returned whose car is the rank of U and
  %whose cdr is the triangularized U. Note that, just as in BAREISS, the
  %lower triangular portion of the returned matrix standard form is only
  %implicitly all nils--the requisite RPLACAs are not performed.  Also,
  %if PERM is non-nil and the rank,r,  of U is less than the order of U,
  %only the first r rows of the upper triangular portion are explicitly
  %set.  The all nil rows are only implicitly all nils.
  %U is a list of lists of standard forms (a matrix standard form)
  %corresponding to the appropriate augmented matrix.
  %If the value of procedure is NIL then U is singular, otherwise the
  %value is the triangularized form of U (in the same form);
  BEGIN SCALAR AA,C0,CI1,CI2,IK1,IJ,KK1,KJ,K1J,K1K1,
	       UI,U1,X,K1COL,KIJ,FLG;
        INTEGER K,K1,COL,MAXCOL;
        %U1 points to K-1th row of U
        %UI points to Ith row of U
        %IJ points to U(I,J)
        %K1J points to U(K-1,J)
        %KJ points to U(K,J)
        %IK1 points to U(I,K-1)
        %KK1 points to U(K,K-1)
        %K1K1 points to U(K-1,K-1)
        %M in comments is number of rows in U
        %N in comments is number of columns in U;
        MAXCOL := LENGTH(U);
        AA:= 1;
        K:= 2;
        K1:=1;
        U1:=U;
        GO TO PIVOT;
   AGN: U1 := CDR U1;
        IF NULL CDR U1 OR NULL CDDR U1 THEN
           IF PERM AND CDR(U1) AND
              NULL(CAR(IJ := PNTH(NTH(U, MAXCOL), MAXCOL))) THEN <<
                 MAPC(CDR(IJ), FUNCTION LAMBDA(X);
                               IF X THEN RETURN NIL);
                 RETURN (MAXCOL-1).U >>
           ELSE
              RETURN MAXCOL.U;
        AA:=NTH(CAR U1,K);              %AA := U(K,K);
        K:=K+2;
        K1:=K-1;
        U1:=CDR U1;
   PIVOT:  %pivot algorithm;
        COL := K1;
        K1J:= K1K1 := PNTH(CAR U1,K1);
  PIV1: K1COL := PNTH(CAR(U1), COL);
        IF CAR K1COL THEN GO TO L2;
        UI:= CDR U1;                    %I := K;
   L:   IF NULL UI THEN
           IF PERM THEN
              IF COL>=MAXCOL THEN
                 RETURN (K1-1).U
              ELSE <<
                 COL := COL+1;
                 GO TO PIV1 >>
           ELSE
              RETURN NIL
        ELSE IF NULL CAR(IJ := PNTH(CAR UI,COL))
          THEN GO TO L1;
   L0:  IF NULL IJ THEN GO TO L2;
        X := CAR IJ;
        RPLACA(IJ,NEGF CAR K1COL);
        RPLACA(K1COL,X);
        IJ:= CDR IJ;
        K1COL:= CDR K1COL;
        GO TO L0;
   L1:  UI:= CDR UI;
        GO TO L;
   L2:  SWAPCOLUMNS(U, K1, COL, PERM);
        COL := K;
  PIV2: UI:= CDR U1;                    %I:= K;
   L21: IF NULL UI THEN
           IF PERM THEN
              IF COL>=MAXCOL THEN <<
                 FLG := T;
                 WHILE FLG AND U1 DO <<
                    IK1 := PNTH(CAR(U1), K1);
                    IJ := PNTH(IK1, MAXCOL-K1+2);
                    KIJ := PNTH(K1K1, MAXCOL-K1+2);
                    WHILE FLG AND IJ DO
                       IF ADDF(MULTF(CAR(K1K1), CAR(IJ)),
                               MULTF(CAR(IK1), NEGF(CAR(KIJ))) )
                       THEN FLG := NIL
                       ELSE IJ := CDR(IJ);
                    U1 := CDR(U1) >>;
                 IF FLG THEN
                    RETURN (K-1).U
                 ELSE
                    RETURN NIL >>
              ELSE <<
                 COL := COL+1;
                 GO TO PIV2 >>
           ELSE
              RETURN NIL;
        IJ:= PNTH(CAR UI,K1);
        C0:= ADDF(MULTF(CAR K1K1,NTH(IJ, COL-K+2)),
                  MULTF(NTH(K1K1, COL-K+2),NEGF CAR IJ));
        IF C0 THEN GO TO L3;
        UI:= CDR UI;                    %I:= I+1;
        GO TO L21;
   L3:  SWAPCOLUMNS(U, K, COL, PERM);
        C0:= QUOTF!*(C0,AA);
        KK1 := KJ := PNTH(CADR U1,K1);  %KK1 := U(K,K-1);
        IF CDR U1 AND NULL CDDR U1 THEN GO TO EV0
         ELSE IF UI EQ CDR U1 THEN GO TO COMP;
   L31: IF NULL IJ THEN GO TO COMP;     %IF I>N THEN GO TO COMP;
        X:= CAR IJ;
        RPLACA(IJ,NEGF CAR KJ);
        RPLACA(KJ,X);
        IJ:= CDR IJ;
        KJ:= CDR KJ;
        GO TO L31;
        %pivoting complete;
    COMP:
        IF NULL CDR U1 THEN GO TO EV;
        UI:= CDDR U1;                   %I:= K+1;
    COMP1:
        IF NULL UI THEN GO TO EV;       %IF I>M THEN GO TO EV;
        IK1:= PNTH(CAR UI,K1);
        CI1:= QUOTF!*(ADDF(MULTF(CADR K1K1,CAR IK1),
                           MULTF(CAR K1K1,NEGF CADR IK1)),
                     AA);
        CI2:= QUOTF!*(ADDF(MULTF(CAR KK1,CADR IK1),
                           MULTF(CADR KK1,NEGF CAR IK1)),
                     AA);
        IF NULL CDDR K1K1 THEN GO TO COMP3;%IF J>N THEN GO TO COMP3;
        IJ:= CDDR IK1;                  %J:= K+1;
        KJ:= CDDR KK1;
        K1J:= CDDR K1K1;
    COMP2:
        IF NULL IJ THEN GO TO COMP3;
        RPLACA(IJ,QUOTF!*(ADDF(MULTF(CAR IJ,C0),
                               ADDF(MULTF(CAR KJ,CI1),
                                  MULTF(CAR K1J,CI2))),
                     AA));
        IJ:= CDR IJ;
        KJ:= CDR KJ;
        K1J:= CDR K1J;
        GO TO COMP2;
    COMP3:
        UI:= CDR UI;
        GO TO COMP1;
    EV0:IF NULL C0 THEN RETURN;
    EV: KJ := CDR KK1;
        X := CDDR K1K1;                 %X := U(K-1,K+1);
        RPLACA(KJ,C0);
    EV1:KJ:= CDR KJ;
        IF NULL KJ THEN GO TO AGN;
        RPLACA(KJ,QUOTF!*(ADDF(MULTF(CAR K1K1,CAR KJ),
                               MULTF(CAR KK1,NEGF CAR X)),
                     AA));
        X := CDR X;
        GO TO EV1
   END;

SYMBOLIC PROCEDURE SWAPCOLUMNS(MATRX, COL1, COL2, PERM);
IF COL1=COL2 THEN
   MATRX
ELSE <<
   SWAPELEMENTS(PERM, COL1, COL2);
   FOR EACH U IN MATRX DO
      SWAPELEMENTS(U, COL1, COL2);
   MATRX >>;

SYMBOLIC PROCEDURE SWAPELEMENTS(LST, I, J);
% Swaps the  Ith and Jth elements of the list LST al la
%  RPLACA and returns nil.;
BEGIN SCALAR TEMP;
   IF I>J THEN <<
      TEMP := I;
      I := J;
      J := TEMP >>;
   LST := PNTH(LST, I);
   I := J-I+1;
   TEMP := NTH(LST, I);
   RPLACA(PNTH(LST, I), CAR(LST));
   RPLACA(LST, TEMP)
END;

SYMBOLIC PROCEDURE INDEXLIS(M, N);
% M,N are integers.  Returns the list (M M+1 M+2 ... N-1 N);
IF M<=N THEN M . INDEXLIS(M+1,N);

SYMBOLIC PROCEDURE INSERTARBCONSTS(M, ZEROROW, ARBFN);
% M is a matrix standard form, representing a
% matrix which has been row reduced.  All elements below
% the principal diagonal are implicitly nil, as are all
% elements in row ZEROROW and below.  It is such a form
% as is returned by SOLVEBAREISS with a non-nil second
% argument.  It inserts approriate arbitrary constants in
% the inhomogeneous portion, and 1's on the main diagonal
% except for the last row, which gets the new determinant
% of the square submatrix.  Calls ARBFN to generate arbitrary
% constants.;
BEGIN SCALAR U, V, NEWDET; INTEGER N;
   N := LENGTH(M);
   IF ZEROROW<=N THEN <<
      NEWDET := 1;
      U := M;
      FOR I := 1:(ZEROROW-1) DO <<
         NEWDET := MULTF(NEWDET, NTH(CAR(U), I));
         U := CDR(U) >>;
      FOR I := ZEROROW:(N-1) DO <<
         V := PNTH(CAR(U), I);
         RPLACA(V, 1);
         V := CDR(V);
         FOR J := I+1:N DO <<
            RPLACA(V, NIL);
            V := CDR(V) >>;
         WHILE V DO <<
	    RPLACA(V, !*K2F EVAL LIST ARBFN);
            V := CDR(V) >>;
         U := CDR(U) >>;
      V := PNTH(CAR(U), N);
      RPLACA(V, NEWDET);
      V := CDR(V);
      WHILE V DO <<
	 RPLACA(V, MULTF(NEWDET, !*K2F EVAL LIST ARBFN));
         V := CDR(V) >> >>;
   RETURN M
END;

SYMBOLIC PROCEDURE PERMUTE(U, V);
% U is a list.  V is a list of the numbers 1,2,...LENGTH(U), permuted;
% Returns a constructed list of the elements of U permuted by V.;
IF V THEN NCONC(LIST(NTH(U,CAR(V))), PERMUTE(U, CDR(V)));
   
SYMBOLIC PROCEDURE MAKEARBCOMPLEX();
BEGIN SCALAR ANS;
   ANS := NUMR(SIMP!*(LIST('ARBCOMPLEX, !!ARBINT)));
   !!ARBINT := !!ARBINT+1;
   RETURN ANS
END;

%******** Algebraic Let Statements and related declarations ********;

PUT('ASIN, 'INVERSE, 'SIN);
PUT('ACOS, 'INVERSE, 'COS);

ALGEBRAIC <<

OPERATOR SOL, INTERVL, ARBCOMPLEX, ARBREAL, ARBINT, LST;

COMMENT Supply missing argument and simplify 1/4 roots of unity;
LET   E**(I*PI/2) = I,
      E**(I*PI) = -1,
      E**(3*I*PI/2)=-I;

FOR ALL N SUCH THAT FIXP N
   LET COS((N*PI)/2)= 0;

LET COS(PI/2)=0;

FOR ALL N SUCH THAT FIXP N
   LET SIN((N*PI)/2)=
	IF REMAINDER(ABS N,4)<2 THEN 1 ELSE -1;

LET SIN(PI/2)=1;

FOR ALL N SUCH THAT FIXP N
   LET COS((N*PI)/3)=
	(IF N=4 OR REMAINDER(ABS N+2,6)>3 THEN -1 ELSE 1)/2;

LET COS(PI/3)=1/2;

FOR ALL N SUCH THAT FIXP N
   LET SIN((N*PI)/3)=
	(IF REMAINDER(ABS N,6)<3 THEN 1 ELSE -1)*SQRT(3)/2;

LET SIN(PI/3)=SQRT(3)/2;

FOR ALL N SUCH THAT FIXP N
   LET COS((N*PI)/4)=
       (IF REMAINDER(ABS N+2,8)<4 THEN 1 ELSE -1)*SQRT(2)/2;

LET COS(PI/4)=SQRT 2/2;

FOR ALL N SUCH THAT FIXP N
   LET SIN((N*PI)/4)=
	(IF REMAINDER(ABS N,8)<4 THEN 1 ELSE -1)*SQRT(2)/2;

LET SIN(PI/4)=SQRT(2)/2;

FOR ALL N SUCH THAT FIXP N
   LET COS((N*PI)/6)=

      (IF REMAINDER(ABS N+2,12)<6 THEN 1 ELSE -1)*SQRT(3)/2;

LET COS(PI/6)=SQRT 3/2;

FOR ALL N SUCH THAT FIXP N
   LET SIN((N*PI)/6)=
	(IF REMAINDER(ABS N,12)<6 THEN 1 ELSE -1)/2;

LET SIN(PI/6)=1/2;

COMMENT Rules for reducing the number of distinct kernels in an
   equation;

FOR ALL A,B,X SUCH THAT RATNUMP C AND RATNUMP D LET
   SOL(A**C-B**D, X) = A**(C*LCMD(C,D)) - B**(D*LCMD(C,D));

FOR ALL A,B,C,D,X SUCH THAT FREEOFKERN(X,A) AND FREEOFKERN(X,C) LET
   SOL(A**B-C**D, X) = E**(B*LOG A - D*LOG C);

FOR ALL A,B,C,D,X SUCH THAT FREEOFKERN(X,A) AND FREEOFKERN(X,C) LET
   SOL(A*LOG B + C*LOG D, X) = B**A*D**C - 1,
   SOL(A*LOG B - C*LOG D, X) = B**A - D**C;

FOR ALL A,B,C,D,F,X SUCH THAT FREEOFKERN(X,A) AND FREEOFKERN(X,C) LET
   SOL(A*LOG B + C*LOG D + F, X) = SOL(LOG(B**A*D**C) + F, X),
   SOL(A*LOG B + C*LOG D - F, X) = SOL(LOG(B**A*D**C) - F, X),
   SOL(A*LOG B - C*LOG D + F, X) = SOL(LOG(B**A/D**C) + F, X),
   SOL(A*LOG B - C*LOG D - F, X) = SOL(LOG(B**A/D**C) - F, X);

FOR ALL A,B,D,F,X SUCH THAT FREEOFKERN(X,A) LET
   SOL(A*LOG B + LOG D + F, X) = SOL(LOG(B**A*D) + F, X),
   SOL(A*LOG B + LOG D - F, X) = SOL(LOG(B**A*D) - F, X),
   SOL(A*LOG B - LOG D + F, X) = SOL(LOG(B**A/D) + F, X),
   SOL(A*LOG B - LOG D - F, X) = SOL(LOG(B**A/D) - F, X),
   SOL(LOG D - A*LOG B + F, X) = SOL(LOG(D/B**A) + F, X),
   SOL(LOG D - A*LOG B - F, X) = SOL(LOG(D/B**A) - F, X);

FOR ALL A,B,D,X SUCH THAT FREEOFKERN(X,A) LET
   SOL(A*LOG B + LOG D, X) = B**A*D - 1,
   SOL(A*LOG B - LOG D, X) = B**A - D,
   SOL(LOG D - A*LOG B, X) = D - B**A;

FOR ALL A,B,C,X LET
   SOL(LOG A + LOG B + C, X) = SOL(LOG(A*B) + C, X),
   SOL(LOG A - LOG B + C, X) = SOL(LOG(A/B) + C, X),
   SOL(LOG A + LOG B - C, X) = SOL(LOG(A*B) - C, X),
   SOL(LOG A - LOG B - C, X) = SOL(LOG(A/B) - C, X);

FOR ALL A,C,X SUCH THAT FREEOFKERN(X,C) LET
   SOL(LOG A + C, X) = A - E**C,
   SOL(LOG A - C, X) = A - E**(-C);

FOR ALL A,B,X LET
   SOL(LOG A + LOG B, X) = A*B - 1,
   SOL(LOG A - LOG B, X) = A - B,
   SOL(COS A - SIN B, X) = SOL(COS A - COS(PI/2-B), X),
   SOL(SIN A + COS B, X) = SOL(SIN A - SIN(B-PI/2), X),
   SOL(SIN A - COS B, X) = SOL(SIN A - SIN(PI/2-B), X),
   SOL(SIN A + SIN B, X) = SOL(SIN A - SIN(-B), X),
   SOL(SIN A - SIN B, X) = IF !*ALLBRANCH THEN SIN((A-B)/2)*
       COS((A+B)/2)  ELSE A-B,
   SOL(COS A + COS B, X) = IF !*ALLBRANCH THEN COS((A+B)/2)*
       COS((A-B)/2)  ELSE A+B,
   SOL(COS A - COS B, X) = IF !*ALLBRANCH THEN SIN((A+B)/2)*
       SIN((A-B)/2)  ELSE A-B,
   SOL(ASIN A - ASIN B, X) = A-B,
   SOL(ASIN A + ASIN B, X) = A+B,
   SOL(ACOS A - ACOS B, X) = A-B,
   SOL(ACOS A + ACOS B, X) = A+B;

LET COS(PI/2)=0>>;


END;