reduce-historical/r37/packages/support/0patches.red

module patches; % Patches to correct problems in REDUCE 3.7.

% Author: Anthony C. Hearn.

% Copyright (c) 2001 Anthony C. Hearn.  All Rights Reserved.

global '(patch!-date!*);

patch!-date!* := "7-Sep-2001";

% Bugs fixed by these patches.

% 28 Jun 99.  Gnuplot handling on the Macintosh was not correct.

%  7 Aug 99.  The evaluation of df(tan((sqrt(1-x^2)*asin acos x
%                + 2*sqrt(1-x^2)*x)/x),x) did not terminate.

% 20 Oct 99.  The sequence a1:=12x^2-16x+3; a2:=3x-4; off mcd;
%             on combineexpt; e^(a1/a2); gave the wrong answer.

%  8 Nov 99.  factorize(2*c*s*u^3*v^5-2*c*s*u^3*v +2*c*s*u*v^5-2*c*s*u*v
%               -s^2*u^4*v^4+s^2*u^4+s^2*u^2*v^6-s^2*u^2*v^4-s^2*u^2*v^2
%               +s^2*u^2 +s^2*v^6-s^2*v^2+u^4*v^4-u^4*v^2 -v^4+v^2) gave
%               a catastrophic error.

%  9 Nov 99.  Patched procedures generated a "redefined" message.

% 16 Nov 99.  Some EXCALC calculations could cause a catastrophic error.

% 18 Dec 99.  Integrations could give catastrophic errors because some
%             kernels were not unique.

% 31 Jan 00.  The sequence weight x=1,y=1; wtlevel 10; factor x; led to
%             the error that x was invalid as a kernel.

%  5 Feb 00.  The sequence x := mat((1,2)); sign sqrt 42; led to a
%             spurious error.

%  6 Feb 00.  The sequence on complex; sqrt(i*sqrt(3)-1); gave a wrong
%             result.

% 10 Feb 00.  Some root evaluations could lead to an error like
%             <equation> invalid as scalar.

% 18 Feb 00.  A sequence like m := mat((a,b),(c,d)); det sub(a=1,m);
%             would cause a type mismatch error.

% 18 Apr 00.  Complaints about the pattern matching limit of 5 terms
%             are resolved by the addition of a variable matchlength!*,
%             whose initial value of 5 can be changed as needed.

% 22 Apr 00.  The RULE mechanism left spurious expressions in various
%             non-local variables.

% 28 Jul 00.  A sum index within a derivative was treated as an identifier
%             (e.g., sum(x^n/factorial n*sub(x=0,df(cos x,x,n)),n,0,5);

%  2 Aug 00.  With complex on, some factorizations seemed to run forever
%             (e.g., factorize (400y^12+400y^10*z+40y^9*z^2+100y^8*z^2
%              +20y^7*z^5+120y^7*z^4+20y^7*z^3+41y^6*z^4+60y^5*z^7
%              +60y^5*z^5+20y^4*z^7+6y^4*z^6+20y^4*z^5
%              +2y^3*z^6+9y^2*z^8+6y*z^8+z^8))

% 29 Aug 00.  The sequence load_package gentran,scope; matrix a(10,10);
%             on gentranopt; gentran a(1,1) ::=: a(1,1); caused a
%             segmentation violation or similar error.

% 19 Sep 00.  Clearing some sqrt rules could lead to a spurious
%             "not found" message.

% 20 Sep 00.  The sequence load_package algint;
%             int(1/sqrt((2*e^c-y)/(e^c*y)),y);
%             caused a catastrophic error.

%  8 Nov 00.  Some sequences did not optimize completely when the SCOPE
%             command "optimize" was used.

% 20 Nov 00.  The sum operator did not always preserve a noncom order
%             (e.g., noncom u,v; sum(u(m)*v(1-m),m,0,1);)

% 12 Dec 00.  int with four arguments did not automatically load the
%             defint package.

% 13 Dec 00.  Some gcd calculations could produce an endless loop. E.g.,
%             in on numval,rounded; y:=x^4+x3*x^3+x2*x^2+x1*x+x0;
%             on fullroots; solve(y,x);

%  9 Jan 01.  SOLVE did not return results in same order as the given
%             variables (e.g., solve({y=x+t^2,x=y+u^2},{x,y,u,t});

% 14 Jan 01.  Some resultants (e.g. resultant(p^3-3p^2-a,3p*(p-2),p))
%             caused an error.

% 19 Jan 01.  Some algebraic integrals could produce a catastrophic
%             error when the algint package was loaded.

% 22 Jan 01.  Some algebraic integrals could produce a spurious zero
%             divisor message when the algint package was loaded (e.g.,
%               int((sqrt((-sqrt(a^4*x^2+4)+a^2*x)/(2*x))
%               *(-sqrt(a^4*x^2+4)*a^2*x-a^4*x^2-4))/(2*(a^4*x^2+4)),x))

% 23 Jan 01.  Inverses of matrices containing non-commuting objects
%             could be incorrect (e.g. noncom q; 1/mat((1,0,0),
%                (x/p*q 1,1,0),(x*y/(2p*(p-1))*q 1*q 1,y/(p-2)*q 1,1))).

%  2 Feb 01.  Some calls of SOLVE could produce a "zero divisor" error
%             error (e.g., solve(sqrt x*sqrt((4x^2*x+1)/x)-1=0,x)).

%  9 Feb 01.  The patched version of combine!-logs included an undefined
%             macro.

% 20 Feb 01.  Even with combineexpt on, expressions like a*a^x and
%             e*e^(2/(2-x)) did not simplify adequately.

%  6 Mar 01.  With algint loaded, some integrals would abort before
%             completion (e.g., int((x^(2/3)*sqrt(sqrt(y)*sqrt(pi) + 2pi
%             *y*x)*sqrt(- sqrt(y)*sqrt(pi)+2pi*y*x))/(4pi*y*x^3 - x),x)).

%  7 Apr 01.  Factorizations with rounded and complex rounded were not
%             supported.

% 23 Apr 01.  A term could be dropped from the solution of a set of
%             linear equations involving surds.

%  7 May 01.  Patches of 19 Sep 00 updated to correct bug.

%  1 Jun 01.  With precise on, sqrt(x^2) returned x rather than abs(x).

% 11 Jun 01.  Expressions involving noncommuting arguments of a
%             commuting operator were not handled correctly.
%             E.g., operator p,q; noncom q; p q a*p q b -p q b*p q a;

% 15 Jun 01.  Patch of 9 Jan 01 gave a catastrophic error with depend
%             statements, e.g., depend u,z; solve(u=y/x,z).

%  7 Aug 01.  Resultants could return an error or the wrong sign
%             (replaces patch of 14 Jan 01).

%  9 Aug 01.  Rules for tan, like let tan(~x)=>sin x/cos x, made some
%             integrals not return a closed form (e.g., int(1/sin(y)^2,y)).

% 24 Aug 01.  Some factorizations with complex on could have the wrong
%             sign.

%  7 Sep 01.  Some solve problems could cause a catastrophic error in
%             SQFRF.


% Alg declarations.

fluid '(matchlength!*);

matchlength!* := 5;

flag('(matchlength!*),'share);

fluid '(!*sqrtrulep);

patch alg;

% 20 Oct 99, 20 Feb 01.

symbolic procedure exptunwind(u,v);
   begin scalar x,x1,x2,y,z,z2;
  a:  if null v then return u;
      x := caar v;
      x1 := cadr x;
      x2 := caddr x;
      y := cdar v;
      v := cdr v;
      if !*combineexpt and length u=1 and null cdr(z2 := kernels u)
	then u := {(({'expt,car z2,ldeg u} . 1) . lc u)};
      while (z := assocp1(x1,v)) and
         (z2 := simp {'plus,{'times,x2,y},{'times,caddar z,cdr z}})
         and (!*combineexpt or (fixp numr z2 and fixp denr z2))
	do <<if fixp numr z2 and fixp denr z2
               then <<x2 := divide(numr z2,denr z2);
                      if car x2>0
			then <<if fixp x1 then u := multf(x1**car x2,u)
				else u := multpf(mksp(x1,car x2),u);
                               z2 := cdr x2 ./ denr z2>>;
                      y := numr z2>>
	      else y := 1;
             x2 := prepsq(quotf(numr z2,y) ./ denr z2);
             v := delete(z,v)>>;
      if !*combineexpt and y=1 and fixp x1 then
         <<while (z := assocp2(x2,v)) and cdr z=1 and fixp cadar z do
              <<x1 := cadar z * x1; v := delete(z,v)>>;
           if eqcar(x2,'quotient) and fixp cadr x2 and fixp caddr x2
		and cadr x2<caddr x2
             then <<z := nrootn(x1**cadr x2,caddr x2);
                    if cdr z = 1 then u := multd(car z,u)
                     else if car z = 1
                      then u := multf(formsf(x1,x2,1),u)
		     else <<u := multd(car z,u);
                            v := (list('expt,cdr z,x2) . 1) . v>>>>
            else u := multf(formsf(x1,x2,y),u)>>
       else u := multf(formsf(x1,x2,y),u);
      go to a
   end;

% 31 Jan 00.

symbolic procedure factor1(u,v,w);
   begin scalar x,y,z,r;
        y := lispeval w;
        for each j in u do
          if (x := getrtype j) and (z := get(x,'factor1fn))
              then apply2(z,u,v)
            else <<while eqcar(j:=reval j,'list) and cdr j do 
                     <<r:=append(r,cddr j); j:=cadr j>>;
		   x := !*a2kwoweight j;
                   if v then y := aconc!*(delete(x,y),x)
	            else if not(x member y)
                     then msgpri(nil,j,"not found",nil,nil)
                    else y := delete(x,y)>>;
        set(w,y);
        if r then return factor1(r,v,w)
   end;

% 5 Feb 00.

algebraic (let sign(sqrt ~a)  => 1 when sign a=1);

% 18 Feb 00.

symbolic procedure getrtype u;
   begin scalar x,y;
    return
    if null u then nil
     else if atom u
      then if not idp u then not numberp u and getrtype1 u
	    else if flagp(u,'share)
	     then if (x := eval u) eq u then nil else getrtype x
	    else if (x := get(u,'avalue)) and
		       not(car x memq '(scalar generic))
		    or (x := get(u,'rtype)) and (x := list x)
                    then if y := get(car x,'rtypefn) then apply1(y,nil)
                          else car x
                  else nil
     else if not idp car u then nil
     else if (x := get(car u,'avalue)) and (x := get(car x,'rtypefn))
      then apply1(x,cdr u)
     else if car u eq 'sub then 'yetunknowntype
     else getrtype2 u
   end;

symbolic procedure let3(u,v,w,b,flgg);
   begin scalar x,y1,y2,z;
        x := u;
        if null x then <<u := 0; return errpri1 u>>
         else if numberp x then return errpri1 u;
       y2 := getrtype v;
       if b and idp x then <<remprop(x,'rtype); remprop(x,'avalue)>>;
	if (y1 := getrtype x)
	  then return if z := get(y1,'typeletfn)
			then lispapply(z,list(x,v,y1,b,getrtype v))
		       else typelet(x,v,y1,b,getrtype v)
	 else if y2 and not(y2 eq 'yetunknowntype)
	  then return if z := get(y2,'typeletfn)
			then lispapply(z,list(x,v,nil,b,y2))
		       else typelet(x,v,nil,b,y2)
         else letscalar(u,v,w,x,b,flgg)
   end;

% 18 Apr 00.

symbolic procedure mcharg1(u,v,w);
   if null u and null v then list nil
    else begin integer m,n;
        m := length u;
        n := length v;
        if flagp(w,'nary) and m>2
	  then if m<=matchlength!* and flagp(w,'symmetric)
                             then return mchcomb(u,v,w)
                else if n=2 then <<u := cdr mkbin(w,u); m := 2>>
		else return nil;
        return if m neq n then nil
                else if flagp(w,'symmetric) then mchsarg(u,v,w)
                else if mtp v then list pair(v,u)
                else mcharg2(u,v,list nil,w)
   end;

% 19 Sep 00, 7 May 01, 15 Jun 01.

symbolic procedure letscalar(u,v,w,x,b,flgg);
   begin
      if not atom x
               then if not idp car x then return errpri2(u,'hold)
                     else if car x eq 'df
                      then if null letdf(u,v,w,x,b) then nil
                            else return nil
                     else if getrtype car x
                      then return let2(reval x,v,w,b)
                     else if not get(car x,'simpfn)
                      then <<redmsg(car x,"operator");
                             mkop car x;
                             return let3(u,v,w,b,flgg)>>
                     else nil
         else if null b and null w
          then <<remprop(x,'avalue);
		 remprop(x,'rtype);
                 remflag(list x,'antisymmetric);
                 remprop(x,'infix);
		 remprop(x,'kvalue);
		 remflag(list x,'linear);
		 remflag(list x,'noncom);
                 remprop(x,'op);
                 remprop(x,'opmtch);
                 remprop(x,'simpfn);
                 remflag(list x,'symmetric);
                 wtl!* := delasc(x,wtl!*);
                 if flagp(x,'opfn)
                   then <<remflag(list x,'opfn); remd x>>;
		 rmsubs();
                 return nil>>;
        if eqcar(x,'expt) and caddr x memq frlis!*
          then letexprn(u,v,w,!*k2q x,b,flgg)
         else if eqcar(x,'sqrt)
	  then <<!*sqrtrulep := t;
		 let2({'expt,cadr x,'(quotient 1 2)},v,w,b)>>;
	x := simp0 x where !*precise = t;
        return if not domainp numr x then letexprn(u,v,w,x,b,flgg)
                else errpri1 u
   end;

symbolic procedure setk1(u,v,b);
   begin scalar x,y,z,!*uncached;
      !*uncached := t;
      if atom u
        then <<if null b
                 then <<if not get(u,'avalue)
                          then msgpri(nil,u,"not found",nil,nil)
                         else remprop(u,'avalue);
                        return nil>>
		else if (x:= get(u,'avalue)) then put!-avalue(u,car x,v)
		else put!-avalue(u,'scalar,v);
               return v>>
       else if not atom car u
	then rerror(alg,25,"Invalid syntax: improper assignment");
      u := car u . revlis cdr u;
      if null b
        then <<z:=assoc(u,wtl!*);
               if not(y := get(car u,'kvalue))
                  or not (x := assoc(u,y))
		 then <<if null z and null !*sqrtrulep then
                            msgpri(nil,u,"not found",nil,nil)>>
                else put(car u,'kvalue,delete(x,y));
		if z then wtl!*:=delasc(u,wtl!*);
               return nil>>
       else if not (y := get(car u,'kvalue))
	then put!-kvalue(car u,nil,u,v)
       else <<if x := assoc(u,y)
		then <<updoldrules(u,v); y := delasc(car x,y)>>;
	      put!-kvalue(car u,y,u,v)>>;
      return v
     end;

% 7 May 01.

symbolic procedure clearrules u;
   begin scalar !*sqrtrulep;
      return rule!-list(u,nil)
   end;

% 2 Feb 01.

symbolic procedure simprad(u,n);
   if !*reduced then multsq(radfa(numr u,n),invsq radfa(denr u,n))
     else begin scalar iflag,x,y,z;
       if !*rationalize then <<
	  y:=list(denr u,1);
          u:=multf(numr u, exptf(denr u,n-1)) ./ 1 >>
         else y := radf(denr u,n);
       if n=2 and minusf numr u
	 then <<iflag := t; x := radf(negf numr u,n)>>
	else x := radf(numr u,n);
       z := simp list('quotient,retimes cdr x, retimes cdr y);
       if domainp numr z and domainp denr z
	 then z := multsq(mkrootsq(prepf numr z,n),
			  invsq mkrootsq(prepf denr z,n))
	else <<if iflag
		 then <<iflag := nil;
			z := negsq z>>;
	       z := mkrootsq(prepsq z,n)>>;
       z := multsq(multsq(if !*precise and evenp n 
			    then car x ./ 1
                           else car x ./ 1, 1 ./ car y), z);
       if iflag then z := multsq(z,mkrootsq(-1,2));
       return z
   end;

symbolic procedure radfa(u,n);
   begin scalar x,y;
      x := fctrf u;
      if numberp car x then x := append(zfactor car x,cdr x)
       else x := (car x ./ 1) . cdr x;
      y := 1 ./ 1;
      for each j in x do y := multsq(y,radfb(car j,cdr j,n));
      return y
   end;

symbolic procedure radfb(u,m,n);
   begin scalar x,y;
      x := radf(u,n);
      y := exptf(car x,m) ./ 1;
      return multsq(exptsq(mkrootlsq(cdr x,n),m),y)
   end;

% 20 Feb 01.

symbolic procedure reval2(u,v);
   if v or null !*combineexpt or dmode!* then !*q2a1(simp!* u,v)
    else !*q2a1((simp!* u where !*mcd = nil),v);

% 1 Jun 01.

symbolic procedure simpexpt2(u,n,flg);
   begin scalar m,n,x,y;
    if u=1 then return 1 ./ 1;
    m:=numr n;
    if pairp u then <<
     if car u eq 'expt
      then <<n:=multsq(m:=simp caddr u,n);
	     if !*precise
	       then u := list('abs,cadr u)
	      else u := cadr u;
	     return simpexpt1(u,n,flg)>>
     else if car u eq 'sqrt and not !*keepsqrts
      then return simpexpt2(cadr u, multsq(1 ./ 2,n),flg)
     else if car u eq 'times and not !*precise
      then <<x := 1 ./ 1;
	     for each z in cdr u do x := multsq(simpexpt1(z,n,flg),x);
	     return x>>
     else if car u eq 'times and (y:=split!-sign cdr u) and car y
      then <<x := simpexpt1(retimes append(cadr y,cddr y),n,flg);
             for each z in car y do x := multsq(simpexpt1(z,n,flg),x);
	     return x>>
     else if car u eq 'quotient
	and (not !*precise
		or posnump caddr u and posnump prepsq n
	    )
      then <<if not flg and !*mcd then
		return simpexpt1(prepsq simp!* u,n,t);
	     n := prepsq n;
	     return quotsq(simpexpt{cadr u,n},simpexpt{caddr u,n})>>
     else if car u eq 'minus and not !*precise and not(cadr u = 1)
      then return (multsq(simpexpt list(-1,expon),
		 simpexpt list(cadr u,expon))) where expon=prepsq n>>;
    if null flg
      then <<if null(dmode!* and idp u and get(u,dmode!*))
	       then u := prepsq simp!* u;
	     return simpexpt1(u,n,t)>>
     else if numberp u and zerop u then return nil ./ 1
     else if not numberp m then m := prepf m;
    n := prepf denr n;
    if m memq frlis!* and n=1 then return list ((u . m) . 1) . 1;
    if !*mcd or not numberp m or n neq 1
      or atom u or denr simp!* u neq 1 then return simpx1(u,m,n)
      else return mksq(u,m)
   end;

endpatch;

% Algint declarations.

fluid '(!*noacn !*structure !*tra !*trmin gaussiani intvar sqrtflag);

fluid '(!*pvar listofallsqrts listofnewsqrts);

global '(modevalcount);

patch algint;

% 20 Sep 00.

symbolic procedure algebraiccase(expression,zlist,varlist);
begin
  scalar rischpart,deriv,w,firstterm;
  scalar sqrtflag,!*structure;
  sqrtflag:=t;
  sqrtsave(listofallsqrts,listofnewsqrts,list(intvar . intvar));
  rischpart:= errorset!*(list('doalggeom,mkquote expression),
			 !*backtrace);
  newplace list (intvar.intvar);
  if atom rischpart
    then <<
      if !*tra then prin2t "Inner integration failed";
      deriv:=nil ./ 1;
      rischpart:=deriv >>
    else
      if atom car rischpart
        then <<
      if !*tra or !*trmin then
	prin2t "The 'logarithmic part' is not elementary";
          return (nil ./ 1) . expression >>
      else <<
        rischpart:=car rischpart;
    deriv:=!*diffsq(rischpart,intvar) where sqrtflag=nil;
    if !*tra or !*trmin then <<
      prin2t "Inner working yields";
          printsq rischpart;
      prin2t "with derivative";
          printsq deriv >> >>;
  deriv:=!*addsq(expression,negsq deriv);
  if null numr deriv
    then return rischpart . (nil ./ 1);
  if null involvesq(deriv,intvar)
    then return !*addsq(rischpart,
                !*multsq(deriv,((mksp(intvar,1) .* 1) .+ nil) ./ 1))
          . (nil ./ 1);
  varlist:=getvariables deriv;
  zlist:=findzvars(varlist,list intvar,intvar,nil);
  varlist:=setdiff(varlist,zlist);
  firstterm:=simp!* car zlist;
  w:=sqrt2top !*multsq(deriv,invsq !*diffsq(firstterm,intvar));
  if null involvesq(w,intvar)
    then return !*addsq(rischpart,!*multsq(w,firstterm)) . (nil ./ 1);
  if !*noacn then interr "Testing only logarithmic code";
  deriv:=transcendentalcase(deriv,intvar,nil,zlist,varlist);
  return !*addsq(car deriv, rischpart) . cdr deriv
  end;

% 22 Jan 01, 9 Feb 01.

symbolic procedure combine!-logs(coef,logarg);
begin
  scalar ans,dencoef,parts,logs,lparts,!*rationalize,trueimag;
  !*rationalize:=t;
  coef:=simp!* coef;
  if null numr logarg then return coef;
  parts:=split!-real!-imag numr coef;
  if null numr cdr parts then return multsq(coef,logarg);
  dencoef:=multf(denr coef,denr logarg);
  if !*tra then <<
     prin2t "attempting to find 'real' form for";
     mathprint list('times,list('plus,prepsq car parts,
                                      list('times,prepsq cdr parts,'i)),
			   prepsq logarg) >>;
  logarg:=numr logarg;
  logs:= 1 ./ 1;
  while pairp logarg do <<
	if ldeg logarg neq 1 then interr "what a log";
	if atom mvar logarg then interr "what a log";
	if car mvar logarg neq 'log then interr "what a log";
        logs:=!*multsq(logs,
		       !*exptsq(simp!* cadr mvar logarg,lc logarg));
	logarg:=red logarg >>;
  logs:=rationalizesq logs;
  ans:=multsq(!*multsq(car parts,logs),1 ./ dencoef);
  lparts:=split!-real!-imag numr logs;
  if numr difff(denr cdr lparts,intvar)
    then interr "unexpected denominator";
  lparts:=!*multsq(denr cdr lparts ./ 1,car lparts) . cdr lparts;
  if not onep denr car lparts then interr "unexpected denominator";
  trueimag:=quotsq(addf(!*exptf(numr car lparts,2),
                        !*exptf(numr cdr lparts,2)) ./ 1,
                   !*exptf(denr logs,2) ./ 1);
  if numr diffsq(trueimag,intvar)
     then ans:=!*addsq(ans,
                     !*multsq(gaussiani ./ multf(2,dencoef),
                              !*multsq(simplogsq trueimag,cdr parts)));
  trueimag:=!*multsq(car lparts,!*invsq(numr cdr lparts ./ 1));
  if numr diffsq(trueimag,intvar)
     then ans:=!*addsq(ans,!*multsq(!*multsq(cdr parts,1 ./ dencoef),
                                  !*k2q list('atan,prepsq!* trueimag)));
  return ans;
  end;

%  6 Mar 01.

symbolic procedure modevalvar v;
   begin scalar w;
      if atom v
	then <<if (w := get(v,'modvalue)) then return w;
	       put(v,'modvalue,modevalcount);
	       modevalcount := modevalcount+1;
	       return modevalcount-1>>
       else if car v neq 'sqrt
	then <<if !*tra then <<princ "Unexpected algebraic:"; print v>>;
	       error1()>>
       else if numberp cadr v then return (mksp(v,1) .* 1) .+ nil;
      w := modeval(!*q2f simp cadr v,!*pvar);
      w := assoc(w,listofallsqrts);
      if w then return cdr w else return 'failed
   end;

endpatch;

% Excalc declarations.

global '(basisforml!* detm!* indxl!* metricd!* metricu!*);

smacro procedure ldpf u;
   caar u;

smacro procedure lowerind u;
   list('minus,u);

patch excalc;

% 16 Nov 99.

symbolic procedure mkmetric u;
   begin scalar x,y,z,okord;
     putform(list(cadr u,nil,nil),0);
     put(cadr u,'indxsymmetries,
         '((lambda (indl) (tot!-sym!-indp 
                             (evlis '((nth indl 1) 
                                      (nth indl 2)))))));
     put(cadr u,'indxsymmetrize,
         '((lambda (indl) (symmetrize!-inds '(1 2) indl))));
     flag(list cadr u,'covariant);
     okord := kord!*;
     kord!* := basisforml!*;
     x := simp!* caddr u;
     y := indxl!*;
     metricu!* := t;
     for each j in indxl!* do
       <<for each k in y do
           setk(list(cadr u,lowerind j,lowerind k),0);
         y := cdr y>>;
     for each j on partitsq(x,'basep) do
       if ldeg ldpf j = 2 then
           setk(list(cadr u,lowerind cadr mvar ldpf j,
                            lowerind cadr mvar ldpf j),
                mk!*sq lc j)
        else
           setk(list(cadr u,lowerind cadr mvar ldpf j,
                            lowerind cadr mvar lc ldpf j),
                mk!*sq multsq(lc j,1 ./ 2));
     kord!* := okord;
     x := for each j in indxl!* collect
            for each k in indxl!* collect
               simpindexvar list(cadr u,lowerind j,lowerind k);
	 z := subfg!*;
     subfg!* := nil;
     y := lnrsolve(x,generateident length indxl!*);
	 subfg!* := z;
     metricd!* := mkasmetric x;
     metricu!* := mkasmetric y;
     detm!* := mk!*sq detq x
   end;

endpatch;


% Ezgcd declarations.

fluid '(image!-set reduced!-degree!-lclst unlucky!-case);

symbolic smacro procedure polyzerop u; null u;

patch ezgcd;

% 8 Nov 99.

symbolic procedure ezgcdf(u,v);
   begin scalar kordx,x;
      kordx := kord!*;
      x := errorset2{'ezgcdf1,mkquote u,mkquote v};
      if null errorp x then return first x;
      setkorder kordx;
      return gcdf1(u,v)
   end;
 
symbolic procedure poly!-gcd(u,v);
   begin scalar !*exp,z;
        if polyzerop u then return poly!-abs v
         else if polyzerop v then return poly!-abs u
         else if u=1 or v=1 then return 1;
        !*exp := t;
        if quotf1(u,v) then z := v
	 else if quotf1(v,u) then z := u
	 else if !*gcd then z := gcdlist list(u,v)
	 else z := 1;
        return poly!-abs z
   end;
 
symbolic procedure gcdlist3(l,onestep,vlist);
  begin
    scalar unlucky!-case,image!-set,gg,gcont,l1,w,w1,w2,
           reduced!-degree!-lclst,p1,p2;
    l1:=for each p in l collect p . ezgcd!-comfac p;
    l:=for each c in l1 collect
        quotfail1(car c,comfac!-to!-poly cdr c,
                  "Content divison in GCDLIST3 failed");
    gcont:=gcdlist for each c in l1 collect cddr c;
    if domainp gcont then if not(gcont=1)
      then errorf "GCONT has numeric part";
    l := sort(for each p in l collect poly!-abs p,function ordp);
    w := nil;
    while l do <<
       w := car l . w;
       repeat l := cdr l until null l or not(car w = car l)>>;
    l := reversip w;
    w := nil;
    if null cdr l then return multf(gcont,car l);
    if domainp (gg:=car (l:=sort(l,function degree!-order))) then
      return gcont;
    if ldeg gg=1 then <<
       if division!-test(gg,l) then return multf(poly!-abs gg,gcont)
       else return gcont >>;
    if onestep then <<
       p1 := poly!-abs car l; p2 := poly!-abs cadr l;
       if p1=p2 then <<
             if division!-test(p1,cddr l) then return multf(p1,gcont) >>
       else <<
       gg := poly!-gcd(lc p1,lc p2);
       w1 := multf(red p1, quotfail1(lc p2, gg,
        "Division failure when just one pseudoremainder step needed"));
       w2 := multf(red p2,negf quotfail1(lc p1, gg,
        "Division failure when just one pseudoremainder step needed"));
       w := ldeg p1 - ldeg p2;
       if w > 0 then w2 := multf(w2, (mksp(mvar p2, w) .* 1) .+ nil)
	else if w < 0
	 then w1 := multf(w1, (mksp(mvar p1, -w) .* 1) .+ nil);
       gg := ezgcd!-pp addf(w1, w2);
       if division!-test(gg,l) then return multf(gg,gcont) >>>>;
      return gcdlist31(l,vlist,gcont,gg,l1)
   end;

endpatch;

% Int declarations.

fluid '(!*failhard !*purerisch !*reverse !*trdint sqrt!-places!-alist
        badpart ccount cmap cmatrix content cval denbad denominator!*
        lhs!* loglist orderofelim rhs!* tanlist !*intflag!* indexlist
        intvar listofnewsqrts listofallsqrts lorder power!-list!*
        sqrtflag sqrtlist sqrt!-intvar basic!-listofnewsqrts
        basic!-listofallsqrts sqfr sillieslist varlist);

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

patch int;

% 18 Dec 99.

symbolic procedure findtrialdivs zl;
   begin scalar dlists1,args1;
      for each z in zl do
	 if exportan z
	   then <<if car z eq 'tan
		    then <<args1 := (mksp(z,2) .* 1) .+ 1;
			   tanlist := (args1 ./ 1) . tanlist>>
		   else args1 := !*kk2f z;
		  dlists1 := (z . args1) . dlists1>>;
      return dlists1
   end;

% 12 Dec 00.

symbolic procedure simpdint u;
   begin scalar low,upp,fn,var,x,y;
      if length u neq 4
	then rerror(int,2,"Improper number of arguments to INT");
      load!-package 'defint;
      fn := car u;
      var := cadr u;
      low := caddr u;
      upp := cadddr u;
      low := reval low;
      upp := reval upp;
      if low = upp then return nil ./ 1
       else if null getd 'new_defint then nil
       else if upp = 'infinity
	then if low = 0
	       then if not smemql('(infinity unknown),
				  x := defint!* {fn,var})
		      then return simp!* x else nil
	      else if low = '(minus infinity)
	       then return mkinfint(fn,var)
	      else if freeof(var,low)
	       then if not smemql('(infinity unknown),
				  x := defint!* {fn,var})
		     and not smemql('(infinity unknown),
				  y := indefint!* {fn,var,low})
		      then return simp!* {'difference,x,y} else nil
	      else nil
       else if upp = '(minus infinity) or low = 'infinity
	then return negsq simpdint {fn,var,upp,low}
       else if low = '(minus infinity)
	then return
	   simpdint{prepsq simp{'sub,{'equal,var,{'minus,var}},fn},
		     var,{'minus,upp},'infinity}
       else if low = 0
	then if freeof(var,upp)
		and not smemql('(infinity unknown),
			       x := indefint!* {fn,var,upp})
	       then return simp!* x else nil
       else if freeof(var,upp) and freeof(var,low)
		 and not smemq('(infinity unknown),
			       x := indefint!* {fn,var,upp})
		 and not smemql('(infinity unknown),
			       y := indefint!* {fn,var,low})
	then return simp!* {'difference,x,y};
      return mkdint(fn,var,low,upp)
   end;

% 1 Jun 01.

symbolic procedure simpint u;
   if atom u or null cdr u or cddr u and (null cdddr u or cddddr u)
     then rerror(int,1,"Improper number of arguments to INT")
    else if cddr u then simpdint u
    else begin scalar ans,dmod,expression,variable,loglist,oldvarstack,
		 !*intflag!*,!*purerisch,cflag,intvar,listofnewsqrts,
		 listofallsqrts,sqrtfn,sqrt!-intvar,sqrt!-places!-alist,
		 basic!-listofallsqrts,basic!-listofnewsqrts,coefft,
		 varchange,w,!*precise;
    !*intflag!* := t;
    variable := !*a2k cadr u;
    if not(idp variable or pairp variable and numlistp cdr variable)
      then <<varchange := variable . intern gensym();
	     if !*trint
	       then prin2t {"Integration kernel", variable,
			  "replaced by simple variable", cdr varchange};
	     variable := cdr varchange>>;
    intvar := variable;
    w := cddr u;
    if w then rerror(int,3,"Too many arguments to INT");
    listofnewsqrts:= list mvar gaussiani;
    listofallsqrts:= list (cadr mvar gaussiani . gaussiani);
    sqrtfn := get('sqrt,'simpfn);
    put('sqrt,'simpfn,'proper!-simpsqrt);
    if dmode!* then
       << if (cflag:=get(dmode!*, 'cmpxfn)) then onoff('complex, nil);
	  if (dmod := get(dmode!*,'dname)) then
	     onoff(dmod,nil)>> where !*msg := nil;
    begin scalar dmode!*,!*exp,!*gcd,!*keepsqrts,!*limitedfactors,!*mcd,
		 !*rationalize,!*structure,!*uncached,kord!*,
		 ans1,denexp,badbit,nexp,oneterm;
       !*keepsqrts := !*limitedfactors := t;
       !*exp := !*gcd := !*mcd := !*structure := !*uncached := t;
       dmode!* := nil;
       if !*algint
	 then <<
	    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))>>;
       coefft := (1 ./ 1);
       expression := int!-simp car u;
       if varchange
	 then <<depend1(car varchange,cdr varchange,t);
		expression := int!-subsq(expression,{varchange})>>;
       denexp := 1 ./ denr expression;
       nexp := numr expression;
       while not atom nexp and null cdr nexp and
	  not depends(mvar nexp,variable) do
	      <<coefft := multsq(coefft,(((caar nexp) . 1) . nil) ./ 1);
		nexp := lc nexp>>;
       ans1 := nil;
       while nexp do begin
	   scalar x,zv,tmp;
	   if atom nexp then <<x := !*f2q nexp; nexp := nil>>
	    else <<x := !*t2q car nexp; nexp := cdr nexp>>;
	   x := multsq(x,denexp);
	   zv := findzvars(getvariables x,list variable,variable,nil);
	   begin scalar oldzlist;
	       while oldzlist neq zv do <<
		    oldzlist := zv;
		    foreach zz in oldzlist do
		    zv:=findzvars(distexp(pseudodiff(zz,variable)),
				  zv,variable,t)>>;
	       % The following line was added to make, for example,
	       % int(df(sin(x)/x),x) return the expected result.
	       zv := sort(zv, function ordp)
	   end;
	   tmp := ans1;
	   while tmp do
	      <<if zv=caar tmp
		  then <<rplacd(car tmp,addsq(cdar tmp,x));
			 tmp := nil; zv := nil>>
		 else tmp := cdr tmp>>;
	   if zv then ans1 := (zv . x) . ans1
       end;
       if length ans1 = 1 then oneterm := t;
       nexp := ans1;
       ans := nil ./ 1;
       badbit:=nil ./ 1;
       while nexp do
	<<u := cdar nexp;
	  if !*trdint
	    then <<princ "Integrate"; printsq u;
		   princ "with Zvars "; print caar nexp>>;
	  ans1 := errorset!*(list('integratesq,mkquote u,
			     mkquote variable,mkquote loglist,
			     mkquote caar nexp),
			     !*backtrace);
	  nexp := cdr nexp;
	  if errorp ans1 then badbit := addsq(badbit,u)
	   else <<ans := addsq(caar ans1, ans);
		  badbit:=addsq(cdar ans1,badbit)>>>>;
       if !*trdint
	 then <<prin2 "Partial answer="; printsq ans;
		prin2 "To do="; printsq badbit>>;
       if badbit neq '(nil . 1)
	 then <<setkorder nil;
		badbit := reordsq badbit;
		ans := reordsq ans;
		coefft := reordsq coefft;
	  if !*trdint then <<princ "Retrying..."; printsq badbit>>;
	  if oneterm and ans = '(nil . 1) then ans1 := nil
	    else ans1 := errorset!*(list('integratesq,mkquote badbit,
				  mkquote variable,mkquote loglist,nil),
				  !*backtrace);
	  if null ans1 or errorp ans1
	    then ans := addsq(ans,simpint1(badbit . variable . w))
	   else <<ans := addsq(ans,caar ans1);
	      if not smemq(variable, ans) then ans := nil ./ 1;
		  if cdar ans1 neq '(nil . 1)
		    then ans := addsq(ans,
				    simpint1(cdar ans1 . variable . w))
		>>>>;
    end;
    ans := multsq(coefft,ans);
    if !*trdint then << prin2t "Resimp and all that"; printsq ans >>;
    put('int,'simpfn,'simpiden);
    put('sqrt,'simpfn,sqrtfn);
    << if dmod then onoff(dmod,t);
       if cflag then onoff('complex,t)>> where !*msg := nil;
    oldvarstack := varstack!*;
    varstack!* := nil;
    ans := errorset!*(list('int!-resub,mkquote ans,mkquote
			   varchange),t);
    put('int,'simpfn,'simpint);
    varstack!* := oldvarstack;
    return if errorp ans then error1() else car ans
   end;

symbolic procedure mkdint(fn,var,low,upp);
   begin scalar x,!*precise;
      if getd 'defint0
	 and not((x := defint0 {fn,var,low,upp}) eq 'failed)
	then return simp x
       else if not smemq('infinity,low) and not smemq('infinity,upp)
	then <<x := prepsq!* simpint {fn,var};
	       if not eqcar(x,'int)
		 then return simp!* {'difference,
				     subeval {{'equal,var,upp},x},
				     subeval {{'equal,var,low},x}}>>;
      return mksq({'int,fn,var,low,upp},1)
   end;

%  9 Aug 01.

symbolic 
   procedure transcendentalcase(integrand,svar,xlogs,zlist,varlist);
   begin scalar divlist,jhd!-content,content,prim,sqfr,dfu,indexlist,
      sillieslist,originalorder,wrongway,power!-list!*,
      sqrtlist,tanlist,loglist,dflogs,eprim,dfun,unintegrand,
      sqrtflag,badpart,rhs!*,lhs!*,gcdq,cmap,cval,orderofelim,cmatrix;
      scalar ccount,denominator!*,result,denbad,temp;
        gensymcount:=0;
      integrand:=sqrt2top integrand;
      if !*trint then << printc "Extension variables z<i> are";
          print zlist>>;
     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 "After unnormalization the integrand is ";
                printsq unintegrand >>
      end;
      divlist := findtrialdivs zlist;
      sqrtlist := findsqrts zlist;
      divlist := trialdiv(denr unintegrand,divlist);
      prim := sqfree(cdr divlist,zlist);
      jhd!-content := content;
      printfactors(prim,nil);
      eprim := sqmerge(countz car divlist,prim,nil);
      printfactors(eprim,t);
      sqfr := for each u in eprim collect multup u;
      if !*reverse then zlist := reverse zlist;
      indexlist := createindices zlist;
      dfu:=dfnumr(svar,car divlist);
      loglist := append(loglist,factorlistlist prim);
      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);
      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;
      rhs!* := multbyarbpowers f2df dfu;
      if checkdffail(rhs!*,svar)
        then <<if !*trint then printsq checkdffail(rhs!*,svar);
      interr "Simplification fails on above expression">>;
      if !*trint then << 
          printc "Distributed Form of Numerator is:";
          printdf rhs!*>>;
      lhs!* := f2df dflogs;
      if !*trint then << printc "Distributed Form of integrand is:";
          printdf lhs!*;
          terpri()>>;
      cval := mkvect(ccount);
      for i := 0:ccount do putv(cval,i,nil ./ 1);
      power!-list!* := tansfrom(rhs!*,zlist,indexlist,0);
      lorder:=maxorder(power!-list!*,zlist,0);
      originalorder := for each x in lorder collect x;
        if !*trint then << 
                printc "Maximum order for variables determined as ";
                print lorder >>;
        if !*statistics then << !*number!*:=0;
                !*spsize!*:=1;
                foreach xx in lorder do
                   !*spsize!*:=!*spsize!* * (xx+1) >>;
      dfun:=solve!-for!-u(rhs!*,lhs!*,nil);
      backsubst4cs(nil,orderofelim,cmatrix);
        if !*statistics then << prin2 !*number!*; prin2 " used out of ";
                printc !*spsize!* >>;
      badpart:=substinulist badpart;
      dfun:=df2q substinulist dfun;
      result:= !*multsq(dfun,!*invsq(denominator!* ./ 1));
      result:= !*multsq(result,!*invsq(jhd!-content ./ 1));
      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  <<
          terpri();
          printc
          "*****************************************************";
          printc
           "************ THE INTEGRAL IS : **********************";
          printc
           "*****************************************************";
          terpri();
          printsq result;
          terpri()>>;
      if badpart then begin scalar n,oorder;
          if !*trint
            then printc "plus a part which has not been integrated";
          lhs!*:=badpart;
          lorder:=maxorder(power!-list!*,zlist,0);
          oorder:=originalorder;
          n:=length lorder;
          while lorder do <<
                if car lorder > car originalorder then
                        wrongway:=t;
                if car lorder=car originalorder then n:= n-1;
                lorder:=cdr lorder;
                originalorder:=cdr originalorder >>;
          if !*trint and wrongway then printc "Went wrong way";
          dfun:=df2q badpart;
          dfun:= !*multsq(dfun,invsq(denbad ./ 1));
	  badpart := dfun;
          if wrongway then <<
                if !*trint then printc "Resetting....";
                result:=nil ./ 1;
                dfun := integrand; badpart:=dfun >>;
          if rootcheckp(unintegrand,svar) then
                return simpint1(integrand . svar.nil) . (nil ./ 1)
          else if !*purerisch or allowedfns zlist then 
              << badpart := dfun;
		 dfun := nil ./ 1 >>
           else << !*purerisch:=t;
                if !*trint
                  then <<printc "   Applying transformations ..."; 
                         printsq dfun>>;
	      temp := get('tan,'opmtch);
	      remprop('tan,'opmtch);
              denbad:=transform(dfun,svar);
              if denbad=dfun
		then <<dfun:=nil ./ 1;
		       badpart:=denbad;
		       put('tan,'opmtch,temp)>>
             else <<denbad:=errorset!*(list('integratesq,mkquote denbad,
                                      mkquote svar,mkquote xlogs, nil),
                                      !*backtrace);
		    put('tan,'opmtch,temp);
                     if not atom denbad then <<
			denbad:=car denbad;
                        dfun:=untan car denbad;

                        if (dfun neq '(nil . 1)) then
                            badpart:=untan cdr denbad;
                        if car badpart and not(badpart=denbad) then <<
                          wrongway:=nil;
                          lhs!*:=f2df car badpart;
                          lorder:=maxorder(power!-list!*,zlist,0);
                          n:=length lorder;
                          while lorder do <<
                                if car lorder > car oorder then
                                        wrongway:=t;
                                if car lorder=car oorder then n:= n-1;
                                lorder:=cdr lorder;
                                oorder:=cdr oorder >>;
                          if wrongway or (n=0) then <<
                               if !*trint then printc "Still backwards";
                               dfun := nil ./ 1;
                               badpart := integrand>>>>>>
                     else <<badpart := dfun; dfun := nil ./ 1 >>>>>>;
          if !*failhard then rerror(int,9,"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) end
         else badpart:=nil ./ 1;
      return (sqrt2top result . badpart)
   end;

endpatch;

unfluid '(indexlist);

patch limits;

% 1 Jun 01.

symbolic procedure simplimit u;
   begin scalar fn,exprn,var,val,old,v,!*precise,!*protfg;
     if length u neq 4
       then rerror(limit,1,
		   "Improper number of arguments to limit operator");
     fn:= car u; exprn := cadr u; var := !*a2k caddr u; val := cadddr u;
     !*protfg := t;
     old := get('cot,'opmtch);
     put('cot,'opmtch,
	 '(((!~x) (nil . t) (quotient (cos !~x) (sin !~x)) nil)));
     v := errorset!*({'apply,mkquote fn,mkquote {exprn,var,val}},nil);
     put('cot,'opmtch,old);
     !*protfg := nil;
     return if errorp v or (v := car v) = aeval 'failed then mksq(u,1)
	     else simp!* v
   end;

endpatch;

% Matrix declarations.

fluid '(!*bezout);

patch matrix;

% 7 Aug 99, 23 Apr 01.

symbolic procedure quotfexf!*1(u,v);
   if null u then nil
    else (if x then x
	   else (if denr y = 1 then numr y
		  else if denr (y := (rationalizesq y
				     where !*rationalize = t))=1
		   then numr y
		  else rerror(matrix,11,
			      "Catastrophic division failure"))
		 where y=rationalizesq(u ./ v))
	  where x=quotf(u,v);

% 23 Jan 01.

symbolic procedure lnrsolve(u,v);
   begin scalar temp,vlhs,vrhs,ok,
                !*exp,!*solvesingular;
   if !*ncmp then return clnrsolve(u,v);
   !*exp := t;
   if asymplis!* or wtl!* then
    <<temp := asymplis!* . wtl!*;
      asymplis!* := wtl!* := nil>>;
   vlhs := for i:=1:length car u collect intern gensym();
   vrhs := for i:=1:length car v collect intern gensym();
   u := car normmat augment(u,v);
   v := append(vlhs,vrhs);
   ok := setkorder v;
   u := foreach r in u collect prsum(v,r);
   v := errorset!*({function solvebareiss, mkquote u,mkquote vlhs},t);
   if caar v memq {'singular,'inconsistent} then 
      <<setkorder ok; rerror(matrix,13,"Singular matrix")>>;
   v := pair(cadr s,car s) where s = cadar v;
   u := foreach j in vlhs collect
	   coeffrow(negf numr q,vrhs,denr q) where q = cdr atsoc(j,v);
   setkorder ok;
   if temp then <<asymplis!* := car temp; wtl!* := cdr temp>>;
   return for each j in u collect
             for each k in j collect
                if temp then resimp k else cancel k;
   end;

% 23 Apr 01.

symbolic procedure extmult(u,v);
   if null u or null v then nil
    else (if x then cdr x .* (if car x then negf c!:subs2multf(lc u,lc v)
			       else c!:subs2multf(lc u,lc v))
			  .+ extadd(extmult(!*t2f lt u,red v),
				    extmult(red u,v))
	   else extadd(extmult(red u,v),extmult(!*t2f lt u,red v)))
	  where x = ordexn(car lpow u,lpow v);

% 14 Jan 01, 7 Aug 01.

symbolic procedure resultant(u,v,var);
   if domainp u and domainp v then 1
    else begin scalar x;
       kord!* := var . kord!*;
       if null domainp u and null(mvar u eq var) then u := reorder u;
       if null domainp v and null(mvar v eq var) then v := reorder v;
       x := if !*bezout then bezout_resultant(u,v,var)
	     else polyresultantf(u,v,var);
       setkorder cdr kord!*;
       return x
   end;

symbolic procedure resultantsq(u,v,var);
   if domainp numr u and domainp numr v and denr u = 1 and denr v = 1
     then 1 ./ 1
    else begin scalar x;
       kord!* := var . kord!*;
       if null domainp numr u and null(mvar numr u eq var)
	 then u := reordsq u;
       if null domainp numr v and null(mvar numr v eq var)
	 then v := reordsq v;
       x := if !*bezout then bezout_resultant(!*q2f u,!*q2f v,var)
	     else polyresultantf(!*q2f u,!*q2f v,var);
       setkorder cdr kord!*;
       return !*f2q x
   end;

symbolic procedure polyresultantf(u,v,var);
   begin scalar beta,cd,cn,delta,gam,r,s,temp,x;
      cd := cn := r := s := 1;
      gam := -1;
      if domainp u or domainp v then return 1
       else if ldeg u<ldeg v
	then <<s := (-1)^(ldeg u*ldeg v); temp := u; u := v; v := temp>>;
      while v do
       <<delta := ldeg u-ldegr(v,var);
	 beta := negf(multf(r,exptf(gam,delta)));
	 r := lcr(v,var);
	 gam := multf(exptf(negf r,delta),exptf(gam,1-delta));
	 temp := u;
	 u := v;
	 if not evenp ldeg temp and not evenp ldegr(u,var) then s := -s;
	 v := quotf(pseudo_remf(temp,v,var),beta);
	 if v
	   then <<cn := multf(cn,exptf(beta,ldeg u));
		  cd := multf(cd,
		     exptf(r,(1+delta)*ldeg u-ldeg temp+ldegr(v,var)));
		  if (x := quotf(cd,cn)) then <<cn := 1; cd := x>>>>>>;
      return if not domainp u and mvar u eq var then nil
	      else if ldeg temp neq 1
	       then quotf(multf(s,multf(cn,exptf(u,ldeg temp))),cd)
	      else u
   end;

symbolic procedure lcr(u,var);
   if domainp u or mvar u neq var then u else lc u;

symbolic procedure ldegr(u,var);
   if domainp u or mvar u neq var then 0 else ldeg u;

symbolic procedure pseudo_remf(u,v,var);
   !*q2f simp pseudo!-remainder {mk!*sq(u ./ 1),mk!*sq(v ./ 1),var};

symbolic procedure bezout_resultant(u,v,w);
   begin integer n,nm; scalar ap,ep,uh,ut,vh,vt;
     if domainp u or null(mvar u eq w)
        then return if not domainp v and mvar v eq w
                       then exptf(u,ldeg v)
                     else 1
      else if domainp v or null(mvar v eq w)
        then return if mvar u eq w then exptf(v,ldeg u) else 1;
     n := ldeg v - ldeg u;
     if n < 0 then return multd((-1)**(ldeg u*ldeg v),
                                bezout_resultant(v,u,w));
     ep := 1;
     nm := ldeg v;
     uh := lc u;
     vh := lc v;
     ut := if n neq 0 then multpf(w to n,red u) else red u;
     vt := red v;
     ap := addf(multf(uh,vt),negf multf(vh,ut));
     ep := b!:extmult(!*sf2exb(ap,w),ep);
     for j := (nm - 1) step -1 until (n + 1) do
        <<if degr(ut,w) = j then
             <<uh := addf(lc ut,multf(!*k2f w,uh));
               ut := red ut>>
           else uh := multf(!*k2f w,uh);
          if degr(vt,w) = j then
             <<vh := addf(lc vt,multf(!*k2f w,vh));
               vt := red vt>>
           else vh := multf(!*k2f w,vh);
          ep := b!:extmult(!*sf2exb(addf(multf(uh,vt),
                                    negf multf(vh,ut)),w),ep)>>;
      if n neq 0
         then <<ep := b!:extmult(!*sf2exb(u,w),ep);
                for j := 1:(n-1) do
                ep := b!:extmult(!*sf2exb(multpf(w to j,u),w),ep)>>; 
      return if null ep then nil else lc ep
   end;

endpatch;

% Ncpoly declarations.

fluid '(!*complex !*trnc dipvars!*);

patch ncpoly;

%  9 Jan 01.

symbolic procedure nc_factsolve(s,vl,all);
  begin scalar v,sb,ns,so,soa,sol,nz,w,q,z,r,abort;
   v:= numr simp car vl;
   ns:=for each e in s collect numr simp e;
   r:=t;
   while r do
   <<r:=nil; s:=ns; ns:=nil;
     for each e in s do if not abort then
     <<e:=absf numr subf(e,sb);
       while(q:=quotf(e,v)) do e:=q;
       if null e then nil else
       if domainp e or not(mvar e member vl) then abort:=t else
       if null red e and domainp lc e then
       <<w:=mvar e; sb:=(w . 0).sb; r:=t;
         vl:=delete(w,vl)>>
       else if not member(e,ns) then ns:=e.ns
     >>;
   >>;
   if abort or null vl then return nil;
   nc_factorize_timecheck();
   if null ns and vl then 
   <<sol:={for each x in vl collect x.1};
     goto done>>;
   s:=for each e in ns collect prepf e;
   if !*trnc then
    <<prin2 "solving ";
      prin2 length s; prin2 " polynomial equations for ";
      prin2 length vl;
      prin2t "variables";
      for each e in s do writepri(mkquote e,'only);>>;
   w:=(cdr solveeval{'list.s,'list.vl} where dipvars!*=nil);
 loop:
   nc_factorize_timecheck();
   if null w then goto done;
   so:= cdr car w; w:=cdr w; soa:=nil;
   if smemq('i,so) and null !*complex then go to loop;
   for each y in vl do if not smember(y,so) then
       <<soa:=(y . 1) . soa; nz:=t>>;
   for each y in so do
   <<z:=nc_factorize_unwrap(reval caddr y,soa); 
     nz:=nz or z neq 0;
     soa:=(cadr y . z).soa;
   >>;
   if not nz then goto loop;
   q:=assoc(car vl,soa);
   if null q or cdr q=0 then go to loop;
   soa := for each j in soa collect (car j . sublis(soa,cdr j));
   sol := soa . sol;
   if all then go to loop;
 done:
   sol:=for each s in sol collect append(sb,s);
   if !*trnc then
    <<prin2t "solutions:";
      for each w in sol do
       writepri(mkquote('list.
         for each s in w collect {'equal,car s,cdr s}),'only);
      prin2t "-------------------------";
    >>;
   return sol
  end;

endpatch;

% Plot declarations.

global '(!*plotpause !*plotusepipe dirchar!* opsys!* plotcleanup!*
	 plotcmds!* plotcommand!* plotdir!* plotdta!* plotheader!*
	 tempdir!*);

patch plot;

% 28 Jun 99.

symbolic procedure init_gnuplot();
 <<
!*plotpause := -1;
plotcleanup!* := {};
tempdir!* := getenv 'tmp;
if null tempdir!* then tempdir!* := getenv 'temp;
dirchar!* := "/";
plotcommand!* := "gnuplot";
opsys!* := assoc('opsys, lispsystem!*);
if null opsys!* then opsys!* := 'unknown
else opsys!* := cdr opsys!*;
if getenv "gnuplot" then plotdir!* := getenv "gnuplot"
 else if null plotdir!* and not (opsys!* = 'unix)
  then plotdir!* := get!-lisp!-directory();
if opsys!* = 'win32 then <<
    plotcommand!* := "wgnuplot";
    plotheader!* := "";
    dirchar!* := "\";
    plotdta!* := for each n in
       {"gnutmp.tm1", "gnutmp.tm2", "gnutmp.tm3", "gnutmp.tm4",
        "gnutmp.tm5", "gnutmp.tm6", "gnutmp.tm7", "gnutmp.tm8"}
       collect gtmpnam n;
    plotcleanup!* := if null tempdir!* then {"erase gnutmp.tm*"}
                     else {bldmsg("erase %w\gnutmp.tm*", tempdir!*)} >>
else if opsys!* = 'msdos then <<
    plotheader!* := "";      % ?? "set term vga";
    dirchar!* := "\";
    plotdta!* := for each n in
       {"gnutmp.tm1", "gnutmp.tm2", "gnutmp.tm3", "gnutmp.tm4",
        "gnutmp.tm5", "gnutmp.tm6", "gnutmp.tm7", "gnutmp.tm8"}
       collect gtmpnam n;
    plotcmds!*:= gtmpnam "gnutmp.tm0";
    plotcleanup!* := if null tempdir!* then {"erase gnutmp.tm*"}
                     else {bldmsg("erase %w\gnutmp.tm*", tempdir!*)} >>
else if opsys!* = 'riscos then <<
    plotheader!* := "";
    dirchar!* := ".";
    plotdta!* := for i:=1:10 collect tmpnam();
    plotcmds!*:= tmpnam();
    plotcleanup!* :=
       bldmsg("remove %w", plotcmds!*) .
          for each f in plotdta!* collect bldmsg("remove %w", f)
    >>
else if opsys!* = 'unix then <<
    plotheader!* := "set term x11";
    plotdta!* := for i:=1:10 collect tmpnam();
    plotcmds!*:= tmpnam();
    plotcleanup!* :=
       bldmsg("rm %w", plotcmds!*) .
          for each f in plotdta!* collect bldmsg("rm %w", f) >>
else if opsys!* = 'finder then <<
    plotcommand!* := "gnuplot";
    plotcmds!*:= "::::gnuplot:reduce.plt";
    plotheader!* := "";
    dirchar!* := ":";
    plotdta!* := for each n in
       {"::::gnuplot:gnutmp.tm1", "::::gnuplot:gnutmp.tm2",
        "::::gnuplot:gnutmp.tm3", "::::gnuplot:gnutmp.tm4",
        "::::gnuplot:gnutmp.tm5", "::::gnuplot:gnutmp.tm6",
        "::::gnuplot:gnutmp.tm7", "::::gnuplot:gnutmp.tm8"}
       collect gtmpnam n;
    plotcleanup!* := nil  >>
else <<
    rederr bldmsg("gnuplot for %w not available yet", opsys!*);
    plotdta!* := for i:=1:10 collect tmpnam();
    plotcmds!*:= tmpnam();
    plotheader!* := "set term dumb" >>;
if 'pipes member lispsystem!* then !*plotusepipe:=t
else plotcommand!* := bldmsg("%w %w", plotcommand!*, plotcmds!*);
if plotdir!* then
    plotcommand!* := bldmsg("%w%w%w",
                            plotdir!*, dirchar!*, plotcommand!*);
   nil >>;

endpatch;

patch poly;

% 7 Aug 99.

symbolic procedure rationalizesq u;
   begin scalar !*structure,!*sub2,v,x;
      if x := get(dmode!*,'rationalizefn) then u := apply1(x,u);
      powlis!* := '(i 2 (nil . t) -1 nil) . powlis!*;
      v := subs2q u;
      powlis!* := cdr powlis!*;
      return if domainp denr v then v
	      else if (x := rationalizef denr v) neq 1
	       then <<v := multf(numr v,x) ./ multf(denr v,x);
		      if null !*algint and null !*rationalize 
                        then v := gcdchk v;
		      subs2q v>>
	     else u
   end;

%  6 Feb 00, 7 Sep 01.

symbolic procedure sqfrf u;
   begin integer n; scalar !*gcd,units,v,w,x,y,z,!*msg,r;
      !*gcd := t;
      if (r := !*rounded) then
         <<on rational; u := numr resimp !*f2q u>>;
      n := 1;
      x := mvar u;
      v := gcdf(u,diff(u,x));
      u := quotf(u,v);
      if flagp(dmode!*,'field) and ((y := lnc u) neq 1)
	then <<u := multd(!:recip y,u); v := multd(y,v)>>;
      while degr(v,x)>0 do
       <<w := gcdf(v,u);
         if u neq w then z := (quotf(u,w) . n) . z;
         v := quotf(v,w);
         u := w;
         n := n + 1>>;
         if r then
            <<on rounded;
	      u := numr resimp !*f2q u;
	      z := for each j in z
		       collect numr resimp !*f2q car j . cdr j>>;
      if v neq 1 and assoc(v,units) then v := 1;
      if v neq 1 then if n=1 then u := multf(v,u)
       else if (w := rassoc(1,z)) then rplaca(w,multf(v,car w))
       else if null z and ((w := rootxf(v,n)) neq 'failed)
	then u := multf(w,u)
       else if not domainp v then z := aconc(z,v . 1)
       else errach {"sqfrf failure",u,n,z};
      return (u . n) . z
   end;

% 2 Aug 00.

symbolic procedure sqfrp u;
   begin scalar !*ezgcd, dmode!*;
     if null getd 'ezgcdf1 then load_package ezgcd;
     !*ezgcd := t;
     return domainp gcdf!*(u,diff(u,mvar u))
   end;

% 13 Dec 00.

symbolic procedure gcdk(u,v);
   begin scalar lclst,var,w,x;
        if u=v then return u
         else if domainp u or degr(v,(var := mvar u))=0 then return 1
         else if ldeg u<ldeg v then <<w := u; u := v; v := w>>;
        if quotf1(u,v) then return v
         else if !*heugcd and (x := heu!-gcd(u,v)) then return x
         else if ldeg v=1
           or getd 'modular!-multicheck and modular!-multicheck(u,v,var)
           or not !*mcd
          then return 1;
    a:  w := remk(u,v);
        if null w then return v
         else if degr(w,var)=0 then return 1;
        lclst := addlc(v,lclst);
        if x := quotf1(w,lc w) then w := x
         else for each y in lclst do
	    if atom y and not flagp(dmode!*,'field)
	      or not
	       (domainp y and (flagp(dmode!*,'field)
		  or ((x := get(car y,'units))
		       and y member (for each z in x collect car z))))
	    then while (x := quotf1(w,y)) do w := x;
        u := v; v := prim!-part w;
        if degr(v,var)=0 then return 1 else go to a
   end;

% 19 Jan 01.

symbolic procedure quarticf pol;
   begin scalar !*sub2,a,a2,a0,b,dsc,p,p1,p2,q,shift,var;
      var := mvar pol;
      p := shift!-pol pol;
      a := coeffs car p;
      shift := caddr p;
      if cadr a then rerror(poly,16,list(pol,"not correctly shifted"))
	else if cadddr a then return list(1,pol);
      a2 := cddr a;
      a0 := caddr a2;
      a2 := car a2;
      a := car a; 
      q := quadraticf1(a,a2,a0);
      if not(q eq 'failed)
	then <<a2 := car q; q := cdr q;
	       a := exptsq(addsq(!*k2q mvar pol,shift),2);
	       b := numr subs2q quotsq(addsq(multsq(!*f2q car q,a),
					     !*f2q cadr q),
				       !*f2q cadr p);
	       a := numr subs2q quotsq(addsq(multsq(!*f2q caddr q,a),
					     !*f2q cadddr q),
				       !*f2q cadr p);
	       a := quadraticf!*(a,var);
	       b := quadraticf!*(b,var);
	       return multf(a2,multf(car a,car b))
			 . nconc!*(cdr a,cdr b)>>
       else if null !*surds or denr shift neq 1
	then return list(1,pol);
      shift := numr shift;
      if knowndiscrimsign eq 'negative then go to complex;
      dsc := powsubsf addf(exptf(a2,2),multd(-4,multf(a,a0)));
      p2 := minusf a0;
      if not p2 and minusf dsc then go to complex;
      p1 := not a2 or minusf a2;
      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 := rootxf(a,2);
      if a eq 'failed then return list(1,pol);
      dsc := rootxf(dsc,2);
      if dsc eq 'failed then return list(1,pol);
      p := invsq !*f2q addf(a,a);
      q := multsq(!*f2q addf(a2,negf dsc),p);
      p := multsq(!*f2q addf(a2,dsc),p);
      b := multf(a,exptf(addf(!*k2f mvar pol,shift),2));
      a := powsubsf addf(b,q);
      b := powsubsf addf(b,p);
      knowndiscrimsign := p1;
      a := quadraticf!*(a,var);
      knowndiscrimsign := p2;
      b := quadraticf!*(b,var);
      knowndiscrimsign := nil;
      return multf(car a,car b) . nconc!*(cdr a,cdr b);
   complex:
      a := rootxf(a,2);
      if a eq 'failed then return list(1,pol);
      a0 := rootxf(a0,2);
      if a0 eq 'failed then return list(1,pol);
      a2 := powsubsf addf(multf(2,multf(a,a0)),negf a2);
      a2 := rootxf(a2,2);
      if a2 eq 'failed then return list(1,pol);
      p := addf(!*k2f mvar pol,shift);
      q := addf(multf(a,exptf(p,2)),a0);
      p := multf(a2,p);
      a := powsubsf addf(q,p);
      b := powsubsf addf(q,negf p);
      knowndiscrimsign := 'negative;
      a := quadraticf!*(a,var);
      b := quadraticf!*(b,var);
      knowndiscrimsign := nil;
      return multf(car a,car b) . nconc!*(cdr a,cdr b)
   end;

%  7 Apr 01.

symbolic procedure factorize u;
   (begin scalar x,y;
      x := simp!* u;
      y := denr x;
      if not domainp y then typerr(u,"polynomial");
      u := numr x;
      if u = 1 then return
	 {'list, if !*nopowers then 1 else {'list,1,1}}
       else if fixp u then !*ifactor := t;
      if !*force!-prime and not primep !*force!-prime
	then typerr(!*force!-prime,"prime");
      u := if dmode!* and not(dmode!* memq '(!:rd!: !:cr!:))
	 then if get(dmode!*,'factorfn)
                 then begin scalar !*factor;
                        !*factor := t;
			 return fctrf u
                      end
	   else rerror(poly,14,
		       list("Factorization not supported over domain",
		       get(dmode!*,'dname)))
       else fctrf u;
      return facform2list(u,y)
   end) where !*ifactor = !*ifactor;

% 7 Apr 01, 24 Aug 01.

symbolic procedure factor!-prim!-sqfree!-f u;
   begin scalar x,y,!*msg,r;
      r := !*rounded;
      if r and univariatep numr u and lc numr u=1 and denr u=1
	then return unifactor u
       else if r or !*complex or !*rational then
	 <<if r then on rational;
	   u := numr resimp !*f2q car u . cdr u>>;
      if null !*limitedfactors
        then <<if null dmode!* then y := 'factorf
                else <<x := get(dmode!*,'sqfrfactorfn);
                       y := get(dmode!*,'factorfn);
                       if x and not(x eq y) then y := 'factorf>>;
               if y
                 then <<y := apply1(y,car u);
			u := (exptf(car y,cdr u) . for each j in cdr y
                                           collect(car j . cdr u));
                        go to ret>>>>;
      u := factor!-prim!-sqfree!-f!-1(car u,cdr u);
 ret: if r then
      <<on rounded;
        u := car u . for each j in cdr u collect
           (numr resimp !*f2q car j . cdr j)>>;
      return u
   end;

% 7 Apr 01.

symbolic procedure unifactor u;
   if not eqcar(u := root_val list mk!*sq u,'list)
     then errach {"unifactor1",u}
    else 1 . for each j in cdr u collect
		 if not eqcar(j,'equal) then errach{"unifactor2",u}
		  else addsq(!*k2q cadr j,negsq simp caddr j);

% 11 Jun 01.

symbolic procedure noncomp u; 
   !*ncmp and noncomp1 u;

symbolic procedure noncomp1 u;
  if null pairp u then nil
   else if pairp car u then noncomfp u
   else flagp(car u,'noncom) or noncomlistp cdr u;

symbolic procedure noncomlistp u;
   pairp u and (noncomp1 car u or noncomlistp cdr u);

endpatch;

% Rlisp declarations.

fluid '(newrules!*);

patch rlisp;

%  9 Nov 99.

symbolic procedure load!-package u;
   begin scalar x,y;
      if stringp u then return load!-package intern compress explode2 u
       else if null idp u then rederr list(u,"is not a package name")
       else if memq(u,loaded!-packages!*)
        then return u
       else if or(atom(x:= errorset(list('evload,list('quote,list u)),
                               nil,!*backtrace)),
                  cdr x)
        then rederr
           list("error in loading package",u,"or package not found");
      loaded!-packages!* := u . loaded!-packages!*;
      x := get(u,'package);
      if x then x := cdr x;
   a: if null x then go to b
       else if null atom get(car x,'package) then load!-package car x
       else if or(atom(y := errorset(list('evload,
                                         list('quote,list car x)),
                                    nil,!*backtrace)),
                  cdr y)
        then rederr list("module",car x,"of package",u,
                         "cannot be loaded");
      x := cdr x;
      go to a;
   b: if (x := get(u,'patchfn))
	then begin scalar !*usermode,!*redefmsg; eval list x end
   end;

% 22 April 00.

symbolic procedure begin11 x;
   begin scalar mode,result,newrule!*;
      if cursym!* eq 'end
	 then if terminalp() and null !*lisp!_hook
		then progn(cursym!* := '!*semicol!*, !*nosave!* := t,
			   return nil)
	       else progn(comm1 'end, return 'end)
       else if eqcar((if !*reduce4 then x else cadr x),'retry)
	then if programl!* then x := programl!*
	      else progn(lprim "No previous expression",return nil);
      if null !*reduce4 then progn(mode := car x,x := cadr x);
      program!* := x;
      if eofcheck() then return 'c else eof!* := 0;
      add2inputbuf(x,if !*reduce4 then nil else mode);
      if null atom x
	  and car x memq '(bye quit)
	then if getd 'bye
	       then progn(lispeval x, !*nosave!* := t, return nil)
	      else progn(!*byeflag!* := t, return nil)
       else if null !*reduce4 and eqcar(x,'ed)
	then progn((if getd 'cedit and terminalp()
		      then cedit cdr x
		     else lprim "ED not supported"),
		   !*nosave!* := t, return nil)
       else if !*defn
	then if erfg!* then return nil
	      else if null flagp(key!*,'ignore)
		and null eqcar(x,'quote)
	       then progn((if x then dfprint x else nil),
			  if null flagp(key!*,'eval) then return nil);
      if !*output and ifl!* and !*echo and null !*lessspace
	then terpri();
      result := errorset!*(x,t);
      if errorp result or erfg!*
	then progn(programl!* := list(mode,x),return 'err2)
       else if !*defn then return nil;
      if null !*reduce4
	then if null(mode eq 'symbolic) then x := getsetvars x else nil
       else progn(result := car result,
		  (if null result then result := mkobject(nil,'noval)),
		  mode := type result,
		  result := value result);
      add2resultbuf((if null !*reduce4 then car result else result),
		    mode);
      if null !*output then return nil
       else if null(semic!* eq '!$)
	then if !*reduce4 then (begin
		   terpri();
		   if mode eq 'noval then return nil
		    else if !*debug then prin2t "Value:";
		   rapply1('print,list list(mode,result))
		 end)
       else if mode eq 'symbolic
	      then if null car result and null(!*mode eq 'symbolic)
		     then nil
	      else begin
		  terpri();
		  result:=
		       errorset!*(list('print,mkquote car result),t)
		    end
       else if car result
	then result := errorset!*(list('assgnpri,mkquote car result,
				       (if x then 'list . x else nil),
				       mkquote 'only),
				  t);
      if null !*reduce4
	then return if errorp result then 'err3 else nil
       else if null(!*mode eq 'noval)
	then progn(terpri(), prin2 "of type: ", print mode);
      return nil
   end;

endpatch;

% Roots declarations.

% fluid '(rootacc!#!# rootacc!#!# !*noeqns);

% patch roots;

% 10 Feb 00.

% Commented out since now solved another way (7 Apr 01).

% symbolic procedure root_val x;
%    roots x
%      where rootacc!#!#=p, iniprec!#=p where p=precision 0, !*msg=nil,
%            !*noeqns=t;

% endpatch;

% Scope declarations.

global '(kvarlst prevlst varlst!*);

patch scope;

% 29 Aug 00.

symbolic procedure maxtype type;
   if atom type then type
    else if pairp cdr type then cadr type else car type;

%  8 Nov 00.

symbolic procedure prepmultmat(preprefixlist);
begin scalar tlcm,var,varexp,kvl,kfound,pvl,pfound,tel,ratval,ratlst,
                                                      newvarlst,hvarlst;
 hvarlst:= nil;
 while not null (varlst!*) do
 <<var := car varlst!*; varlst!* := cdr varlst!*;
   if flagp(var,'ratexp)
    then
     <<tlcm:=1;
       remflag(list var,'ratexp);
       foreach elem in get(var,'varlst!*) do
        if pairp cdr elem then tlcm := lcm2(tlcm,cddr elem);
       varexp:=fnewsym();
       tel:=(varexp.(if tlcm = 2
                 then list('sqrt,var)
                 else list('expt,var,
                         if onep cdr(tel:=simpquot list(1,tlcm)) then
                            car tel
                         else
                            list('quotient,car tel,cdr tel))));
       if assoc(var,kvarlst)
        then
         <<kvl:=kfound:=nil;
           while kvarlst and not(kfound) do
            if caar(kvarlst) eq var
             then
              << kvl:=tel.kvl; kfound:=t;
                 pvl:=pfound:=nil; prevlst:=reverse(prevlst);
                 while prevlst and not(pfound) do
                  if cdar(prevlst) eq var
                   then << pvl:=cons(caar prevlst,varexp).pvl;
                           pfound:=t
                        >>
                   else << pvl:=car(prevlst).pvl;
                           prevlst:=cdr(prevlst)
                        >>;
                 if pvl then
                  if prevlst then prevlst:=append(reverse prevlst,pvl)
                             else prevlst:=pvl
              >>
             else
              << kvl:=car(kvarlst).kvl; kvarlst:=cdr kvarlst>>;
           if kvl then
             if kvarlst then kvarlst:=append(reverse kvl,kvarlst)
                      else kvarlst:=reverse kvl
         >>
        else preprefixlist:=tel.preprefixlist;
       ratlst:=newvarlst:=nil;
       foreach elem in get(var,'varlst!*) do
        if pairp cdr elem
         then 
          << ratval:=divide((tlcm * cadr elem)/(cddr elem),tlcm);
             ratlst:=cons(car elem,cdr ratval).ratlst;
             if car(ratval)>0
              then newvarlst:=cons(car elem,car ratval).newvarlst
          >>
         else newvarlst:=elem.newvarlst;
       if ratlst
        then << put(varexp,'varlst!*,reverse ratlst);
                hvarlst:=varexp.hvarlst
             >>;
       if newvarlst
        then << put(var,'varlst!*,reverse newvarlst);
                hvarlst:=var.hvarlst
             >>
        else remprop(var,'varlst!*)
     >>
    else hvarlst:=var.hvarlst
 >>;
 varlst!* := hvarlst;
 return preprefixlist
   end;

endpatch;

% Solve declarations.

fluid '(!*multiplicities vars!*);

global '(multiplicities!*);

patch solve;

%  9 Jan 01.

symbolic procedure !*solvelist2solveeqlist u;
   begin scalar x,y,z;
      u := for each j in u collect solveorder j;
      for each j in u do
         <<if caddr j=0 then rerror(solve,2,"zero multiplicity")
            else if null cadr j
             then  x := for each k in car j collect
                                               list('equal,!*q2a k,0)
            else x := for each k in pair(cadr j,car j)
                          collect list('equal,car k,!*q2a cdr k);
           if length vars!* > 1 then x := 'list . x else x := car x;
           z := (caddr j . x) . z>>;
      z := sort(z,function ordp);
      x := nil;
      if !*multiplicities
	 then <<for each k in z do for i := 1:car k do x := cdr k . x;
		multiplicities!* := nil>>
       else <<for each k in z do << x := cdr k . x; y := car k . y>>;
	      multiplicities!* := 'list . reversip y>>;
      return 'list . reversip x
   end;

% 9 Jan 01, 15 Jun 01.

symbolic procedure solveorder u;
   begin scalar v,w,x,y,z;
      v := vars!*;
      x := cadr u;
      if length x<length v then v := setdiff(v,setdiff(v,x));
      if null x or x=v then return u;
      y := car u;
      while x do <<z := (car x . car y) . z; x := cdr x; y := cdr y>>;
      w := v;
  a:  if null w then return reversip x . v . cddr u
       else if null(y := depassoc(car w,z)) then return u
       else x := cdr y . x;
      w := cdr w;
      go to a
   end;

symbolic procedure depassoc(u,v);
   if null v then nil
    else if u = caar v then car v
    else if depends(caar v,u) then nil
    else depassoc(u,cdr v);

% 2 Feb 01.

symbolic procedure check!-solns(z,ex,var);
   begin scalar x;
      if errorp (x :=
	    errorset2 {'check!-solns1,mkquote z,mkquote ex,mkquote var})
	then return
	   check!-solns1(z,(numr simp!* prepf ex where !*reduced=t),var)
       else return car x
   end;

symbolic procedure check!-solns1(z,ex,var);
   begin scalar x,y,fv,sx,vs;
      fv := freevarl(ex,var);
      for each z1 in z do
        fv := union(fv,union(freevarl(numr caar z1,var),
                             freevarl(denr caar z1,var)));
      fv := delete('i,fv);
      if fv then for each v in fv do
        if not flagp(v,'constant) then
           vs := (v . list('quotient,1+random 999,1000)) . vs;
      sx := if vs then numr subf(ex,vs) else ex;
      while z do
         if null cadar z then <<z := nil; x := 'unsolved>>
         else if
           <<y := numr subf(ex,list(caadar z . mk!*sq caaar z));
             null y
             or fv and null(y := numr subf(sx,list(caadar z .
                   mk!*sq subsq(caaar z,vs))))
             or null numvalue y>>
           then <<x := car z . x; z := cdr z>>
          else z := cdr z;
      return if null x then 'unsolved else x
   end;

%  7 Apr 01.

symbolic procedure solvequadratic(a2,a1,a0);
   if !*rounded and numcoef a0 and numcoef a1 and numcoef a2
      then for each z in cdr root_val list mkpolyexp2(a2,a1,a0)
	 collect simp!* (if eqcar(z,'equal) then caddr z
			  else errach {"Quadratic confusion",z})
    else begin scalar d;
      d := sqrtq subtrsq(quotsqf(exptsq(a1,2),4),multsq(a2,a0));
      a1 := quotsqf(negsq a1,2);
      return list(subs2!* quotsq(addsq(a1,d),a2),
                  subs2!* quotsq(subtrsq(a1,d),a2))
   end;

endpatch;

% Sum declarations.

fluid '(sum_last_attempt_rules!* !*zeilberg);

patch sum;

% 28 Jul 00.

symbolic procedure freeof!-df(u, v);
   if atom u then t
    else if car(u) eq 'df
     then freeof!-df(cadr u, v) and not smember(v,cddr u)
    else freeof!-dfl(cdr u, v);

symbolic procedure freeof!-dfl(u, v);
   if null u then t else freeof!-df(car u,v) and freeof!-dfl(cdr u,v);

symbolic procedure simp!-sum u;
   begin scalar y;
      y := cdr u;
      u := car u;
      if not atom y and not freeof!-df(u, car y) then
	if atom y
	       then return !*p2f(car fkern(list('sum,u)) .* 1) ./ 1
	 else return sum!-df(u, y);
      u := simp!* u;
      return if null numr u then u
	      else if atom y
               then !*p2f(car fkern(list('sum,prepsq u)) .* 1) ./ 1
	      else if !*zeilberg then gosper!*(mk!*sq u,y)
	      else simp!-sum0(u,y)
   end;

symbolic procedure sum!-subst(u,x,a);
   if u = x then a
    else if atom u then u
    else sum!-subst(car u, x,a) . sum!-subst(cdr u,x,a);

symbolic procedure sum!-df(u,y);
   begin scalar w,z,upper,lower,dif;
      if length(y) = 3 then  <<
	 lower := cadr y;
	 upper := caddr y;
	 dif := addsq(simp!* upper, negsq simp!* lower);
	 if denr dif = 1 then
	    if null numr dif
	       then return simp!* sum!-subst(u, car y, upper)
	       else if fixp numr dif then dif := numr dif
	       else dif := nil
	   else dif := nil;
	 if dif and dif <= 0 then return nil ./ 1 >>;
     if null dif then <<
	  z := 'sum . (u . y);
	  let sum_last_attempt_rules!*;
	  w:= opmtch z;
	  rule!-list (list sum_last_attempt_rules!*,nil);
	  return if w then simp w else mksq(z,1)>>;
     z := nil ./ 1;
  a: if dif < 0 then return z;
     z := addsq(z,simp!* sum!-subst(u, car y, list('plus,lower,dif)));
     dif := dif - 1;
     go to a
   end;

% 20 Nov 00.

symbolic procedure termlst(u,v,klst);
   begin scalar x,kern,lst;
      if null u then return nil
       else if null klst or domainp u
	then return list multsq(v,!*f2q u);
      kern := car klst;
      klst := cdr klst;
      x := setkorder list kern;
      u := reorder u;
      v := reorder(numr v) ./ reorder(denr v);
      while not domainp u and mvar u eq kern do <<
         lst := nconc(termlst(lc u, multsq(!*p2q lpow u, v),klst),lst);
         u := red u>>;
      if u then lst := nconc(termlst(u,v,klst),lst);
      setkorder x;
      return lst
   end;

endpatch;

% Taylor declarations.

fluid '(!*taylorautocombine);

patch taylor;

% 1 Jun 01.

symbolic procedure simptaylor u;
  if remainder(length u,3) neq 1
    then Taylor!-error('wrong!-no!-args,'taylor)
   else if null subfg!* then mksq('taylor . u,1)
   else begin scalar !*precise,arglist,degree,f,ll,result,var,var0;
     if !*taylorautocombine and not ('taysimpsq memq mul!*)
       then mul!* := aconc!*(mul!*,'taysimpsq);
     f := simp!* car u;
     u := revlis cdr u;
     arglist := u;
     while not null arglist do
       << var := car arglist;
          var := if eqcar(var,'list) then cdr var else {var};
          for each el in var collect begin
            el := simp!* el;
            if kernp el then return mvar numr el
             else typerr(prepsq el,'kernel)
            end;
          var0 := cadr arglist;
          degree := caddr arglist;
          if not fixp degree
            then typerr(degree,"order of Taylor expansion");
          arglist := cdddr arglist;
          ll := {var,var0,degree,degree + 1} . ll>>;
     result := taylorexpand(f,reversip ll);
     return if smember('Taylor!*,result) then result
             else mksq('taylor . prepsq f . u,1)
   end;

endpatch;

endmodule;

end;