reduce-historical/r30/int.red

COMMENT REDUCE INTEGRATION PACKAGE WITHOUT ALGEBRAIC EXTENSIONS;

COMMENT Messages look better if one does OFF RAISE;

OFF ECHO;

SYMBOLIC;
 
FLAG('(INTERR),'TRANSFER);   %For the compiler;

COMMENT SMACRO's needed to support Cambridge LISP constructs;

SMACRO PROCEDURE EVENP X; REMAINDER(X,2)=0;

SMACRO PROCEDURE GCD(U,V); GCDN(U,V);

INFIX IEQUAL;

SYMBOLIC SMACRO PROCEDURE U IEQUAL V; EQN(U,V);

SMACRO PROCEDURE READCLOCK; TIME();

SMACRO PROCEDURE REVERSEWOC U; REVERSIP U;

SMACRO PROCEDURE SUPERPRINT U; PRETTYPRINT U;

%the next two are needed since arguments may not be numbers;

SMACRO PROCEDURE ONEP U; U=1;

SMACRO PROCEDURE ZEROP U; U=0;

COMMENT The following three smacros can be used if there is a reason
for not using actual vectors;

%SMACRO PROCEDURE MKVECT N; %MKNILL(N+1);

%SMACRO PROCEDURE PUTV(U,N,V); %CAR RPLACA(PNTH(U,N+1),V);

%SMACRO PROCEDURE GETV(U,N); %NTH(U,N+1);

COMMENT End of Cambridge LISP compatibility section;

FLUID '(LORDER SILLIESLIST VARLIST);

GLOBAL '(GENSYMCOUNT);

SYMBOLIC SMACRO PROCEDURE !*F2POL U;
   %U is a standard form;
   %Value is a polynomial form after power substitutions made;
   %If a quotient results from substitutions, an error occurs;
   !*Q2F SUBS2F U;

SYMBOLIC SMACRO PROCEDURE !*MULTF!*(U,V); MULTF(U,V);

SYMBOLIC PROCEDURE FLATTEN U;
   IF NULL U THEN NIL
    ELSE IF ATOM U THEN LIST U
    ELSE IF ATOM CAR U THEN CAR U . FLATTEN CDR U
    ELSE NCONC(FLATTEN CAR U,FLATTEN CDR U);

SYMBOLIC PROCEDURE GENSYM1 U;
    << GENSYMCOUNT:=GENSYMCOUNT+1;
       COMPRESS APPEND(EXPLODE U,EXPLODE GENSYMCOUNT) >>;
SYMBOLIC SMACRO PROCEDURE PRINTC X; PRIN2T X;

SYMBOLIC PROCEDURE MKNILL N;
   IF N=0 THEN NIL ELSE NIL . MKNILL(N-1);

SYMBOLIC PROCEDURE SQRT N;
% return sqrt of n if same is exact, or something non-numeric
% otherwise;
    IF NOT NUMBERP N THEN 'NONNUMERIC
    ELSE IF N<0 THEN 'NEGATIVE
    ELSE IF FLOATP N THEN SQRT!-FLOAT N
    ELSE IF N<2 THEN N
    ELSE NR(N,(N+1)/2);

SYMBOLIC PROCEDURE NR(N,ROOT);
% root is an overestimate here. nr moves downwards to root;
 BEGIN
    SCALAR W;
    W:=ROOT*ROOT;
    IF N=W THEN RETURN ROOT;
    W:=(ROOT+N/ROOT)/2;
    IF W>=ROOT THEN RETURN !*P2F MKSP(MKSQRT N,1);
    RETURN NR(N,W)
 END;

GLOBAL '(SQRT!-FLOAT!-TOLERANCE);

SQRT!-FLOAT!-TOLERANCE := 0.00001;

SYMBOLIC PROCEDURE SQRT!-FLOAT N;
% Simple Newton-Raphson floating point square root calculator.
% Not warranted against truncation errors, etc;
BEGIN INTEGER SCALE; SCALAR ANS;
  IF N<0.0 THEN REDERR "SQRT!-FLOAT GIVEN NEGATIVE ARGUMENT";
  % Scale argument to within 1e-10 to 1e+10;
  SCALE := 0;
  WHILE N > 1E+10 DO <<
    SCALE := SCALE + 1;
    N := N/1E+10 >>;
  WHILE N < 1E-10 DO <<
    SCALE := SCALE - 1;
    N := N*1E-10 >>;
  ANS := IF N>2.0 THEN (N+1)/2
         ELSE IF N<0.5 THEN 2/(N+1)
         ELSE N;
  WHILE ABS(ANS**2/N - 1.0) > SQRT!-FLOAT!-TOLERANCE DO
    ANS := 0.5*(ANS+N/ANS);
  RETURN ANS*10**(5*SCALE)
END;

COMMENT Kludge to define derivative of an integral;

SYMBOLIC PUT('DF,'OPMTCH,'(((INT !&Y !&X) !&X) (NIL . T)
			   (EVL!* !&Y) NIL) . GET('DF,'OPMTCH));

GLOBAL '(FRLIS!*);

SYMBOLIC FRLIS!* := '!&X . '!&Y . FRLIS!*;

SYMBOLIC IF NOT GETD 'MODBIND
   THEN <<PUT('EVL!*,'OPMTCH,'(((!&X) (NIL . T) !&X NIL)));
	  PUT('EVL!*,'SIMPFN,'SIMPIDEN)>>;
%	  MKOP 'SQRT>>;
   %distinguish between mode and non-mode system;

ALGEBRAIC;

%FOR ALL X LET SQRT X**2=X;

SYMBOLIC;


COMMENT support for module use;

GLOBAL '(EXPORTSLIST!* IMPORTSLIST!* !*MODULEP);

DEFLIST('((EXPORTS RLIS) (IMPORTS RLIS) (MODULE RLIS)
	  (ENDMODULE ENDSTAT)),'STAT);

SYMBOLIC PROCEDURE EXPORTS U;
   BEGIN
      EXPORTSLIST!* := UNION(U,EXPORTSLIST!*);
   END;

SYMBOLIC PROCEDURE IMPORTS U;
   BEGIN
      IMPORTSLIST!* := UNION(U,IMPORTSLIST!*);
   END;

SYMBOLIC PROCEDURE MODULE U;
   %Sets up a module definition;
   BEGIN
      !*MODULEP := T;
   END;

SYMBOLIC PROCEDURE ENDMODULE;
   BEGIN
      EXPORTSLIST!* := NIL;
      IMPORTSLIST!* := NIL;
      !*MODULEP := NIL
   END;




%**********************************************************************;
% SET REDUCE AND LISP OPTIONS ONCE AND FOR ALL;

%ON COMP;



% ALL FLUID VARIABLES ARE DECLARED HERE;

FLUID '(CONTENT SQFR ZLIST INDEXLIST SQRTLIST )$
FLUID '(!*MCD !*GCD !*EXP !*SQRT !*STRUCTURE);
FLUID '(	PT ULIST
	REDUCTIONEQ LOGLIST CLIST CCOUNT CVAL CMAP TANLIST LHS
	BADPART CUBEROOTFLAG VARLIST CLOGFLAG EXPRESSION RESIDUE
	VARIABLE ORDEROFELIM CMATRIX DENOMINATOR TAYLORVARIABLE
	!*PURERISCH !*NOLNR);

%FLAGS TO BE SET USING 'ON' AND 'OFF' STATEMENTS;

GLOBAL '(!*RATINTSPECIAL !*TRINT !*SEPLOGS !*FAILHARD !*TRDIV
	!*STATISTICS !*NUMBER!* !*SPSIZE!*
    BTRLEVEL !*GENSYMLIST!*);

BTRLEVEL:=5; %DEFAULT TO A REASONABLY FULL BACKTRACE;
ON SEPLOGS;%,OVERLAYMODE;
%TOPLEVELCODE:='(COMPILER RLISP APROC);

%**********************************************************************;

SMACRO PROCEDURE FIRSTSUBS U;
CAR U;
% THE FIRST SUBSTITUTION IN A SUBSTITUTION LIST;

SMACRO PROCEDURE RSUBS U;
CDR U;

SMACRO PROCEDURE LSUBS U;
CAR U;

% THE ABOVE TWO FUNCTIONS DEFINE LEFT AND RIGHT HALVES OF A
% SUBSTITUTION RULE;


SMACRO PROCEDURE LFIRSTSUBS U;
CAAR U;

SMACRO PROCEDURE RFIRSTSUBS U;
CDAR U;

% SOME COMBINATIONS OF THE ABOVE;

SMACRO PROCEDURE ARGOF U;
CADR U;

% THE ARGUMENT OF A UNARY FUNCTION;


FLAG ('(ATAN DILOG ERF EXPINT EXPT LOG TAN),'TRANSCENDENTAL);
ALGEBRAIC;
%Patterns for integration of various logarithmic cases;
%FOR ALL X,A,B,C,D LET INT(LOG(A*X+B)/(C*X+D),X)=
%	LOG(C*X+D)*LOG(B*C-A*D)/C - LOG C*LOG(C*X+D)/C 
%	- DILOG((A*C*X+B*C)/(B*C-A*D))/C;
%% A=1;
%FOR ALL X,B,C,D LET INT(LOG(X+B)/(C*X+D),X)=
%	LOG(C*X+D)*(LOG(B*C-D)-LOG C)/C -DILOG((C*X+B*C)/(B*C-D))/C;
%% B=0;
%FOR ALL X,A,C,D LET INT(LOG(A*X)/(C*X+D),X)=
%	LOG(C*X+D)*(LOG(-1)+LOG(A)+LOG(D)-LOG C)/C - DILOG(-C*X/D)/C;
%% C=1;
%FOR ALL X,A,B,D LET INT(LOG(A*X+B)/(X+D),X)=
%	LOG(X+D)*LOG(B-A*D)-DILOG((A*X+B)/(B-A*D));
%% D=0;
%FOR ALL X,A,B,C LET INT(LOG(A*X+B)/(C*X),X)=
%	LOG(C*X)*LOG(B)/C - DILOG((A*X+B)/B)/C;
%% A=1, B=0;
%FOR ALL X,C,D LET INT(LOG(X)/(C*X+D),X)=
%	LOG(C*X+D)*(LOG(-1)+LOG(D)-LOG(C))/C - DILOG(-C*X/D)/C;
%% A=1, C=1;
%FOR ALL X,B,D LET INT(LOG(X+B)/(X+D),X)=
%	LOG(X+D)*LOG(B-D) - DILOG((X+B)/(B-D));
%% A=1, D=0;
%FOR ALL X,B,C LET INT(LOG(X+B)/(C*X),X)=
%	LOG(C*X)*LOG(B)/C - DILOG((X+B)/B)/C;
%% B=0, C=1;
%FOR ALL X,A,D LET INT(LOG(A*X)/(X+D),X)=
%	LOG(X+D)*(LOG(-1)+LOG(A)+LOG(D)) - DILOG(-X/D);
%% C=1, D=0;
%FOR ALL X,A,B LET INT(LOG(A*X+B)/X,X)=
%	LOG(X+D)*(LOG(-1)+LOG(D)) - DILOG(-X/D);
%% A=1, C=1, D=0;
%FOR ALL X,B LET INT(LOG(X+B)/X,X)=
%	LOG(X)*LOG(B) - DILOG((X+B)/B);
%% A=1, B=0, C=1;
%FOR ALL X,D LET INT(LOG(X)/(X+D),X)=
%	LOG(X+D)*(LOG(-1)+LOG(D)) - DILOG(-X/D);
%
LISP;
!*NOLNR:=NIL;
MODULE CONTENTS;

EXPORTS CONTENTS,CONTENTSMV,DFNUMR,DIFFLOGS,FACTORLISTLIST,MULTSQFREE,
	MULTUP,SQFREE,SQMERGE;

IMPORTS INT!-FAC,FQUOTF,GCDF,INTERR,!*MULTF!*,PARTIALDIFF,QUOTF,ORDOP,
	ADDF,NEGF,DOMAINP,DIFFF,MKSP,NEGSQ,INVSQ,ADDSQ,MULTSQ,DIFFSQ;


COMMENT we assume that no power substitution is necessary in
	this module;

SYMBOLIC PROCEDURE CONTENTS(P,V);
% FIND THE CONTENTS OF THE POLYNOMIAL P WRT VARIABLE V;
% NOTE THAT V MAY NOT BE THE MAIN VARIABLE OF P;
    IF DOMAINP(P) THEN P
    ELSE IF V=MVAR P THEN CONTENTSMV(P,V,NIL)
    ELSE IF ORDOP(V,MVAR P) THEN P
    ELSE CONTENTSMV(MAKEMAINVAR(P,V),V,NIL);

SYMBOLIC PROCEDURE CONTENTSMV(P,V,SOFAR);
% FIND CONTENTS OF POLYNOMIAL P;
% V IS MAIN VARIABLE OF P;
% SOFAR IS PARTIAL RESULT;
    IF SOFAR=1 THEN 1
    ELSE IF DOMAINP P THEN GCDF(P,SOFAR)
    ELSE IF NOT V=MVAR P THEN GCDF(P,SOFAR)
    ELSE CONTENTSMV(RED P,V,GCDF(LC P,SOFAR));



SYMBOLIC PROCEDURE MAKEMAINVAR(P,V);
% BRING V UP TO BE THE MAIN VARIABLE IN POLYNOMIAL P;
% NOTE THAT THE RECONSTRUCTED P MUST BE USED WITH CARE SINCE;
% IT DOES NOT CONFORM TO THE NORMAL REDUCE ORDERING RULES;
    IF DOMAINP P THEN P
    ELSE IF V=MVAR P THEN P
    ELSE MERGEADD(MULCOEFFSBY(MAKEMAINVAR(LC P,V),LPOW P,V),
      MAKEMAINVAR(RED P,V),V);

SYMBOLIC PROCEDURE MULCOEFFSBY(P,POW,V);
% MULTIPLY EACH COEFFICIENT IN P BY THE STANDARD POWER POW;
    IF NULL P THEN NIL
    ELSE IF DOMAINP P OR NOT V=MVAR P THEN ((POW .* P) .+ NIL)
    ELSE (LPOW P .* ((POW .* LC P) .+ NIL)) .+ MULCOEFFSBY(RED P,POW,V);

SYMBOLIC PROCEDURE MERGEADD(A,B,V);
% ADD POLYNOMIALS A AND B GIVEN THAT THEY HAVE SAME MAIN VARIABLE V;
    IF DOMAINP A OR NOT V=MVAR A THEN
      IF DOMAINP B OR NOT V=MVAR B THEN ADDF(A,B)
      ELSE LT B .+ MERGEADD(A,RED B,V)
    ELSE IF DOMAINP B OR NOT V=MVAR B THEN
      LT A .+ MERGEADD(RED A,B,V)
    ELSE (LAMBDA XC;
      IF XC=0 THEN (LPOW A .* ADDF(LC A,LC B)) .+
	    MERGEADD(RED A,RED B,V)
      ELSE IF XC>0 THEN LT A .+ MERGEADD(RED A,B,V)
      ELSE LT B .+ MERGEADD(A,RED B,V))
	(TDEG LT A-TDEG LT B);



SYMBOLIC PROCEDURE SQFREE(P,VL);
    IF (NULL VL) OR (DOMAINP P) THEN
	<<CONTENT:=P; NIL>>
    ELSE BEGIN    SCALAR W,V,DP,GG,PG,DPG,P1,W1;
	W:=CONTENTS(P,CAR VL); % CONTENT OF P ;
	P:=QUOTF(P,W); % MAKE P PRIMITIVE;
	W:=SQFREE(W,CDR VL); % PROCESS CONTENT BY RECURSION;
	IF P=1 THEN RETURN W;
	V:=CAR VL; % PICK OUT VARIABLE FROM LIST;
	WHILE NOT (P=1) DO <<
	    DP:=PARTIALDIFF(P,V);
	    GG:=GCDF(P,DP);
	    PG:=QUOTF(P,GG);
	    DPG:=NEGF PARTIALDIFF(PG,V);
	    P1:=GCDF(PG,ADDF(QUOTF(DP,GG),DPG));
	    W1:=P1.W1;
	    P:=GG>>;
	RETURN SQMERGE(REVERSE W1,W,T)
	END;

SYMBOLIC PROCEDURE SQMERGE(W1,W,SIMPLEW1);
% W AND W1 ARE LISTS OF FACTORS OF EACH POWER. IF SIMPLEW1 IS TRUE
% THEN W1 CONTAINS ONLY SINGLE FACTORS FOR EACH POWER. ;
    IF NULL W1 THEN W
    ELSE IF NULL W THEN IF CAR W1=1 THEN NIL.SQMERGE(CDR W1,W,SIMPLEW1)
	  ELSE (IF SIMPLEW1 THEN LIST CAR W1 ELSE CAR W1).
SQMERGE(CDR W1,W,SIMPLEW1)
    ELSE IF CAR W1=1 THEN (CAR W).SQMERGE(CDR W1,CDR W,SIMPLEW1) ELSE
	APPEND(IF SIMPLEW1 THEN LIST CAR W1 ELSE CAR W1,CAR W).
	SQMERGE(CDR W1,CDR W,SIMPLEW1);

SYMBOLIC PROCEDURE MULTUP L;
% L IS A LIST OF S.F.'S. RESULT IS S.Q. FOR PRODUCT OF ELEMENTS OF L;
   BEGIN	 SCALAR RES;
      RES:=1 ./ 1;
      WHILE NOT NULL L DO <<
	 RES:=MULTSQ(RES,(CAR L) ./ 1);
	 L:=CDR L >>;
      RETURN RES
   END;

SYMBOLIC PROCEDURE DIFLIST(L,CL,X,RL);
% DIFFERENTIATES L (LIST OF S.F.'S) WRT X TO PRODUCE THE SUM OF;
% TERMS FOR THE DERIVATIVE OF NUMR OF 1ST PART OF ANSWER.  CL IS;
% COEFFICIENT LIST (S.F.'S) & RL IS LIST OF DERIVATIVES WE HAVE;
% DEALT WITH SO FAR;
% RESULT IS S.Q.;
   IF NULL L THEN NIL ./ 1
   ELSE BEGIN    SCALAR TEMP;
      TEMP:=MULTSQ(MULTUP RL,MULTUP CDR L);
      TEMP:=MULTSQ(DIFFF(CAR L,X),TEMP);
      TEMP:=MULTSQ(TEMP,(CAR CL) ./ 1);
      RETURN ADDSQ(TEMP,DIFLIST(CDR L,CDR CL,X,(CAR L).RL))
   END;

SYMBOLIC PROCEDURE MULTSQFREE W;
% W IS LIST OF SQFREE FACTORS. RESULT IS PRODUCT OF EACH LIST IN W
% TO GIVE ONE POLYNOMIAL FOR EACH SQFREE POWER;
   IF NULL W THEN NIL
   ELSE (!*Q2F MULTUP CAR W).MULTSQFREE CDR W;

SYMBOLIC PROCEDURE L2LSF L;
% L IS A LIST OF KERNELS. RESULT IS A LIST OF SAME MEMBERS AS S.F.'S;
   IF NULL L THEN NIL
   ELSE ((MKSP(CAR L,1) .* 1) .+ NIL).L2LSF CDR L;

SYMBOLIC PROCEDURE DFNUMR(X,DL);
% GIVES THE DERIVATIVE OF THE NUMR OF THE 1ST PART OF ANSWER.;
% DL IS LIST OF ANY EXPONENTIAL OR 1+TAN**2 THAT OCCUR IN INTEGRAND;
% DENR. THESE ARE DIVIDED OUT FROM RESULT BEFORE HANDING IT BACK.;
% RESULT IS S.Q., READY FOR PRINTING;
   BEGIN	 SCALAR TEMP1,TEMP2,COEFLIST,QLIST,COUNT;
      IF NOT NULL SQFR THEN <<
      COUNT:=0;
      QLIST:=CDR SQFR;
      COEFLIST:=NIL;
      WHILE NOT NULL QLIST DO <<
	 COUNT:=COUNT+1;
	 COEFLIST:=COUNT.COEFLIST;
	 QLIST:=CDR QLIST >>;
      COEFLIST:=REVERSE COEFLIST >>;
      TEMP1:=MULTSQ(DIFLIST(L2LSF ZLIST,L2LSF INDEXLIST,X,NIL),
		    MULTUP SQFR);
      IF NOT NULL SQFR AND NOT NULL CDR SQFR THEN <<
      TEMP2:=MULTSQ(DIFLIST(CDR SQFR,COEFLIST,X,NIL),
	    MULTUP L2LSF ZLIST);
      TEMP2:=MULTSQ(TEMP2,(CAR SQFR) ./ 1) >>
      ELSE TEMP2:=NIL ./ 1;
      TEMP1:=ADDSQ(TEMP1,NEGSQ TEMP2);
      TEMP2:=CDR TEMP1;
      TEMP1:=CAR TEMP1;
      QLIST:=NIL;
      WHILE NOT NULL DL DO <<
         IF NOT CAR DL MEMBER QLIST THEN QLIST:=(CAR DL).QLIST;
         DL:=CDR DL >>;
      WHILE NOT NULL QLIST DO <<
	 TEMP1:=QUOTF(TEMP1,CAR QLIST);
	 QLIST:=CDR QLIST >>;
      RETURN TEMP1 ./ TEMP2
   END;

SYMBOLIC PROCEDURE DIFFLOGS(LL,DENM1,X);
% LL IS LIST OF LOG TERMS (WITH COEFFTS), DEN IS COMMON DENOMINATOR;
% OVER WHICH THEY ARE TO BE PUT.  RESULT IS S.Q. FOR DERIVATIVE OF ALL;
% THESE WRT X;
   IF NULL LL THEN NIL ./ 1
   ELSE BEGIN    SCALAR TEMP,QU,CVAR,LOGORATAN,ARG;
      LOGORATAN:=CAAR LL;
      CVAR:=CADAR LL;
      ARG:=CDDAR LL;
      TEMP:=MULTSQ(CVAR ./ 1,DIFFSQ(ARG,X));
      IF LOGORATAN='IDEN THEN QU:=1 ./ 1
	ELSE IF LOGORATAN='LOG THEN QU:=ARG
	ELSE IF LOGORATAN='ATAN THEN QU:=ADDSQ(1 ./ 1,MULTSQ(ARG,ARG))
	ELSE INTERR "LOGORATAN=? IN DIFFLOGS";
%NOTE CALL TO SPECIAL DIVISION ROUTINE;
      QU:=FQUOTF(!*F2POL !*MULTF!*(!*MULTF!*(DENM1,NUMR TEMP),
		DENR QU),NUMR QU);
			%*MUST* GO EXACTLY;
     TEMP:=MULTSQ(INVSQ (DENR TEMP ./ 1),QU);
		 %RESULT OF FQUOTF IS A S.Q;
      RETURN SUBS2Q ADDSQ(TEMP,DIFFLOGS(CDR LL,DENM1,X))
   END;

SYMBOLIC PROCEDURE FACTORLISTLIST (W,CLOGFLAG);
% W IS LIST OF LISTS OF SQFREE FACTORS IN S.F.	RESULT IS LIST OF LOG;
% TERMS REQUIRED FOR INTEGRAL ANSWER. THE ARGUMENTS FOR EACH LOG FN;
% ARE IN S.Q.;
    BEGIN SCALAR RES,X,Y;
	WHILE NOT NULL W DO <<
	    X:=CAR W;
	    WHILE NOT NULL X DO <<
		Y:=FACBYPP(CAR X,VARLIST);
		WHILE NOT NULL Y DO <<
		    RES:=APPEND(INT!-FAC CAR Y,RES);
		    Y:=CDR Y >>;
		X:=CDR X >>;
	    W:=CDR W >>;
	RETURN RES
    END;

SYMBOLIC PROCEDURE FACBYPP(P,VL);
%USE CONTENTS/PRIMITIVE PARTS TO TRY TO FACTOR P;
    IF NULL VL THEN LIST P
    ELSE BEGIN SCALAR PRINCILAP!-PART,CO;
	CO:=CONTENTS(P,CAR VL);
	VL:=CDR VL;
	IF CO=1 THEN RETURN FACBYPP(P,VL); %THIS VAR NO HELP;
	PRINCILAP!-PART:=QUOTF(P,CO); %PRIMITIVE PART;
	IF PRINCILAP!-PART=1 THEN RETURN FACBYPP(P,VL); %AGAIN NO HELP;
	RETURN NCONC(FACBYPP(PRINCILAP!-PART,VL),FACBYPP(CO,VL))
    END;


ENDMODULE;


MODULE CSOLVE;

EXPORTS BACKSUBST4CS,CREATECMAP,FINDPIVOT,PRINTSPREADC,PRINTVECSQ,
   SPREADC,SUBST4ELIMINATEDS;

IMPORTS NTH,INTERR,!*MULTF!*,PRINTSF,PRINTSQ,QUOTF,PUTV,NEGF,INVSQ,
   NEGSQ,ADDSQ,MULTSQ,MKSP,ADDF,DOMAINP,PNTH;


% routines to do with the C constants;

SYMBOLIC PROCEDURE FINDPIVOT CVEC;
% Finds first non-zero element in CVEC and returns its cell number.;
% If no such element exists, result is nil.;
   BEGIN	 SCALAR I,X;
      I:=1;
      X:=GETV(CVEC,I);
      WHILE I<CCOUNT AND NULL X DO
      << I:=I+1;
	 X:=GETV(CVEC,I) >>;
      IF NULL X THEN RETURN NIL;
      RETURN I
   END;

SYMBOLIC PROCEDURE SUBST4ELIMINATEDCS(NEWEQN,SUBSTORDER,CEQNS);
% Substitutes into NEWEQN for all the C's that have been eliminated so;
% far. These are given by CEQNS. SUBSTORDER gives the order of;
% substitution as well as the constant multipliers. Result is the;
% transformed NEWEQN.;
   IF NULL SUBSTORDER THEN NEWEQN
   ELSE BEGIN    SCALAR NXT,ROW,CVAR,TEMP;
      ROW:=CAR CEQNS;
      NXT:=CAR SUBSTORDER;
      IF NULL (CVAR:=GETV(NEWEQN,NXT)) THEN
	 RETURN SUBST4ELIMINATEDCS(NEWEQN,CDR SUBSTORDER,CDR CEQNS);
      NXT:=GETV(ROW,NXT);
      FOR I:=0 : CCOUNT DO
      << TEMP:=!*MULTF!*(NXT,GETV(NEWEQN,I));
	 TEMP:=ADDF(TEMP,NEGF !*MULTF!*(CVAR,GETV(ROW,I)));
	 PUTV(NEWEQN,I,!*F2POL TEMP) >>;
      RETURN SUBST4ELIMINATEDCS(NEWEQN,CDR SUBSTORDER,CDR CEQNS)
   END;


SYMBOLIC PROCEDURE BACKSUBST4CS(CS2SUBST,CS2SOLVE,CMATRIX);
% Solves the C-eqns and sets vector CVAL to the C-constant values;
% CMATRIX is a list of matrix rows for C-eqns after Gaussian ;
% elimination has been performed. CS2SOLVE is a list of the remaining;
% C's to evaluate and CS2SUBST are the C's we have evaluated already.;
   IF NULL CMATRIX THEN NIL
   ELSE BEGIN    SCALAR EQNN,CVAR,ALREADY,SUBSTLIST,TEMP,TEMP2;
      EQNN:=CAR CMATRIX;
      CVAR:=CAR CS2SOLVE;
      ALREADY:=NIL ./ 1; % The S.Q. nil ;
      SUBSTLIST:=CS2SUBST;
% NOW SUBSTITUTE FOR PREVIOUSLY EVALUATED C'S:;
      WHILE NOT NULL SUBSTLIST DO
      << TEMP:=CAR SUBSTLIST;
	 IF NOT NULL GETV(EQNN,TEMP) THEN
	    ALREADY:=ADDSQ(ALREADY,MULTSQ(GETV(EQNN,TEMP) ./ 1,
				 GETV(CVAL,TEMP)));
	 SUBSTLIST:=CDR SUBSTLIST >>;
% NOW SOLVE FOR THE C GIVEN BY CVAR (ANY REMAINING C'S ASSUMED ZERO);
      TEMP:=NEGSQ ADDSQ(GETV(EQNN,0) ./ 1,ALREADY);
      IF NOT NULL (TEMP2:=QUOTF(NUMR TEMP,GETV(EQNN,CVAR))) THEN
				       TEMP:=TEMP2 ./ DENR TEMP
      ELSE TEMP:=MULTSQ(TEMP,INVSQ(GETV(EQNN,CVAR) ./ 1));
      IF NOT NULL NUMR TEMP THEN PUTV(CVAL,CVAR,
		RESIMP ROOTEXTRACTSQ SUBS2Q TEMP);
      BACKSUBST4CS(REVERSEWOC(CVAR . REVERSEWOC CS2SUBST),
	    CDR CS2SOLVE,CDR CMATRIX)
   END;

%**********************************************************************;
% Routines to deal with linear equations for the constants C;
%**********************************************************************;

SYMBOLIC PROCEDURE CREATECMAP;
%Sets LOGLIST to list of things of form (LOG C-constant f), where f is;
% function linear in one of the z-variables and C-constant is in S.F.;
% When creating these C-constant names, the CMAP is also set up and ;
% returned as the result.;
   BEGIN	 SCALAR I,L,C;
      L:=LOGLIST;
      I:=1;
      WHILE NOT NULL L DO <<
	 C:=(GENSYM1('C) . I) . C;
	 I:=I+1;
	 RPLACD(CAR L,((MKSP(CAAR C,1) .* 1) .+ NIL) . CDAR L);
	 L:=CDR L >>;
      IF !*TRINT THEN PRINTC ("Constants Map" . C);
      RETURN C
   END;


SYMBOLIC PROCEDURE SPREADC(EQNN,CVEC1,W);
%SETS A VECTOR 'CVEC1' TO COEFFICIENTS OF C<I> IN EQNN;
    IF DOMAINP EQNN THEN PUTV(CVEC1,0,ADDF(GETV(CVEC1,0),
				!*F2POL !*MULTF!*(EQNN,W)))
    ELSE BEGIN    SCALAR MV,T1,T2;
	SPREADC(RED EQNN,CVEC1,W);
	MV:=MVAR EQNN;
	T1:=ASSOC(MV,CMAP); %TESTS IF IT IS A C VAR;
	IF NOT NULL T1 THEN RETURN <<
	    T1:=CDR T1; %LOC IN VECTOR FOR THIS C;
	    IF NOT (TDEG LT EQNN=1) THEN INTERR "NOT LINEAR IN C EQN";
	    T2:=ADDF(GETV(CVEC1,T1),!*MULTF!*(W,LC EQNN));
	    PUTV(CVEC1,T1,!*F2POL T2) >>;
	T1:=((LPOW EQNN) .* 1) .+ NIL; %THIS MAIN VAR AS SF;
	SPREADC(LC EQNN,CVEC1,!*F2POL !*MULTF!*(W,T1))
    END;

SYMBOLIC PROCEDURE PRINTSPREADC CVEC1;
    BEGIN
	FOR I:=0 : CCOUNT DO <<
	   PRIN2 I;
	   PRINTC ":";
	   PRINTSF(GETV(CVEC1,I)) >>;
	PRINTC "END OF PRINTSPREADC OUTPUT"
    END;

%SYMBOLIC PROCEDURE PRINTVECSQ CVEC;
%% PRINT CONTENTS OF CVEC WHICH CONTAINS S.Q.'S (NOT S.F.'S);
%% STARTS FROM CELL 1 NOT 0 AS ABOVE ROUTINE (PRINTSPREADC);
%   BEGIN
%      FOR I:=1 : CCOUNT DO <<
%	 PRIN2 I;
%	 PRINTC ":";
%	 IF NULL GETV(CVEC,I) THEN PRINTC "0"
%	 ELSE PRINTSQ(GETV(CVEC,I)) >>;
%      PRINTC "END OF PRINTVECSQ OUTPUT"
%   END;


ENDMODULE;


MODULE CUBEROOT;

EXPORTS CUBEROOTDF;

IMPORTS CONTENTSMV,GCDF,!*MULTF!*,NROOTN,PARTIALDIFF,PRINTDF,QUOTF,VP2,
   MKSP,MK!*SQ,DOMAINP;

%CUBE-ROOT OF STANDARD FORMS;





SYMBOLIC PROCEDURE CUBEROOTSQ A;
    CUBEROOTF NUMR A ./ CUBEROOTF DENR A;

SYMBOLIC PROCEDURE CUBEROOTF P;
    BEGIN	SCALAR IP,QP;
	IF NULL P THEN RETURN NIL;
	IP:=CUBEROOTF1 P;
	QP:=CDR IP;
	IP:=CAR IP; %RESPECTABLE AND NASTY PARTS OF THE CUBEROOT;
	IF ONEP QP THEN RETURN IP; %EXACT ROOT FOUND;
	QP:=LIST('EXPT,PREPF QP,'(QUOTIENT 1 3));
	CUBEROOTFLAG:=T; %SYMBOLIC CUBE-ROOT INTRODUCED;
	QP:=(MKSP(QP,1).* 1) .+ NIL;
	RETURN !*F2POL !*MULTF!*(IP,QP)
    END;

SYMBOLIC PROCEDURE CUBEROOTF1 P;
	
%RETURNS A . B WITH P=A**2*B;
	%does this need power reduction??;
    IF DOMAINP P THEN NROOTN(P,3)
    ELSE BEGIN SCALAR CO,PPP,G,PG;
	CO:=CONTENTSMV(P,MVAR P,NIL); %CONTENTS OF P;
	PPP:=QUOTF(P,CO); %PRIMITIVE PART;
%NOW CONSIDER PPP=P1*P2**2*P3**3*P4**4*...;
	CO:=CUBEROOTF1(CO); %PROCESS CONTENTS VIA RECURSION;
	G:=GCDF(PPP,PARTIALDIFF(PPP,MVAR PPP));
%G=P2*P3**2*P4**3*...;
	IF NOT DOMAINP G THEN <<
	    PG:=QUOTF(PPP,G);
    %PG=P1*P2*P3*P4*...;
	    G:=GCDF(G,PARTIALDIFF(G,MVAR G));
    % G=G3*G4**2*G5**3*...;
	    G:=GCDF(G,PG)>>; %A TRIPLE FACTOR OF PPP;
	IF DOMAINP G THEN PG:=1 . PPP
	ELSE <<
	    PG:=QUOTF(PPP,!*MULTF!*(G,!*MULTF!*(G,G))); %WHAT'S LEFT;
	    PG:=CUBEROOTF1(!*F2POL PG); %SPLIT THAT UP;
	    RPLACA(PG,!*MULTF!*(CAR PG,G))>>;
		 %PUT IN THE THING FOUND HERE;
	RPLACA(PG,!*F2POL !*MULTF!*(CAR PG,CAR CO));
	RPLACD(PG,!*F2POL !*MULTF!*(CDR PG,CDR CO));
	RETURN PG
    END;

ENDMODULE;


MODULE DEPEND;

EXPORTS DEPENDSPL,DEPENDSP,INVOLVESQ,INVOLVSF;

IMPORTS TAYLORP,DOMAINP;


SYMBOLIC PROCEDURE DEPENDSP(X,V);
    IF NULL V THEN T
     ELSE IF ATOM X THEN IF X EQ V THEN X ELSE NIL
    ELSE IF CAR X = '!*SQ
      THEN INVOLVESQ(CADR X,V)
    ELSE IF TAYLORP X
     THEN IF V EQ TAYLORVARIABLE THEN TAYLORVARIABLE ELSE NIL
    ELSE BEGIN
     SCALAR W;
    IF X=V THEN RETURN V;
% CHECK IF A PREFIX FORM EXPRESSION DEPENDS ON THE VARIABLE V;
% NOTE THAT THIS ASSUMES THE FORM X IS IN NORMAL PREFIX NOTATION;
      W := X; % preserve the dependency;
      X:=CDR X; % READY TO RECURSIVELY CHECK ARGUMENTS;
SCAN: IF NULL X THEN RETURN NIL; % NO DEPENDENCY FOUND;
	IF DEPENDSP(CAR X,V) THEN RETURN W;
	X:=CDR X;
	GO TO SCAN
    END;

SYMBOLIC PROCEDURE TAYLORP U; NIL;  %dummy for now;


SYMBOLIC PROCEDURE INVOLVESQ(SQ,TERM);
INVOLVESF(NUMR SQ,TERM) OR INVOLVESF(DENR SQ,TERM);


SYMBOLIC PROCEDURE INVOLVESF(SF,TERM);
IF DOMAINP SF OR NULL SF
  THEN NIL
  ELSE IF DEPENDSP(MVAR SF,TERM)
    THEN T
    ELSE INVOLVESF(LC SF,TERM) OR
	 INVOLVESF(RED SF,TERM);

ENDMODULE;


MODULE DF2Q;

EXPORTS DF2Q;

IMPORTS ADDF,GCDF,MKSP,!*MULTF!*,QUOTF;

COMMENT This module converts distributed forms to standard forms.
	We assume that results already have reduced powers, so
	that no power substitution is necessary;

%TRIAL REPLACEMENT FOR DF2Q;
SYMBOLIC PROCEDURE DF2Q P;
% Converts distributed form P to standard quotient;
    BEGIN	SCALAR N,D,GG,W;
	IF NULL P THEN RETURN NIL ./ 1;
	D:=DENR LC P;
	W:=RED P;
	WHILE NOT NULL W DO <<
	    GG:=GCDF(D,DENR LC W); %GET DENOMINATOR OF ANSWER...;
	    D:=!*MULTF!*(D,QUOTF(DENR LC W,GG));
		 %..AS LCM OF DENOMS IN INPUT;
	    W:=RED W >>;
	N:=NIL; %PLACE TO BUILD NUMERATOR OF ANSWER;
	WHILE NOT NULL P DO <<
	    N:=ADDF(N,!*MULTF!*(XL2F(LPOW P,ZLIST,INDEXLIST),
		!*MULTF!*(NUMR LC P,QUOTF(D,DENR LC P))));
	    P:=RED P >>;
	RETURN N ./ D
    END;

SYMBOLIC PROCEDURE XL2F(L,Z,IL);
% L is an exponent list from a D.F., Z is the Z-list,
% IL is the list of indices.
% Value is L converted to standard form. ;
    IF NULL Z THEN 1
	ELSE IF CAR L=0 THEN XL2F(CDR L,CDR Z,CDR IL)
	ELSE IF NOT ATOM CAR L THEN
	    BEGIN	SCALAR TEMP;
		IF CAAR L=0 THEN TEMP:= CAR IL
		ELSE TEMP:=LIST('PLUS,CAR IL,CAAR L);
		TEMP:=MKSP(LIST('EXPT,CAR Z,TEMP),1);
		RETURN !*MULTF!*(((TEMP .* 1) .+ NIL),
			       XL2F(CDR L,CDR Z,CDR IL))
	    END
%	ELSE IF MINUSP CAR L THEN				      ;
%	     MULTSQ(INVSQ (((MKSP(CAR Z,-CAR L) .* 1) .+ NIL)),       ;
%		   XL2F(CDR L,CDR Z,CDR IL))			      ;
	ELSE !*MULTF!*((MKSP(CAR Z,CAR L) .* 1) .+ NIL,
		    XL2F(CDR L,CDR Z,CDR IL));


ENDMODULE;


MODULE DISTRIB;

EXPORTS DFPRINTFORM,MULTBYARBPOWERS,NEGDF,QUOTDFCONST,SUB1IND,VP1,
   VP2,PLUSDF,MULTDF,MULTDFCONST,ORDDF;

IMPORTS INTERR,ADDSQ,NEGSQ,EXPTSQ,SIMP,DOMAINP,MK!*SQ,ADDF,
   MULTSQ,INVSQ,MINUSP,MKSP,SUB1;

%***********************************************************************
%  ROUTINES FOR MANIPULATING DISTRIBUTED FORMS.
%	NOTE:
%	    THE EXPRESSIONS LT,RED,LC,LPOW HAVE BEEN USED ON DISTRIBUTED
%	    FORMS AS THE LATTER'S STRUCTURE IS SUFFICIENTLY SIMILAR TO
%	    S.F.'S.  HOWEVER LC DF IS A S.Q. NOT A S.F. AND LPOW DF IS A
%	    LIST OF THE EXPONENTS OF THE VARIABLES.  THIS ALSO MAKES
%	    LT DF DIFFERENT.  RED DF IS D.F. AS EXPECTED.
%**********************************************************************;

SYMBOLIC PROCEDURE PLUSDF(U,V);
% U and V are D.F.'s. Value is D.F. for U+V;
    IF NULL U THEN V
	ELSE IF NULL V THEN U
	ELSE IF LPOW U=LPOW V THEN
	    (LAMBDA(X,Y); IF NULL NUMR X THEN Y ELSE (LPOW U .* X) .+ Y)
	    (ADDSQ(LC U,LC V),PLUSDF(RED U,RED V))
	ELSE IF ORDDF(LPOW U,LPOW V) THEN LT U .+ PLUSDF(RED U,V)
	ELSE (LT V) .+ PLUSDF(U,RED V);

SYMBOLIC PROCEDURE ORDDF(U,V);
% U and V are the LPOW of a D.F. - i.e. the list of exponents ;
% Value is true if LPOW U '>' LPOW V and false otherwise ;
    IF NULL U THEN IF NULL V THEN INTERR "ORDDF = CASE"
	ELSE INTERR "ORDDF V LONGER THAN U"
	ELSE IF NULL V THEN INTERR "ORDDF U LONGER THAN V"
	ELSE IF EXPTCOMPARE(CAR U,CAR V) THEN T
	ELSE IF EXPTCOMPARE(CAR V,CAR U) THEN NIL
	ELSE ORDDF(CDR U,CDR V);

SYMBOLIC PROCEDURE EXPTCOMPARE(X,Y);
    IF ATOM X THEN IF ATOM Y THEN X>Y ELSE NIL
	ELSE IF ATOM Y THEN T
	ELSE CAR X > CAR Y;

SYMBOLIC PROCEDURE NEGDF U;
    IF NULL U THEN NIL
	ELSE (LPOW U .* NEGSQ LC U) .+ NEGDF RED U;

SYMBOLIC PROCEDURE MULTDF(U,V);
% U and V are D.F.'s. Value is D.F. for U*V;
% reduces squares of square-roots as it goes;
    IF NULL U OR NULL V THEN NIL
    ELSE BEGIN SCALAR Y;
%use (a+b)*(c+d) = (a*c) + a*(c+d) + b*(c+d);
	Y:=MULTERM(LT U,LT V); %leading terms;
	Y:=PLUSDF(Y,MULTDF(RED U,V));
	Y:=PLUSDF(Y,MULTDF((LT U) .+ NIL,RED V));
	RETURN Y
    END;

SYMBOLIC PROCEDURE MULTERM(U,V);
%multiply two terms to give a D.F.;
    BEGIN SCALAR COEF;
       COEF:= SUBS2Q MULTSQ(CDR U,CDR V); %coefficient part;
       RETURN MULTDFCONST(COEF,MULPOWER(CAR U,CAR V))
    END;

SYMBOLIC PROCEDURE MULPOWER(U,V);
% u and v are exponent lists. multiply corresponding forms;
    BEGIN SCALAR R,S;
       R:=ADDEXPTSDF(U,V);
	IF NOT NULL SQRTLIST THEN S:=REDUCEROOTS(R,ZLIST);
       R:=(R .* (1 ./ 1)) .+ NIL;
       IF NOT (S=NIL) THEN R:=MULTDF(R,S);
       RETURN R
    END;

SYMBOLIC PROCEDURE REDUCEROOTS(R,ZL); 
    BEGIN SCALAR S; 
       WHILE NOT NULL R DO << 
          IF EQCAR(CAR ZL,'SQRT) THEN 
              S:=TRYREDUCTION(R,CAR ZL,S); 
          R:=CDR R; ZL:=CDR ZL >>; 
       RETURN S 
    END; 

SYMBOLIC PROCEDURE TRYREDUCTION(R,VAR,S);
   BEGIN SCALAR X;
      X:=CAR R; %CURRENT EXPONENT;
      IF NOT ATOM X THEN << R:=X; X:=CAR R >>; %NUMERIC PART;
      IF (X=0) OR (X=1) THEN RETURN S; %NO REDUCTION POSSIBLE;
      X:=DIVIDE(X,2);
      RPLACA(R,CDR X); %REDUCE EXPONENT AS REDORDED;
      X:=CAR X;
      VAR:=SIMP CADR VAR; %SQRT ARG AS A S Q;
      VAR:=EXPTSQ(VAR,X);
      X:=MULTDFCONST(1 ./ DENR VAR,F2DF NUMR VAR); %DISTRIBUTE;
      IF S=NIL THEN S:=X
      ELSE S:=MULTDF(S,X);
      RETURN S
   END;



SYMBOLIC PROCEDURE ADDEXPTSDF(X,Y);
% X and Y are LPOW's of D.F. Value is list of sum of exponents;
    IF NULL X THEN IF NULL Y THEN NIL ELSE INTERR "X TOO LONG"
	ELSE IF NULL Y THEN INTERR "Y TOO LONG"
	ELSE EXPTPLUS(CAR X,CAR Y).ADDEXPTSDF(CDR X,CDR Y);

SYMBOLIC PROCEDURE EXPTPLUS(X,Y);
    IF ATOM X THEN IF ATOM Y THEN X+Y ELSE LIST (X+CAR Y)
	ELSE IF ATOM Y THEN LIST (CAR X +Y)
	ELSE INTERR "BAD EXPONENT SUM";

SYMBOLIC PROCEDURE MULTDFCONST(X,U);
% X is S.Q. not involving Z variables of D.F. U. Value is D.F.;
% for X*U;
    IF (NULL U) OR (NULL NUMR X) THEN NIL
	ELSE LPOW U .* SUBS2Q MULTSQ(X,LC U) .+ MULTDFCONST(X,RED U);

SYMBOLIC PROCEDURE F2DF P;
% P is standard form. Value is P in D.F.;
    IF DOMAINP P THEN DFCONST(P ./ 1)
	ELSE IF MVAR P MEMBER ZLIST THEN
	     PLUSDF(MULTDF(VP2DF(MVAR P,TDEG LT P,ZLIST),F2DF LC P),
		    F2DF RED P)
	ELSE PLUSDF(MULTDFCONST(((LPOW P .* 1) .+ NIL) ./ 1,F2DF LC P),
		    F2DF RED P);

SYMBOLIC PROCEDURE VP1(VAR,DEGG,Z);
% Takes VAR and finds it in Z (=list), raises it to power DEGG and puts;
% the result in exponent list form for use in a distributed form.;
    IF NULL Z THEN INTERR "VAR NOT IN Z-LIST AFTER ALL"
	ELSE IF VAR=CAR Z THEN DEGG.VP2 CDR Z
	ELSE 0 . VP1(VAR,DEGG,CDR Z);

SYMBOLIC PROCEDURE VP2 Z;
% Makes exponent list of zeroes;
    IF NULL Z THEN NIL
	ELSE 0 . VP2 CDR Z;

SYMBOLIC PROCEDURE VP2DF(VAR,EXPRN,Z);
% Makes VAR**EXPRN into exponent list and then converts the resulting
% power into a distributed form.
% special care with square-roots;
IF EQCAR(VAR,'SQRT) AND EXPRN>1 THEN 
	MULPOWER(VP1(VAR,EXPRN,Z),VP2 Z)
   ELSE (VP1(VAR,EXPRN,Z) .* (1 ./ 1)) .+ NIL;

SYMBOLIC PROCEDURE DFCONST Q;
% Makes a distributed form from standard quotient constant Q;
    IF NUMR Q=NIL THEN NIL
	ELSE ((VP2 ZLIST) .* Q) .+ NIL;

%DF2Q MOVED TO A SECTION OF ITS OWN;
SYMBOLIC PROCEDURE DF2PRINTFORM P;
%CONVERT TO A STANDARD FORM GOOD ENOUGH FOR PRINTING;
    IF NULL P THEN NIL
    ELSE BEGIN
	SCALAR MV,CO;
	MV:=XL2Q(LPOW P,ZLIST,INDEXLIST);
	IF MV=(1 ./ 1) THEN <<
	    CO:=LC P;
	    IF DENR CO=1 THEN RETURN ADDF(NUMR CO,
		DF2PRINTFORM RED P);
	    CO:=MKSP(MK!*SQ CO,1);
	    RETURN (CO .* 1) .+ DF2PRINTFORM RED P >>;
	CO:=LC P;
	IF NOT (DENR CO=1) THEN MV:=MULTSQ(MV,1 ./ DENR CO);
	MV:=MKSP(MK!*SQ MV,1) .* NUMR CO;
	RETURN MV .+ DF2PRINTFORM RED P
    END;


SYMBOLIC PROCEDURE XL2Q(L,Z,IL);
% L is an exponent list from a D.F., Z is the Z-list,
% IL is the list of indices.
% Value is L converted to standard quotient. ;
    IF NULL Z THEN 1 ./ 1
	ELSE IF CAR L=0 THEN XL2Q(CDR L,CDR Z,CDR IL)
	ELSE IF NOT ATOM CAR L THEN
	    BEGIN	  SCALAR TEMP;
		IF CAAR L=0 THEN TEMP:= CAR IL
		ELSE TEMP:=LIST('PLUS,CAR IL,CAAR L);
		TEMP:=MKSP(LIST('EXPT,CAR Z,TEMP),1);
		RETURN MULTSQ(((TEMP .* 1) .+ NIL) ./ 1,
			       XL2Q(CDR L,CDR Z,CDR IL))
	    END
	ELSE IF MINUSP CAR L THEN
	     MULTSQ(INVSQ (((MKSP(CAR Z,-CAR L) .* 1) .+ NIL) ./ 1),
		   XL2Q(CDR L,CDR Z,CDR IL))
	ELSE MULTSQ(((MKSP(CAR Z,CAR L) .* 1) .+ NIL) ./ 1,
		    XL2Q(CDR L,CDR Z,CDR IL));


SYMBOLIC PROCEDURE MULTBYARBPOWERS U;
% Multiplies the ordinary D.F., U, by arbitrary powers
% of the z-variables;
%	i-1  j-1  k-1
% i.e. x    z	 z    ... so result is D.F. with the exponent list
%	     1	  2
% appropriately altered to contain list elements instead of numeric
% ones;
   IF NULL U THEN NIL
   ELSE ((ADDARBEXPTSDF LPOW U) .* LC U) .+ MULTBYARBPOWERS RED U;

SYMBOLIC PROCEDURE ADDARBEXPTSDF X;
% Adds the arbitrary powers to powers in exponent list, X, to produce
% new exponent list. e.g. 3 -> (2) to represent x**3 now becoming:
%	   3	i-1    i+2
%	  x  * x    = x      . ;
   IF NULL X THEN NIL
   ELSE LIST EXPTPLUS(CAR X,-1) . ADDARBEXPTSDF CDR X;


ENDMODULE;


MODULE DIVIDE;

EXPORTS FQUOTF,TESTDIVDF,DFQUOTDF;

IMPORTS DF2Q,F2DF,GCDF,INTERR,MULTDF,NEGDF,PLUSDF,PRINTDF,PRINTSF,
   QUOTF,MULTSQ,INVSQ,NEGSQ;

%EXACT DIVISION OF STANDARD FORMS TO GIVE A STANDARD QUOTIENT;
%INTENDED FOR DIVIDING OUT KNOWN FACTORS AS PRODUCED BY THE;
%INTEGRATION PROGRAM. HORRIBLE AND SLOW, I EXPECT!!;

SYMBOLIC PROCEDURE DFQUOTDF(A,B);
    BEGIN	SCALAR RESIDUE;
	IF (!*TRINT OR !*TRDIV) THEN <<
	    PRINTC "DFQUOTDF CALLED ON ";
	    PRINTDF A; PRINTDF B>>;
	A:=DFQUOTDF1(A,B);
	IF (!*TRINT OR !*TRDIV) THEN << PRINTC "QUOTIENT GIVEN AS ";
	    PRINTDF A >>;
	IF NOT NULL RESIDUE THEN BEGIN
	    SCALAR GRES,W;
	    IF !*TRINT OR !*TRDIV THEN <<
	    PRINTC "RESIDUE IN DFQUOTDF =";
	    PRINTDF RESIDUE;
	    PRINTC "WHICH SHOULD BE ZERO";
	    W:=RESIDUE;
	    GRES:=NUMR LC W; W:=RED W;
	    WHILE NOT NULL W DO <<
		GRES:=GCDF(GRES,NUMR LC W);
		W:=RED W >>;
	    PRINTC "I.E. THE FOLLOWING VANISHES";
	    PRINTSF GRES>>;
	    INTERR "NON-EXACT DIVISION DUE TO A LOG TERM"
	    END;
	RETURN A
   END;

SYMBOLIC PROCEDURE FQUOTF(A,B);
% INPUT: A AND B STANDARD QUOTIENTS WITH (A/B) AN EXACT;
% DIVISION WITH RESPECT TO THE VARIABLES IN ZLIST, ;
% BUT NOT NECESSARILY OBVIOUSLY SO. THE 'NON-OBVIOUS' PROBLEMS;
% WILL BE BECAUSE OF (E.G.) SQUARE-ROOT SYMBOLS IN B;
% OUTPUT: STANDARD QUOTIENT FOR (A/B);
% (PRINTS MESSAGE IF REMAINDER IS NOT 'CLEARLY' ZERO;
% A MUST NOT BE ZERO;
    BEGIN	  SCALAR T1;
	IF NULL A THEN INTERR "A=0 IN FQUOTF";
	T1:=QUOTF(A,B); %TRY IT THE EASY WAY;
	IF NOT NULL T1 THEN RETURN T1 ./ 1; %OK;
	RETURN DF2Q DFQUOTDF(F2DF A,F2DF B)
    END;

SYMBOLIC PROCEDURE DFQUOTDF1(A,B);
    BEGIN	SCALAR Q;
	IF NULL B THEN INTERR "ATTEMPT TO DIVIDE BY ZERO";
        Q:=SQRTLIST; %REMOVE SQRTS FROM DENOMINATOR, MAYBE; 
        WHILE NOT NULL Q DO BEGIN 
            SCALAR CONJ; 
            CONJ:=CONJSQRT(B,CAR Q); %CONJUGATE WRT GIVEN SQRT; 
            IF NOT (B=CONJ) THEN << 
                A:=MULTDF(A,CONJ); 
                B:=MULTDF(B,CONJ) >>; 
            Q:=CDR Q END; 
        Q:=DFQUOTDF2(A,B);
	RESIDUE:=REVERSEWOC RESIDUE;
	RETURN Q
    END;

SYMBOLIC PROCEDURE DFQUOTDF2(A,B);
%AS ABOVE BUT A AND B ARE DISTRIBUTED FORMS, AS IS THE RESULT;
    IF NULL A THEN NIL
    ELSE BEGIN SCALAR XD,LCD;
	XD:=XPDIFF(LPOW A,LPOW B);
	IF XD='FAILED THEN <<
	    XD:=LT A; A:=RED A;
	    RESIDUE:=XD .+ RESIDUE;
	    RETURN DFQUOTDF2(A,B) >>;
	LCD:=SUBS2Q MULTSQ(LC A,INVSQ LC B);
	IF NULL NUMR LCD THEN RETURN DFQUOTDF2(RED A,B);
	LCD := XD .* LCD;
	XD:=PLUSDF(A,MULTDF(NEGDF (LCD .+ NIL),B));
	IF XD AND (LPOW XD = LPOW A 
		   OR XPDIFF(LPOW XD,LPOW B) = 'FAILED)
	  THEN <<IF !*TRINT OR !*TRDIV
		   THEN <<PRINTC "DFQUOTDF TROUBLE:"; PRINTDF XD>>;
	         XD := ROOTEXTRACTDF XD;
		 IF !*TRINT OR !*TRDIV THEN PRINTDF XD>>;
	RETURN LCD .+ DFQUOTDF2(XD,B)
    END;

SYMBOLIC PROCEDURE ROOTEXTRACTDF U;
   IF NULL U THEN NIL
    ELSE BEGIN SCALAR V;
      V := RESIMP ROOTEXTRACTSQ LC U;
      RETURN IF NULL NUMR V THEN ROOTEXTRACTDF RED U
	      ELSE (LPOW U .* V) .+ ROOTEXTRACTDF RED U
    END;

SYMBOLIC PROCEDURE ROOTEXTRACTSQ U;
   IF NULL NUMR U THEN U
    ELSE ROOTEXTRACTF NUMR U ./ ROOTEXTRACTF DENR U;

SYMBOLIC PROCEDURE ROOTEXTRACTF V;
   IF DOMAINP V THEN V
    ELSE BEGIN SCALAR U,R,C,X,P;
      U := MVAR V;  P := LDEG V;
      R := ROOTEXTRACTF RED V;
      C := ROOTEXTRACTF LC V;
      IF NULL C THEN RETURN R
       ELSE IF ATOM U THEN RETURN (LPOW V .* C) .+ R
       ELSE IF CAR U EQ 'SQRT
	OR CAR U EQ 'EXPT AND EQCAR(CADDR U,'QUOTIENT)
	   AND CAR CDADDR U = 1 AND NUMBERP CADR CDADDR U
	THEN <<P := DIVIDE(P,IF CAR U EQ 'SQRT THEN 2
			      ELSE CADR CDADDR U);
      IF CAR P = 0 
        THEN RETURN IF NULL C THEN R ELSE (LPOW V .* C) .+ R
       ELSE IF NUMBERP CADR U
	THEN <<C := MULTD(CADR U ** CAR P,C); P := CDR P>>
       ELSE <<X := SIMPEXPT LIST(CADR U,CAR P);
	      IF DENR X = 1
		THEN <<C := MULTF(NUMR X,C); P := CDR P>>>>>>;
      RETURN IF P=0 THEN ADDF(C,R)
	      ELSE IF NULL C THEN R
	      ELSE ((U TO P) .* C) .+ R
   END;

PUT('DF,'SIMPFN,'SIMPDF!*);

SYMBOLIC PROCEDURE SIMPDF!* U;
  BEGIN SCALAR V,V1;
	V:=SIMPDF U;
	V1:=ROOTEXTRACTSQ V;
	IF NOT(V1=V) THEN RETURN RESIMP V1
	ELSE RETURN V
END;

SYMBOLIC PROCEDURE XPDIFF(A,B);
%RESULT IS LIST A-B, OR 'FAILED' IF A MEMBER OF THIS WOULD BE NEGATIVE;
    IF NULL A THEN IF NULL B THEN NIL
	ELSE INTERR "B TOO LONG IN XPDIFF"
    ELSE IF NULL B THEN INTERR "A TOO LONG IN XPDIFF"
    ELSE IF CAR B>CAR A THEN 'FAILED
    ELSE (LAMBDA R;
	IF R='FAILED THEN 'FAILED
	ELSE (CAR A-CAR B) . R) (XPDIFF(CDR A,CDR B));


SYMBOLIC PROCEDURE CONJSQRT(B,VAR); 
%SUBST(VAR=-VAR,B); 
    IF NULL B THEN NIL 
    ELSE CONJTERM(LPOW B,LC B,VAR) .+ CONJSQRT(RED B,VAR); 
 
SYMBOLIC PROCEDURE CONJTERM(XL,COEF,VAR); 
%DITTO BUT WORKING ON A TERM; 
    IF INVOLVESP(XL,VAR,ZLIST) THEN XL .* NEGSQ COEF 
    ELSE XL .* COEF; 
 
SYMBOLIC PROCEDURE INVOLVESP(XL,VAR,ZL); 
%CHECK IF EXPONENT LIST HAS NON-ZERO POWER FOR VARIABLE; 
    IF NULL XL THEN INTERR "VAR NOT FOUND IN INVOLVESP" 
    ELSE IF CAR ZL=VAR THEN (NOT ZEROP CAR XL) 
    ELSE INVOLVESP(CDR XL,VAR,CDR ZL); 


ENDMODULE;


MODULE DRIVER;

EXPORTS INTEGRATESQ,SIMPINT,PURGE,SIMPINT1;

IMPORTS ALGEBRAICCASE,ALGFNPL,FINDZVARS,GETVARIABLES,INTERR,PRINTSQ,
  TRANSCENDENTALCASE,VARSINLIST,KERNP,SIMPCAR,PREPSQ,MKSQ,SIMP,
   OPMTCH,FORMLNR;


%FORM IS   INT(EXPR,VAR,X1,X2,...);
%MEANING IS INTEGRATE EXPR WRT VAR, GIVEN THAT THE RESULT MAY;
%CONTAIN LOGS OF X1,X2,...;
% X1, ETC ARE INTENDED FOR USE WHEN THE SYSTEM HAS TO BE HELPED;
% IN THE CASE THAT EXPR IS ALGEBRAIC;
SYMBOLIC PROCEDURE SIMPINT U;
% Simplify an integral, links up with general prefix mode system;
    BEGIN SCALAR EXPRESSION,VARIABLE,TT,LOGLIST,W,!*GCD,!*MCD,!*EXP,
		 !*PURERISCH,!*SQRT,!*STRUCTURE;
% ARGUMENT IS A LIST OF TWO ELEMENTS, WHICH ARE PREFIX FORMS;
% OF THE INTEGRAND AND VARIABLE OF INTEGRATION;
    !*GCD:=T;
    !*MCD:=T;
    !*EXP:=T;
    !*SQRT:=T;
    !*STRUCTURE := T;
    VARIABLE:=CDR U;
    EXPRESSION:=SIMPP CAR U; %CONVERT INTEGRAND INTO A SQ;
    IF NULL VARIABLE THEN GO TO NOTENOUGHARGS;
    W:=CDR VARIABLE;
    VARIABLE:= !*Q2K SIMPP CAR VARIABLE; %CONVERT VARIABLE;
%NOW ARGUMENTS HAVE BEEN CHECKED. START WORK;
    LOGLIST:=MAPCAR(W,FUNCTION SIMPP);
    U:=ERRORSET('(INTEGRATESQ EXPRESSION VARIABLE LOGLIST),
		 NIL,!*BACKTRACE);
    IF NOT ATOM U THEN RETURN CAR U; %INTEGRATION OK;
    RETURN SIMPINT1(EXPRESSION . VARIABLE.W);
    % LEAVE IT FORMAL & LINEARISED;
NOTENOUGHARGS:	INTERR "NOT ENOUGH ARGS FOR INT";
TOOMANYARGS: INTERR "TOO MANY ARGS FOR INT"
    END;

SYMBOLIC PROCEDURE SIMPP U;
   %converts U to canonical form. Resimplifies if U is a *sq form;
   IF EQCAR(U,'!*SQ) THEN RESIMP CADR U ELSE SIMP U;

PUT('INT,'SIMPFN,'SIMPINT);


SYMBOLIC PROCEDURE INTEGRATESQ(INTEGRAND,VAR,XLOGS);
 BEGIN SCALAR VARLIST,ZLIST;
    IF !*TRINT THEN <<
	PRINTC "INTEGRAND IS...";
	PRINTSQ INTEGRAND >>;
    VARLIST:=GETVARIABLES INTEGRAND;
    VARLIST:=VARSINLIST(XLOGS,VARLIST); %IN CASE MORE EXIST IN XLOGS;
    ZLIST:=FINDZVARS(VARLIST,LIST VAR,VAR,NIL); %%IMPORTSANT KERNELS;
%the next section causes problems with nested exponentials or logs;
    BEGIN SCALAR OLDZLIST;
        WHILE OLDZLIST NEQ ZLIST DO <<
            OLDZLIST:=ZLIST;
	    FOREACH ZZ IN OLDZLIST DO
		ZLIST:=FINDZVARS(PSEUDODIFF(ZZ,VAR),ZLIST,VAR,T) >>
    END;
    IF !*TRINT  THEN <<
      PRINTC "WITH 'NEW' FUNCTIONS :";
      PRINT ZLIST >>;
    IF !*PURERISCH AND NOT ALLOWEDFNS ZLIST
      THEN RETURN SIMPINT1 (INTEGRAND . VAR.NIL);
      % IF IT IS NOT SUITABLE FOR RISCH;
    VARLIST:=PURGE(ZLIST,VARLIST);
% NOW ZLIST IS LIST OF THINGS THAT DEPEND ON X, AND VARLIST IS LIST;
% OF CONSTANT KERNELS IN INTEGRAND;
    RETURN TRANSCENDENTALCASE(INTEGRAND,VAR,XLOGS,ZLIST,VARLIST)
 END;

SYMBOLIC PROCEDURE PSEUDODIFF(A,VAR);
    IF ATOM A THEN NIL
    ELSE IF CAR A MEMQ '(EXPT PLUS TIMES QUOTIENT LOG SQRT)
	THEN BEGIN SCALAR AA,BB;
	    FOREACH ZZ IN CDR A DO <<
		BB:=PSEUDODIFF(ZZ,VAR);
		IF AA THEN AA:=BB . AA ELSE BB >>;
	    RETURN AA
	END
    ELSE LIST PREPSQ SIMPDF(LIST(A,VAR));

MKOP 'INT!*;

SYMBOLIC PROCEDURE SIMPINT1 U;
   BEGIN SCALAR V,!*SQRT;
      U := 'INT . PREPSQ CAR U . CDR U;
      IF (V := FORMLNR U) NEQ U
	THEN IF !*NOLNR THEN <<
		V:= SIMP SUBST('INT!*,'INT,V);
		RETURN REMAKESF NUMR V ./ REMAKESF DENR V>>
	      ELSE <<!*NOLNR:= NIL . !*NOLNR;
		     U:=ERRORSET(LIST('SIMP,MKQUOTE V),NIL,!*BACKTRACE);
		     IF PAIRP U THEN V:=CAR U;
		     !*NOLNR:= CDR !*NOLNR;
		     RETURN V>>;
      RETURN IF (V := OPMTCH U) THEN SIMP V ELSE MKSQ(U,1)
   END;

SYMBOLIC PROCEDURE REMAKESF U;
   %remakes standard form U, substituting operator INT for INT!*;
   IF DOMAINP U THEN U
    ELSE ADDF(MULTPF(IF EQCAR(MVAR U,'INT!*)
		       THEN MKSP('INT . CDR MVAR U,LDEG U)
		      ELSE LPOW U,REMAKESF LC U),
	       REMAKESF RED U);

SYMBOLIC PROCEDURE ALLOWEDFNS U;
IF NULL U
  THEN T
  ELSE IF ATOM CAR U OR
      FLAGP(CAAR U,'TRANSCENDENTAL)
    THEN ALLOWEDFNS CDR U
    ELSE NIL;


SYMBOLIC PROCEDURE PURGE(A,B);
    IF NULL A THEN B
    ELSE IF NULL B THEN NIL
    ELSE PURGE(CDR A,DELETE(CAR A,B));


ENDMODULE;


MODULE D3D4;

EXPORTS CUBIC,QUARTIC;

IMPORTS COVECDF,CUBEROOTF,NTH,FORCEAZERO,MAKEPOLYDF,MULTDF,MULTDFCONST,
   !*MULTF!*,NEGDF,PLUSDF,PRINTDF,PRINTSF,QUADRATIC,SQRTF,VP1,VP2,ADDF,
   NEGF;

%SPLITTING OF CUBICS AND QUARTICS;

SYMBOLIC PROCEDURE CUBIC(POL,VAR,RES);
%SPLIT THE UNIVARIATE (WRT Z-VARS) CUBIC POL, AT LEAST IF A;
%CHANGE OF ORIGIN PUTS IT IN THE FORM (X-A)**3-B=0;
    BEGIN	SCALAR A,B,C,D,V,SHIFT,P,Q,DSC;
	V:=COVECDF(POL,VAR,3);
	SHIFT:=FORCEAZERO(V,3); %MAKE COEFF X**2 VANISH;
				%ALSO CHECKS UNIVARIATE;
%	IF SHIFT='FAILED THEN GO TO PRIME;
	A:=GETV(V,3); B:=GETV(V,2); %=0, I HOPE!;
	C:=GETV(V,1); D:=GETV(V,0);
	IF !*TRINT THEN << PRINTC "CUBIC HAS COEFFICIENTS";
	    PRINTSF A; PRINTSF B;
	    PRINTSF C; PRINTSF D >>;
	IF NOT NULL C THEN <<
	    PRINTC "CUBIC TOO HARD TO SPLIT";
	    GO TO EXIT >>;
	A:=CUBEROOTF(A); %CAN'T EVER FAIL;
	D:=CUBEROOTF(D);
	IF !*TRINT THEN << PRINTC "CUBE ROOTS OF A AND D ARE";
	    PRINTSF A; PRINTSF D>>;
	%NOW A*(X+SHIFT)+D IS A FACTOR OF POL;
	%CREATE X+SHIFT IN P;
	P:=(VP2 ZLIST .* SHIFT) .+ NIL;
	P:=(VP1(VAR,1,ZLIST) .* (1 ./ 1)) .+ P; %(X+SHIFT);
	B:=NIL;
	B:=(VP2 ZLIST .* (D ./ 1)) .+ B;
	B:=PLUSDF(B,MULTDFCONST(A ./ 1,P));
	B:=MAKEPOLYDF B; %GET RID OF DENOMINATOR;
	IF !*TRINT THEN << PRINTC "ONE FACTOR OF THE CUBIC IS";
	    PRINTDF B >>;
	RES:=('LOG . B) . RES;
	%NOW FORM THE (QUADRATIC) COFACTOR;
	B:=(VP2 ZLIST .* (!*F2POL !*MULTF!*(D,D) ./ 1)) .+ NIL;
	B:=PLUSDF(B,MULTDFCONST(NEGF !*F2POL !*MULTF!*(A,D) ./ 1,P));
	B:=PLUSDF(B,MULTDFCONST(!*F2POL !*MULTF!*(A,A) ./ 1,
				MULTDF(P,P)));
	RETURN QUADRATIC(MAKEPOLYDF B,VAR,RES); %DEAL WITH WHAT IS LEFT;
   PRIME:
	PRINTC "THE FOLLOWING CUBIC DOES NOT SPLIT";
  EXIT:
	PRINTDF POL;
	RETURN ('LOG . POL) . RES
    END;

FLUID '(KNOWNDISCRIMSIGN);

SYMBOLIC PROCEDURE QUARTIC(POL,VAR,RES);
%SPLITS UNIVARIATE (WRT Z-VARS) QUARTICS THAT CAN BE WRITTEN;
%IN THE FORM (X-A)**4+B*(X-A)**2+C;
    BEGIN	SCALAR A,B,C,D,E,V,SHIFT,P,Q,P1,P2,DSC;
	V:=COVECDF(POL,VAR,4);
	SHIFT:=FORCEAZERO(V,4); %MAKE COEFF X**3 VANISH;
%	IF SHIFT='FAILED THEN GO TO PRIME;
	A:=GETV(V,4); B:=GETV(V,3); %=0, I HOPE!;
	C:=GETV(V,2); D:=GETV(V,1);
	E:=GETV(V,0);
	IF !*TRINT THEN << PRINTC "QUARTIC HAS COEFFICIENTS";
	    PRINTSF A; PRINTSF B;
	    PRINTSF C; PRINTSF D;
	    PRINTSF E >>;
	IF NOT NULL D THEN << PRINTC "QUARTIC TOO HARD TO SPLIT";
	    GO TO EXIT >>;
	B:=C; C:=E; %SQUASH UP THE NOTATION;
	IF KNOWNDISCRIMSIGN EQ 'NEGATIVE THEN GO TO COMPLEX;
	DSC := !*F2POL ADDF(MULTF(B,B),MULTF(-4,MULTF(A,C)));
	P2 := MINUSF C;
	IF NOT P2 AND MINUSF DSC THEN GO TO COMPLEX;
	P1 := NULL B OR MINUSF B;
	IF NOT P1 THEN IF P2 THEN P1 := T ELSE P2 := T;
	P1 := IF P1 THEN 'POSITIVE ELSE 'NEGATIVE;
	P2 := IF P2 THEN 'NEGATIVE ELSE 'POSITIVE;
	A := SQRTF A;
	DSC := SQRTF DSC;
	E := INVSQ(ADDF(A,A) ./ 1);
	D := MULTSQ(ADDF(B,NEGF DSC) ./ 1,E);
	E := MULTSQ(ADDF(B,DSC) ./ 1,E);
	IF !*TRINT
	  THEN <<PRINTC "QUADRATIC FACTORS WILL HAVE COEFFICIENTS";
		 PRINTSF A; PRINT 0; PRINTSQ D;
		 PRINTC "OR"; PRINTSQ E>>;
	P := (VP2 ZLIST .* SHIFT) .+ NIL;
	P := (VP1(VAR,1,ZLIST) .* (1 ./ 1)) .+ P; %(X+SHIFT);
	Q := MULTDF(P,P);   %SQUARE OF SAME;
	Q := MULTDFCONST(A ./ 1,Q);
	P := PLUSDF(Q,(VP2 ZLIST .* D) .+ NIL);
	Q := PLUSDF(Q,(VP2 ZLIST .* E) .+ NIL);
	IF !*TRINT
	  THEN <<PRINTC "ALLOWING FOR CHANGE OF ORIGIN:";
		 PRINTDF P; PRINTDF Q>>;
	KNOWNDISCRIMSIGN := P1;
	RES := QUADRATIC(P,VAR,RES);
	KNOWNDISCRIMSIGN := P2;
	RES := QUADRATIC(Q,VAR,RES);
	GO TO QUARTICDONE;
 COMPLEX:
	A:=SQRTF(A);
	C:=SQRTF(C);
	B:=ADDF(!*F2POL !*MULTF!*(2,!*MULTF!*(A,C)),NEGF B);
	B:=SQRTF B;
%NOW A*(X+SHIFT)**2 (+/-) B*(X+SHIFT) + C IS A FACTOR;
	IF !*TRINT
	  THEN << PRINTC "QUADRATIC FACTORS WILL HAVE COEFFICIENTS";
	    PRINTSF A; PRINTSF B; PRINTSF C>>;
	P:=(VP2 ZLIST .* SHIFT) .+ NIL;
	P:=(VP1(VAR,1,ZLIST) .* (1 ./ 1)) .+ P; %(X+SHIFT);
	Q:=MULTDF(P,P); %SQUARE OF SAME;
	P:=MULTDFCONST(B ./ 1,P);
	Q:=MULTDFCONST(A ./ 1,Q);
	Q:=PLUSDF(Q,(VP2 ZLIST .* (C ./ 1)) .+ NIL);
	IF !*TRINT THEN <<
	    PRINTC "ALLOWING FOR CHANGE OF ORIGIN, P (+/-) Q WITH P,Q=";
	    PRINTDF P; PRINTDF Q>>;
%NOW P+Q AND P-Q ARE THE FACTORS OF THE QUARTIC;
	KNOWNDISCRIMSIGN := 'NEGATIVE;
	RES:=QUADRATIC(PLUSDF(Q,P),VAR,RES);
	RES:=QUADRATIC(PLUSDF(Q,NEGDF P),VAR,RES);
 QUARTICDONE:
	KNOWNDISCRIMSIGN := NIL;
	IF !*TRINT THEN PRINTC "QUARTIC DONE";
	RETURN RES;
    PRIME:
	PRINTC "THE FOLLOWING QUARTIC DOES NOT SPLIT";
   EXIT:
	PRINTDF POL;
	RETURN ('LOG . POL) . RES
    END;


ENDMODULE;


MODULE FACTR;

EXPORTS INT!-FAC,VAR2DF;

IMPORTS CUBIC,DF2Q,F2DF,INTERR,MULTDF,PRINTDF,QUADRATIC,QUARTIC,UNIFAC,
   UNIFORM,VP1,VP2,SUB1;


SYMBOLIC PROCEDURE INT!-FAC X;
%INPUT: PRIMITIVE, SQUARE-FREE POLYNOMIAL (S.FORM);
%OUTPUT:
% LIST OF 'FACTORS' WRT ZLIST;
% EACH ITEM IN THIS LIST IS EITHER;
%     LOG . SQ;
% OR  ATAN . SQ;
% AND THESE LOGS AND ARCTANS ARE ALL THAT IS NEEDED IN THE;
% INTEGRATION OF 1/(ARGUMENT);
    BEGIN	  SCALAR RES,POL,DSET,VAR,DEGREE,VARS;
	POL:=F2DF X; %CONVERT TO DISTRIBUTED FORM;
	DSET:=DEGREESET(POL);
%NOW EXTRACT FACTORS OF THE FORM 'X' OR 'LOG(X)' ETC;
%THESE CORRESPOND TO ITEMS IN DSET WITH A NON-ZERO CDR;
	BEGIN    SCALAR ZL,DS;
	   ZL:=ZLIST; DS:=DSET;
	   WHILE NOT NULL DS DO <<
	       IF ONEP CDAR DS THEN <<
		   RES:=('LOG . VAR2DF(CAR ZL,1,ZLIST)) . RES;
			%RECORD IN ANSWER;
		   POL:=MULTDF(VAR2DF(CAR ZL,-1,ZLIST),POL);
			 %DIVIDE OUT;
		   IF !*TRINT THEN << PRINTC "TRIVIAL FACTOR FOUND";
		       PRINTDF CDAR RES>>;
		   RPLACA(DS,SUB1 CAAR DS . CDAR DS) >>
	       ELSE IF NULL ZEROP CDAR DS THEN
		  INTERR "REPEATED TRIVIAL FACTOR IN ARG TO FACTOR";
	       ZL:=CDR ZL; DS:=CDR DS >>;
	END; %SINGLE TERM FACTORS ALL REMOVED NOW;
	DSET:=MAPCAR(DSET,FUNCTION CAR); %GET LOWER BOUNDS;
	IF !*TRINT
	  THEN PRINTC ("UPPER BOUNDS OF REMAINING FACTORS ARE NOW: " .
			 DSET);
	IF DSET=VP2 ZLIST THEN GO TO FINISHED; %THING LEFT IS CONSTANT;
	BEGIN    SCALAR DS,ZL;
	    VAR:=CAR ZLIST; DEGREE:=CAR DSET;
	    IF NOT ZEROP DEGREE THEN VARS:=VAR . VARS;
	    DS:=CDR DSET; ZL:=CDR ZLIST;
	    WHILE NOT NULL DS DO <<
		IF NOT ZEROP CAR DS THEN <<
		    VARS:=CAR ZL . VARS;
		    IF ZEROP DEGREE OR DEGREE>CAR DS THEN <<
			VAR:=CAR ZL; DEGREE:=CAR DS >> >>;
		ZL:=CDR ZL; DS:=CDR DS >>
	END;
% NOW VAR IS VARIABLE THAT THIS POLY INVOLVES TO LOWEST DEGREE;
% DEGREE IS THE DEGREE OF THE POLY IN SAME VARIABLE;
	IF !*TRINT
	  THEN PRINTC ("BEST VAR IS " . VAR . "WITH EXPONENT " .
			 DEGREE);
	IF ONEP DEGREE THEN <<
	    RES:=('LOG . POL) . RES; %CERTAINLY IRREDUCIBLE;
	    IF !*TRINT
	      THEN << PRINTC "THE FOLLOWING IS CERTAINLY IRREDUCIBLE";
		PRINTDF POL>>;
	    GO TO FINISHED >>;
	IF DEGREE=2 THEN <<
	    IF !*TRINT THEN << PRINTC "QUADRATIC";
		PRINTDF POL>>;
	    RES:=QUADRATIC(POL,VAR,RES);
	    GO TO FINISHED >>;
	DSET:=UNIFORM(POL,VAR);
	IF NOT (DSET='FAILED) THEN <<
	    IF !*TRINT THEN << PRINTC "UNIVARIATE POLYNOMIAL";
		PRINTDF POL >>;
	    RES:=UNIFAC(DSET,VAR,DEGREE,RES);
	    GO TO FINISHED >>;
	IF NOT NULL CDR VARS THEN GO TO NASTY; %ONLY TRY UNIVARIATE NOW;
	IF DEGREE=3 THEN <<
	    IF !*TRINT THEN << PRINTC "CUBIC";
		PRINTDF POL>>;
	    RES:=CUBIC(POL,VAR,RES);
%	    IF !*OVERLAYMODE
%	      THEN EXCISE 'D3D4;
	    GO TO FINISHED >>;
	IF DEGREE=4 THEN <<
	    IF !*TRINT THEN << PRINTC "QUARTIC";
		PRINTDF POL>>;
	    RES:=QUARTIC(POL,VAR,RES);
%	    IF !*OVERLAYMODE
%	      THEN EXCISE 'D3D4;
	    GO TO FINISHED>>;
%ELSE ABANDON HOPE AND HAND BACK SOME RUBBISH.;
NASTY:
	RES:=('LOG . POL) . RES;
	PRINTC
	  "THE FOLLOWING POLYNOMIAL HAS NOT BEEN PROPERLY FACTORED";
	PRINTDF POL;
	GO TO FINISHED;


   FINISHED: %RES IS A LIST OF D.F. S AS REQUIRED;
	POL:=NIL; %CONVERT BACK TO STANDARD FORMS;
	WHILE NOT NULL RES DO
	    BEGIN	  SCALAR TYPE,ARG;
	    TYPE:=CAAR RES; ARG:=CDAR RES;
	    ARG:=DF2Q ARG;
	    IF TYPE='LOG THEN RPLACD(ARG,1);
	    POL:=(TYPE . ARG) . POL;
	    RES:=CDR RES END;
	RETURN POL
    END;


SYMBOLIC PROCEDURE VAR2DF(VAR,N,ZLIST);
    ((VP1(VAR,N,ZLIST) .* (1 ./ 1)) .+ NIL);

SYMBOLIC PROCEDURE DEGREESET POL;
%FINDS DEGREE BOUNDS FOR ALL VARS IN DISTRIBTED FORM POLY;
    DEGREESUB(DBL LPOW POL,RED POL);

SYMBOLIC PROCEDURE DBL X;
% CONVERTS LIST OF X INTO LIST OF (X . X);
    IF NULL X THEN NIL
    ELSE (CAR X . CAR X) . DBL CDR X;

SYMBOLIC PROCEDURE DEGREESUB(CUR,POL);
% UPDATE DEGREE BOUNDS 'CUR' TO INCLUDE INFO ABOUT POL;
    <<
	WHILE NOT NULL POL DO <<
	    CUR:=DEGREESUB1(CUR,LPOW POL);
	    POL:=RED POL >>;
	CUR >>;

SYMBOLIC PROCEDURE DEGREESUB1(CUR,NXT);
%MERGE INFORMATION FROM EXPONENT SET NEXT INTO CUR;
    IF NULL CUR THEN NIL
    ELSE DEGREESUB2(CAR CUR,CAR NXT) . DEGREESUB1(CDR CUR,CDR NXT);

SYMBOLIC PROCEDURE DEGREESUB2(TWO,ONE);
    MAX(CAR TWO,ONE) . MIN(CDR TWO,ONE);


ENDMODULE;


MODULE IBASICS;

EXPORTS PARTIALDIFF,PRINTDF,PRINTSQ,RATIONALINTEGRATE,PRINTSF,INTERR;

IMPORTS DF2PRINTFORM,SQPRINT,VARSINSF,TERPRI!*,ADDSQ,MULTSQ,MULTD,MKSP;


%PRINT STANDARD QUOTIENT (RATIONAL FUNCTION);
% CRUDE EQUIVALENT TO PRINTSF NUMR U: "/": PRINTSF DENO U;

SYMBOLIC PROCEDURE PRINTSQ U;
   BEGIN
      TERPRI!*(T); %START ON A NEW LINE;
      SQPRINT U; %LOGICAL PRINT ROUTINE;
      TERPRI!*(T)
   END;

% PRINT STANDARD FORM (POLYNOMIAL);
FLUID '(U!*); %NEEDED BECAUSE OF THE ERRORSET;

SYMBOLIC PROCEDURE PRINTSF U!*;
    IF NULL U!* THEN PRINT 0
   ELSE BEGIN    SCALAR W;
    W:=ERRORSET('(PROG NIL (TERPRI!* T)
	    (XPRINF U!* NIL NIL) (TERPRI!* T)),2,!*BACKTRACE);
    IF NOT ATOM W THEN RETURN CAR W;
    PRINTC "REDUCE PRINTING FAILED ON STANDARD FORM";
    PRINT U!*;
    TERPRI!*(T);
    RETURN U!*
   END;
UNFLUID '(U!*);

SYMBOLIC PROCEDURE PRINTDF U;
% PRINT DISTRIBUTED FORM VIA CHEAP CONVERSION TO REDUCE STRUCTURE;
    BEGIN SCALAR !*GCD;
       PRINTSF DF2PRINTFORM U;
    END;


SYMBOLIC PROCEDURE INTERR MESS;
   BEGIN
     PRINTC "INTEGRATION PACKAGE ERROR";
     PRINTC MESS;
     ERROR1()
   END;


SYMBOLIC PROCEDURE RATIONALINTEGRATE(X,VAR);
    BEGIN	  SCALAR N,D;
      N:=NUMR X; D:=DENR X;
      IF NOT VAR MEMBER VARSINSF(D,NIL) THEN
	    RETURN SUBS2Q MULTSQ(POLYNOMIALINTEGRATE(N,VAR),1 ./ D);
      INTERR "RATIONAL INTEGRATION NOT CODED YET"
    END;


% INTEGRATE STANDARD FORM. RESULT IS STANDARD QUOTIENT;
SYMBOLIC PROCEDURE POLYNOMIALINTEGRATE(X,V);
    IF NULL X THEN NIL ./ 1
    ELSE IF ATOM X THEN ((MKSP(V,1) .* 1) .+ NIL) ./ 1
    ELSE BEGIN    SCALAR R;
      R:=POLYNOMIALINTEGRATE(RED X,V); % DEAL WITH REDUCTUM;
      IF V=MVAR X THEN BEGIN    SCALAR DEGREE,NEWLT;
	 DEGREE:=1+TDEG LT X;
	 NEWLT:=((MKSP(V,DEGREE) .* LC X) .+ NIL) ./ 1; % UP EXPONENT;
	 R:=ADDSQ(MULTSQ(NEWLT,1 ./ DEGREE),R)
	 END
      ELSE BEGIN	 SCALAR NEWTERM;
	NEWTERM:=(((LPOW X) .* 1) .+ NIL) ./ 1;
	NEWTERM:=MULTSQ(NEWTERM,POLYNOMIALINTEGRATE(LC X,V));
	R:=ADDSQ(R,NEWTERM)
	END;
      RETURN SUBS2Q R
    END;

% PARTIAL DIFFERENTIATION OF P WRT V - P IS S.F. AS IS RESULT;
SYMBOLIC PROCEDURE PARTIALDIFF(P,V);
    IF ATOM P THEN NIL
    ELSE
	IF V=MVAR P THEN
	    (LAMBDA X; IF X=1 THEN LC P
	     ELSE ((MKSP(V,X-1) .* MULTD(X,LC P))
			 .+ PARTIALDIFF(RED P,V)))
	    (TDEG LT P)
	ELSE
	    (LAMBDA X; IF NULL X THEN PARTIALDIFF(RED P,V)
	     ELSE ((LPOW P .* X) .+ PARTIALDIFF(RED P,V)))
	    (PARTIALDIFF(LC P,V));

PUT('PDIFF,'SIMPFN,'SIMPPDIFF);


ENDMODULE;


MODULE JPATCHES;

EXPORTS !*MULTF!*;

IMPORTS !*MULTF!*SQRT,SIMPSQRTI,RETIMES,MULTSQ,SIMPEXPT,INVSQ,MKSQ,XN,
   FLATTEN,MKSPM,MKSP,EXPTF,SIMP,GCDN,ADDF,ORDOP,NONCOMP,MKSFPF,
   MULTD,DOMAINP;


%SYMBOLIC PROCEDURE SIMPX1(U,M,N);
%   %U,M AND N ARE PREFIX EXPRESSIONS;
%   %VALUE IS THE STANDARD QUOTIENT EXPRESSION FOR U**(M/N);
%   BEGIN SCALAR FLG,Z;
%	IF NULL FRLIS!* OR NULL XN(FRLIS!*,FLATTEN (M . N))
%	  THEN GO TO A;
%	EXPTP!* := T;
%	RETURN !*K2Q LIST('EXPT,U,IF N=1 THEN M
%				   ELSE LIST('QUOTIENT,M,N));
%    A:	IF NUMBERP M AND FIXP M THEN GO TO E
%	 ELSE IF ATOM M THEN GO TO B
%	 ELSE IF CAR M EQ 'MINUS THEN GO TO MNS
%	 ELSE IF CAR M EQ 'PLUS THEN GO TO PLS
%	 ELSE IF CAR M EQ 'TIMES AND NUMBERP CADR M AND FIXP CADR M
%		AND NUMBERP N
%	  THEN GO TO TMS;
%    B:	Z := 1;
%    C:	IF ATOM U AND NOT NUMBERP U THEN FLAG(LIST U,'USED!*);
%	U := LIST('EXPT,U,IF N=1 THEN M ELSE LIST('QUOTIENT,M,N));
%	IF NOT U MEMBER EXPTL!* THEN EXPTL!* := U . EXPTL!*;
%    D:	RETURN MKSQ(U,IF FLG THEN -Z ELSE Z); %U IS ALREADY IN LOWEST
%	%TERMS;
%    E:	IF NUMBERP N AND FIXP N THEN GO TO INT;
%	Z := M;
%	M := 1;
%	GO TO C;
%    MNS: M := CADR M;
%	IF !*MCD THEN RETURN INVSQ SIMPX1(U,M,N);
%	FLG := NOT FLG;
%	GO TO A;
%    PLS: Z := 1 ./ 1;
%    PL1: M := CDR M;
%	IF NULL M THEN RETURN Z;
%	Z := MULTSQ(SIMPEXPT LIST(U,
%			LIST('QUOTIENT,IF FLG THEN LIST('MINUS,CAR M)
%					ELSE CAR M,N)),
%		    Z);
%	GO TO PL1;
%    TMS: Z := GCDN(N,CADR M);
%	N := N/Z;
%	Z := CADR M/Z;
%	M := RETIMES CDDR M;
%	GO TO C;
%    INT:Z := DIVIDE(M,N);
%	IF CDR Z<0 THEN Z:= (CAR Z - 1) . (CDR Z+N);
%	IF CDR Z=0
%	  THEN RETURN SIMPEXPT LIST(U,CAR Z);
%	IF N=2 AND !*SQRT
%	  THEN RETURN MULTSQ(SIMPEXPT LIST(U,CAR Z),
%			     SIMPSQRTI U);
%	RETURN MULTSQ(SIMPEXPT LIST(U,CAR Z),
%			MKSQ(LIST('EXPT,U,LIST('QUOTIENT,1,N)),CDR Z))
%   END;


ENDMODULE;


MODULE KRON;

EXPORTS LINFAC,QUADFAC;

IMPORTS EVALAT,LINETHROUGH,QUADTHROUGH,TESTDIV;

%KRONEKER FACTORIZATION FOR UNIVARIATE POLYS OVER THE INTEGERS;
%ONLY LINEAR AND QUADRATIC FACTORS ARE FOUND HERE;



SYMBOLIC PROCEDURE LINFAC(W);
    TRYKR(W,'(0 1));

SYMBOLIC PROCEDURE QUADFAC(W);
    TRYKR(W,'(-1 0 1));


SYMBOLIC PROCEDURE TRYKR(W,POINTS);
%LOOK FOR FACTOR OF W BY EVALUATION AT (POINTS) AND USE OF;
% INTERPOLATE. RETURN (FAC . COFAC) WITH FAC=NIL IF NONE;
%FOUND AND COFAC=NIL IF NOTHING WORTHWHILE IS LEFT;
    BEGIN	  SCALAR VALUES,ATTEMPT;
	IF NULL W THEN RETURN NIL . NIL;
	IF  (LENGTH POINTS > CAR W) THEN RETURN W . NIL;
%THAT SAYS IF W IS ALREADY TINY, IT IS ALREADY FACTORED;
	VALUES:=MAPCAR(POINTS,FUNCTION (LAMBDA X;
	   EVALAT(W,X)));
	IF !*TRINT THEN << PRINTC ("AT X= " . POINTS);
	    PRINTC ("P(X)= " . VALUES)>>;
	IF 0 MEMBER VALUES THEN GO TO LUCKY; %(X-1) IS A FACTOR!;
	VALUES:=MAPCAR(VALUES,FUNCTION ZFACTORS);
	RPLACD(VALUES,MAPCAR(CDR VALUES,FUNCTION (LAMBDA Y;
	    APPEND(Y,MAPCAR(Y,FUNCTION MINUS)))));
	IF !*TRINT THEN <<PRINTC "POSSIBLE FACTORS GO THROUGH SOME OF";
	    PRINT VALUES>>;
	ATTEMPT:=SEARCH4FAC(W,VALUES,NIL);
	IF NULL ATTEMPT THEN ATTEMPT:=NIL . W;
	RETURN ATTEMPT;
  LUCKY: %HERE (X-1) IS A FACTOR BECAUSE P(0) OR P(1) OR P(-1);
	 %VANISHED AND CASES P(0), P(-1) WILL HAVE BEEN REMOVED;
	 %ELSEWHERE;
	ATTEMPT:='(1 1 -1); %THE FACTOR;
	RETURN ATTEMPT . TESTDIV(W,ATTEMPT)
    END;

SYMBOLIC PROCEDURE SEARCH4FAC(W,VALUES,CV);
%COMBINATORIAL SEARCH. CV GETS CURRENT SELECTED VALUE-SET;
%RETURNS NIL IF FAILS, ELSE FACTOR . COFACTOR;
    IF NULL VALUES THEN TRYFACTOR(W,CV)
    ELSE BEGIN    SCALAR FF,Q;
	FF:=CAR VALUES; %TRY ALL VALUES HERE;
 LOOP:	IF NULL FF THEN RETURN NIL; %NO FACTOR FOUND;
	Q:=SEARCH4FAC(W,CDR VALUES,(CAR FF) . CV);
	IF NULL Q THEN << FF:=CDR FF; GO TO LOOP>>;
	RETURN Q
    END;

SYMBOLIC PROCEDURE TRYFACTOR(W,CV);
%TESTS IF CV REPRESENTS A FACTOR OF W;
    BEGIN	  SCALAR FF,Q;
	IF NULL CDDR CV THEN FF:=LINETHROUGH(CADR CV,CAR CV)
	ELSE FF:=QUADTHROUGH(CADDR CV,CADR CV,CAR CV);
	IF FF='FAILED THEN RETURN NIL; %IT DOES NOT INTERPOLATE;
	Q:=TESTDIV(W,FF);
	IF Q='FAILED THEN RETURN NIL; %NOT A FACTOR;
	RETURN FF . Q
    END;


ENDMODULE;


MODULE LOWDEG;

EXPORTS FORCEAZERO,MAKEPOLYDF,QUADRATIC,COVECDF,EXPONENTDF;

IMPORTS DFQUOTDF,GCDF,INTERR,MINUSDFP,MULTDF,MULTDFCONST,!*MULTF!*,
   NEGSQ,MINUSP,PRINTSQ,MULTSQ,INVSQ,PNTH,NTH,MKNILL,
   NEGDF,PLUSDF,PRINTDF,PRINTSQ,QUOTF,SQRTDF,VAR2DF,VP2,ADDSQ,SUB1;

%SPLITTING OF LOW DEGREE POLYNOMIALS;

SYMBOLIC PROCEDURE COVECDF(POL,VAR,DEGREE);
%EXTRACT COEFFICIENTS OF POLYNOMIAL WRT VAR, GIVEN A DEGREE-BOUND
% DEGREE;
%RESUL IS A LISP VECTOR;
    BEGIN	  SCALAR I,V,X,W;
	W:=POL;
	V:=MKVECT(DEGREE);
	WHILE NOT NULL W DO <<
	    X:=EXPONENTOF(VAR,LPOW W,ZLIST);
	    IF (X<0) OR (X>DEGREE) THEN INTERR "BAD DEGREE IN COVECDF";
	    PUTV(V,X,LT W . GETV(V,X));
	    W:=RED W >>;
	FOR I:=0:DEGREE DO PUTV(V,I,MULTDF(REVERSEWOC GETV(V,I),
	    VAR2DF(VAR,-I,ZLIST)));
	RETURN V
    END;

SYMBOLIC PROCEDURE QUADRATIC(POL,VAR,RES);
%ADD IN TO RES LOGS OR ARCTANS CORRESPONDING TO SPLITTING THE
% POLYNOMIAL;
% POL GIVEN THAT IT IS QUADRATIC WRT VAR;
%;
%DOES NOT ASSUME POL IS UNIVARIATE;
    BEGIN	SCALAR A,B,C,W,DISCRIM;
	 W:=COVECDF(POL,VAR,2);
	 A:=GETV(W,2); B:=GETV(W,1); C:=GETV(W,0);
% THAT SPLIT THE QUADRATIC UP TO FIND THE COEFFICIENTS A,B,C;
	IF !*TRINT THEN << PRINTC "A="; PRINTDF A;
	    PRINTC "B="; PRINTDF B;
	    PRINTC "C="; PRINTDF C>>;
	DISCRIM:=PLUSDF(MULTDF(B,B),
	    MULTDFCONST((-4) . 1,MULTDF(A,C)));
	IF !*TRINT THEN << PRINTC "DISCRIMINANT IS";
	    PRINTDF DISCRIM>>;
	IF NULL DISCRIM THEN INTERR "DISCRIM=0 IN QUADRATIC";
	IF KNOWNDISCRIMSIGN
	  THEN <<IF KNOWNDISCRIMSIGN EQ 'NEGATIVE THEN GO TO ATANCASE>>
	 ELSE IF (NOT CLOGFLAG) AND (MINUSDFP DISCRIM)
	  THEN GO TO ATANCASE;
	DISCRIM:=SQRTDF(DISCRIM);
	IF DISCRIM='FAILED THEN GO TO NOFACTORS;
	IF !*TRINT THEN << PRINTC "SQUARE-ROOT IS";
	    PRINTDF DISCRIM>>;
	W:=VAR2DF(VAR,1,ZLIST);
	W:=MULTDF(W,A);
	B:=MULTDFCONST(1 ./ 2,B);
	DISCRIM:=MULTDFCONST(1 ./ 2,DISCRIM);
	W:=PLUSDF(W,B); %A*X+B/2;
	A:=PLUSDF(W,DISCRIM); B:=PLUSDF(W,NEGDF(DISCRIM));
	IF !*TRINT THEN << PRINTC "FACTORS ARE";
	    PRINTDF A; PRINTDF B>>;
	RETURN ('LOG . A) . ('LOG . B) . RES;
ATANCASE:
	DISCRIM:=SQRTDF NEGDF DISCRIM; %SQRT(4*A*C-B**2) THIS TIME!;
	IF DISCRIM='FAILED THEN GO TO NOFACTORS; %SQRT DID NOT EXIST?;
	RES := ('LOG . POL) . RES; %ONE PART OF THE ANSWER;
	A:=MULTDF(A,VAR2DF(VAR,1,ZLIST));
	A:=PLUSDF(B,MULTDFCONST(2 ./ 1,A));
	A:=DFQUOTDF(A,DISCRIM); %ASSUMES DIVISION IS EXACT;
	RETURN ('ATAN . A) . RES;
NOFACTORS:
	PRINTC "THE FOLLOWING QUADRATIC DOES NOT SEEM TO FACTOR";
	PRINTDF POL;
	RETURN ('LOG . POL) . RES
    END;

SYMBOLIC PROCEDURE EXPONENTOF(VAR,L,ZL);
    IF NULL ZL THEN INTERR "VAR NOT FOUND IN EXPONENTOF"
    ELSE IF VAR=CAR ZL THEN CAR L
    ELSE EXPONENTOF(VAR,CDR L,CDR ZL);


SYMBOLIC PROCEDURE DF2SF A;
    IF NULL A THEN NIL
    ELSE IF ((NULL RED A) AND
	(ONEP DENR LC A) AND
	(LPOW A=VP2 ZLIST)) THEN NUMR LC A
    ELSE INTERR "NASTY CUBIC OR QUARTIC";



SYMBOLIC PROCEDURE MAKEPOLYDF P;
%MULTIPLY DF BY LCM OF DENOMINATORS OF ALL COEFFICIENT DENOMINATORS;
    BEGIN	SCALAR H,W;
	IF NULL(W:=P) THEN RETURN NIL; %POLY IS ZERO ALREADY;
	H:=DENR LC W; %A GOOD START;
	W:=RED W;
	WHILE NOT NULL W DO <<
	    H:=QUOTF(!*MULTF!*(H,DENR LC W),GCDF(H,DENR LC W));
	    W:=RED W >>;
	%H IS NOW LCM OF DENOMINATORS;
	RETURN MULTDFCONST(!*F2POL H ./ 1,P)
    END;


SYMBOLIC PROCEDURE FORCEAZERO(P,N);
%SHIFT POLYNOMIAL P SO THAT COEFF OF X**(N-1) VANISHES;
%RETURN THE AMOUNT OF THE SHIFT, UPDATE (VECTOR) P;
    BEGIN	SCALAR R,I,W;
	FOR I:=0:N DO PUTV(P,I,DF2SF GETV(P,I)); %CONVERT TO POLYS;
	R:=GETV(P,N-1);
	IF NULL R THEN RETURN NIL ./ 1; %ALREADY ZERO;
	R:= SUBS2Q MULTSQ(R ./ 1,INVSQ(!*MULTF!*(N,GETV(P,N)) ./ 1));
			%THE SHIFT AMOUNT;
%NOW I HAVE TO SET P:=SUBST(X-R,X,P) AND THEN REDUCE TO SF AGAIN;
	IF !*TRINT THEN << PRINTC "SHIFT IS BY ";
	    PRINTSQ R>>;
	W:=MKVECT(N); %WORKSPACE VECTOR;
	FOR I:=0:N DO PUTV(W,I,NIL ./ 1); %ZERO IT;
	I:=N;
	WHILE NOT MINUSP I DO <<
	    MULVECBYXR(W,NEGSQ R,N); %W:=(X-R)*W;
	    PUTV(W,0,ADDSQ(GETV(W,0),GETV(P,I) ./ 1));
	    I:=I-1 >>;
	IF !*TRINT THEN << PRINTC "SQ SHIFTED POLY IS";
	    PRINT W>>;
	FOR I:=0:N DO PUTV(P,I,GETV(W,I));
	W:=DENR GETV(P,0);
	FOR I:=1:N DO W:=QUOTF(!*MULTF!*(W,DENR GETV(P,I)),
	    GCDF(W,DENR GETV(P,I)));
	FOR I:=0:N DO PUTV(P,I,NUMR SUBS2Q MULTSQ(GETV(P,I),W ./ 1));
	W:=GETV(P,0);
	FOR I:=1:N DO W:=GCDF(W,GETV(P,I));
	IF NOT (W=1) THEN
	    FOR I:=0:N DO PUTV(P,I,QUOTF(GETV(P,I),W));
	IF !*TRINT THEN << PRINTC "FINAL SHIFTED POLY IS ";
	    PRINT P>>;
	RETURN R
    END;

SYMBOLIC PROCEDURE MULVECBYXR(W,R,N);
%W IS A VECTOR REPRESENTING A POLY OF DEGREE N;
%MULTIPLY IT BY (X+R);
    BEGIN	SCALAR I,IM1;
	I:=N;
	IM1:=SUB1 I;
	WHILE NOT MINUSP IM1 DO <<
	    PUTV(W,I,SUBS2Q ADDSQ(GETV(W,IM1),MULTSQ(R,GETV(W,I))));
	    I:=IM1; IM1:=SUB1 I >>;
	PUTV(W,0,SUBS2Q MULTSQ(GETV(W,0),R));
	RETURN W
    END;




ENDMODULE;


MODULE REFORM;

EXPORTS LOGSTOSQ,SUBSTINULIST;

IMPORTS PREPSQ,MKSP,NTH,MULTSQ,ADDSQ,DOMAINP,INVSQ,PLUSDF;

SYMBOLIC PROCEDURE SUBSTINULIST ULIST;
% Substitutes for the C-constants in the values of the U's given in;
% ULIST. Result is a D.F.;
   IF NULL ULIST THEN NIL
   ELSE BEGIN SCALAR TEMP,LCU;
      LCU:=LC ULIST;
      TEMP:=EVALUATEUCONST NUMR LCU;
      IF NULL NUMR TEMP THEN TEMP:=NIL
      ELSE TEMP:=((LPOW ULIST) .*
	SUBS2Q MULTSQ(TEMP,INVSQ(DENR LCU ./ 1))) .+ NIL;
      RETURN PLUSDF(TEMP,SUBSTINULIST RED ULIST)
   END;

SYMBOLIC PROCEDURE EVALUATEUCONST COEFFT;
% Substitutes for the C-constants into COEFFT (=S.F.). Result is S.Q.;
    IF NULL COEFFT OR DOMAINP COEFFT THEN COEFFT ./ 1
    ELSE BEGIN SCALAR TEMP;
      IF NULL(TEMP:=ASSOC(MVAR COEFFT,CMAP)) THEN
	    TEMP:=(!*P2F LPOW COEFFT) ./ 1
      ELSE TEMP:=GETV(CVAL,CDR TEMP);
      TEMP:=MULTSQ(TEMP,EVALUATEUCONST(LC COEFFT));
      RETURN SUBS2Q ADDSQ(TEMP,EVALUATEUCONST(RED COEFFT))
    END;

SYMBOLIC PROCEDURE LOGSTOSQ;
% Converts LOGLIST to sum of the log terms as a S.Q.;
   BEGIN SCALAR LGLST,LOGSQ,I,TEMP;
      I:=1;
      LGLST:=LOGLIST;
      LOGSQ:=NIL ./ 1;
LOOP: IF NULL LGLST THEN RETURN LOGSQ;
      TEMP:=CDDR CAR LGLST;
      IF !*TRINT
	THEN << PRINTC "Standard Form ARG FOR ADDITIONAL LOG ETC =";
	  PRINT TEMP >>;
      IF NOT (CAAR LGLST='IDEN) THEN <<
	  TEMP:=PREPSQ TEMP; %CONVERT TO PREFIX FORM;
	  TEMP:=LIST(CAAR LGLST,TEMP); %FUNCTION NAME;
	  TEMP:=((MKSP(TEMP,1) .* 1) .+ NIL) ./ 1 >>;
      TEMP:=MULTSQ(TEMP,GETV(CVAL,I));
      LOGSQ:= SUBS2Q ADDSQ(TEMP,LOGSQ);
      LGLST:=CDR LGLST;
      I:=I+1;
      GO TO LOOP
   END;

ENDMODULE;


MODULE SIMPLOG;

EXPORTS SIMPLOG,SIMPLOGSQ;

IMPORTS QUOTF,PREPF,MKSP,SIMP!*,MULTSQ,SIMPTIMES,ADDSQ,MINUSF,NEGF,
   ADDF,COMFAC,NEGSQ,MK!*SQ,CARX;

SYMBOLIC PROCEDURE SIMPLOG(EXXPR);
 SIMPLOGI(CARX(EXXPR,'LOG));


SYMBOLIC PROCEDURE SIMPLOGI(SQ);
BEGIN
   IF ATOM SQ
     THEN GO TO SIMPLIFY;
   IF CAR SQ EQ 'TIMES
     THEN RETURN ADDSQ(SIMPLOGI CADR SQ,SIMPLOGI CADDR SQ);
   IF CAR SQ EQ 'QUOTIENT
     THEN RETURN ADDSQ(SIMPLOGI CADR SQ,
		       NEGSQ SIMPLOGI CADDR SQ);
   IF CAR SQ EQ 'EXPT
     THEN RETURN SIMPTIMES LIST(CADDR SQ,
				MK!*SQ SIMPLOGI CADR SQ);
  IF CAR SQ = '!*SQ
    THEN RETURN SIMPLOGSQ CADR SQ;
 SIMPLIFY:
   SQ:=SIMP!* SQ;
  RETURN SIMPLOGSQ SQ
  END;


SYMBOLIC PROCEDURE SIMPLOGSQ SQ;
ADDSQ((SIMPLOG2 NUMR SQ),NEGSQ(SIMPLOG2 DENR SQ));


 SYMBOLIC PROCEDURE SIMPLOG2(SF);
 IF ATOM SF
   THEN IF NULL SF
     THEN REDERR "LOG 0 FORMED"
     ELSE IF NUMBERP SF
       THEN IF SF IEQUAL 1
	 THEN NIL ./ 1
	 ELSE IF SF IEQUAL 0
	   THEN REDERR "LOG 0 FORMED"
	   ELSE((MKSP(LIST('LOG,SF),1) .* 1) .+ NIL) ./ 1
       ELSE FORMLOG(SF)
   ELSE BEGIN
     SCALAR FORM;
     FORM:=COMFAC SF;
     IF NOT NULL CAR FORM
       THEN RETURN ADDSQ(FORMLOG(FORM .+ NIL),
			 SIMPLOG2 QUOTF(SF,FORM .+ NIL));
     % WE HAVE KILLED COMMON POWERS;
     FORM:=CDR FORM;
     IF FORM NEQ 1
       THEN RETURN ADDSQ(SIMPLOG2 FORM,
			  SIMPLOG2 QUOTF(SF,FORM));
     % REMOVE A COMMON FACTOR FROM THE SF;
     RETURN (FORMLOG SF)
     END;


 SYMBOLIC PROCEDURE FORMLOG(SF);
 IF (NULL RED SF)
   THEN  IF EQCAR(MVAR SF,'EXPT)
     THEN ADDSQ(SIMPLOG2 LC SF,
		SUBS2Q MULTSQ(SIMPLOGI MVAR SF,SIMP!* LDEG SF))
     ELSE IF (LC SF IEQUAL 1) AND (LDEG SF IEQUAL 1)
       THEN ((MKSP(LIST('LOG,MVAR SF),1) .* 1) .+ NIL) ./ 1
       ELSE ADDSQ(SIMPTIMES LIST(LIST('LOG,MVAR SF),LDEG SF),
		SIMPLOG2 LC SF)
   ELSE IF MINUSF SF
     THEN ADDF((MKSP(LIST('LOG,-1),1) .* 1) .+ NIL,
	       FORMLOG2 NEGF SF) ./ 1
     ELSE (FORMLOG2 SF) ./ 1;


SYMBOLIC PROCEDURE FORMLOG2 SF;
((MKSP(LIST('LOG,PREPF SF),1) .* 1) .+ NIL);


ENDMODULE;


MODULE SIMPSQRT;

SYMBOLIC PROCEDURE SIMPSQRTSQ SQ;
(SIMPSQRT2 NUMR SQ) ./ (SIMPSQRT2 DENR SQ);

 SYMBOLIC PROCEDURE SIMPSQRT2(SF);
 IF ATOM SF
   THEN IF NULL SF
     THEN NIL
     ELSE IF NUMBERP SF
       THEN IF MINUSP SF
	 THEN !*F2POL !*MULTF!*(SIMPSQRT2 (-SF),
		    (MKSP(MKSQRT(-1),1) .* 1) .+ NIL)
	 ELSE BEGIN
	   SCALAR N;
	   N:=SQRT SF;
	   IF IDP N
	     THEN RETURN (MKSP(MKSQRT!* SF,1) .* 1) .+ NIL
	     ELSE RETURN N
	   END
     ELSE FORMSQRT(SF)
   ELSE BEGIN
     SCALAR FORM;
     FORM:=COMFAC SF;
     IF NOT NULL CAR FORM
       THEN RETURN !*F2POL !*MULTF!*(FORMSQRT(FORM .+ NIL),
			 SIMPSQRT2 QUOTF(SF,FORM .+ NIL));
     % WE HAVE KILLED COMMON POWERS;
     FORM:=CDR FORM;
     IF FORM NEQ 1
       THEN RETURN !*F2POL !*MULTF!*(SIMPSQRT2 FORM,
			  SIMPSQRT2 QUOTF(SF,FORM));
     % REMOVE A COMMON FACTOR FROM THE SF;
     RETURN FORMSQRT SF
     END;


 SYMBOLIC PROCEDURE FORMSQRT(SF);
	%Is *F2POL really necessary here??;
 IF (NULL RED SF)
   THEN IF (LC SF IEQUAL 1) AND (LDEG SF IEQUAL 1)
     THEN (MKSP(MKSQRT!* MVAR SF,1) .* 1) .+ NIL
    ELSE !*F2POL
      !*MULTF!*(NUMR SIMPEXPT(LIST(MKSQRT!* MVAR SF,LDEG SF)),
		SIMPSQRT2 LC SF)
   ELSE (MKSP(MKSQRT!* SF,1) .* 1) .+ NIL;

SYMBOLIC PROCEDURE MKSQRT!* U;
   IF SFP U THEN MKSQRT !*F2A U ELSE MKSQRT U;

ALGEBRAIC;
% OPERATOR SQRT;
 SYMBOLIC;
% DEFLIST ('((SQRT (((X) QUOTIENT (SQRT X) (TIMES 2 X))))),'DFN);

SYMBOLIC PROCEDURE SIMPSQRTI SQ;
BEGIN
   IF ATOM SQ
     THEN IF NUMBERP SQ
       THEN RETURN (SIMPSQRT2 SQ) ./ 1
       ELSE RETURN ((MKSP(MKSQRT SQ,1) .* 1) .+ NIL) ./ 1;
   IF CAR SQ EQ 'TIMES
     THEN RETURN SUBS2Q MULTSQ(SIMPSQRTI CADR SQ,SIMPSQRTI CADDR SQ);
   IF CAR SQ EQ 'QUOTIENT
     THEN RETURN SUBS2Q MULTSQ(SIMPSQRTI CADR SQ,
		       INVSQ SIMPSQRTI CADDR SQ);
   IF CAR SQ EQ 'EXPT
     THEN RETURN SIMPEXPT
		   LIST(MK!*SQ SIMPSQRTI CADR SQ,CADDR SQ);
  IF CAR SQ = '!*SQ
    THEN RETURN SIMPSQRTSQ CADR SQ;
  RETURN SIMPSQRTSQ SIMP!* SQ
  END;

ENDMODULE;


MODULE SOLVE;

EXPORTS SOLVE!-FOR!-U;

IMPORTS NTH,FINDPIVOT,GCDF,GENSYM1,MKVECT,INTERR,MULTDFCONST,
   !*MULTF!*,NEGDF,ORDDF,PLUSDF,PRINTDF,PRINTSF,PRINTSPREADC,PRINTSQ,
   QUOTF,PUTV,SPREADC,SUBST4ELIMINATEDCS,MKNILL,PNTH,DOMAINP,ADDF,
   INVSQ,MULTSQ;


%***********************************************************************
% ROUTINES FOR SOLVING THE FINAL REDUCTION EQUATION:
%**********************************************************************;


SYMBOLIC PROCEDURE UTERM(POWU,RHS);
% Finds the contribution from RHS of reduction equation, of the;
% U-coefficient given by POWU. Result is in D.F.;
   IF NULL RHS THEN NIL
   ELSE BEGIN    SCALAR COEF,POWER;
      POWER:=ADDINDS(POWU,LPOW RHS);
      COEF:=EVALUATECOEFFTS(NUMR LC RHS,POWU);
      IF NULL COEF THEN RETURN UTERM(POWU,RED RHS);
      COEF:=COEF ./ DENR LC RHS;
      RETURN PLUSDF((POWER .* COEF) .+ NIL,UTERM(POWU,RED RHS))
   END;

SYMBOLIC PROCEDURE SOLVE!-FOR!-U(RHS,LHS,ULIST);
% Solves the reduction eqn LHS = RHS. Returns list of U-coefficients;
% and their values (ULIST are those we have so far), and a list of;
% C-equations to be solved (CLIST are the eqns we have so far);
   IF NULL LHS THEN ULIST
   ELSE BEGIN    SCALAR U,LPOWLHS;
      LPOWLHS:=LPOW LHS;
      BEGIN SCALAR LL,MM,CHGE; LL:=MAXORDER(RHS,ZLIST,0);
	MM:=LORDER;
	WHILE MM DO << IF CAR LL < CAR MM THEN
		<< CHGE:=T; RPLACA(MM,CAR LL) >>;
	    LL:=CDR LL; MM:=CDR MM >>;
	IF !*TRINT AND CHGE THEN << PRINT ("Maxorder now ".LORDER) >>
      END;
      U:=PICKUPU(RHS,LPOW LHS,T);
      IF NULL U THEN
      << IF !*TRINT THEN << PRINTC "****** C-EQUATION TO SOLVE:";
	     PRINTSF NUMR LC LHS;
	     PRINTC "	 = 0";
	     PRINTC " ">>;
          % Remove a zero constant from the lhs, rather than use
	  % Gauss Elim;
	IF GAUSSELIMN(NUMR LC LHS,LT LHS) THEN
		 LHS:=SQUASHCONSTANTS(RED LHS)
        ELSE LHS:=RED LHS >>
      ELSE
      << ULIST:=(CAR U .
	   SUBS2Q MULTSQ(COEFDF(LHS,LPOWLHS),INVSQ CDR U)).ULIST;
	IF !*STATISTICS THEN !*NUMBER!*:=!*NUMBER!*+1;
	 IF !*TRINT THEN << PRINTC ("**** U(".CAR U);
	     PRINTC "	 =";
	     PRINTSQ MULTSQ(COEFDF(LHS,LPOWLHS),INVSQ CDR U);
	     PRINTC " ">>;
	 LHS:=PLUSDF(LHS,
		NEGDF MULTDFCONST(CDAR ULIST,UTERM(CAR U,RHS))) >>;
      IF !*TRINT THEN << PRINTC ".... LHS is now:";
	  PRINTDF LHS;
	  PRINTC " ">>;
      RETURN SOLVE!-FOR!-U(RHS,LHS,ULIST)
   END;

SYMBOLIC PROCEDURE SQUASHCONSTANTS(EXPRESS);
BEGIN SCALAR CONSTLST,II,XP,CL,SUBBY,CMT,XX;
	CONSTLST:=REVERSE CMAP;
	CMT:=CMATRIX;
XXX:	XX:=CAR CMT;		% Look at next row of Cmatrix;
	CL:=CONSTLST;		% and list of the names;
	II:=1;		% will become index of removed constant;
	WHILE NOT GETV(XX,II) DO
		<< II:=II+1; CL:=CDR CL >>;
	SUBBY:=CAAR CL;		%II is now index, and SUBBY the name;
	IF MEMBER(SUBBY,SILLIESLIST) THEN
		<<CMT:=CDR CMT; GO TO XXX>>; %This loop must terminate;
			% This is because at least one constant remains;
	XP:=PREPSQ !*F2Q GETV(XX,0);	% start to build up the answer;
	CL:=CDR CL;
	IF NOT (CCOUNT=II) THEN FOR JJ=II+1:CCOUNT DO <<
		IF GETV(XX,JJ) THEN
			XP:=LIST('PLUS,XP,
				LIST('TIMES,CAAR CL,
					PREPSQ !*F2Q GETV(XX,JJ)));
		CL:=CDR CL >>;
	XP:=LIST('QUOTIENT,LIST('MINUS,XP),
			PREPSQ !*F2Q GETV(XX,II));
	IF !*TRINT THEN << PRIN2 "Replace "; PRIN2 SUBBY;
		PRIN2 " by "; PRINTSQ SIMP XP >>;
	SILLIESLIST:=SUBBY . SILLIESLIST;
	RETURN SUBDF(EXPRESS,XP,SUBBY)
END;

SYMBOLIC PROCEDURE CHECKU(ULIST,U);
% Checks that U is not already in ULIST - ie. that this u-coefficient;
% has not already been given a value;
   IF NULL ULIST THEN NIL
   ELSE IF (CAR U) = CAAR ULIST THEN T
   ELSE CHECKU(CDR ULIST,U);

SYMBOLIC PROCEDURE CHECKU1(POWU,RHS);
%Checks that use of a particular U-term will not cause trouble;
%by introducing negative exponents into lhs when it is used;
    BEGIN
    TOP:
	IF NULL RHS THEN RETURN NIL;
	IF NEGIND(POWU,LPOW RHS) THEN
	  IF NOT NULL EVALUATECOEFFTS(NUMR LC RHS,POWU) THEN RETURN T;
	RHS:=RED RHS;
	GO TO TOP
    END;

SYMBOLIC PROCEDURE NEGIND(PU,PR);
%check if substituting index values in power gives rise to -ve
% exponents;
    IF NULL PU THEN NIL
    ELSE IF (CAR PU+CAAR PR)<0 THEN T
    ELSE NEGIND(CDR PU,CDR PR);


SYMBOLIC PROCEDURE EVALUATECOEFFTS(COEFFT,INDLIST);
% Substitutes the values of the i,j,k,...'s that appear in the S.F. ;
% COEFFT (=coefficient of r.h.s. of reduction equation). Result is S.F.;
   IF NULL COEFFT OR DOMAINP COEFFT THEN
      IF ZEROP COEFFT THEN NIL ELSE COEFFT
   ELSE BEGIN    SCALAR TEMP;
      IF MVAR COEFFT MEMBER INDEXLIST THEN
	 TEMP:=VALUECOEFFT(MVAR COEFFT,INDLIST,INDEXLIST)
      ELSE TEMP:=!*P2F LPOW COEFFT;
      TEMP:=!*MULTF!*(TEMP,EVALUATECOEFFTS(LC COEFFT,INDLIST));
      RETURN ADDF(!*F2POL TEMP,EVALUATECOEFFTS(RED COEFFT,INDLIST))
   END;

SYMBOLIC PROCEDURE VALUECOEFFT(VAR,INDVALUES,INDLIST);
% Finds the value of VAR, which should be in INDLIST, given INDVALUES;
% - the corresponding values of INDLIST variables;
   IF NULL INDLIST THEN INTERR "VALUECOEFFT - NO VALUE"
   ELSE IF VAR EQ CAR INDLIST THEN
      IF ZEROP CAR INDVALUES THEN NIL
      ELSE CAR INDVALUES
   ELSE VALUECOEFFT(VAR,CDR INDVALUES,CDR INDLIST);

SYMBOLIC PROCEDURE ADDINDS(POWU,POWRHS);
% Adds indices in POWU to those in POWRHS. Result is LPOW of D.F.;
   IF NULL POWU THEN IF NULL POWRHS THEN NIL
      ELSE INTERR "POWRHS TOO LONG"
   ELSE IF NULL POWRHS THEN INTERR "POWU TOO LONG"
   ELSE (CAR POWU + CAAR POWRHS).ADDINDS(CDR POWU,CDR POWRHS);


SYMBOLIC PROCEDURE PICKUPU(RHS,POWLHS,FLG);
% Picks up the 'lowest' U coefficient from RHS if it exists and returns;
% it in the form of LT of D.F.;
% returns NIL if no legal term in RHS can be found;
% POWLHS is the power we want to match (LPOW of D.F);
% and COEFFU is the list of previous coefficients that must be zero;
 BEGIN SCALAR COEFFU,U;
    PT:=RHS;
TOP:
    IF NULL PT THEN RETURN NIL; %no term found - failed;
    U:=NEXTU(LT PT,POWLHS); %check this term...;
    IF NULL U THEN GO TO NOTTHISONE;
    IF NOT TESTORD(CAR U,LORDER) THEN GO TO NEVERTHISONE;
    IF NOT CHECKCOEFFTS(COEFFU,CAR U) THEN GO TO NOTTHISONE;
    %that inhibited clobbering things already passed over;
    IF CHECKU(ULIST,U) THEN GO TO NOTTHISONE;
    %that avoided redefining a u value;
    IF CHECKU1(CAR U,RHS) THEN GO TO NEVERTHISONE;
    %avoid introduction of negative exponents;
    IF FLG THEN
	U:=PATCHUPTAN(LIST U,POWLHS,RED PT,RHS);
    RETURN U;
NEVERTHISONE:
    COEFFU:=(LC PT) . COEFFU;
NOTTHISONE:
    PT:=RED PT;
    GO TO TOP
 END;

SYMBOLIC PROCEDURE PATCHUPTAN(U,POWLHS,RPT,RHS);
	BEGIN
	    SCALAR UU,CC,DD,TANLIST,REDU,REDU1;
	    PT:=RPT;
	    WHILE PT DO <<
		IF (UU:=PICKUPU(PT,POWLHS,NIL)) 
			AND TESTORD(CAR UU,LORDER) THEN <<
				% Nasty found, patch it up;
		    CC:=(GENSYM1('!C).CAAR U).CC;
				% CC is an alist of constants;
		    IF !*TRINT THEN << PRINTC ("****** U(".CAAR U);
			PRINTC "     =";
			PRINT CAAR CC >>;
		    REDU:=PLUSDF(REDU,
			MULTDFCONST(!*K2Q CAAR CC,UTERM(CAAR U,RHS)));
		    U:=UU.U
		>>;
		IF PT THEN PT:=RED PT >>;
	    REDU1:=REDU;
	    WHILE REDU1 DO BEGIN SCALAR XX; XX:=CAR REDU1;
IF !*TRINT THEN << PRIN2 "Introduced RESIDUE "; PRINT XX >>;
		IF (NOT TESTORD(CAR XX,LORDER)) THEN <<
		    IF !*TRINT THEN <<
			PRINTSQ CDR XX; PRINTC "  =  0" >>;
		    IF DD:=KILLSINGLES(CADR XX,CC) THEN <<
			REDU:=SUBDF(REDU,0,CAR DD);
			REDU1:=SUBDF(REDU1,0,CAR DD);
			ULIST:=((CDR DD).(NIL ./ 1)).ULIST;
			U:=RMVE(U,CDR DD);
			CC:=PURGECONST(CC,DD) >>
		    ELSE REDU1:=CDR REDU1  >>
		ELSE REDU1:=CDR REDU1  END;
	    FOREACH XX IN REDU DO <<
		IF (NOT TESTORD(CAR XX,LORDER)) THEN <<
		    WHILE CC DO << 
				ADDCTOMAP(CAAR CC);
				ULIST:=((CDAR CC).(!*K2Q CAAR CC))
					  . ULIST;
				IF !*STATISTICS
				  THEN !*NUMBER!*:=!*NUMBER!*+1;
				CC:=CDR CC >>;
			GAUSSELIMN(NUMR LC REDU,LT REDU)>> >>;
	    IF REDU THEN << WHILE CC DO << ADDCTOMAP(CAAR CC);
			ULIST:=((CDAR CC).(!*K2Q CAAR CC)).ULIST;
			IF !*STATISTICS THEN !*NUMBER!*:=!*NUMBER!*+1;
			CC:=CDR CC >>;
		LHS:=PLUSDF(LHS,NEGDF REDU) >>;
    RETURN CAR U
END;

SYMBOLIC PROCEDURE KILLSINGLES(XX,CC);
  IF ATOM XX THEN NIL
  ELSE IF NOT (CDR XX EQ NIL) THEN NIL
  ELSE BEGIN SCALAR DD;
    DD:=ASSOC(CAAAR XX,CC);
    IF DD THEN RETURN DD;
    RETURN KILLSINGLES(CDAR XX,CC)
END;

SYMBOLIC PROCEDURE RMVE(L,X);
   IF CAAR L=X THEN CDR L ELSE CONS(CAR L,RMVE(CDR L,X));

SYMBOLIC PROCEDURE SUBDF(A,B,C);
% SUBSTITUTE B FOR C INTO THE DF A;
% Used to get rid of silly constants introduced;
IF A=NIL THEN NIL ELSE
  BEGIN SCALAR X;
    X:=SUBF(NUMR LC A,LIST (C . B)) ;
    IF X=(NIL . 1) THEN RETURN SUBDF(RED A,B,C)
	ELSE RETURN PLUSDF(
		LIST ((LPOW A).((CAR X).MULTF(CDR X,DENR LC A))),
		SUBDF(RED A,B,C))
END;

SYMBOLIC PROCEDURE TESTORD(A,B);
% Test order of two DF's in recursive fashion;
  IF NULL A THEN T
    ELSE IF CAR A LEQ CAR B THEN TESTORD(CDR A,CDR B)
    ELSE NIL;

SYMBOLIC PROCEDURE TANFROM(RHS,Z,NN);
% We notice that in all bad cases we have (j-num)tan**j...;
% Extract the num;
BEGIN SCALAR N,ZZ,R,RR;
    R:=RHS;
    N:=0; ZZ:=ZLIST;
    WHILE CAR ZZ NEQ Z DO << N:=N+1; ZZ:=CDR ZZ >>;
    WHILE R DO <<
	RR:=CAAR R;  % The list of powers;
	FOR I=1:N DO RR:=CDR RR;
	IF FIXP CAAR RR THEN IF CAAR RR>0 THEN <<
		RR:=NUMR CDAR R;
		IF NULL RED RR THEN RR:=NIL ./ 1
         ELSE IF FIXP (RR:=QUOTF(RED RR,LC RR))
		THEN RR:=-RR ELSE RR:=0>>;
	IF ATOM RR THEN RETURN RR;
	R:=CDR R >>;
    IF NULL R THEN RETURN MAXFROM(LHS,NN)+1;
   RETURN MAX(RR,MAXFROM(LHS,NN)+1)
END;


SYMBOLIC PROCEDURE COEFDF(Y,U);
  IF Y=NIL THEN NIL
  ELSE IF LPOW Y=U THEN LC Y
  ELSE COEFDF(RED Y,U);


SYMBOLIC PROCEDURE PURGECONST(A,B);
% Remove a const from and expression. May be the same as DELETE?;
  IF NULL A THEN NIL
  ELSE IF CAR A=B THEN PURGECONST(CDR A,B)
  ELSE CONS(CAR A,PURGECONST(CDR A,B));

SYMBOLIC PROCEDURE MAXORDER(RHS,Z,N);
% Find a limit on the order of terms, theis is ad hoc;
  IF NULL Z THEN NIL
    ELSE IF EQCAR(CAR Z,'SQRT) THEN
	CONS(1,MAXORDER(RHS,CDR Z,N+1))
    ELSE IF (ATOM CAR Z) OR (CAAR Z NEQ 'TAN) THEN
	CONS(MAXFROM(LHS,N)+1,MAXORDER(RHS,CDR Z,N+1))
    ELSE CONS(TANFROM(RHS,CAR Z,N),MAXORDER(RHS,CDR Z,N+1));

SYMBOLIC PROCEDURE MAXFROM(L,N);
% Largest order in the nth varable;
  IF NULL L THEN 0
  ELSE MAX(NTH(CAAR L,N+1),MAXFROM(CDR L,N));


SYMBOLIC PROCEDURE COPY U;
  IF ATOM U THEN U
    ELSE CONS(COPY CAR U,COPY CDR U);


SYMBOLIC PROCEDURE ADDCTOMAP CC;
BEGIN
    SCALAR NCVAL;
    CCOUNT:=CCOUNT+1;
    NCVAL:=MKVECT(CCOUNT);
    FOR I=0:(CCOUNT-1) DO PUTV(NCVAL,I,GETV(CVAL,I));
    PUTV(NCVAL,CCOUNT,NIL ./ 1);
    CVAL:=NCVAL;
    CMAP:=(CC . CCOUNT).CMAP;
    IF !*TRINT THEN << PRIN2 "Constant Map CHANGED TO "; PRINT CMAP >>;
    CMATRIX:=MAPCAR(CMATRIX,FUNCTION ADDTOVECTOR);
END;

SYMBOLIC PROCEDURE ADDTOVECTOR V;
    BEGIN SCALAR VV;
	VV:=MKVECT(CCOUNT);
	FOR I=0:(CCOUNT-1) DO PUTV(VV,I,GETV(V,I));
	PUTV(VV,CCOUNT,NIL);
	RETURN VV
    END;

SYMBOLIC PROCEDURE CHECKCOEFFTS(CL,INDV);
% checks to see that the coefficients in CL (coefficient list - S.Q.s);
% are zero when the i,j,k,... are given values in INDV (LPOW of;
% D.F.). if so the result is true else NIL=false;
    IF NULL CL THEN T
    ELSE BEGIN    SCALAR RES;
	RES:=EVALUATECOEFFTS(NUMR CAR CL,INDV);
	IF NOT(NULL RES OR RES=0) THEN RETURN NIL
	ELSE RETURN CHECKCOEFFTS(CDR CL,INDV)
    END;

SYMBOLIC PROCEDURE NEXTU(LTRHS,POWLHS);
% picks out the appropriate U coefficients for term: LTRHS to match the;
% powers of the z-variables given in POWLHS (= exponent list of D.F.). ;
% return this coefficient in form LT of D.F. If U coefficient does;
% not exist then result is NIL. If it is multiplied by a zero then;
% result is NIL;
   IF NULL LTRHS THEN NIL
   ELSE BEGIN    SCALAR INDLIST,UCOEFFT;
      INDLIST:=SUBTRACTINDS(POWLHS,CAR LTRHS,NIL);
      IF NULL INDLIST THEN RETURN NIL;
      UCOEFFT:=EVALUATECOEFFTS(NUMR CDR LTRHS,INDLIST);
      IF NULL UCOEFFT OR UCOEFFT=0 THEN RETURN NIL;
      RETURN INDLIST .* (UCOEFFT ./ DENR CDR LTRHS)
   END;

SYMBOLIC PROCEDURE SUBTRACTINDS(POWLHS,L,SOFAR);
% subtract the indices in list L from those in POWLHS to find;
% appropriate values for i,j,k,... when equating coefficients of terms;
% on lhs of reduction eqn. SOFAR is the resulting value list we;
% have constructed so far. if any i,j,k,... value is -ve then result;
% is NIL;
    IF NULL L THEN REVERSEWOC SOFAR
    ELSE IF ((CAR POWLHS)-(CAAR L))<0 THEN NIL
    ELSE SUBTRACTINDS(CDR POWLHS,CDR L,
	((CAR POWLHS)-(CAAR L)) . SOFAR);

SYMBOLIC PROCEDURE GAUSSELIMN(EQUATION,TOKILL);
% Performs Gaussian elimination on the matrix for the c-equations;
% as each c-equation is found. EQUATION is the next one to deal with;
   BEGIN	 SCALAR NEWROW,PIVOT;
      IF ZEROP CCOUNT THEN GO TO NOWAY; %FAILURE;
      NEWROW:=MKVECT(CCOUNT);
      SPREADC(EQUATION,NEWROW,1);
      SUBST4ELIMINATEDCS(NEWROW,REVERSE ORDEROFELIM,REVERSE CMATRIX);
      PIVOT:=FINDPIVOT NEWROW;
      IF NULL PIVOT THEN GO TO NOPIVOTFOUND;
      ORDEROFELIM:=PIVOT . ORDEROFELIM;
      NEWROW:=MAKEPRIM NEWROW; %REMOVE HCF FROM NEW EQUATION;
      CMATRIX:=NEWROW . CMATRIX;
%      IF !*TRINT THEN PRINTSPREADC NEWROW;
      RETURN T;
 NOPIVOTFOUND:
      IF NULL GETV(NEWROW,0) THEN <<
	IF !*TRINT THEN PRINTC "Already included";
	RETURN NIL>>; %EQUATION WAS 0=0;
 NOWAY:
      BADPART:=TOKILL . BADPART; %NON-INTEGRABLE TERM;
      IF !*TRINT THEN PRINTC "Inconsistent";
      RETURN NIL
   END;

SYMBOLIC PROCEDURE MAKEPRIM ROW;
    BEGIN	  SCALAR I,G;
	G:=GETV(ROW,0);
	FOR I:=1:CCOUNT DO G:=GCDF(G,GETV(ROW,I));
	IF G NEQ 1 THEN 
	   FOR I:=0:CCOUNT DO PUTV(ROW,I,QUOTF(GETV(ROW,I),G));
	FOR I := 0:CCOUNT DO
	  <<G := GETV(ROW,I);
	    IF G AND NOT DOMAINP G
	      THEN PUTV(ROW,I,NUMR RESIMP((ROOTEXTRACTF G) ./ 1))>>;
	RETURN ROW
    END;




ENDMODULE;


MODULE SQRTF;

EXPORTS MINUSDFP,SQRTDF,NROOTN,DOMAINP,MINUSF;

IMPORTS CONTENTSMV,GCDF,INTERR,!*MULTF!*,PARTIALDIFF,PRINTDF,QUOTF,
   SIMPSQRT2,VP2;

%SQUARE-ROOT OF STANDARD FORMS;





SYMBOLIC PROCEDURE MINUSDFP A;
%TEST SIGN OF LEADING COEDD OF D.F;
    IF NULL A THEN INTERR "MINUSDFP 0 ILLEGAL"
    ELSE MINUSF NUMR LC A;

SYMBOLIC PROCEDURE SQRTDF L;
%TAKES SQUARE ROOT OF D.F.;
    IF NULL L THEN NIL
    ELSE IF NOT NULL RED L THEN 'FAILED
    ELSE BEGIN SCALAR C;
	IF LPOW L=VP2 ZLIST THEN GO TO OK;
	PRINTC "SQRTDF NOT COMPLETE";
	PRINTDF L;
	RETURN 'FAILED;
    OK: RETURN (LPOW L .* SQRTSQ LC L) .+ NIL
    END;

SYMBOLIC PROCEDURE SQRTSQ A;
    SQRTF NUMR A ./ SQRTF DENR A;

SYMBOLIC PROCEDURE SQRTF P;
    BEGIN	SCALAR IP,QP;
	IF NULL P THEN RETURN NIL;
	IP:=SQRTF1 P;
	QP:=CDR IP;
	IP:=CAR IP; %RESPECTABLE AND NASTY PARTS OF THE SQRT;
	IF ONEP QP THEN RETURN IP; %EXACT ROOT FOUND;
        QP:=SIMPSQRT2 QP;
	RETURN !*F2POL !*MULTF!*(IP,QP)
    END;

SYMBOLIC PROCEDURE SQRTF1 P;
%RETURNS A . B WITH P=A**2*B;
    IF DOMAINP P THEN NROOTN(P,2)
    ELSE BEGIN SCALAR CO,PP,G,PG;
	CO:=CONTENTSMV(P,MVAR P,NIL); %CONTENTS OF P;
	PP:=QUOTF(P,CO); %PRIMITIVE PART;
	CO:=SQRTF1(CO); %PROCESS CONTENTS VIA RECURSION;
	G:=GCDF(PP,PARTIALDIFF(PP,MVAR PP));
	PG:=QUOTF(PP,G);
	G:=GCDF(G,PG); %A REPEATED FACTOR OF PP;
	IF G=1 THEN PG:=1 . PP
	ELSE <<
	    PG:= !*F2POL QUOTF(PP,!*MULTF!*(G,G)); %WHAT IS STILL LEFT;
	    PG:=SQRTF1(PG); %SPLIT THAT UP;
	    RPLACA(PG,!*MULTF!*(CAR PG,G))>>;
		 %PUT IN THE THING FOUND HERE;
	RPLACA(PG,!*F2POL !*MULTF!*(CAR PG,CAR CO));
	RPLACD(PG,!*F2POL !*MULTF!*(CDR PG,CDR CO));
	RETURN PG
    END;

% NROOTN removed as in REDUCE base;

ENDMODULE;


MODULE TDIFF;

EXPORTS !-!-SIMPDF;

IMPORTS SIMPCAR,KERNP,DIFFSQ,PREPSQ,MSGPRI;

FLAG('(!-!-SIMPDF),'LOSE);

%TDF(EXPR,VAR) DIFFERENTIATES BUT WITH TIMING SERVICE;

SYMBOLIC PROCEDURE !-!-SIMPDF U;
   %U IS A LIST OF FORMS, THE FIRST AN EXPRESSION AND THE REMAINDER
   %KERNELS AND NUMBERS.
   %VALUE IS DERIVATIVE OF FIRST FORM WRT REST OF LIST;
   BEGIN    SCALAR V,X,Y,TT;
	TT := TIME(); %start the clock;
	V := CDR U;
	U := SIMPCAR U;
    A:	IF NULL V OR NULL NUMR U THEN GO TO EXIT;
	X := IF NULL Y OR Y=0 THEN SIMPCAR V ELSE Y;
	IF NULL KERNP X THEN GO TO E;
	X := CAAAAR X;
	V := CDR V;
	IF NULL V THEN GO TO C;
	Y := SIMPCAR V;
	IF NULL NUMR Y THEN GO TO D
	 ELSE IF NOT DENR Y=1 OR NOT NUMBERP NUMR Y THEN GO TO C;
	Y := CAR Y;
	V := CDR V;
    B:	IF Y=0 THEN GO TO A;
	U := DIFFSQ(U,X);
	Y := Y-1;
	GO TO B;
    C:	U := DIFFSQ(U,X);
	GO TO A;
    D:	Y := NIL;
	V := CDR V;
	GO TO A;
    EXIT:
       PRINT LIST('TIME,TIME()-TT);
       RETURN U;
    E:  MSGPRI("DIFFERENTIATION WRT",PREPSQ X,"NOT ALLOWED",NIL,T)
   END;

PUT('TDF,'SIMPFN,'!-!-SIMPDF);


ENDMODULE;


MODULE TIDYSQRT; 
 
EXPORTS SQRT2TOP;

%GENERAL TIDYING UP ABOUT SQUARE ROOTS; 
 
%SYMBOLIC PROCEDURE TIDYSQRTDF A; 
%    IF NULL A THEN NIL 
%    ELSE BEGIN    SCALAR TT,R; 
%        TT:=TIDYSQRT LC A; 
%        R:=TIDYSQRTDF RED A; 
%        IF NULL NUMR TT THEN RETURN R; 
%        RETURN ((LPOW A) .* TT) .+ R 
%    END; 
% 
%SYMBOLIC PROCEDURE TIDYSQRT Q; 
%    BEGIN    SCALAR NN,DD; 
%        NN:=TIDYSQRTF NUMR Q; 
%        IF NULL NN THEN NIL ./ 1; %ANSWER IS ZERO; 
%        DD:=TIDYSQRTF DENR Q; 
%        RETURN MULTSQ(NN,INVSQ DD) 
%    END; 
% 
% 
%SYMBOLIC PROCEDURE TIDYSQRTF P; 
%%INPUT - STANDARD FORM; 
%%OUTPUT - STANDARD QUOTIENT; 
%% SIMPLIFIES SQRT(A)**N WITH N>1; 
%    IF DOMAINP P THEN P ./ 1 
%    ELSE BEGIN    SCALAR V,W; 
%        V:=LPOW P; 
%        IF CAR V='I THEN V:=MKSP('(SQRT -1),CDR V); %I->SQRT(-1); 
%        IF EQCAR(CAR V,'SQRT) AND NOT ONEP CDR V THEN BEGIN SCALAR X; 
% %HERE WE HAVE A REDUCTION TO APPLY; 
%            X:=DIVIDE(CDR V,2); %HALVE EXPONENT; 
%            W:=EXPTSQ(SIMP CADAR V,CAR X); %RATIONAL PART OF ANSWER; 
%            IF NOT ZEROP CDR X THEN W:=MULTSQ(W, 
%                ((MKSP(CAR V,1) .* 1) .+ NIL) ./ 1); 
%            %THE NEXT LINE ALLOWS FOR THE HORRORS OF NESTED SQRTS; 
%            W:=TIDYSQRT W 
%            END 
%        ELSE W:=((V .* 1) .+ NIL) ./ 1; 
%        V:=MULTSQ(W,TIDYSQRTF LC P); 
%        RETURN ADDSQ(V,TIDYSQRTF RED P) 
%    END; 
% 
%
%MOVE SQRTS IN A SQ TO THE NUMERATOR; 
 
SYMBOLIC PROCEDURE MULTOUTDENR Q; 
    BEGIN  SCALAR N,D,ROOT,CONJ; 
        N:=NUMR Q; 
        D:=DENR Q; 
   LOOP:ROOT:=FINDSQUAREROOT D; %SEARCH DENOM; 
        IF NULL ROOT THEN RETURN (N . D);
	%NOTHING TO BE DONE; 
        CONJ:=CONJUGATEWRT(D,ROOT); 
        N:=!*F2POL !*MULTF!*(N,CONJ); 
        D:=!*F2POL !*MULTF!*(D,CONJ); 
        GO TO LOOP 
        END; 
 
 
SYMBOLIC PROCEDURE SQRT2TOP Q; 
BEGIN 
  SCALAR N,D; 
  N:=MULTOUTDENR Q; 
  D:=DENR N; 
  N:=NUMR N; 
  IF D EQ DENR Q 
    THEN RETURN Q;%NO CHANGE; 
  IF D IEQUAL 1 
    THEN RETURN (N ./ 1); 
  Q:=GCDCOEFFSOFSQRTS N; 
  IF Q IEQUAL 1 
    THEN IF MINUSF D 
      THEN RETURN (NEGF N ./ NEGF D) 
      ELSE RETURN (N ./ D); 
  Q:=GCDF(Q,D); 
  N:=QUOTF(N,Q); 
  D:=QUOTF(D,Q); 
  IF MINUSF D 
    THEN RETURN (NEGF N ./ NEGF D) 
    ELSE RETURN (N ./ D) 
    END; 
 
 
%SYMBOLIC PROCEDURE DENRSQRT2TOP Q; 
%BEGIN 
%  SCALAR N,D; 
%  N:=MULTOUTDENR Q; 
%  D:=DENR N; 
%  N:=NUMR N; 
%  IF D EQ DENR Q 
%    THEN RETURN D;  %NO CHANGES; 
%  IF D IEQUAL 1 
%    THEN RETURN 1; 
%  Q:=GCDCOEFFSOFSQRTS N; 
%  IF Q IEQUAL 1 
%    THEN RETURN D; 
%  Q:=GCDF(Q,D); 
%  IF Q IEQUAL 1 
%    THEN RETURN D 
%    ELSE RETURN QUOTF(D,Q) 
%  END; 
 
SYMBOLIC PROCEDURE FINDSQUAREROOT P; 
%LOCATE A SQRT SYMBOL IN POLY P; 
    IF DOMAINP P THEN NIL 
    ELSE BEGIN SCALAR W; 
        W:=MVAR P; %CHECK MAIN VAR FIRST; 
        IF ATOM W 
          THEN RETURN NIL; %WE HAVE PASSED ALL SQRTS; 
        IF EQCAR(W,'SQRT) THEN RETURN W; 
        W:=FINDSQUAREROOT LC P; 
        IF NULL W THEN W:=FINDSQUAREROOT RED P; 
        RETURN W 
    END; 
 
SYMBOLIC PROCEDURE CONJUGATEWRT(P,VAR); 
% VAR -> -VAR IN FORM P; 
    IF DOMAINP P THEN P 
    ELSE IF MVAR P=VAR THEN BEGIN 
        SCALAR X,C,R; 
        X:=TDEG LT P; %DEGREE; 
        C:=LC P; %COEFFICIENT; 
        R:=RED P; %REDUCTUM; 
        X:=REMAINDER(X,2); %NOW JUST 0 OR 1; 
        IF X=1 THEN C:=NEGF C; %-COEFFICIENT; 
        RETURN (LPOW P .* C) .+ CONJUGATEWRT(R,VAR) END 
    ELSE IF ORDOP(VAR,MVAR P) THEN P 
    ELSE (LPOW P .* CONJUGATEWRT(LC P,VAR)) .+ 
        CONJUGATEWRT(RED P,VAR); 
 
SYMBOLIC PROCEDURE GCDCOEFFSOFSQRTS U; 
IF ATOM U 
  THEN IF NUMBERP U AND MINUSP U 
    THEN -U 
    ELSE U 
  ELSE IF EQCAR(MVAR U,'SQRT) 
    THEN BEGIN 
      SCALAR V; 
      V:=GCDCOEFFSOFSQRTS LC U; 
      IF V IEQUAL 1 
        THEN RETURN V 
        ELSE RETURN GCDF(V,GCDCOEFFSOFSQRTS RED U) 
      END 
    ELSE BEGIN 
      SCALAR ROOT; 
      ROOT:=FINDSQUAREROOT U; 
      IF NULL ROOT 
        THEN RETURN U; 
      U:=MAKEMAINVAR(U,ROOT); 
      ROOT:=GCDCOEFFSOFSQRTS LC U; 
      IF ROOT IEQUAL 1 
        THEN RETURN 1 
        ELSE RETURN GCDF(ROOT,GCDCOEFFSOFSQRTS RED U) 
      END; 

ENDMODULE;


MODULE TRCASE;

EXPORTS TRANSCENDENTALCASE;

IMPORTS BACKSUBST4CS,COUNTZ,CREATECMAP,CREATEINDICES,DF2Q,DFNUMR,
  DIFFLOGS,FSDF,FACTORLISTLIST,FINDSQRTS,FINDTRIALDIVS,GCDF,MKVECT,
  INTERR,LOGSTOSQ,MERGIN,MULTBYARBPOWERS,!*MULTF!*,MULTSQFREE,
  PRINTDF,PRINTFACTORS,PRINTSQ,QUOTF,RATIONALINTEGRATE,PUTV,
  SIMPINT1,SOLVE!-FOR!-U,SQFREE,SQMERGE,SQRT2TOP,SUBSTINULIST,TRIALDIV,
  MERGEIN,NEGSQ,ADDSQ,F2DF,MKNILL,PNTH,INVSQ,MULTSQ,DOMAINP,MK!*SQ,
  MKSP,PRETTYPRINT,PREPSQ;

FLUID '(DENBAD VAR XLOGS);      % For the ERRORSET below;

SYMBOLIC 
   PROCEDURE TRANSCENDENTALCASE(INTEGRAND,VAR,XLOGS,ZLIST,VARLIST);
   BEGIN SCALAR DIVLIST,W,JHD!-CONTENT,CONTENT,PRIM,SQFR,DFU,INDEXLIST,
%      JHD!-CONTENT is local, while CONTENT is free (set in SQFREE);
	SILLIESLIST,ORIGINALORDER,ORIGINALLHS,WRONGWAY,
      SQRTLIST,TANLIST,LOGLIST,DFLOGS,EPRIM,DFUN,UNINTEGRAND,
      SQRTFLAG,BADPART,RHS,LHS,GCDQ,CMAP,CVAL,ORDEROFELIM,CMATRIX;
      SCALAR CUBEROOTFLAG,CCOUNT,DENOMINATOR,RESULT,DENBAD;
	GENSYMCOUNT:=0;
      INTEGRAND:=SQRT2TOP INTEGRAND; % Move the sqrts to the numerator;
      IF !*TRINT THEN << PRINTC "EXTENSION VARIABLES Z<I> ARE";
	  PRINT ZLIST>>;
      IF !*RATINTSPECIAL AND NULL CDR ZLIST THEN
	    RETURN RATIONALINTEGRATE(INTEGRAND,VAR);
% *** NOW UNNORMALIZE INTEGRAND, MAYBE *** ; 
     BEGIN SCALAR W,Z,GG; 
	GG:=1; 
	FOREACH Z IN ZLIST DO <<
	    W:=DIFFSQ(SIMP Z,VAR); 
	    GG:=MULTF(GG,QUOTF(DENR W,GCDF(DENR W,GG))) >>; 
	GG:=QUOTF(GG,GCDF(GG,DENR INTEGRAND)); 
	UNINTEGRAND:=(MULTF(GG,NUMR INTEGRAND) 
			./ MULTF(GG,DENR INTEGRAND)); 
	IF !*TRINT THEN <<
		PRINTC "UNNORMALIZED INTEGRAND ="; 
		PRINTSQ UNINTEGRAND >> END; 
      DIVLIST:=FINDTRIALDIVS ZLIST;
		 %ALSO PUTS SOME THINGS ON LOGLIST SOMETIMES;
%     IF !*TRINT THEN << PRINTC "EXPONENTIALS AND TANS TO TRY DIVIDING:";
%	  PRINT DIVLIST>>;
	SQRTLIST:=FINDSQRTS ZLIST;
%     IF !*TRINT THEN << PRINTC "SQUARE-ROOT Z-VARIABLES";
%	  PRINT SQRTLIST >>;
      DIVLIST:=TRIALDIV(DENR UNINTEGRAND,DIVLIST);
%     IF !*TRINT THEN << PRINTC "DIVISORS:";
%	  PRINT CAR DIVLIST;
%	  PRINT CDR DIVLIST>>;
%N.B. THE NEXT LINE ALSO SETS 'CONTENT' AS A FREE VARIABLE;
% Since SQFREE may be used later, we copy it into JHD!-CONTENT;
      PRIM:=SQFREE(CDR DIVLIST,ZLIST);
      JHD!-CONTENT:=CONTENT;
      PRINTFACTORS(PRIM,NIL);
      EPRIM:=SQMERGE(COUNTZ CAR DIVLIST,PRIM,NIL);
      PRINTFACTORS(EPRIM,T);
%     IF !*TRINT THEN << TERPRI();
%	  PRINTSF DENOMINATOR;
%	  TERPRI();
%	  PRINTC "...CONTENT IS:";
%	  PRINTSF JHD!-CONTENT>>;
      SQFR:=MULTSQFREE EPRIM;
%     IF !*TRINT THEN << PRINTC "...SQFR IS:";
%	  SUPERPRINT SQFR>>;
      INDEXLIST:=CREATEINDICES ZLIST;
%     IF !*TRINT THEN << PRINTC "...INDICES ARE:";
%	  SUPERPRINT INDEXLIST>>;
      DFU:=DFNUMR(VAR,CAR DIVLIST);
%     IF !*TRINT THEN << TERPRI();
%	  PRINTC "************ DERIVATIVE OF U IS:";
%	  PRINTSQ DFU>>;
      LOGLIST:=APPEND(LOGLIST,FACTORLISTLIST (PRIM,NIL));
      LOGLIST:=MERGEIN(XLOGS,LOGLIST);
      LOGLIST:=MERGEIN(TANLIST,LOGLIST);
      CMAP:=CREATECMAP();
      CCOUNT:=LENGTH CMAP;
      IF !*TRINT THEN << PRINTC "LOGLIST ";
	   PRINT LOGLIST >>;
      DFLOGS:=DIFFLOGS(LOGLIST,DENR UNINTEGRAND,VAR);
      IF !*TRINT THEN << PRINTC "************ 'DERIVATIVE' OF LOGS IS:";
	  PRINTSQ DFLOGS>>;
      DFLOGS:=ADDSQ((NUMR UNINTEGRAND) ./ 1,NEGSQ DFLOGS);
      % Put everything in reduction eqn over common denominator: ;
      GCDQ:=GCDF(DENR DFLOGS,DENR DFU);
      DFUN:= !*F2POL !*MULTF!*(NUMR DFU,
				DENBAD:=QUOTF(DENR DFLOGS,GCDQ));
      DENBAD:=!*MULTF!*(DENR DFU,DENBAD);
      DENBAD:= !*F2POL !*MULTF!*(DENR UNINTEGRAND,DENBAD);
      DFLOGS:= !*F2POL !*MULTF!*(NUMR DFLOGS,QUOTF(DENR DFU,GCDQ));
      DFU:=DFUN;
      % Now DFU and DFLOGS are S.F.s;
      RHS:=MULTBYARBPOWERS F2DF DFU;
      IF !*TRINT THEN << PRINTC "Distributed Form of U is:";
	  PRINTDF RHS>>;
      LHS:=F2DF DFLOGS;
      IF !*TRINT THEN << PRINTC "Distributed Form of l.h.s. is:";
	  PRINTDF LHS;
	  TERPRI()>>;
      CVAL:=MKVECT(CCOUNT);
      FOR I:=0 : CCOUNT DO PUTV(CVAL,I,NIL ./ 1);
      LORDER:=MAXORDER(RHS,ZLIST,0);
	ORIGINALORDER:=LORDER;
	ORIGINALLHS:=LHS;
	IF !*TRINT THEN << PRINTC "Maximum order determined as ";
		PRINT LORDER >>;
	IF !*STATISTICS THEN << !*NUMBER!*:=0;
		!*SPSIZE!*:=1;
		FOREACH XX IN LORDER DO
		   !*SPSIZE!*:=!*SPSIZE!* * (XX+1) >>;
		% That calculates the largest U that can appear;
      DFUN:=SOLVE!-FOR!-U(RHS,LHS,NIL);
      BACKSUBST4CS(NIL,ORDEROFELIM,CMATRIX);
%      IF !*TRINT THEN IF NOT (CCOUNT=0) THEN PRINTVECSQ CVAL;
	IF !*STATISTICS THEN << PRIN2 !*NUMBER!*; PRIN2 " used out of ";
		PRINTC !*SPSIZE!* >>;
      BADPART:=SUBSTINULIST BADPART;
		 %SUBSTITUTE FOR C<I> STILL IN BADPART;
      DFUN:=DF2Q SUBSTINULIST DFUN;
%     IF !*TRINT THEN SUPERPRINT DFUN;
      RESULT:= SUBS2Q MULTSQ(DFUN,INVSQ(DENOMINATOR ./ 1));
      RESULT:= SUBS2Q MULTSQ(RESULT,INVSQ(JHD!-CONTENT ./ 1));
%     IF !*TRINT THEN SUPERPRINT RESULT;
      DFLOGS:=LOGSTOSQ();
      IF NOT NULL NUMR DFLOGS
	THEN RESULT:=ADDSQ(RESULT,DFLOGS);
      IF !*TRINT THEN << SUPERPRINT RESULT;
	  TERPRI();
	  PRINTC
	  "*****************************************************";
	  PRINTC
	   "************ THE INTEGRAL IS : **********************";
	  PRINTC
	   "*****************************************************";
	  TERPRI();
	  PRINTSQ RESULT;
	  TERPRI()>>;
      IF NOT NULL BADPART THEN <<
	  IF !*TRINT THEN PRINTC "PLUS A BAD PART";
	  LHS:=BADPART;
	  LORDER:=MAXORDER(RHS,ZLIST,0);
	  WHILE LORDER DO <<
		IF CAR LORDER > CAR ORIGINALORDER THEN
			WRONGWAY:=T;
		LORDER:=CDR LORDER;
		ORIGINALORDER:=CDR ORIGINALORDER >>;
	  DFUN:=DF2Q BADPART;
	  IF !*TRINT
	    THEN <<PRINTSQ DFUN; PRINTC "DENBAD = "; PRINTSF DENBAD>>;
	  DFUN:= SUBS2Q MULTSQ(DFUN,INVSQ(DENBAD ./ 1));
	  IF WRONGWAY THEN << RESULT:= NIL ./ 1; DFUN:=INTEGRAND >>;
	  IF ROOTCHECKP(UNINTEGRAND,VAR) THEN
		RETURN SIMPINT1(INTEGRAND . VAR.NIL)
          ELSE IF !*PURERISCH OR ALLOWEDFNS ZLIST THEN 
    	      DFUN:=SIMPINT1 (DFUN . VAR.NIL)
           ELSE << !*PURERISCH:=T;
		IF !*TRINT
		  THEN <<PRINTC "   [Transforming ..."; PRINTSQ DFUN>>;
              DENBAD:=TRANSFORM(DFUN,VAR);
	      IF DENBAD=DFUN
		THEN DFUN:=SIMPINT1(DFUN . VAR.NIL)
              ELSE <<DENBAD:=ERRORSET('(INTEGRATESQ DENBAD VAR XLOGS),
				      NIL,!*BACKTRACE);
		IF NOT ATOM DENBAD THEN DFUN:=UNTAN CAR DENBAD
                ELSE DFUN:=SIMPINT1(DFUN . VAR.NIL) >> >>;
	      IF !*TRINT THEN PRINTSQ DFUN;
	      IF !*FAILHARD THEN INTERR "FAILHARD SWITCH SET";
	  RESULT:=ADDSQ(RESULT,DFUN) >>;
%      IF !*OVERLAYMODE
%	THEN EXCISE TRANSCODE;
      RETURN SQRT2TOP RESULT
   END;

%UNFLUID '(DFUN VAR XLOGS);

ENDMODULE;


MODULE HALFANGLE;

EXPORTS HALFANGLE,UNTAN;

SYMBOLIC PROCEDURE TRANSFORM(U,X);
% Transform the SQ U to remove the 'bad' functions sin, cos, cot etc
% in favor of half angles;
    HALFANGLE(U,X);


% Rest of this page is due to Harrington;

%PROCEDURES FOR CONVERSION TO HALF ANGLE TANGENTS;


% SOME NEWRED PROCEDURES THAT IM USED TO;

SYMBOLIC PROCEDURE QUOTQQ(U1,V1);
MULTSQ(U1, INVSQ(V1));

SYMBOLIC PROCEDURE !*SUBTRQ(U1,V1);
ADDSQ(U1, NEGSQ(V1));


SYMBOLIC PROCEDURE !*INT2QM(U1);
IF U1=0 THEN NIL . 1 ELSE U1 . 1;

SYMBOLIC PROCEDURE HALFANGLE(R,X);
% TOP LEVEL PROCEDURE FOR CONVERTING;
% R IS A RATIONAL EXPRESSION TO BE CONVERTED,
% X THE INTEGRATION VARIABLE;
% A RATIONAL EXPRESSION IS RETURNED;
QUOTQQ(HFAGLF(NUMR(R),X), HFAGLF(DENR(R),X));

SYMBOLIC PROCEDURE HFAGLF(P,X);
% CONVERTING POLYNOMIALS,  A RATIONAL EXPRESSION IS RETURNED;
IF DOMAINP(P) THEN !*F2Q(P)
ELSE SUBS2Q ADDSQ(MULTSQ(EXPTSQ(HFAGLK(MVAR(P),X), LDEG(P)),
		  HFAGLF(LC(P),X)),
  HFAGLF(RED(P),X));

SYMBOLIC PROCEDURE HFAGLK(K,X);
% CONVERTING KERNELS,  A RATIONAL EXPRESSION IS RETURNED;
BEGIN
   SCALAR KT;
   IF ATOM K OR NOT MEMBER(X,FLATTEN(CDR(K))) THEN RETURN !*K2Q K;
   K := CAR(K) . HFAGLARGS(CDR(K), X);
   KT := SIMP LIST('TAN, LIST('QUOTIENT, CADR(K), 2));
   RETURN IF CAR(K) = 'SIN
    THEN QUOTQQ(MULTSQ(!*INT2QM(2),KT), ADDSQ(!*INT2QM(1),
			 EXPTSQ(KT,2)))
   ELSE IF CAR(K) = 'COS
    THEN QUOTQQ(!*SUBTRQ(!*INT2QM(1), EXPTSQ(KT,2)), ADDSQ(!*INT2QM(1),
      EXPTSQ(KT,2)))
   ELSE IF CAR(K) = 'TAN
    THEN QUOTQQ(MULTSQ(!*INT2QM(2),KT), !*SUBTRQ(!*INT2QM(1),
			 EXPTSQ(KT,2)))
   ELSE IF CAR(K) = 'SINH THEN
     QUOTQQ(!*SUBTRQ(EXPTSQ(!*K2Q('EXPT.('E. CDR K)),2),
     !*INT2QM(1)), MULTSQ(!*INT2QM(2), !*K2Q('EXPT . ('E . CDR(K)))))
   ELSE IF CAR(K) = 'COSH THEN
     QUOTQQ(ADDSQ(EXPTSQ(!*K2Q('EXPT.('E. CDR K)),2),
     !*INT2QM(1)), MULTSQ(!*INT2QM(2), !*K2Q('EXPT . ('E . CDR(K)))))
   ELSE IF CAR(K) = 'TANH THEN
     QUOTQQ(!*SUBTRQ(EXPTSQ(!*K2Q('EXPT.('E. CDR K)),2),
     !*INT2QM(1)), ADDSQ(EXPTSQ(!*K2Q ('EXPT.('E.CDR(K))),2),
     !*INT2QM(1)))
   ELSE !*K2Q(K);  % ADDITIONAL TRANSFORMATION MIGHT BE ADDED HERE;
END;


SYMBOLIC PROCEDURE HFAGLARGS(L,X);
%CONVERSION OF ARGUMENT LIST;
IF NULL L THEN NIL
ELSE PREPSQ(HFAGLK(CAR(L),X)) . HFAGLARGS(CDR(L), X);

SYMBOLIC PROCEDURE UNTANF X; 
   BEGIN SCALAR Y,Z,W; 
      IF DOMAINP X THEN RETURN X . 1; 
      Y := MVAR X; 
      IF EQCAR(Y,'INT) THEN ERROR(99,NIL);  %assume all is hopeless;
      Z := LDEG X; 
      W := 1 . 1; 
      Y := 
       IF ATOM Y THEN !*K2Q Y
	ELSE IF CAR Y EQ 'TAN
         THEN IF REMAINDER(Z,2)=0
                THEN <<Z := Z/2; 
                       SIMP LIST('QUOTIENT,
                                 LIST('PLUS,
                                      LIST('MINUS,
                                           LIST('COS,
                                                'TIMES
                                                  . (2 . CDR Y))),
                                      1),LIST('PLUS,
                                              LIST('COS,
                                                   'TIMES
                                                     . (2 . CDR Y)),
                                              1))>>
               ELSE IF Z=1
                THEN SIMP LIST('QUOTIENT,
                               LIST('PLUS,
                                    LIST('MINUS,
                                         LIST('COS,
                                              'TIMES . (2 . CDR Y))),
                                    1),LIST('SIN,
                                            'TIMES . (2 . CDR Y)))
               ELSE <<Z := (Z - 1)/2; 
                      W := 
                       SIMP LIST('QUOTIENT,
                                 LIST('PLUS,
                                      LIST('MINUS,
                                           LIST('COS,
                                                'TIMES
                                                  . (2 . CDR Y))),
                                      1),LIST('SIN,
                                              'TIMES
                                                . (2 . CDR Y))); 
                      SIMP LIST('QUOTIENT,
                                LIST('PLUS,
                                     LIST('MINUS,
                                          LIST('COS,
                                               'TIMES
                                                 . (2 . CDR Y))),
                                     1),LIST('PLUS,
                                             LIST('COS,
                                                  'TIMES
                                                    . (2 . CDR Y)),
                                             1))>>
	ELSE SIMP Y;
      RETURN ADDSQ(MULTSQ(MULTSQ(EXPTSQ(Y,Z),UNTANF LC X),W),
                   UNTANF RED X)
   END;

SYMBOLIC PROCEDURE UNTANLIST(Y);
IF NULL Y THEN NIL ELSE (PREPSQ (UNTAN(SIMP CAR Y)) . UNTANLIST(CDR Y));

SYMBOLIC PROCEDURE UNTAN(X);
COMMENT EXPECTS X TO BE CANONICAL QUOTIENT;
BEGIN SCALAR Y;
Y:=COSSQCHK SINSQRDCHK MULTSQ(UNTANF(NUMR X), INVSQ  UNTANF(DENR X));
RETURN IF LENGTH FLATTEN Y>LENGTH FLATTEN X THEN X ELSE Y
END;

SYMBOLIC PROCEDURE SINSQRDCHK(X);
MULTSQ(SINSQCHKF(NUMR X), INVSQ SINSQCHKF(DENR X));

SYMBOLIC PROCEDURE SINSQCHKF(X);
BEGIN
   SCALAR Y,Z,W;
   IF DOMAINP X THEN RETURN X . 1;
   Y := MVAR X;
   Z := LDEG X;
   W := 1 . 1;
   Y := IF EQCAR(Y,'SIN) THEN IF REMAINDER(Z,2) = 0
    THEN <<Z := QUOTIENT(Z,2);
	   SIMP LIST('PLUS,1,LIST('MINUS,
				  LIST('EXPT,('COS . CDR(Y)),2)))>>
   ELSE IF Z = 1 THEN !*K2Q Y
   ELSE  << Z := QUOTIENT(DIFFERENCE(Z,1),2); W := !*K2Q Y;
          SIMP LIST('PLUS,1,LIST('MINUS,
				 LIST('EXPT,('COS . CDR(Y)),2)))>>
    ELSE !*K2Q Y;
   RETURN ADDSQ(MULTSQ(MULTSQ(EXPTSQ(Y,Z),SINSQCHKF(LC X)),W),
		SINSQCHKF(RED X));
END;


SYMBOLIC PROCEDURE COSSQCHKF(X);
BEGIN
   SCALAR Y,Z,W,X1,X2;
   IF DOMAINP X THEN RETURN X . 1;
   Y := MVAR X;
   Z := LDEG X;
   W := 1 . 1;
   X1 := COSSQCHKF(LC X);
   X2 := COSSQCHKF(RED X);
   X := ADDSQ(MULTSQ(!*P2Q LPOW X,X1),X2);
   Y := IF EQCAR(Y,'COS) THEN IF REMAINDER(Z,2) = 0
    THEN <<Z := QUOTIENT(Z,2);
	   SIMP LIST('PLUS,1,LIST('MINUS,
				  LIST('EXPT,('SIN . CDR(Y)),2)))>>
   ELSE IF Z = 1 THEN !*K2Q Y
   ELSE  << Z := QUOTIENT(DIFFERENCE(Z,1),2); W := !*K2Q Y;
          SIMP LIST('PLUS,1,LIST('MINUS,
				 LIST('EXPT,('SIN . CDR(Y)),2)))>>
    ELSE !*K2Q Y;
   Y := ADDSQ(MULTSQ(MULTSQ(EXPTSQ(Y,Z),W),X1),X2);
   RETURN IF LENGTH(Y) > LENGTH(X) THEN X ELSE Y;
END;

SYMBOLIC PROCEDURE COSSQCHK(X);
BEGIN
   SCALAR GCD1;
   GCD1 := !*GCD;
   !*GCD := T;
   X := MULTSQ(COSSQCHKF(NUMR X), INVSQ COSSQCHKF(DENR X));
   !*GCD := GCD1;
   RETURN X;
END;


SYMBOLIC PROCEDURE LROOTCHK(L,X);
% CHECKS EACH MEMBER OF LIST L FOR A ROOT;
IF NULL L THEN NIL ELSE KROOTCHK(CAR L, X) OR LROOTCHK(CDR L, X);

SYMBOLIC PROCEDURE KROOTCHK(F,X);
% CHECKS A KERNEL TO SEE IF IT IS A ROOT;
IF ATOM F THEN NIL
ELSE IF CAR(F) = 'SQRT
     AND MEMBER(X, FLATTEN CDR F)  THEN T
ELSE IF CAR(F) = 'EXPT
     AND NOT ATOM CADDR(F)
     AND CAADDR(F) = 'QUOTIENT
     AND MEMBER(X, FLATTEN CADR F)  THEN T
ELSE LROOTCHK(CDR F, X);

SYMBOLIC PROCEDURE ROOTCHK1P(F,X);
% CHECKS POLYNOMIAL FOR A ROOT;
IF DOMAINP F THEN NIL
ELSE KROOTCHK(MVAR F,X) OR ROOTCHK1P(LC F, X) OR ROOTCHK1P(RED F, X);

SYMBOLIC PROCEDURE ROOTCHECKP(F,X);
% CHECKS RATIONAL (STANDARD QUOTIENT) FOR A ROOT;
ROOTCHK1P(NUMR F, X) OR ROOTCHK1P(DENR F, X);

ENDMODULE;


MODULE TRIALDIV;

EXPORTS COUNTZ,FINDSQRTS,FINDTRIALDIVS,PRINTFACTORS,TRIALDIV,SIMP,MKSP;

IMPORTS !*MULTF!*,PRINTSF,QUOTF;


SYMBOLIC PROCEDURE COUNTZ DL;
% DL is a list of S.F.s;
    BEGIN	  SCALAR S,N,RL;
LOOP2:	IF NULL DL THEN RETURN ARRANGELISTZ RL;
	N:=1;
LOOP1:	N:=N+1;
	S:=CAR DL;
	DL:=CDR DL;
	IF NOT NULL DL AND (S EQ CAR DL) THEN
	    GO TO LOOP1
	ELSE RL:=(S.N).RL;
	GO TO LOOP2
    END;

SYMBOLIC PROCEDURE ARRANGELISTZ D;
    BEGIN	  SCALAR N,S,RL,R;
	N:=1;
	IF NULL D THEN RETURN RL;
LOOPD:	IF (CDAR D)=N THEN S:=(CAAR D).S
	ELSE R:=(CAR D).R;
	D:=CDR D;
	IF NOT NULL D THEN GO TO LOOPD;
	D:=R;
	RL:=S.RL;
	S:=NIL;
	R:=NIL;
	N:=N+1;
	IF NOT NULL D THEN GO TO LOOPD;
	RETURN REVERSEWOC RL
    END;

SYMBOLIC PROCEDURE PRINTFACTORS(W,PRDENOM);
    % W is a list of factors to each power. If PRDENOM is true ;
    % this prints denominator of answer, else prints square-free ;
    % decomposition. ;
    BEGIN	  SCALAR I,WX;
	I:=1;
	IF PRDENOM THEN <<
	    DENOMINATOR:=1;
	    IF !*TRINT
	      THEN PRINTC "DENOMINATOR OF 1ST PART OF ANSWER IS:";
	    IF NOT NULL W THEN W:=CDR W >>;
LOOPX:	IF W=NIL THEN RETURN;
	IF !*TRINT THEN PRINTC ("FACTORS OF MULTIPLICITY".I);
	WX:=CAR W;
	WHILE NOT NULL WX DO <<
	    IF !*TRINT THEN PRINTSF CAR WX;
	    FOR J:=1 : I DO 
		DENOMINATOR:= !*F2POL !*MULTF!*(CAR WX,DENOMINATOR);
		%this call of F2POL is probably not necessary??;
	    WX:=CDR WX >>;
	I:=I+1;
	W:=CDR W;
	GO TO LOOPX
    END;

SYMBOLIC PROCEDURE FINDTRIALDIVS ZL;
%ZL IS LIST OF KERNELS FOUND IN INTEGRAND. RESULT IS A LIST;
%GIVING THINGS TO BE TREATED SPECIALLY IN THE INTEGRATION;
%VIZ: EXPS AND TANS;
%RESULT IS LIST OF FORM ((A . B) ...);
% WITH A A KERNEL AND CAR A=EXPT OR TAN;
% AND B A STANDARD FORM FOR EITHER EXPT OR (1+TAN**2);
    BEGIN	  SCALAR DLISTS1,ARGS1;
	WHILE NOT NULL ZL DO <<
	    IF EXPORTAN CAR ZL THEN <<
		IF CAAR ZL='TAN
		  THEN << ARGS1:=(MKSP(CAR ZL,2) .* 1) .+ 1;
		    TANLIST:=(ARGS1 ./ 1) . TANLIST>>
		ELSE ARGS1:=!*K2F CAR ZL;
		DLISTS1:=(CAR ZL . ARGS1) . DLISTS1>>;
	    ZL:=CDR ZL >>;
	RETURN DLISTS1
    END;

SYMBOLIC PROCEDURE EXPORTAN DL;
    IF ATOM DL THEN NIL
    ELSE BEGIN
    % EXTRACT EXP OR TAN FNS FROM THE Z-LIST;
    IF EQ(CAR DL,'TAN) THEN RETURN T;
NXT:	IF NOT EQ(CAR DL,'EXPT) THEN RETURN NIL;
	DL:=CADR DL;
        IF ATOM DL THEN RETURN T;
	GO TO NXT
    END;


SYMBOLIC PROCEDURE FINDSQRTS Z; 
    BEGIN  SCALAR R; 
        WHILE NOT NULL Z DO << 
            IF EQCAR(CAR Z,'SQRT) THEN R:=(CAR Z) . R; 
            Z:=CDR Z >>; 
        RETURN R 
    END; 
SYMBOLIC PROCEDURE TRIALDIV(X,DL);
    BEGIN	  SCALAR QLIST,Q;
    WHILE NOT NULL DL DO
	IF NOT NULL(Q:=QUOTF(X,CDAR DL)) THEN <<
	    IF (CAAAR DL='TAN) AND NOT EQCAR(QLIST,CDAR DL) THEN
		LOGLIST:=('IDEN . SIMP CADR CAAR DL) . LOGLIST;
			 %TAN FIDDLE!;
	    QLIST:=(CDAR DL).QLIST;
	    X:=Q >>
	ELSE DL:=CDR DL;
    RETURN QLIST.X
    END;


ENDMODULE;


MODULE UNIFAC;

EXPORTS EVALAT,LINETHROUGH,QUADTHROUGH,TESTDIV,UNIFAC,ZFACTORS;

IMPORTS CUBIC,LINFAC,PRINTDF,QUADFAC,QUADRATIC,QUARTIC,VP1,ZFACTOR,
   GCD,MINUSP,PRETTYPRINT;

%UNIVARIATE FACTORIZATION FOR INTEGRATION;

SYMBOLIC PROCEDURE ZFACTORS N;
%PRODUCES A LIST OF ALL (POSITIVE) INTEGER FACTORS OF THE ;
%INTEGER N;
    IF N=0 THEN LIST 0
    ELSE IF (N:=ABS N)=1 THEN LIST 1
    ELSE COMBINATIONTIMES ZFACTOR N;

SYMBOLIC PROCEDURE ZFACTOR N;
% INPUT N A POSITIVE INTEGER;
% OUTPUT A LIST ((PRIME . EXPONENT) ...) GIVING FACTORS OF N;
    BEGIN	  SCALAR FL,Q,W,C;
	C:=0; %MULTIPLICITY;
 TRY2:	Q:=DIVIDE(N,2); %PULL OUT FACTORS OF 2;
	IF ZEROP CDR Q THEN <<
	    C:=C+1;
	    N:=CAR Q;
	    GO TO TRY2 >>;
	IF NOT ZEROP C THEN FL:=(2 . C) . FL;
	W:=3; C:=0;
 TRYW:	Q:=DIVIDE(N,W);
	IF ZEROP CDR Q THEN <<
	    C:=C+1;
	    N:=CAR Q;
	    GO TO TRYW >>;
	IF NOT ZEROP C THEN FL:=(W . C) . FL;
	IF REMAINDER(W,3)=1 THEN W:=W+4
	    ELSE W:=W+2;
	C:=0;
	IF NOT ((W*W)>N) THEN GO TO TRYW;
	IF NOT ONEP N THEN FL:=(N . 1) . FL;
	RETURN FL
    END;

SYMBOLIC PROCEDURE COMBINATIONTIMES FL;
    IF NULL FL THEN LIST 1
    ELSE BEGIN    SCALAR N,C,RES,PR;
	N:=CAAR FL; C:=CDAR FL;
	PR:=COMBINATIONTIMES CDR FL;
	WHILE NOT MINUSP C DO <<
	    RES:=PUTIN(EXPT(N,C),PR,RES);
	    C:=C-1 >>;
	RETURN RES
    END;

SYMBOLIC PROCEDURE PUTIN(N,L,W);
    IF NULL L THEN W
    ELSE PUTIN(N,CDR L,(N*CAR L) . W);


SYMBOLIC PROCEDURE UNIFAC(POL,VAR,DEGREE,RES);
    BEGIN	  SCALAR W,Q,C;
	W:=POL;
	IF !*TRINT THEN SUPERPRINT W;
%NOW TRY LOOKING FOR LINEAR FACTORS;
TRYLIN: Q:=LINFAC(W);
	IF NULL CAR Q THEN GO TO NOMORELIN;
	RES := ('LOG . BACK2DF(CAR Q,VAR)) . RES;
	W:=CDR Q;
	GO TO TRYLIN;
NOMORELIN:
	Q:=QUADFAC(W);
	IF NULL CAR Q THEN GO TO NOMOREQUAD;
	RES := QUADRATIC(BACK2DF(CAR Q,VAR),VAR,RES);
	W:=CDR Q;
	GO TO NOMORELIN;
NOMOREQUAD:
	IF NULL W THEN RETURN RES; %ALL DONE;
	DEGREE:=CAR W; %DEGREE OF WHAT IS LEFT;
	C:=BACK2DF(W,VAR);
	IF DEGREE=3 THEN RES:=CUBIC(C,VAR,RES)
	ELSE IF DEGREE=4 THEN RES:=QUARTIC(C,VAR,RES)
	ELSE IF ZEROP REMAINDER(DEGREE,2) AND
		PAIRP (Q := HALFPOWER CDDR W)
	 THEN <<W := (DEGREE/2) . (CADR W . Q);
	        W := UNIFAC(W,VAR,CAR W,NIL);
		RES := PLUCKFACTORS(W,VAR,RES)>>
	ELSE <<
	    PRINTC "THE FOLLOWING HAS NOT BEEN SPLIT";
	    PRINTDF C;
	    RES:=('LOG . C) . RES>>;
	RETURN RES
    END;

SYMBOLIC PROCEDURE HALFPOWER W;
   IF NULL W THEN NIL
    ELSE IF CAR W=0 
     THEN (LAMBDA R;
	   IF R EQ 'FAILED THEN R ELSE CADR W . R) HALFPOWER CDDR W
    ELSE 'FAILED;

SYMBOLIC PROCEDURE PLUCKFACTORS(W,VAR,RES);
   BEGIN SCALAR S,P,Q,R,KNOWNDISCRIMSIGN;
      WHILE W DO
	<<P := CAR W;
	  IF CAR P EQ 'ATAN THEN NIL
	   ELSE IF CAR P EQ 'LOG
	    THEN <<Q := DOUBLEPOWER CDR P . Q;
		   %PRIN2 "Q="; %PRINTDF CAR Q;
		  >>
	   ELSE INTERR "BAD FORM";
	  W := CDR W>>;
      WHILE Q DO
       <<P := CAR Q;
	 IF CAAAR P=4 
	   THEN <<KNOWNDISCRIMSIGN := 'NEGATIVE;
		  RES := QUARTIC(P,VAR,RES);
	          KNOWNDISCRIMSIGN := NIL>>
	   ELSE IF CAAAR P=2 
	    THEN RES := QUADRATIC(P,VAR,RES)
	   ELSE RES := ('LOG . P) . RES;
	  Q := CDR Q>>;
      RETURN RES
   END;

SYMBOLIC PROCEDURE DOUBLEPOWER R;
   IF NULL R THEN NIL
    ELSE (LIST(2*CAAAR R) . CDAR R) . DOUBLEPOWER CDR R;

SYMBOLIC PROCEDURE BACK2DF(P,V);
%UNDO THE EFFECT OF UNIFORM;
    BEGIN	  SCALAR R,N;
	N:=CAR P;
	P:=CDR P;
	WHILE NOT MINUSP N DO <<
	    IF NOT ZEROP CAR P THEN R:=
		(VP1(V,N,ZLIST) .* (CAR P ./ 1)) .+ R;
	    P:=CDR P;
	    N:=N-1 >>;
	RETURN REVERSEWOC R
    END;

SYMBOLIC PROCEDURE EVALAT(P,N);
%EVALUATE POLYNOMIAL AT INTEGER POINT N;
    BEGIN	  SCALAR R;
	R:=0;
	P:=CDR P;
	WHILE NOT NULL P DO <<
	    R:=N*R+CAR P;
	    P:=CDR P >>;
	RETURN R
    END;

SYMBOLIC PROCEDURE TESTDIV(A,B);
% QUOTIENT A/B OR FAILED;
    BEGIN	  SCALAR Q;
	Q:=TESTDIV1(CDR A,CAR A,CDR B,CAR B);
	IF Q='FAILED THEN RETURN Q;
	RETURN (CAR A-CAR B) . Q
    END;

SYMBOLIC PROCEDURE TESTDIV1(A,DA,B,DB);
    IF DA<DB THEN BEGIN
    CHECK0: IF NULL A THEN RETURN NIL
	    ELSE IF NOT ZEROP CAR A THEN RETURN 'FAILED;
	    A:=CDR A;
	    GO TO CHECK0
	END
    ELSE BEGIN    SCALAR Q;
	Q:=DIVIDE(CAR A,CAR B);
	IF ZEROP CDR Q THEN Q:=CAR Q
	ELSE RETURN 'FAILED;
	A:=TESTDIV1(AMBQ(CDR A,CDR B,Q),DA-1,B,DB);
	IF A='FAILED THEN RETURN A;
	RETURN Q . A
    END;

SYMBOLIC PROCEDURE AMBQ(A,B,Q);
% A-B*Q WITH Q AN INTEGER;
    IF NULL B THEN A
    ELSE ((CAR A)-(CAR B)*Q) . AMBQ(CDR A,CDR B,Q);


SYMBOLIC PROCEDURE LINETHROUGH(Y0,Y1);
    BEGIN	  SCALAR A;
	A:=Y1-Y0;
	IF ZEROP A THEN RETURN 'FAILED;
	IF A<0 THEN <<A:=-A; Y0:=-Y0 >>;
	IF ONEP GCDN(A,Y0) THEN RETURN LIST(1,A,Y0);
	RETURN 'FAILED
    END;


SYMBOLIC PROCEDURE QUADTHROUGH(YM1,Y0,Y1);
    BEGIN	  SCALAR A,B,C;
	A:=DIVIDE(YM1+Y1,2);
	IF ZEROP CDR A THEN A:=(CAR A)-Y0
	ELSE RETURN 'FAILED;
	IF ZEROP A THEN RETURN 'FAILED; %LINEAR THINGS ALREADY DONE;
	C:=Y0;
	B:=DIVIDE(Y1-YM1,2);
	IF ZEROP CDR B THEN B:=CAR B
	ELSE RETURN 'FAILED;
	IF NOT ONEP GCDN(A,GCD(B,C)) THEN RETURN 'FAILED;
	IF A<0 THEN <<A:=-A; B:=-B; C:=-C>>;
	RETURN LIST(2,A,B,C)
    END;


ENDMODULE;


MODULE UNIFORM;

EXPORTS UNIFORM;

IMPORTS EXPONENTOF;


SYMBOLIC PROCEDURE UNIFORM(P,V);
%CONVERT FROM D.F. IN ONE VARIABLE (V) TO A SIMPLE LIST OF;
%COEFFS (WITH DEGREE CONSED ONTO FRONT);
%FAILS IF COEFFICIENTS ARE NOT ALL SIMPLE INTEGERS;
    IF NULL P THEN 0 . (0 . NIL)
    ELSE BEGIN    SCALAR A,B,C,D;
	A:=EXPONENTOF(V,LPOW P,ZLIST);
	B:=LC P;
	IF NOT ONEP DENR B THEN RETURN 'FAILED;
	B:=NUMR B;
	IF NULL B THEN B:=0
	ELSE IF NOT NUMBERP B THEN RETURN 'FAILED;
	IF A=0 THEN RETURN A . (B . NIL); %CONSTANT TERM;
	C:=UNIFORM(RED P,V);
	IF C='FAILED THEN RETURN 'FAILED;
	D:=CAR C;
	C:=CDR C;
	D:=D+1;
	WHILE NOT (A=D) DO <<
	    C:=0 . C;
	    D:=D+1>>;
	RETURN A . (B . C)
    END;


ENDMODULE;


MODULE MAKEVARS;

EXPORTS GETVARIABLES,VARSINLIST,VARSINSQ,VARSINSF,FINDZVARS,
	CREATEINDICES,MERGEIN;

IMPORTS DEPENDSP,UNION;


% Note that 'i' is already maybe committed for sqrt(-1);
%also 'l' and 'o' are not used as the print badly on certain;
%terminals etc and may lead to confusion;

!*GENSYMLIST!* := '(! j ! k ! l ! m ! n ! o ! p ! q ! r ! s
		    ! t ! u ! v ! w ! x ! y ! z);

%MAPC(!*GENSYMLIST!*,FUNCTION REMOB); %REMOB protection;


SYMBOLIC PROCEDURE VARSINLIST(L,VL);
%L IS A LIST OF S.Q. - FIND ALL VARIABLES MENTIONED;
%GIVEN THAL VL IS A LIST ALREADY KNOWN ABOUT;
    BEGIN	WHILE NOT NULL L DO <<
	    VL:=VARSINSF(NUMR CAR L,VARSINSF(DENR CAR L,VL));
	    L:=CDR L >>;
	RETURN VL
    END;

SYMBOLIC PROCEDURE GETVARIABLES SQ;
    VARSINSF(NUMR SQ,VARSINSF(DENR SQ,NIL));

SYMBOLIC PROCEDURE VARSINSF(FORM,L);
   IF ATOM FORM THEN L
   ELSE BEGIN
     WHILE NOT ATOM FORM DO <<
	L:=VARSINSF(LC FORM,UNION(L,LIST MVAR FORM));
	FORM:=RED FORM >>;
     RETURN L
   END;

SYMBOLIC PROCEDURE FINDZVARS(VL,ZL,VAR,FLG);
    BEGIN	  SCALAR V;
% VL is the crude list of variables found in the original integrand;
% ZL must have merged into it all EXP, LOG etc terms from this;
% If FLG is true then ignore DF as a function;
SCAN: IF NULL VL THEN RETURN ZL;
	 V:=CAR VL; % NEXT VARIABLE;
	 VL:=CDR VL;
% at present items get put onto ZL if they are non-atomic;
% and they depend on the main variable. The arguments of;
% functions are decomposed by recursive calls to findzvar;
	%give up if V has been declared dependent on other things;
	IF ASSOC(V,DEPL!*) THEN ERROR1()
	 ELSE IF NOT ATOM V AND (NOT V MEMBER ZL) AND DEPENDSP(V,VAR)
	 THEN IF CAR V MEMQ '(TIMES QUOTIENT PLUS MINUS DIFFERENCE INT)
		 OR (((CAR V) EQ 'EXPT) AND FIXP CADDR V)
	     THEN
		 ZL:=FINDZVARS(CDR V,ZL,VAR,FLG)
	    ELSE IF FLG AND CAR V='DF THEN
		<< !*PURERISCH:=T; RETURN ZL >>   % TRY AND STOP IT;
	     ELSE ZL:=V.FINDZVARS(CDR V,ZL,VAR,FLG);
		 % SCAN ARGUMENTS OF FN;
	GO TO SCAN
   END;

SYMBOLIC PROCEDURE CREATEINDICES ZL; 
% Produces a list of unique indices, each associated with a ; 
% different Z-variable; 
     REVERSEWOC CRINDEX1(ZL,!*GENSYMLIST!*); 
 
SYMBOLIC PROCEDURE CRINDEX1(ZL,GL); 
 BEGIN IF NULL ZL THEN RETURN NIL; 
    IF NULL GL THEN << GL:=LIST GENSYM1 'i; %new symbol needed; 
        NCONC(!*GENSYMLIST!*,GL) >>; 
    RETURN (CAR GL) . CRINDEX1(CDR ZL,CDR GL) END; 

SYMBOLIC PROCEDURE RMEMBER(A,B);
    IF NULL B THEN NIL
    ELSE IF A=CDAR B THEN CAR B
    ELSE RMEMBER(A,CDR B);

SYMBOLIC PROCEDURE MERGEIN(DL,LL);
%ADJOIN LOGS OF THINGS IN DL TO EXISTING LIST LL;
    IF NULL DL THEN LL
    ELSE IF RMEMBER(CAR DL,LL) THEN MERGEIN(CDR DL,LL)
    ELSE MERGEIN(CDR DL,('LOG . CAR DL) . LL);


ENDMODULE;


MODULE VECTOR;

EXPORTS MKIDENM,MKVEC2,MKVEC;

IMPORTS MKNILL,PNTH;


SYMBOLIC PROCEDURE MKVEC(L);
BEGIN
  SCALAR V,I;
  V:=MKVECT(-1+LENGTH L);
  I:=0;
  WHILE L DO <<
    PUTV(V,I,(CAR L) ./ 1);
    I:=I+1;
    L:=CDR L >>;
  RETURN V
  END;

ENDMODULE;


END;