reduce-historical/r33/int.red

module int!-intro; % General support for REDUCE integrator.

% Authors: A. C. Norman and P. M. A. Moore.
% Modified by: J. Davenport, J. P. Fitch, A. C. Hearn.

% Note that at one point, INT had been flagged SIMP0FN.  However, that
% lead to problems when the arguments of INT contained pattern
% variables.

fluid '(!*conscount !*noextend !*pvar gaussiani);

global '(btrlevel frlis!* gensymcount initl!*);

!*conscount:=10000; % default maximum number of conses in certain
                    % operations.

!*pvar:='!_a;

btrlevel := 5; %default to a reasonably full backtrace.

% The next smacro is needed at this point to define gaussiani.

symbolic smacro procedure !*kk2f u; !*p2f mksp(u,1);

gaussiani := !*kk2f '(sqrt -1);

gensymcount := 0;

initl!* := append('(!*noextend), initl!*);

flag('(interr),'transfer);   %For the compiler;

flag ('(atan dilog erf expint expt log tan),'transcendental);

comment Kludge to define derivative of an integral and integral of
        a derivative;

frlis!* := union('(!=x !=y),frlis!*);

put('df,'opmtch,'(((int !=y !=x) !=x) (nil . t)
                  (evl!* !=y) nil) . get('df,'opmtch));

put('int,'opmtch,'(((df !=y !=x) !=x) (nil . t)
                  (evl!* !=y) nil) . get('int,'opmtch));

put('evl!*,'opmtch,'(((!=x) (nil . t) !=x nil)));

put('evl!*,'simpfn,'simpiden);


%Various functions used throughout the integrator.

smacro procedure !*kk2q a; ((mksp(a,1) .* 1) .+ nil) ./ 1;

symbolic smacro procedure divsf(u,v); sqrt2top(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 int!-gensym1 u;
    << gensymcount:=gensymcount+1;
       compress append(explode u,explode gensymcount) >>;

symbolic smacro procedure maninp(u,v,w);
   interr "MANINP called -- not implemented";

symbolic procedure mknill n; if n=0 then nil else nil . mknill(n-1);


% Various selectors written as macros.

smacro procedure argof u;
   % Argument of a unary function.
   cadr u;

smacro procedure firstsubs u;
  % The first substitution in a substitution list.
  car u;

smacro procedure lsubs u; car u;

smacro procedure rsubs u; cdr u;

smacro procedure lfirstsubs u; caar u;

smacro procedure rfirstsubs u; cdar u;

put('nthroot,'simpfn,'simpiden);
% The binary n-th root operator nthroot(x,2)=sqrt(x)
% no simplification is used here.
% Hope is that pbuild introduces it, and simplog removes it.


% Selectors for the taylor series structure.

% Format is:
%function.((first.last computed so far) . assoc list of computed terms).

% ***store-hack-1***:
% remove this macro if more store is available.

smacro procedure tayshorten u;nil;

smacro procedure taylordefn u; car u;

symbolic smacro procedure taylorfunction u; caar u;

smacro procedure taylornumbers u; cadr u;

smacro procedure taylorfirst u; caadr u;

smacro procedure taylorlast u; cdadr u;

smacro procedure taylorlist u; cddr u;

smacro procedure taylormake(fn,nums,alist); fn.(nums.alist);

endmodule;


module contents;

% Authors: Mary Ann Moore and Arthur C. Norman

fluid '(clogflag content indexlist sqfr varlist zlist);

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 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.f. for product of elements of l;
   begin         scalar res;
      res:=1;
      while not null l do <<
         res:=multf(res,car l);
         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:=!*multf(multup rl,multup cdr l);
      temp:=!*multsq(difff(car l,x),!*f2q 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 (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),
                      !*f2q multup sqfr);
      if not null sqfr and not null cdr sqfr then <<
          temp2:=!*multsq(diflist(cdr sqfr,coeflist,x,nil),
                          !*f2q 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(!*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 !*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;   % routines to do with the C constants.

% Author: John P. Fitch.

fluid '(ccount cmap cmatrix cval loglist neweqn);

global '(!*trint);

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;

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,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:=(int!-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),
                                !*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,t2) >>;
        t1:=((lpow eqnn) .* 1) .+ nil; %this main var as sf.
        spreadc(lc eqnn,cvec1,!*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;  % Cube roots of standard forms.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(cuberootflag);

exports cuberootdf;

imports contentsmv,gcdf,!*multf,nrootn,partialdiff,printdf,quotf,vp2,
   mksp,mk!*sq,domainp;

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 numberp qp and 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 !*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 pg; %split that up.
            rplaca(pg,!*multf(car pg,g))>>;
                 %put in the thing found here.
        rplaca(pg,!*multf(car pg,car co));
        rplacd(pg,!*multf(cdr pg,cdr co));
        return pg
    end;

endmodule;


module idepend;  % Routines for considering dependency among variables.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(taylorvariable);

exports dependspl,dependsp,involvesq,involvsf;

imports taylorp,domainp;

symbolic procedure dependsp(x,v);
   if null v then t
    else if depends(x,v) then x
    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 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 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 dependsp(mvar sf,term)
       or involvesf(lc sf,term)
       or involvesf(red sf,term);

symbolic procedure dependspl(dep!-list,var);
   % True if any member of deplist (a list of prefix forms) depends on
   % var.
   dep!-list
      and (dependsp(car dep!-list,var) or dependspl(cdr dep!-list,var));

symbolic procedure taylorp exxpr;
   % Sees if a random entity is a taylor expression.
   not atom exxpr
       and not atom car exxpr
       and flagp(taylorfunction exxpr,'taylor);

endmodule;


module df2q;   % Conversion from distributive to standard forms.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(indexlist zlist);

exports df2q;

imports addf,gcdf,mksp,!*multf,quotf;

comment We assume that results already have reduced powers, so
        that no power substitution is necessary;

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;  % Routines for manipulating distributed forms.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(indexlist sqrtlist zlist);

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;

%***********************************************************************
% 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:=!*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 .* !*multsq(x,lc u) .+ multdfconst(x,red u);

%symbolic procedure quotdfconst(x,u);
%    multdfconst(!*invsq x,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 sub1ind power;
%     if atom power then power-1
%     else list sub1 car power;

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;  % Exact division of standard forms to give a S Q.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(residue sqrtlist zlist);

global '(!*trdiv !*trint);

exports fquotf,testdivdf,dfquotdf;

imports df2q,f2df,gcdf,interr,multdf,negdf,plusdf,printdf,printsf,
   quotf,multsq,invsq,negsq;

% 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:= !*multsq(lc a,!*invsq lc b);
        if null numr lcd then return dfquotdf2(red a,b);
           % Should not be necessary;
        lcd := xd .* lcd;
        xd:=plusdf(a,multdf(negdf (lcd .+ nil),b));
        if xd and (lpow xd = lpow a % Again, should not be necessary;
                   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;

% The following hack makes sure that the results of differentiation
% gets passed through ROOTEXTRACT
% a) This should not be done this way, since the effect is global
% b) Should this be done via TIDYSQRT?

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;  % Driving routines for integration program.

% Author: Mary Ann Moore and Arthur C. Norman.
% Modifications by: John P. Fitch.

fluid '(!*backtrace
        !*exp
        !*gcd
        !*keepsqrts
        !*mcd
        !*nolnr
        !*purerisch
        !*rationalize
        !*sqrt
        !*structure
        !*uncached
        basic!-listofnewsqrts
        basic!-listofallsqrts
        expression
        gaussiani
        intvar
        listofnewsqrts
        listofallsqrts
        loglist
        sqrt!-intvar
        sqrt!-places!-alist
        variable
        varlist
        xlogs
        zlist);

global '(!*algint !*failhard !*trint);

exports integratesq,simpint,purge,simpint1;

imports algebraiccase,algfnpl,findzvars,getvariables,interr,printsq,
  transcendentalcase,varsinlist,kernp,simpcar,prepsq,mksq,simp,
   opmtch,formlnr;

switch algint,nolnr,trint;

% 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.
% Extended arguments x1, x2, etc., are not currently supported.

symbolic procedure simpint u;
% Simplifies an integral. First two components of U are the integrand
% and integration variable respectively. Optional succeeding components
% are log forms for the final integral;
    begin scalar ans,expression,variable,loglist,w,
                 !*purerisch,intvar,listofnewsqrts,listofallsqrts,
                 sqrtfn,sqrt!-intvar,sqrt!-places!-alist,
                 basic!-listofallsqrts,basic!-listofnewsqrts;
    if atom u or null cdr u then rederr "Not enough arguments for INT";
    variable := !*a2k cadr u;
    w := cddr u;
    if w then rederr "Too many arguments to INT";
    listofnewsqrts:= list mvar gaussiani; % Initialize for SIMPSQRT.
    listofallsqrts:= list (argof mvar gaussiani . gaussiani);
    sqrtfn := get('sqrt,'simpfn);
    put('sqrt,'simpfn,'proper!-simpsqrt);
    % We need explicit settings of several switches during integral
    % evaluation.  In addition, the current code cannot handle domains
    % like floating point, so we suppress it while the integral is
    % calculated.  UNCACHED is turned on since integrator does its own
    % caching.
    begin scalar dmode!*,!*exp,!*gcd,!*keepsqrts,!*mcd,!*sqrt,
                 !*rationalize,!*structure,!*uncached;
       !*keepsqrts := !*sqrt := t;
       !*exp := !*gcd := !*mcd := !*structure := !*uncached := t;
       dmode!* := nil;
       if !*algint
         then <<intvar:=variable;           % until fix JHD code
            % Start a clean slate (in terms of SQRTSAVE) for this integral
            sqrt!-intvar:=!*q2f simpsqrti variable;
            if (red sqrt!-intvar) or (lc sqrt!-intvar neq 1)
                or (ldeg sqrt!-intvar neq 1)
              then interr "Sqrt(x) not properly formed"
              else sqrt!-intvar:=mvar sqrt!-intvar;
            basic!-listofallsqrts:=listofallsqrts;
            basic!-listofnewsqrts:=listofnewsqrts;
            sqrtsave(basic!-listofallsqrts,basic!-listofnewsqrts,
                         list(variable . variable))>>;
       expression := int!-simp car u;
   %   loglist := for each x in w collect int!-simp x;
       ans := errorset('(integratesq expression variable loglist),
                       !*backtrace,!*backtrace);
    end;
    if errorp ans
      then return <<put('sqrt,'simpfn,sqrtfn);
                    if !*failhard then error1();
                    simpint1(expression . variable . w)>>
     else ans := car ans;
    expression := sqrtchk numr ans ./ sqrtchk denr ans;
    % We now need to check that all simplifications have been done
    % but we have to make sure INT is not resimplified.
    put('int,'simpfn,'simpiden);
    ans := errorset('(resimp expression),t,!*backtrace);
    put('int,'simpfn,'simpint);
    put('sqrt,'simpfn,sqrtfn);
    return if errorp ans then error1() else car ans
   end;

symbolic procedure sqrtchk u;
   % U is a standard form. Result is another standard form with square
   % roots replaced by half powers.
   if domainp u then u
    else if not eqcar(mvar u,'sqrt)
     then addf(multpf(lpow u,sqrtchk lc u),sqrtchk red u)
    else addf(multpf(mksp(list('expt,cadr mvar u,'(quotient 1 2)),
                          ldeg u),
                     sqrtchk lc u),
              sqrtchk red u);

symbolic procedure int!-simp u;
   %converts U to canonical form, including the resimplification of
   % *sq forms;
   subs2 resimp 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); %important 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(distexp(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;
    if !*algint and cdr zlist and algfnpl(zlist,var)
      then return algebraiccase(integrand,zlist,varlist)
     else return transcendentalcase(integrand,var,xlogs,zlist,varlist)
 end;

symbolic procedure distexp(l);
    if null l then nil
    else if atom car l then car l . distexp cdr l
    else if (caar l = 'expt) and (cadar l = 'e) then 
        begin scalar ll;
            ll:=caddr car l;
            if eqcar(ll,'plus) then <<
                ll:=foreach x in cdr ll collect list('expt,'e,x);
                return ('times . ll) . distexp cdr l >>
            else return car l . distexp cdr l
        end
    else distexp car l . distexp cdr l;

symbolic procedure pseudodiff(a,var);
    if atom a then nil
    else if car a memq '(atan equal log plus quotient sqrt times)
        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 if car a eq 'expt
        then if depends(cadr a,var) then
            prepsq simp list('log,cadr a) .
            cadr a .
            caddr a .
            append(pseudodiff(cadr a,var),pseudodiff(caddr a,var))
        else caddr a . pseudodiff(caddr a,var)
    else list prepsq simpdf(list(a,var));

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),
                                 !*backtrace,!*backtrace);
                     if pairp u then v:=car u;
                     !*nolnr:= cdr !*nolnr;
                     return v>>;
      return if (v := opmtch u) then simp v
              else if !*failhard then rederr "FAILHARD switch set"
              else mksq(u,1)
   end;

mkop 'int!*;

put('int!*,'simpfn,'simpint!*);

symbolic procedure simpint!* u;
   begin scalar x;
      return if (x := opmtch('int . u)) then simp x
              else simpiden('int!* . u)
   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;   % Splitting of cubics and quartics.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(knowndiscrimsign zlist);

global '(!*trint);

exports cubic,quartic;

imports covecdf,cuberootf,nth,forceazero,makepolydf,multdf,multdfconst,
   !*multf,negdf,plusdf,printdf,printsf,quadratic,sqrtf,vp1,vp2,addf,
   negf;

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;
        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 <<
            if !*trint 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 .* (!*multf(d,d) ./ 1)) .+ nil;
        b:=plusdf(b,multdfconst(negf !*multf(a,d) ./ 1,p));
        b:=plusdf(b,multdfconst(!*multf(a,a) ./ 1,
                                multdf(p,p)));
        return quadratic(makepolydf b,var,res); %deal with what is left;
   prime:
        if !*trint then printc "The following cubic does not split";
  exit:
        if !*trint then printdf pol;
        return ('log . pol) . res
    end;

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,ee,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);
        ee:=getv(v,0);
        if !*trint then << printc "Quartic has coefficients";
            printsf a; printsf b;
            printsf c; printsf d;
            printsf ee >>;
        if d
          then <<if !*trint then printc "Quartic too hard to split";
                 go to exit >>;
        b:=c; c:=ee; %squash up the notation;
        if knowndiscrimsign eq 'negative then go to complex;
        dsc := 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;
        if a eq 'failed or dsc eq 'failed then go to prime;
        ee := invsq(addf(a,a) ./ 1);
        d := multsq(addf(b,negf dsc) ./ 1,ee);
        ee := multsq(addf(b,dsc) ./ 1,ee);
        if !*trint
          then <<printc "Quadratic factors will have coefficients";
                 printsf a; print 0; printsq d;
                 printc "or"; printsq ee>>;
        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 .* ee) .+ 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);
        if a eq 'failed or c eq 'failed then go to prime;
        b:=addf(!*multf(2,!*multf(a,c)),negf b);
        b:=sqrtf b;
        if b eq 'failed then go to prime;
%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:
        if !*trint then printc "The following quartic does not split";
   exit:
        if !*trint then printdf pol;
        return ('log . pol) . res
    end;

endmodule;


module factr;   % Crude factorization routine for integrator.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(zlist);

global '(!*trint);

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 distributed 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;   % Some basic support routines for integrator.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(!*backtrace !*gcd !*sqfree indexlist sqrtflag sqrtlist 
        varlist zlist);

global '(!*trint);

exports partialdiff,printdf,rationalintegrate,interr;

imports df2printform,printsf,varsinsf,addsq,multsq,multd,mksp;

symbolic procedure printdf u;
% Print distributed form via cheap conversion to reduce structure.
    begin scalar !*gcd;
       printsf df2printform u;
    end;


symbolic procedure !*n2sq(u1);
if u1=0 then nil . 1 else u1 . 1;


symbolic procedure indx(n);
if n<2 then (list 1) else(n . indx(isub1 n));


symbolic procedure interr mess;
   <<(!*trint or !*backtrace)
       and <<prin2 "***** INTEGRATION PACKAGE ERROR:  "; printc mess>>;
     error1()>>;


symbolic procedure rationalintegrate(x,var);
    begin         scalar n,d;
      n:=numr x; d:=denr x;
      if not var member varsinsf(d,nil) then
            return !*multsq(polynomialintegrate(n,var),1 ./ d);
      rederr "Rational integration not coded yet"
    end;


symbolic procedure polynomialintegrate(x,v);
% Integrate standard form. result is standard quotient.
    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 r
    end;

symbolic procedure partialdiff(p,v);
% Partial differentiation of p wrt v - p is s.f. as is result.
    if domainp 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);

symbolic procedure ncdr(l,n);
  % we can use small integer arithmetic here.
  if n=0 then l else ncdr(cdr l,isub1 n);


symbolic procedure mkilist(old,term);
   if null old then nil
    else term.mkilist(cdr old,term);


%symbolic procedure addin(lista,first,listb);
%if null lista
% then nil
% else ((first.car listb).car lista).addin(cdr lista,first,cdr listb);



symbolic procedure removeduplicates(u);
  % Purges duplicates from the list passed to it.
if null u then nil
  else if (atom u) then u.nil
    else if member(car u,cdr u)
      then removeduplicates cdr u
      else (car u).removeduplicates cdr u;


symbolic procedure jsqfree(sf,var);
begin
  varlist:=getvariables(sf ./ 1);
  zlist:=findzvars(varlist,list var,var,nil);
  sqrtlist:=findsqrts varlist; % before the purge
  sqrtflag:=not null sqrtlist;
  varlist:=purge(zlist,varlist);
  if sf eq !*sqfree
    then return list list sf
    else return sqfree(sf,zlist)
  end;

symbolic procedure jfactor(sf,var);
begin
  scalar varlist,zlist,indexlist,sqrtlist,sqrtflag;
  scalar prim,answer,u,v; % scalar var2
  prim:=jsqfree(sf,var);
  indexlist:=createindices zlist;
  prim:=factorlistlist (prim,t);
  while prim do <<
    if caar prim eq 'log then <<
        u:=cdar prim;
        u:=!*multsq(numr u ./ 1,1 ./ cdr stt(numr u,var));
        v:=denr u;
        while not atom v do v:=lc v;
        if  (numberp v) and (0> v)
          then u:=(negf numr u) ./ (negf denr u);
        answer:=u.answer >>
      else if caar prim eq 'atan then <<
% We can ignore this, since we also get LOG (U**2+1) as another term.
       >>
      else interr "Unexpected term in jfactor";
    prim:=cdr prim >>;
  return answer
  end;


  symbolic procedure stt(u,x);
    if domainp u
      then if u eq nil
        then ((-1) . nil)
        else (0 . u)
      else if mvar u eq x
        then ldeg  u . lc u
        else if ordop(x,mvar u)
          then (0 . u)
          else begin
            scalar ltlc,ltrest;
            ltlc:=stt(lc u,x);
            ltrest:= stt(red u,x);
            if car ltlc = car ltrest then go to merge;
            if car ltlc > car ltrest
              then return car ltlc .
                               !*multf(cdr ltlc,(lpow u .* 1) .+ nil)
              else return ltrest;
          merge:
            return car ltlc.addf(cdr ltrest,
                                 !*multf(cdr ltlc,(lpow u .* 1) .+ nil))
            end;


symbolic procedure gcdinonevar(u,v,x);
% Gcd of u and v, regarded as polynnmials in x.
if null u
  then v
  else if null v
    then u
    else begin
      scalar u1,v1,z,w;
      u1:=stt(u,x);
      v1:=stt(v,x);
    loop:
      if (car u1) > (car v1)
        then go to ok;
      z:=u1;u1:=v1;v1:=z;
      z:=u;u:=v;v:=z;
    ok:
      if car v1 iequal 0
        then interr "Coprimeness in gcd";
      z:=gcdf(cdr u1,cdr v1);
      w:=quotf(cdr u1,z);
      if (car u1) neq (car v1)
        then w:=!*multf(w,!*p2f mksp(x,(car u1)-(car v1)));
      u:=addf(!*multf(v,w),
              !*multf(u,negf quotf(cdr v1,z)));
      if null u
        then return v;
      u1:=stt(u,x);
      go to loop
      end;


symbolic procedure mapply(funct,l);
   if null l then rederr "Empty list to mapply"
    else if null cdr l then car l
    else apply(funct,list(car l,mapply(funct,cdr l)));


symbolic procedure !*lcm!*(u,v);
!*multf(u,quotf(v,gcdf(u,v)));


symbolic procedure multsql(u,l);
% Multiplies (!*multsq) each element of l by u.
if null l
  then nil
  else !*multsq(u,car l).multsql(u,cdr l);


symbolic procedure intersect(x,y);
if null x then nil else if member(car x,y) then
    car(x) . intersect(cdr x,y) else
          intersect(cdr x,y);


symbolic procedure mapvec(v,f);
begin
  scalar n;
  n:=upbv v;
  for i:=0:n do
    apply(f,list getv(v,i))
  end;

endmodule;


module jpatches;   % Routines for manipulating sf's with power folding.

% Author: James H. Davenport.

fluid '(sqrtflag);

exports !*addsq,!*multsq,!*invsq,!*multf,squashsqrtsq,!*exptsq,!*exptf;

% !*MULTF(A,B) multiplies the polynomials (standard forms) U and V
% in much the same way as MULTF(U,V) would, EXCEPT...
%     (1) !*MULTF inhibits the action of OFF EXP and of non-commutative
%         multiplications
%     (2) Within !*MULTF powers of square roots, and powers of
%         exponential kernels are reduced as if substitution rules
%         such as FOR ALL X LET SQRT(X)**2=X were being applied;

% Note that !*MULTF comes between MULTF and !*Q2F SUBS2F MULTF in its
% behaviour, and that it is the responsibility of the user to call it
% in sensible places where its services are needed;

% similarly for the other functions defined here;


%symbolic procedure !*addsq(u,v);
   %U and V are standard quotients.
%  %Value is canonical sum of U and V;
%   if null numr u then v
%    else if null numr v then u
%    else if denr u=1 and denr v=1 then addf(numr u,numr v) ./ 1
%    else begin scalar nu,du,nv,dv,x;
%        x := gcdf(du:=denr u,dv:=denr v);
%        du:=quotf(du,x); dv:=quotf(dv,x);
%        nu:=numr u; nv:=numr v;
%        u:=addf(!*multf(nu,dv),!*multf(nv,du));
%        if u=nil then return nil ./ 1;
%        v:=!*multf(du,denr v);
%        return !*ff2sq(u,v)
%    end;

%symbolic procedure !*multsq(a,b);
%  begin
%    scalar n,d;
%    n:=!*multf(numr a,numr b);
%    d:=!*multf(denr a,denr b);
%    return !*ff2sq(n,d)
%  end;

%symbolic procedure !*ff2sq(a,b);
%  begin
%    scalar gg;
%    if null a then return nil ./ 1;
%    gg:=gcdf(a,b);
%    if not (gg=1) then <<
%        a:=quotf(a,gg);
%        b:=quotf(b,gg) >>;
%    if minusf b then <<
%        a:=negf a;
%        b:=negf b >>;
%    return a ./ b
%  end;

symbolic procedure !*addsq(u,v);
   %U and V are standard quotients.
   %Value is canonical sum of U and V;
   if null numr u then v
    else if null numr v then u
    else if denr u=1 and denr v=1 then addf(numr u,numr v) ./ 1
    else begin scalar du,dv,x,y,z;
        x := gcdf(du:=denr u,dv:=denr v);
        du:=quotf(du,x); dv:=quotf(dv,x);
        y:=addf(!*multf(dv,numr u),!*multf(du,numr v));
        if null y then return nil ./ 1;
        z:=!*multf(denr u,dv);
        if minusf z then <<y := negf y; z := negf z>>;
        if x=1 then return y ./ z;
        x := gcdf(y,x);
        return if x=1 then y ./ z else quotf(y,x) ./ quotf(z,x)
    end;

symbolic procedure !*multsq(u,v);
   %U and V are standard quotients. Result is the canonical product of
   %U and V with surd powers suitably reduced.
   if null numr u or null numr v then nil ./ 1
    else if denr u=1 and denr v=1 then !*multf(numr u,numr v) ./ 1
    else begin scalar w,x,y;
     x := gcdf(numr u,denr v);
     y := gcdf(numr v,denr u);
     w := !*multf(quotf(numr u,x),quotf(numr v,y));
     x := !*multf(quotf(denr u,y),quotf(denr v,x));
     if minusf x then <<w := negf w; x := negf x>>;
     y := gcdf(w,x);  % another factor may have been generated.
     return if y=1 then w ./ x else quotf(w,y) ./ quotf(x,y)
    end;

symbolic procedure !*invsq a;
   % Note that several examples (e.g., int(1/(x**8+1),x)) give a more
   % compact result when SQRTFLAG is true if SQRT2TOP is not called.
   if sqrtflag then sqrt2top invsq a else invsq a;

symbolic procedure !*multf(u,v);
  % U and V are standard forms
  % Value is SF for U*V;
  begin
    scalar x,y;
      if null u or null v
      then return nil
      else if u = 1
        then return squashsqrt v
        else if v = 1
        then return squashsqrt u
        else if domainp u
          then return multd(u,squashsqrt v)
          else if domainp v
            then return multd(v,squashsqrt u);
      x:=mvar u;
      y:=mvar v;
      if x eq y
        then go to c
        else if ordop(x,y)
          then go to b;
      x:=!*multf(u,lc v);
      y:=!*multf(u,red v);
      return if null x then y
              else if not domainp lc v
                and mvar u eq mvar lc v
                and not atom mvar u
                and car mvar u memq '(expt sqrt)
               then addf(!*multf(x,!*p2f lpow v),y) % what about noncom?
              else makeupsf(lpow v,x,y);
  b:  x:=!*multf(lc u,v);
      y:=!*multf(red u,v);
      return if null x then y
              else if not domainp lc u
                and mvar lc u eq mvar v
                and not atom mvar v
                and car mvar v memq '(expt sqrt)
               then addf(!*multf(!*p2f lpow u,x),y)
              else makeupsf(lpow u,x,y);
  c:  y:=addf(!*multf(list lt u,red v),!*multf(red u,v));
      if eqcar(x,'sqrt)
        then return addf(squashsqrt y,!*multfsqrt(x,
                        !*multf(lc u,lc v),ldeg u + ldeg v))
        else if eqcar(x,'expt) and prefix!-rational!-numberp caddr x
          then return addf(squashsqrt y,!*multfexpt(x,
                        !*multf(lc u,lc v),ldeg u + ldeg v));
      x:=mkspm(x,ldeg u + ldeg v);
      return if null x or null (u:=!*multf(lc u,lc v))
               then y
               else x .* u .+ y
    end;

symbolic procedure makeupsf(u,x,y);
% Makes u .* x .+ y  except when u is not a valid lpow (because of
% sqrts).
   if atom car u or cdr u = 1 then u .* x .+ y
     else if caar u eq 'sqrt then addf(!*multfsqrt(car u,x,cdr u),y)
       else if <<begin scalar v;
                    v:=car u;
                    if car v neq 'expt then return nil;
                    v:=caddr v;
                    if atom v then return nil;
                    return (car v eq 'quotient 
                            and numberp cadr v
                            and numberp caddr v)
                 end >>
         then addf(!*multfexpt(car u,x,cdr u),y)
         else u .* x .+ y;


symbolic procedure !*multfsqrt(x,u,w);
   % This code (Due to Norman a& Davenport) squashes SQRT(...)**2.
   begin scalar v;
     w:=divide(w,2);
     v:=!*q2f simp cadr x;
     u:=!*multf(u,exptf(v,car w));
     if not zerop cdr w then u:=!*multf(u,!*p2f mksp(x,1));
     return u
   end;


symbolic procedure !*multfexpt(x,u,w);
  begin scalar expon,v;
    expon:=caddr x;
    x:=cadr x;
    w:=w * cadr expon;
    expon:=caddr expon;
    v:=gcdn(w,expon);
    w:=w/v;
    v:=expon/v;
    if not (w > 0) then rederr "Invalid exponent"
     else if v = 1
      then return !*multf(u,exptf(if numberp x then x
                                    else if atom x then !*k2f x
                                    else !*q2f if car x eq '!*sq
                                                 then argof x
                                                else simp x,
                          w));
    expon:=0;
    while not (w < v) do <<expon:=expon + 1; w:=w-v>>;
    if expon>0 then u:=!*multf(u,exptf(!*q2f simp x,expon));
    if w = 0 then return u;
    x:=list('expt,x,list('quotient,1,v));
    return multf(squashsqrt u,!*p2f mksp(x,w))
  end;


symbolic procedure prefix!-rational!-numberp u;
  % Tests for m/n in prefix representation.
    eqcar(u,'quotient) and numberp cadr u and numberp caddr u;


symbolic procedure squashsqrtsq sq;
!*multsq(squashsqrt numr sq ./ 1,
         1 ./ squashsqrt denr sq);


symbolic procedure squashsqrt sf;
if (not sqrtflag) or atom sf or atom mvar sf
  then sf
  else if car mvar sf eq 'sqrt and ldeg sf > 1
    then addf(squashsqrt red sf,!*multfsqrt(mvar sf,lc sf,ldeg sf))
    else if car mvar sf eq 'expt
       and prefix!-rational!-numberp caddr mvar sf
       and ldeg sf >= caddr caddr mvar sf
      then addf(squashsqrt red sf,!*multfexpt(mvar sf,lc sf,ldeg sf))
      else (lpow sf .* squashsqrt lc sf) .+ squashsqrt red sf;




%remd 'simpx1;

%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 0 = cdr z
%          then return simpexpt list(u,car z);
%        if n = 2
%          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;

symbolic procedure !*exptsq(a,n);
% raise A to the power N using !*MULTSQ;
    if n=0 then 1 ./ 1
    else if n=1 then a
    else if n<0 then !*exptsq(invsq a,-n)
    else begin
      scalar q,r;
      q:=divide(n,2);
      r:=cdr q; q:=car q;
      q:=!*exptsq(!*multsq(a,a),q);
      if r=0 then return q
      else return !*multsq(a,q)
    end;


symbolic procedure !*exptf(a,n);
% raise A to the power N using !*MULTF;
    if n=0 then 1
    else if n=1 then a
    else begin
      scalar q,r;
      q:=divide(n,2);
      r:=cdr q; q:=car q;
      q:=!*exptf(!*multf(a,a),q);
      if r=0 then return q
      else return !*multf(a,q)
    end;

endmodule;


module hacksqrt;  % Routines for manipulation of sqrt expressions.

% Author: James H. Davenport.

fluid '(nestedsqrts thisplace);

exports sqrtsintree,sqrtsinsq,sqrtsinsql,sqrtsinsf,sqrtsign;
exports degreenest,sortsqrts;

imports mkvect,interr,getv,dependsp,union;

symbolic procedure sqrtsintree(u,var,slist);
% Adds to slist all the sqrts in the prefix-type tree u.
if atom u
  then slist
  else if car u eq '!*sq
    then union(slist,sqrtsinsq(cadr u,var))
    else if car u eq 'sqrt
      then if dependsp(argof u,var)
        then <<
          slist:=sqrtsintree(argof u,var,slist);
          %  nested square roots
          if member(u,slist)
            then slist
            else u.slist >>
        else slist
      else sqrtsintree(car u,var,sqrtsintree(cdr u,var,slist));


symbolic procedure sqrtsinsq(u,var);
   % Returns list of all sqrts in sq.
   sqrtsinsf(denr u,sqrtsinsf(numr u,nil,var),var);


symbolic procedure sqrtsinsql(u,var);
% Returns list of all sqrts in sq list.
if null u
  then nil
  else sqrtsinsf(denr car u,
      sqrtsinsf(numr car u,sqrtsinsql(cdr u,var),var),var);


symbolic procedure sqrtsinsf(u,slist,var);
% Adds to slist all the sqrts in sf.
if domainp u or null u
  then slist
  else <<
    if  eqcar(mvar u,'sqrt) and
        dependsp(argof mvar u,var) and
        not member(mvar u,slist)
      then begin
        scalar slist2;
        slist2:=sqrtsintree(argof mvar u,var,nil);
        if slist2
          then <<
            nestedsqrts:=t;
            slist:=union(slist2,slist) >>;
        slist:=(mvar u).slist
        end;
    sqrtsinsf(lc u,sqrtsinsf(red u,slist,var),var) >>;


symbolic procedure easysqrtsign(slist,things);
% This procedure builds a list of all substitutions for all possible
% combinations of square roots in list.
if null slist
  then things
  else easysqrtsign(cdr slist,
                    nconc(mapcons(things,(car slist).(car slist)),
                          mapcons(things,
                                  list(car slist,'minus,car slist))));

symbolic procedure hardsqrtsign(slist,things);
% This procedure fulfils the same role for nested sqrts
% ***assumption: the simpler sqrts come further up the list.
if null slist
  then things
  else begin
    scalar thisplace,answers,pos,neg;
    thisplace:=car slist;
    answers:=mapcar(things,function (lambda u;sublis(u,thisplace).u));
    pos:=mapcar(answers,function (lambda u;(thisplace.car u).(cdr u)));
    % pos is sqrt(f) -> sqrt(innersubst f)
    neg:=mapcar(answers,
        function (lambda u;list(thisplace,'minus,car u).(cdr u)));
    % neg is sqrt(f) -> -sqrt(innersubst f)
    return hardsqrtsign(cdr slist,nconc(pos,neg))
    end;


symbolic procedure degreenest(pf,var);
% Returns the maximum degree of nesting of var
% inside sqrts in the prefix form pf.
if atom pf
  then 0
  else if car pf eq 'sqrt
    then if dependsp(cadr pf,var)
      then iadd1 degreenest(cadr pf,var)
      else 0
    else if car pf eq 'expt
      then if dependsp(cadr pf,var)
        then if eqcar(caddr pf,'quotient)
          then iadd1 degreenest(cadr pf,var)
          else degreenest(cadr pf,var)
        else 0
      else degreenestl(cdr pf,var);

symbolic procedure degreenestl(u,var);
%Returns max degreenest from list of pfs u.
if null u
  then 0
  else max(degreenest(car u,var),
           degreenestl(cdr u,var));


symbolic procedure sortsqrts(u,var);
% Sorts list of sqrts into order required by hardsqrtsign
% (and many other parts of the package).
begin
  scalar i,v;
  v:=mkvect(10); %should be good enough!
  while u do <<
    i:=degreenest(car u,var);
    if i iequal 0
      then interr "Non-dependent sqrt found";
    if i > 10
      then interr
         "Degree of nesting exceeds 10 (recompile with 10 increased)";
    putv(v,i,(car u).getv(v,i));
    u:=cdr u >>;
  u:=getv(v,10);
  for i :=9 step -1 until 1 do
    u:=nconc(getv(v,i),u);
  return u
  end;


symbolic procedure sqrtsign(sqrts,x);
   if nestedsqrts then hardsqrtsign(sortsqrts(sqrts,x),list nil)
    else easysqrtsign(sqrts,list nil);

endmodule;


module kron;   % Kronecker factorization of univ. polys over integers.

% Authors: Mary Ann Moore and Arthur C. Norman.

global '(!*trint);

exports linfac,quadfac;

imports zfactor;

% Only linear and quadratic factors are found.

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 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 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;

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 lowdeg;   % Splitting of low degree polynomials.

% Author: To be determined.

fluid '(clogflag knowndiscrimsign zlist);

global '(!*trint);

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;

symbolic procedure covecdf(pol,var,degree);
% Extract coefficients of polynomial wrt var, given a degree-bound
% degree. Result is a lisp vector.
    begin         scalar 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:
        if !*trint
         then <<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
        (denr lc a = 1) 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(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:= !*multsq(r ./ 1,invsq(!*multf(n,getv(p,n)) ./ 1));
           % Used to be subs2q multsq
                        %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 !*multsq(getv(p,i),w ./ 1));
           % Used to be subs2q multsq
        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,!*addsq(getv(w,im1),!*multsq(r,getv(w,i))));
           % Used to be subs2q addsq/multsq
            i:=im1; im1:=sub1 i >>;
        putv(w,0,!*multsq(getv(w,0),r));
           % Used to be subs2q multsq
        return w
    end;

endmodule;


module reform; % Reformulate expressions using C-constant substitution.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(cmap cval loglist ulist);

global '(!*trint);

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) .*
        !*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));
   % Next line had addsq previously
      return !*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 "SF arg for log etc ="; printc 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));
      % Next line had addsq previously
      logsq:=!*addsq(temp,logsq);
      lglst:=cdr lglst;
      i:=i+1;
      go to loop
   end;

endmodule;


module simplog;  % Simplify logarithms.

% Authors: Mary Ann Moore and Arthur C. Norman.

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,'simplog);

symbolic procedure simplogi(sq);
begin
   if atom sq then go to simplify;
   if car sq eq 'times
     then return simpplus(for each u in cdr sq
                            collect mk!*sq simplogi u);
   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 eq 'nthroot
     then return !*multsq(1 ./ caddr sq,simplogi cadr sq);
     % we had (nthroot of n).
   if car sq eq 'sqrt
      then return !*multsq(1 ./ 2,simplogi argof 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 formlogterm(sf);
begin
  scalar u;
  u:=mvar sf;
  if not atom u and (car u member '(times sqrt expt nthroot))
    then u:=addsq(simplog2 lc sf,
                    !*multsq(simplogi u,simp!* ldeg sf))
    else if (lc sf iequal 1) and (ldeg sf iequal 1)
      then u:=((mksp(list('log,u),1) .* 1) .+ nil) ./ 1
      else u:=addsq(simptimes list(list('log,u),ldeg sf),
                    simplog2 lc sf);
  return u
  end;

symbolic procedure formlog sf;
if null red sf
  then formlogterm 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;   % Simplify square roots.

% Authors: Mary Ann Moore and Arthur C. Norman.
% Heavily modified J.H. Davenport for algebraic functions.

fluid '(!*backtrace !*conscount !*galois !*pvar basic!-listofallsqrts
        gaussiani basic!-listofnewsqrts kord!* knowntobeindep
        listofallsqrts listofnewsqrts sqrtflag sqrtlist
        sqrt!-places!-alist variable varlist zlist);

global '(!*tra !*trint);

% This module should be rewritten in terms of the REDUCE function
% SIMPSQRT;

% remd 'simpsqrt;

exports proper!-simpsqrt,simpsqrti,simpsqrtsq,simpsqrt2,sqrtsave,
        newplace,actualsimpsqrt,formsqrt;

symbolic procedure proper!-simpsqrt(exprn);
   simpsqrti carx(exprn,'proper!-simpsqrt);


symbolic procedure simpsqrti sq;
begin
   scalar u;
   if atom sq
     then if numberp sq
       then return (simpsqrt2 sq) ./ 1
       else if (u:=get(sq,'avalue))
         then return simpsqrti cadr u
           % BEWARE!!! This is VERY system dependent.
         else return simpsqrt2((mksp(sq,1) .* 1) .+ nil) ./ 1;
           % If it doesnt have an AVALUE then it is itself;
   if car sq eq 'times
     then return mapply(function multsq,
                        mapcar(cdr sq,function simpsqrti));
   if car sq eq 'quotient
     then return multsq(simpsqrti cadr sq,
                        invsq simpsqrti caddr sq);
   if car sq eq 'expt and numberp caddr sq
     then if evenp caddr sq
       then return simpexpt list(cadr sq,caddr sq / 2)
       else return simpexpt
                     list(mk!*sq simpsqrti cadr sq,caddr sq);
  if car sq = '!*sq
    then return simpsqrtsq cadr sq;
  return simpsqrtsq tidysqrt simp!* sq
  end;


symbolic procedure simpsqrtsq sq;
(simpsqrt2 numr sq) ./ (simpsqrt2 denr sq);


symbolic procedure simpsqrt2 sf;
if minusf sf
  then if sf iequal -1
    then gaussiani
    else begin
      scalar u;
      u:=negf sf;
      if numberp u
        then return multf(gaussiani,simpsqrt3 u);
      % we cannot negate general expressions for the following reason:
%       (%%%thesis remark%%%)
%       sqrt(x*x-1) under x->1/x gives sqrt(1-x*x)/x=i*sqrt(x*x-1)/x
%                 under x->1/x gives x*i*sqrt(-1+1/x*x)=i**2*sqrt(x*x-1)
%       hence an abysmal catastrophe;
      return simpsqrt3 sf
      end
  else simpsqrt3 sf;


symbolic procedure simpsqrt3 sf;
begin
  scalar u;
  u:=assoc(sf,listofallsqrts);
  if u
    then return cdr u;
  % now see if 'knowntobeindep'can help.
  u:=atsoc(listofnewsqrts,knowntobeindep);
  if null u
    then go to no;
  u:=assoc(sf,cdr u);
  if u
    then <<
      listofallsqrts:=u.listofallsqrts;
      return cdr u >>;
no:
  u:=actualsimpsqrt sf;
  listofallsqrts:=(sf.u).listofallsqrts;
  return u
  end;


symbolic procedure sqrtsave(u,v,place);
begin
  scalar a;
  %u is new value of listofallsqrts, v of new.
  a:=assoc(place,sqrt!-places!-alist);
  if null a
    then sqrt!-places!-alist:=(place.(listofnewsqrts.listofallsqrts))
           .sqrt!-places!-alist
    else rplacd(a,listofnewsqrts.listofallsqrts);
      listofnewsqrts:=v;
      % throw away things we are not going to need in future.
      if not !*galois
        then listofallsqrts:=u;
        % we cannot guarantee the validity of our calculations.
      if listofallsqrts eq u
        then return nil;
      v:=listofallsqrts;
      while not (cdr v eq u) do
        v:=cdr v;
      rplacd(v,nil);
      % listofallsqrts is all those added since routine was entered.
      v:=atsoc(listofnewsqrts,knowntobeindep);
      if null v
        then knowntobeindep:=(listofnewsqrts.listofallsqrts)
                              . knowntobeindep
        else rplacd(v,union(cdr v,listofallsqrts));
      listofallsqrts:=u;
      return nil
  end;


symbolic procedure newplace(u);
% Says to restart algebraic bases at a new place u.
begin
  scalar v;
  v:=assoc(u,sqrt!-places!-alist);
  if null v
    then <<
      listofallsqrts:=basic!-listofallsqrts;
      listofnewsqrts:=basic!-listofnewsqrts >>
    else <<
      v:=cdr v;
      listofnewsqrts:=car v;
      listofallsqrts:=cdr v >>;
  return if v then v
              else listofnewsqrts.listofallsqrts
  end;


symbolic procedure mknewsqrt u;
% U is prefix form.
begin
  scalar v,w;
  if not !*galois then go to new;
    % no checking required.
  v:=addf(!*p2f mksp(!*pvar,2),negf !*q2f tidysqrt simp u);
% count !*conscount;
  w:=errorset(list('afactor,mkquote v,mkquote !*pvar),t,!*backtrace);
% if !*tra then <<
%   princ "*** That took ";
%   princ (!*conscount - count nil);
%   printc " conses" >>;
% count 0;
  if atom w
    then go to new
    else w:=car w; %the actual result of afactor.
  if cdr w
    then go to notnew;
new:
  w:=sqrtify u;
  listofnewsqrts:=w . listofnewsqrts;
  return !*kk2f w;
notnew:
  w:=car w;
  v:=stt(w,!*pvar);
  if car v neq 1
    then rederr "HELP ...";
  w:=addf(w,multf(cdr v,(mksp(!*pvar,car v) .* -1) .+nil));
  w:=sqrt2top(w ./ cdr v);
  w:=quotf(numr w,denr w);
  if null w
    then rederr "Division failure";
  return w
  end;


symbolic procedure actualsimpsqrt(sf);
if sf iequal -1
  then gaussiani
  else actualsqrtinner(sf,listofnewsqrts);


symbolic procedure actualsqrtinner(sf,l);
if null l
  then actualsimpsqrt2 sf
  else begin
    scalar z;
%   z:=simp argof car l;
%   if denr z neq 1 or (z := numr z) iequal -1
    z:=!*q2f simp argof car l;
    if z iequal -1
      then return actualsqrtinner(sf,cdr l);
    z:=quotf(sf,z);
    if null z
      then return actualsqrtinner(sf,cdr l);
    return !*multf(!*kk2f car l,actualsimpsqrt z)
    end;


symbolic procedure actualsimpsqrt2(sf);
 if atom sf
   then if null sf
     then nil
     else if numberp sf
       then if sf < 0
         then multf(gaussiani,actualsimpsqrt2(- sf))
           %Above 2 lines inserted JHD 4 Sept 80;
           % test case: SQRT(B*X**2-C)/SQRT(X);
         else begin
           scalar n;
           n:=int!-sqrt sf;
           % Changed for conformity with DEC20 LISP JHD July 1982;
           if not fixp n
             then return mknewsqrt sf
             else return n
           end
     else mknewsqrt(sf)
   else begin
     scalar form;
     form:=comfac sf;
     if car form
       then return multf((if null cdr sf and (car sf = form)
                            then formsqrt(form .+ nil)
                            else simpsqrt2(form .+ nil)),
                            %The above 2 lines changed by JHD;
                            %(following suggestions of Morrison);
                            %to conform to Standard LISP 4 Sept 80;
                         simpsqrt2 quotf(sf,form .+ nil));
     % we have killed common powers.
     form:=cdr form;
     if form neq 1
       then return multf(simpsqrt2 form,
                          simpsqrt2 quotf(sf,form));
     % remove a common factor from the sf.
     return formsqrt sf
     end;


symbolic procedure int!-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 int!-nr(n,(n+1)/2);


symbolic procedure int!-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 !*q2f simpsqrt list n;
    return int!-nr(n,w)
 end;


 symbolic procedure formsqrt(sf);
 if (null red sf)
   then if (lc sf iequal 1) and (ldeg sf iequal 1)
     then mknewsqrt mvar sf
     else multf(if evenp ldeg sf
                  then !*p2f mksp(mvar sf,ldeg sf / 2)
                  else exptf(mknewsqrt mvar sf,ldeg sf),simpsqrt2 lc sf)
   else begin
     scalar varlist,zlist,sqrtlist,sqrtflag;
     scalar v,l,n,w;
     % This returns a list, the i-th member of which is
%      a list of the factors of multiplicity i (as s.f's);
     v:=jsqfree(sf,if variable and involvesf(sf,variable)
                     then variable
                     else findatom mvar sf);
                     % VARIABLE is the best thing to do square-free
%                      decompositions with respect to, but anything
%                      else will do if VARIABLE is not set;
     if null cdr v and null cdar v
       then return mknewsqrt prepf sf;
       % The JSQFREE did nothing;
     l:=nil;
     n:=0;
     while v do <<
       n:=n+1;
       w:=car v;
       while w do <<
         l:=list('expt,mk!*sq !*f2q car w,n) . l;
         w:=cdr w >>;
       v:=cdr v >>;
     if null cdr l
       then l:=car l
       else l:='times.l;
       % makes L into a valid prefix form;
     return !*q2f simpsqrti l
     end;


symbolic procedure findatom pf;
if atom pf
  then pf
  else findatom argof pf;


symbolic procedure sqrtify u;
% Actually creates the sqrt.
begin
  scalar v,n,prinlist;
  n:=degreenest(u,nil);
% v:=list('sqrt,u);
  v := mksqrt u;   % To ensure sqrt**2 rule set.
  prinlist:=member(v,kord!*);
  if prinlist then return car prinlist;
  % This might improve sharing.
  % anyway, we mustn't re-add this object to KORD!*;
  while kord!* and
        eqcar(car kord!*,'sqrt) and
        (degreenest(argof car kord!*,nil) > n) do <<
    prinlist:=(car kord!*) . prinlist;
    kord!*:=cdr kord!* >>;
  kord!*:=v . kord!*;
  while prinlist do <<
    kord!*:=(car prinlist) . kord!*;
    prinlist:=cdr prinlist >>;
  return v
  end;


% deflist ('((sqrt (((x) quotient (sqrt x) (times 2 x))))),'dfn);

% put('sqrt,'simpfn,'proper!-simpsqrt); % Now in driver.

endmodule;


module isolve;   % Routines for solving the final reduction equation.

% Author: Mary Ann Moore and Arthur C. Norman.
% Modifications by: John P. Fitch.

fluid '(badpart
        ccount
        cmap
        cmatrix
        cval
        indexlist
        lhs!*
        lorder
        nnn
        orderofelim
        pt
        rhs!*
        sillieslist
        tanlist
        ulist
        zlist);

global '(!*number!* !*statistics !*trint);

exports solve!-for!-u;

imports nth,findpivot,gcdf,int!-gensym1,mkvect,interr,multdfconst,
   !*multf!*,negdf,orddf,plusdf,printdf,printsf,printspreadc,printsq,
   quotf,putv,spreadc,subst4eliminatedcs,mknill,pnth,domainp,addf,
   invsq,multsq;

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 .
           !*multsq(coefdf(lhs!*,lpowlhs),invsq cdr u)).ulist;
                % used to be subs2q multsq
        if !*statistics then !*number!*:=!*number!*+1;
         if !*trint then <<prin2 "***** U"; prin2 car u; prin2t " =";
                      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 coefft=0 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(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 car indvalues=0 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:=(int!-gensym1('!C).caar u).cc;
                                % CC is an alist of constants;
                    if !*trint then <<prin2 "***** U";
                                      prin2 caar u;
                                      prin2t " =";
                                      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;
            for each 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:=subs2q 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,nnn);
% 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 >>;
 a: if null r then go to b;
        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 go to b;
        r:=cdr r;
    go to a;
 b: if null r then return maxfrom(lhs!*,nnn)+1;
   return max(rr,maxfrom(lhs!*,nnn)+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 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;   % Square root of standard forms.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(!*noextend zlist);

exports minusdfp,sqrtdf,nrootn,domainp,minusf;

imports contentsmv,gcdf,interr,!*multf,partialdiff,printdf,quotf,
   simpsqrt2,vp2;

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 distributive form. "Failed" usually means
   % that the square root is not among already existing objects.
    if null l then nil
    else begin scalar c;
        if lpow l=vp2 zlist then go to ok;
        c:=sqrtsq df2q l;
        if numr c eq 'failed
          then return 'failed;
        if denr c eq 'failed
          then return 'failed;
        return for each u in f2df numr c
                 collect (car u).!*multsq(cdr u,1 ./ denr c);
    ok: c:=sqrtsq lc l;
        if  numr c eq 'failed or
            denr c eq 'failed
          then return 'failed
          else return (lpow l .* c) .+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 qp=1 then return ip; %exact root found.
         if !*noextend then return 'failed;
            % We cannot add new square roots in this case, since it is
            % then impossible to determine if one square root depends
            % on another if new ones are being added all the time.
         if zlistp qp then return 'failed;
            % Liouville's theorem tells you that you never need to add
            % new algebraics depending on the variable of integration.
        qp:=simpsqrt2 qp;
        return !*multf(ip,qp)
    end;

symbolic procedure zlistp qp;
if atom qp then member(qp,zlist)
  else or(member(mvar qp,zlist),zlistp lc qp,zlistp red qp);

symbolic procedure sqrtf1 p;
  % Returns a . b with p=a**2*b.
    if domainp p
     then if fixp p then nrootn(p,2)
           else !*q2f simpsqrt list prepf p . 1
    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:= 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, !*multf(car pg,car co));
        rplacd(pg, !*multf(cdr pg,cdr co));
        return pg
    end;

% NROOTN removed as in REDUCE base.

endmodule;


module tdiff;   % Differentiation routines for integrator.

% Authors: Mary Ann Moore and Arthur C. Norman.

exports !-!-simpdf;

imports simpcar,kernp,diffsq,prepsq,msgpri;

flag('(!-!-simpdf),'lose);

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;  % General tidying up of square roots.

% Authors: Mary Ann Moore and Arthur C. Norman.
% Modifications by J.H. Davenport.

exports sqrt2top,tidysqrt;

%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 return 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;


symbolic procedure multoutdenr q;
  % Move sqrts in a sq to the numerator.
    begin  scalar n,d,root,conj;
        n:=numr q;
        d:=denr q;
        while (root:=findsquareroot d) do <<
          conj:=conjugatewrt(d,root);
          n:=!*multf(n,conj);
          d:=!*multf(d,conj) >>;
        while (root:=findnthroot d) do <<
          conj:=conjugateexpt(d,root,kord!*);
          n:=!*multf(n,conj);
          d:=!*multf(d,conj) >>;
        return (n . d);
        end;

symbolic procedure conjugateexpt(d,root,kord!*);
  begin scalar ord,ans,repl,xi;
  ord:=caddr caddr root; % the denominator of the exponent;
  ans:=1;
  kord!*:= (xi:=gensym()) . kord!*;
  % XI is an ORD'th root of unity;
  for i:=1:ord-1 do <<
    ans:=!*multf(ans,numr subf(d,
                   list(root . list('times,root,list('explt,xi,i)))));
    while (mvar ans eq xi) and ldeg ans > ord do
      ans:=addf(red ans,(xi) to (ldeg ans - ord) .* lc ans .+ nil);
    if (mvar ans eq xi) and ldeg ans = ord then
      ans:=addf(red ans,lc ans) >>;
  if (mvar ans eq xi) and ldeg ans = ord-1 then <<
    repl:=-1;
    for i:=1:ord-2 do
      repl:=(xi) to i .* -1 .+ repl;
    ans:=addf(red ans,!*multf(lc ans,repl)) >>;
  if not domainp ans and mvar ans eq xi
    then interr "Conjugation failure";
  return ans;
  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 findnthroot p;
   nil;   % Until corrected.

symbolic procedure x!-findnthroot p;
    % Locate an n-th root 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,'expt) and eqcar(caddr w,'quotient) then return w;
        w:=findnthroot lc p;
        if null w then w:=findnthroot 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;  % Driving routine for integration of transcendental fns.

% Authors: Mary Ann Moore and Arthur C. Norman.
% Modifications by: John P. Fitch.

fluid '(!*backtrace
        !*nowarnings
        !*purerisch
        !*reverse
        badpart
        ccount
        cmap
        cmatrix
        content
        cuberootflag
        cval
        denbad
        denominator
        indexlist
        lhs!*
        loglist
        lorder
        orderofelim
        rhs!*
        sillieslist
        sqfr
        sqrtflag
        sqrtlist
        tanlist
        svar
        varlist
        xlogs
        zlist);

% !*reverse:       flag to re-order zlist.
% !*nowarnings:    flag to lose messages.

global '(!*failhard
         !*number!*
         !*ratintspecial
         !*seplogs
         !*spsize!*
         !*statistics
         !*trint
         gensymcount);

switch failhard;

exports transcendentalcase;

imports backsubst4cs,countz,createcmap,createindices,df2q,dfnumr,
  difflogs,fsdf,factorlistlist,findsqrts,findtrialdivs,gcdf,mkvect,
  interr,logstosq,mergin,multbyarbpowers,!*multf,multsqfree,
  printdf,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;

% Note that SEPLOGS keeps logarithmic part of result together as a
% kernel form, but this can lead to quite messy results.

symbolic 
   procedure transcendentalcase(integrand,svar,xlogs,zlist,varlist);
   begin scalar divlist,jhd!-content,content,prim,sqfr,dfu,indexlist,
%      JHD!-CONTENT is local, while CONTENT is free (set in SQFREE).
        sillieslist,originalorder,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,svar);
% *** now unnormalize integrand, maybe ***.
     begin scalar w,gg;
        gg:=1; 
        foreach z in zlist do <<
            w:=subs2 diffsq(simp z,svar);
            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:=for each u in eprim collect multup u;
%      sqfr:=multsqfree eprim;
%     if !*trint then << printc "...sqfr is:";
%         superprint sqfr>>;
      if !*reverse then zlist:=reverse zlist; %ALTER ORDER FUNCTION;
      indexlist:=createindices zlist;
%     if !*trint then << printc "...indices are:";
%         superprint indexlist>>;
      dfu:=dfnumr(svar,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,svar);
      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:= !*multf(numr dfu,
                                denbad:=quotf(denr dflogs,gcdq));
      denbad:=!*multf(denr dfu,denbad);
      denbad:= !*multf(denr unintegrand,denbad);
      dflogs:= !*multf(numr dflogs,quotf(denr dfu,gcdq));
      dfu:=dfun;
      % Now DFU and DFLOGS are S.F.s.
      rhs!*:=multbyarbpowers f2df dfu;
      if checkdffail(rhs!*,svar) then interr "Simplification failure";
      if !*trint then << printc "Distributed Form of U is:";
          printdf rhs!*>>;
      lhs!*:=f2df dflogs;
      if checkdffail(lhs!*,svar) then interr "Simplification failure";
      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;
        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:= !*multsq(dfun,!*invsq(denominator ./ 1));
      result:= !*multsq(result,!*invsq(jhd!-content ./ 1));
%     if !*trint then superprint result;
      dflogs:=logstosq();
      if not null numr dflogs then <<
          if !*seplogs and (not domainp numr result) then <<
              result:=mk!*sq result;
              result:=(mksp(result,1) .* 1) .+ nil;
              result:=result ./ 1 >>;
          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:= !*multsq(dfun,invsq(denbad ./ 1));
          if wrongway then << result:= nil ./ 1; dfun:=integrand >>;
          if rootcheckp(unintegrand,svar) then
                return simpint1(integrand . svar.nil)
          else if !*purerisch or allowedfns zlist then 
              dfun:=simpint1 (dfun . svar.nil)
           else << !*purerisch:=t;
                if !*trint
                  then <<printc "   [Transforming ..."; printsq dfun>>;
              denbad:=transform(dfun,svar);
              if denbad=dfun
                then dfun:=simpint1(dfun . svar.nil)
              else <<denbad:=errorset('(integratesq denbad svar xlogs),
                                      !*backtrace,!*backtrace);
                if not atom denbad then dfun:=untan car denbad
                else dfun:=simpint1(dfun . svar.nil) >> >>;
          if !*trint then printsq dfun;
          if !*failhard then rederr "FAILHARD switch set";
          if !*seplogs and not domainp result then <<
                result:=mk!*sq result;
                if not numberp result
                  then result:=(mksp(result,1) .* 1) .+ nil;
                result:=result ./ 1>>;
          result:=addsq(result,dfun) >>;
%      if !*overlaymode then excise transcode;
      return sqrt2top result
   end;

symbolic procedure checkdffail(u,v);
   u and (depends(lc u,v) or checkdffail(red u,v));

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 <<prin2 "Factors of multiplicity "; print i>>;
        wx:=car w;
        while not null wx do <<
            if !*trint then printsf car wx;
            for j:=1 : i do 
                denominator:= !*multf(car wx,denominator);
            wx:=cdr wx >>;
        i:=i+1;
        w:=cdr w;
        go to loopx
    end;

% unfluid '(dfun svar xlogs);

endmodule;


module halfangle;  % Routines for conversion to half angle tangents.

% Author: Steve Harrington.
% Modifications by: John P. Fitch.

fluid '(!*gcd);

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);

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; 
   % This should be done by a table.
   begin scalar y,z,w;
      if domainp x then return x . 1; 
      y := mvar x; 
      if eqcar(y,'int) then error1();  % 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 evenp z
                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);
   % 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 evenp z
       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 evenp z
       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 !*gcd;
   !*gcd := t;
   return multsq(cossqchkf(numr x), invsq cossqchkf(denr 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;  % Trial division routines.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(denominator loglist tanlist);

global '(!*trint);

exports countz,findsqrts,findtrialdivs,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 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;
% Make sure we find nested exponentials?
        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;  % Univariate factorization for integration.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(knowndiscrimsign zlist);

global '(!*trint);

exports unifac;

imports cubic,linfac,printdf,quadfac,quadratic,quartic,vp1,
        gcd,minusp,prettyprint;

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 evenp degree 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 p,q,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 ((for each j in caar r collect 2*j) . 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;

endmodule;


module uniform;  % Convert from D.F. to list of coefficients.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(zlist);

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(denr b=1) 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; % Make dummy variables for integration process.

% Authors: Mary Ann Moore and Arthur C. Norman.

fluid '(!*gensymlist!* !*purerisch);

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 they print badly on certain
% terminals etc and may lead to confusion.

!*gensymlist!* := '(! j ! k ! m ! n ! 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 varsinsq(sq,vl);
  varsinsf(numr sq,varsinsf(denr sq,vl));

symbolic procedure varsinsf(form,l);
   if domainp form then l
   else begin
     while not domainp 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 atom v and v neq var and depends(v,var) then 
           rederr "Can't integrate in the presence of side-relations"
         else if not atom v and (not v member zl) and dependsp(v,var)
         then if car v='!*sq then zl:=findzvarssq(cadr v,zl,var)
         else 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 eq 'df
          then <<!*purerisch := t; return zl>> % try and stop it
             else zl:=v . findzvars(cdr v,zl,var,flg);
                 % scan arguments of fn.
             %ACH: old code used to look only at CADR if a DF involved.
        go to scan
   end;

symbolic procedure findzvarssq(sq,zl,var);
    findzvarsf(numr sq,findzvarsf(denr sq,zl,var),var);

symbolic procedure findzvarsf(sf,zl,var);
    if domainp sf then zl
    else findzvarsf(lc sf,
                    findzvarsf(red sf,
                               findzvars(list mvar sf,zl,var,nil),
                               var),
                  var);

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 int!-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 vect;  % Vector support routines.

% Authors: Mary Ann Moore and Arthur C. Norman.
% Modified by: James H. Davenport.

exports mkuniquevect,mkvec,mkvecf2q,mkidenm,copyvec,vecsort,swap,
        non!-null!-vec,mkvect2;

symbolic procedure mkuniquevect v;
begin scalar u,n;
  n:=upbv v;
  for i:=0:n do begin
    scalar uu;
    uu:=getv(v,i);
    if not (uu member u)
      then u:=uu.u
    end;
  return mkvec u
  end;

symbolic procedure mkvec(l);
begin scalar v,i;
  v:=mkvect(isub1 length l);
  i:=0;
  while l do <<putv(v,i,car l); i:=iadd1 i; l:=cdr l>>;
  return v
  end;

symbolic procedure mkvecf2q(l);
begin
  scalar v,i,ll;
  v:=mkvect(isub1 length l);
  i:=0;
  while l do <<
    ll:=car l;
    if ll = 0 then ll:=nil;
    putv(v,i,!*f2q ll);
    i:=iadd1 i;
    l:=cdr l >>;
  return v
  end;

symbolic procedure mkidenm n;
begin
  scalar ans,u;
  scalar c0,c1;
  c0:=nil ./ 1;
  c1:= 1 ./ 1;
  % constants.
  ans:=mkvect(n);
  for i:=0 step 1 until n do <<
    u:=mkvect n;
    for j:=0 step 1 until n do
      if i iequal j
        then putv(u,j,c1)
        else putv(u,j,c0);
    putv(ans,i,u) >>;
  return ans
  end;

symbolic procedure copyvec(v,n);
   begin scalar new;
    new:=mkvect(n);
    for i:=0:n do putv(new,i,getv(v,i));
    return new
   end;

symbolic procedure vecsort(u,l);
% Sorts vector v of numbers into decreasing order.
% Performs same interchanges of all vectors in the list l.
begin
  scalar j,k,n,v,w;
  n:=upbv u;% elements 0...n exist.
  % algorithm used is a bubble sort.
  for i:=1:n do begin
    v:=getv(u,i);
    k:=i;
  loop:
    j:=k;
    k:=isub1 k;
    w:=getv(u,k);
    if v<=w
      then goto ordered;
    putv(u,k,v);
    putv(u,j,w);
    mapc(l,function (lambda u;swap(u,j,k)));
    if k>0
      then goto loop;
  ordered:
    end;
  return nil
  end;

symbolic procedure swap(u,j,k);
if null u
  then nil
  else begin
    scalar v;
    %swaps elements i,j of vector u.
    v:=getv(u,j);
    putv(u,j,getv(u,k));
    putv(u,k,v)
    end;

symbolic procedure non!-null!-vec v;
begin
  scalar cnt;
  cnt := 0;
  for i:=0:upbv v do
    if getv(v,i)
      then cnt:=iadd1 cnt;
  return cnt
  end;

symbolic procedure mkvect2(n,initial);
begin
  scalar u;
  u:=mkvect n;
  for i:=0:n do
    putv(u,i,initial);
  return u
  end;

endmodule;


end;