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

% Author: Anthony C. Hearn.

% Copyright (c) 2004, 2005, 2006, 2007, 2008 Anthony C. Hearn.  All Rights Reserved.

global '(patch!-date!* patch!-url!-list!*);

patch!-date!* := "2-Aug-2008";

patch!-url!-list!* :=
   '("http://reduce-algebra.com/support/patches/patches.fsl");


% Bugs fixed by these patches.

% 26 Jun 04.  With rounded arithmetic, solving some linear equation
%             problems could lead to a catastrophic error.

%  8 Jul 04.  Some non-zero integrals (e.g., int(e^(a^(1/3)*x)*sin x,x))
%             returned zero.

%  5 Aug 04.  Using RLFI with latex on could lead to invalid operator errors.

%  2 Sep 04.  In rare circumstances, floating point conversion could give
%             an extraneous error.

%  6 Sep 04.  With rational on, some non-zero factorizations could produce
%             a zero coefficient (e.g., on rational;
%             factorize(r^((1/4*n^2 - 1/4*n + 1)/(n - 1)));).

% 28 Sep 04.  Some integrals would not return a closed form solution
%             with algint on that would with algint off
%             (e.g., int(sqrt(x-1)/(sqrt x*(x-1)),x)).

% 10 Dec 04.  With dfprin on, some products and sums printed incorrectly.

% 31 Jan 05.  Some integrals involving square roots could run forever.

% 12 Feb 05.  SOLVE could produce a spurious recursive loop (e.g.,
%             solve((4*e^(y^3/3)*cte+2x^2+y^3+3)/e^(y^3/3),y)).

% 20 Apr 05. int(e^(-a^(1/4)*(-1)^(1/4)*x),x); terminated with an error.

%  2 May 05. int(e^(-a^(1/4)*(-1)^(1/4)*i*x)*b+(1/4)*e^(-a^(1/4)
%                       *(-1)^(1/4)*i*x)*x,x); terminated with an error.

% 22 May 05. Some integrals, e.g., int(e^((3sqrt 5+1)*x)*(sqrt 5+1)
%             +e^((3sqrt 5-1)*x)*(sqrt 5+1),x), never completed.

% 30 May 05. SOLVE could produce a spurious "Zero divisor" error
%            (e.g., solve({log tan(y/2),y+1/x},{x,y})).

%  4 Oct 05. DEG did not work with rational coefficients (e.g.,
%            deg(x**3/a-x/5+1/4,x)).

%  5 Oct 05. Some SOLVE calculations could give a spurious "Zero Divisor"
%            error (e.g., ex0:= sqrt(a^2-y^2); solve((-log(( - x + a + y)
%               /ex0) + log((x + a + y)/ex0) + x - (a^2 - y^2)/ex0),y);

% 16 Nov 05. System errors could occur with rounded and combineexpt on.
%            (E.g., on rounded,combineexpt; 0.183*e^x*t^4.39;).

% 22 Nov 05. Some definite integrals with variables other than x could
%            give a wrong answer, e.g., int(e^(-y),y,0,x).

%  9 Dec 05. With combineexpt on, expressions could be dropped (e.g.,
%            on combineexpt; 4*e^(-3*h/2) - 3*h*e^(-h) + 2*e^(-h)).

%  4 Feb 06. Setcrackflags() was not set in crack, but needed to be.

% 20 Feb 06. The rule for df(Jacobidn(~u,~m),~u) was wrong.

% 21 Feb 06. Evaluating some integrals could suppress the printing of
%            the results.

% 22 Feb 06. Some sub evaluations could include superfluous terms like
%            x = (x^(1/7))^7.

% 23 May 06. Derivatives and integrals of matrices were not computed.

% 18 Aug 06. After nospur, some traces were still evaluated.

% 29 Sep 06. With dfprint on, derivatives of integrals would print in a
%            truncated form.

% 11 Jan 07. With rounded arithmetic and factor on, a non-numeric
%            argument error could occur.

% 18 Oct 07. With algint on, some calculations could give a catastrophic
%            error in mknewsqrt.

% 19 Oct 07. With algint on, some calculations could give a catastrophic
%            error in gcdf.

% 26 Dec 07. In rare cases, factorization could lead to an error
%            "FINAL DELTA IS -VE IN DISTRIBUTE!.LC."

% 31 Jan 08. With precise on, some roots handled "abs" inconsistently
%            (e.g., (x^3)^(1/3), sqrt(1/x^2)).

% 27 Jul 08. When solving quartics over complex integers, a non-numeric
%            error could occur.

%  2 Aug 08. Non-numeric exponents could cause an error in some Groebner
%            calculations.

% Alg declarations.

fluid '(sublist!*);

patch alg;

% 16 Nov 05, 9 Dec 05.

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 null domainp u and null red u
	  and (z2 := kernels u) and null cdr z2
        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 if domainp numr z2 then y := 1
	     else <<y := lcoeffgcd cdr comfac numr z2;
		    if not fixp y then 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;

% 22 Feb 06.

symbolic procedure subeval0 u;
   begin scalar x,y,z,ns;
   while cdr u do <<if not eqcar(car u,'equal) then x := car u . x
		     else if not(cadar u = (y := reval caddar u))
		      then x := {caar u,cadar u,y} . x;
		    u := cdr u>>;
   if null x then return car u else u := nconc(reversip x,u);
      if u member sublist!* then return mk!*sq !*p2q mksp('sub . u,1)
       else sublist!* := u . sublist!*;
      if null(u and cdr u)
       then rederr "SUB requires at least 2 arguments";
      (while cdr u do
	  <<x := reval car u;
	    if getrtype x eq 'list then u := append(cdr x,cdr u)
	     else <<if not eqexpr x then errpri2(car u,t);
		    y := cadr x;
		    if null getrtype y then y := !*a2kwoweight y;
		    if getrtype caddr x then ns := (y . caddr x) . ns
		     else z := (y . caddr x) . z;
		    u := cdr u>>>>) where !*evallhseqp=nil;
      x := aeval car u;
      return subeval1(append(ns,z),x)
   end;

symbolic procedure subsubf(l,expn);
   begin scalar x,y;
      for each j in l do if car j neq (y := prepsq!* simp!* cdr j)
			   then x := (car j . y) . x;
      l := reversip x;
      if null l then return expn;
      y := nil;
      for each j in cddr expn do
         if (x := assoc(j,l)) then <<y := x . y; l := delete(x,l)>>;
      expn := sublis(l,car expn)
                 . for each j in cdr expn collect subsublis(l,j);
      if null y then return expn;
      expn := aconc!*(for each j in reversip!* y
                     collect list('equal,car j,aeval cdr j),expn);
      return if l then subeval expn
	      else mk!*sq !*p2q mksp('sub . expn,1)
   end;

% 23 May 06.

symbolic procedure reval1(u,v);
   (begin scalar x,y;
      if null u then return nil
       else if stringp u then return u
       else if fixp u
        then return if flagp(dmode!*,'convert) then reval2(u,v) else u
       else if atom u
	then if null subfg!* then return u
	      else if idp u and (x := get(u,'avalue))
               then if u memq varstack!* then recursiveerror u
                     else <<varstack!* := u . varstack!*;
			    return if y := get(car x,'evfn)
				     then apply2(y,u,v)
				    else reval1(cadr x,v)>>
              else nil
       else if not idp car u
        then errpri2(u,t)
       else if car u eq '!*sq
	then return if caddr u and null !*resimp
                      then if null v then u else prepsqxx cadr u
                     else reval2(u,v)
       else if flagp(car u,'remember) then return rmmbreval(u,v)
       else if flagp(car u,'opfn) then return reval1(opfneval u,v)
       else if x := get(car u,'psopfn)
	 then <<u := apply1(x,cdr u);
		if x := get(x,'cleanupfn) then u := apply2(x,u,v);
		return u>>
       else if arrayp car u then return reval1(getelv u,v);
       return if x := getrtype u then
		  if y := get(x,'evfn) then apply2(y,u,v)
		    else rerror(alg,101,
				list("Missing evaluation for type",x))
	       else if not atom u
		     and not atom cdr u
		     and (y := getrtype cadr u)
		     and null(y eq 'list and cddr u)
		     and (x := get(y,'aggregatefn))
		     and (not(x eq 'matrixmap) or flagp(car u,'matmapfn))
		     and not flagp(car u,'boolean)
		     and not !*listargs and not flagp(car u,'listargp)
                then apply2(x,u,v)
               else reval2(u,v)
   end) where varstack!* := varstack!*;

symbolic procedure getrtype2 u;
   begin scalar x;
      return if (x := get(car u,'rtype)) and (x := get(x,'rtypefn))
	       then apply1(x,cdr u)
	      else if x := get(car u,'rtypefn) then apply1(x,cdr u)
	      else if flagp(car u,'matmapfn) and cdr u
		 and getrtype cadr u eq 'matrix
	       then 'matrix
	      else nil
   end;

% 31 Jan 08.

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
               and numberp numr m and evenp numr m
	       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;

symbolic procedure radf(u,n);
   begin scalar ipart,rpart,x,y,z,!*gcd,!*mcd;
      if null u then return list u;
      !*gcd := !*mcd := t;
      ipart := 1;
      z := 1;
      while not domainp u do
     <<y := comfac u;
       if car y
         then <<x := divide(pdeg car y,n);
            if car x neq 0
	      then ipart := multf(
		   if evenp car x
		      then !*p2f(mvar u .** car x)
                    else check!-radf!-sign(!*p2f(mvar u .** pdeg car y),
                                           !*p2f(
                            (if !*precise and evenp n then {'abs,mvar u}
                              else mvar u)
                                                  .** car x),
                                           n),
                    ipart);
            if cdr x neq 0
	      then rpart := mkexpt(sfchk mvar u,cdr x) . rpart>>;
       x := quotf(u,comfac!-to!-poly y);
       u := cdr y;
       if !*reduced and minusf x
         then <<x := negf x; u := negf u>>;
       if flagp(dmode!*,'field) then
          <<y := lnc x;
		if y neq 1 then <<x := quotf(x,y); z := multd(y,z)>>>>;
       if x neq 1
         then <<x := radf1(sqfrf x,n);
		y := car x;
		if y neq 1 then
		   <<
		     ipart := multf(y,ipart)>>;
		rpart := append(rpart,cdr x)>>>>;
      if u neq 1
    then <<x := radd(u,n);
           ipart := multf(car x,ipart);
           rpart := append(cdr x,rpart)>>;
      if z neq 1
    then if !*numval
        and (y := domainvalchk('expt,
                       list(!*f2q z,!*f2q !:recip n)))
           then ipart := multd(!*q2f y,ipart)
	  else rpart := prepf z . rpart;
      return ipart . rpart
   end;

endpatch;

% Algint declarations.

fluid '(!*noacn !*tra !*trmin listofallsqrts listofnewsqrts);

patch algint;

% 19 Oct 07.

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 printc "Inner integration failed";
      deriv:=nil ./ 1;
      rischpart:=deriv >>
    else
      if atom car rischpart
        then <<
      if !*tra or !*trmin then
        printc "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 <<
      printc "Inner working yields";
          printsq rischpart;
      printc "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;
  deriv := resimp deriv;
  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;

endpatch;

patch arith;

% 2 Sep 04.

symbolic procedure read!:num(n);
   if fixp n then make!:ibf(n, 0)
    else if not(numberp n or stringp n) then bflerrmsg 'read!:num
    else begin integer j,m,sign;  scalar ch,u,v,l,appear!.,appear!/;
          j := m := 0;
          sign := 1;
          u := v := appear!. := appear!/ := nil;
          l := explode n;
    loop: ch := car l;
          if digit ch then << u := ch . u; j := j + 1 >>
           else if ch eq '!. then << appear!. := t; j := 0 >>
           else if ch eq '!/ then << appear!/ := t; v := u; u := nil >>
           else if ch eq '!- then sign := -1
	   else if ch memq '(!E !D !B !e !d !b) then go to jump;
           if l := cdr l then goto loop else goto make;
    jump: while l := cdr l do
	    <<if digit(ch := car l) or ch eq '!-
		 then v := ch . v >>;
	  l := reverse v;
	  if car l eq '!- then m := - compress cdr l
			  else m:= compress l;
    make: u := reverse u;
          v := reverse v;
          if appear!/ then
            return conv!:r2bf(make!:ratnum(sign*compress v,compress u),
                              if !:bprec!: then !:bprec!: else 170);
          if appear!. then j := - j else j := 0;
          if sign = 1 then u := compress u else u := - compress u;
          return round!:mt (decimal2internal (u, j + m), !:bprec!:)
                   where !:bprec!: := if !:bprec!: then !:bprec!:
                                       else msd!: abs u
    end;

endpatch;

patch crack;

setcrackflags();

endpatch;

% Defint declarations.

symbolic smacro procedure listsq(u);
   for each uu in u collect simp!* uu;

patch defint;

symbolic procedure new_meijer(u);
begin scalar f,y,mellin,new_mellin,m,n,p,q,old_num,old_denom,temp,a1,
b1,a2,b2,alpha,num,denom,n1,temp1,temp2,coeff,v,var,new_var,new_y,
new_v,k;
f := prepsq simp car u;
y := caddr u;
mellin := bastab(car f,cddr f);
temp := car cddddr mellin;
var := cadr f;
if not idp VAR then RETURN error(99,'FAIL);
temp := reval algebraic(sub(x=var,temp));
mellin := {car mellin,cadr mellin,caddr mellin,cadddr mellin,temp};
temp := reduce_var(cadr u,mellin,var);
alpha := simp!* car temp;
new_mellin := cdr temp;
if car cddddr new_mellin neq car cddddr mellin then
        << k := car cddddr mellin;
           y := reval algebraic(sub(var=y,k));
           new_y := simp y>>
else
<< new_var := car cddddr new_mellin;
   new_y := simp reval algebraic(sub(x=y,new_var))>>;
n1 := addsq(alpha,'(1 . 1));
temp1 := {'expt,y,prepsq n1};
temp2  := cadddr new_mellin;
coeff := simp!* reval algebraic(temp1*temp2);
m := caar new_mellin;
n := cadar new_mellin;
p := caddar new_mellin;
q := car cdddar new_mellin;
old_num := cadr new_mellin;
old_denom := caddr new_mellin;
for i:=1 :n do
<< if old_num = nil then a1 := append(a1,{simp!* old_num })
   else <<  a1 := append(a1,{simp!* car old_num});
            old_num := cdr old_num>>;
>>;
for j:=1 :m do
<< if old_denom = nil then b1 := append(b1,{simp!*  old_denom })
   else <<  b1 := append(b1,{simp!* car old_denom});
            old_denom := cdr old_denom>>;
>>;
a2 := listsq old_num;
b2 := listsq old_denom;
if a1 = nil and a2 = nil then
    num := list({negsq alpha})
else if a2 = nil then num := list(append(a1,{negsq alpha}))
else
<< num := append(a1,{negsq alpha}); num := append({num},a2)>>;
if b1 = nil and b2 = nil then
    denom := list({subtrsq(negsq alpha,'(1 . 1))})
else if b2 = nil then
    denom := list(b1,subtrsq(negsq alpha,'(1 . 1)))
else
<< denom := list(b1,subtrsq(negsq alpha,'(1 . 1)));
   denom := append(denom,b2)>>;
v := gfmsq(num,denom,new_y);
if v = 'fail then return simp 'fail
else v := prepsq subsq(v,list(prepsq new_y . y));
if eqcar(v,'meijerg) then new_v := v else new_v := simp v;
return multsq(new_v,coeff);
end;

endpatch;

% Dipoly declarations.

smacro procedure dipfmon(a,e); e . a . dipzero;

smacro procedure vevzero!? u;
 null u or(car u=0 and vevzero!?1 cdr u);

patch dipoly;

% 2 Aug 08.

symbolic procedure a2dip u;
if atom u then a2dipatom u
   else if not atom car u or not idp car u
    then typerr(car u,"dipoly operator")
   else if car u='expt
           then if vevzero!? car a2dip cadr u and
                      vevzero!? car a2dip caddr u
                   then dipfmon(simp!* u,evzero())
                 else dipfnpow(a2dip cadr u,caddr u)
   else (if x then apply(x,list for each y in cdr u collect a2dip y)
        else a2dipatom u)
         where x=get(car u,'dipfn);

endpatch;

patch hephys;

symbolic procedure nospur u; <<rmsubs(); !*nospurp := t; flag(u,'nospur)>>;

endpatch;

% Factor declarations.

fluid '(bad!-case factor!-level factor!-trace!-list);

symbolic smacro procedure get!-image!-content s; cadr cdr cddr s;

symbolic smacro procedure get!-f!-numvec s; cadr cddr cdddr s;

smacro procedure printstr l; << prin2!* l; terpri!*(nil) >>;

symbolic smacro procedure factor!-trace action;
   begin scalar stream;
      if !*trallfac or (!*trfac and factor!-level = 1)
	then stream := nil . nil
       else stream := assoc(factor!-level,factor!-trace!-list);
      if stream then <<stream := wrs cdr stream; action; wrs stream>>
   end;
 
symbolic smacro procedure get!-image!-set s; car s;

patch factor;

% 26 Dec 07.

symbolic procedure distribute!.lc(r,im!.factors,s,v);
  (lambda factor!-level;
  begin scalar k,delta,div!.count,q,uf,i,d,max!.mult,f,numvec,
               dvec,wvec,dtwid,w;
    delta:=get!-image!-content s;
    dist!.lc!.msg1(delta,im!.factors,r,s,v);
    v:=cdr v;
    k:=length v;
    numvec:=get!-f!-numvec s;
    dvec:=mkvect r;
    wvec:=mkvect r;
    for j:=1:r do <<
      putv(dvec,j,1);
      putv(wvec,j,delta*lc getv(im!.factors,j)) >>;
    v:=reverse v;
    for j:=k step -1 until 1 do
    <<
      f:=caar v;
      max!.mult:=cdar v;
      v:=cdr v;
      d:=getv(numvec,j);
      i:=1; % we trial divide d into lc of each image
      div!.count:=0;
      factor!-trace <<
        prin2!* "f("; prin2!* j; prin2!* ")= "; printsf f;
        prin2!* "There are "; prin2!* max!.mult;
        printstr " of these in the leading coefficient.";
        prin2!* "The absolute value of the image of f("; prin2!* j;
        prin2!* ")= "; printstr d >>;
      while ilessp(div!.count,max!.mult)
        and not igreaterp(i,r) do
      << q:=divide(getv(wvec,i),d);
        factor!-trace <<
          prin2!* "  Trial divide into ";
          prin2!* getv(wvec,i); printstr " :" >>;
        while (zerop cdr q) and ilessp(div!.count,max!.mult) do
        << putv(dvec,i,multf(getv(dvec,i),f));
          factor!-trace <<
            prin2!* "    It goes so an f("; prin2!* j;
            prin2!* ") belongs in ";
            printsf getv(im!.factors,i);
            printstr "  Try again..." >>;
          div!.count:=iadd1 div!.count;
          putv(wvec,i,car q);
          q:=divide(car q,d);
        >>;
        i:=iadd1 i;
        factor!-trace
           <<printstr "    no good so try another factor ..." >>>>;
      if ilessp(div!.count,max!.mult) then
        <<bad!-case:=t; div!.count := max!.mult>>
    >>;
    if bad!-case then return;
    dist!.lc!.msg2(dvec,im!.factors,r);
    if onep delta then
    << for j:=1:r do <<
         w:=lc getv(im!.factors,j) /
          evaluate!-in!-order(getv(dvec,j),get!-image!-set s);
         if w<0 then begin
           scalar oldpoly;
           delta:= -delta;
           oldpoly:=getv(im!.factors,j);
           putv(im!.factors,j,negf oldpoly);
           multiply!-alphas(-1,oldpoly,getv(im!.factors,j));
         end;
         putv(dvec,j,multf(abs w,getv(dvec,j))) >>;
      dist!.lc!.msg3(dvec,im!.factors,r);
      return (delta . dvec)
    >>;
    dist!.lc!.msg4 delta;
    for j:=1:r do
    << dtwid:=evaluate!-in!-order(getv(dvec,j),get!-image!-set s);
       uf:=getv(im!.factors,j);
       d:=gcddd(lc uf,dtwid);
       putv(dvec,j,multf(lc uf/d,getv(dvec,j)));
       putv(im!.factors,j,multf(dtwid/d,uf));
       multiply!-alphas!-recip(dtwid/d,uf,getv(im!.factors,j));
       delta:=delta/(dtwid/d)
    >>;
    if delta<=0 then <<
      factor!-trace <<
	 prin2!* "FINAL DELTA IS -VE IN DISTRIBUTE!.LC";
	 printstr delta >>;
      delta := 1 >>;
    factor!-trace <<
      printstr "     Finally we have:";
      for j:=1:r do <<
	prinsf getv(im!.factors,j);
	prin2!* " with l.c. ";
	printsf getv(dvec,j) >> >>;
    return (delta . dvec)
  end) (factor!-level * 10);

endpatch;

% Int declarations.

fluid '(!*purerisch !*trdint gaussiani indexlist intvar sqrt!-places!-alist
        loglist !*intflag!* listofnewsqrts listofallsqrts sqrt!-intvar
        basic!-listofallsqrts basic!-listofnewsqrts !*precise dmode!*
        !*exp !*gcd !*keepsqrts !*limitedfactors !*mcd !*rationalize
	!*structure !*uncached kord!* !*galois !*pvar listofnewsqrts);

smacro procedure argof u; cadr u;

patch int;

% 8 Jul 04, 31 Jan 05, 20 Apr 05, 2 May 05.

symbolic procedure df2q p;
   begin scalar n,d,w,x,y,z;
      if null p then return nil ./ 1;
      d:=denr lc p;
      w:=red p;
      while w do
       <<d := multf(d,quotf(denr lc w,gcdf(d,denr lc w)));
	 w := red w>>;
      while p do begin
	   w := sqrt2top lc p;
	   x := multf(xl2f(lpow p,zlist,indexlist),multf(numr w,d));
	   if null x then return (p := red p);
	   y := denr w;
	   z := quotf(x,y);
	   if null z
	     then <<z := rationalizesq(x ./ y);
		    if denr z neq 1
		      then <<d := multf(denr z,d); n := multf(denr z,n)>>;
		    z := numr z>>;
	   n := addf(n,z);
	   p := red p
	end;
      return tidy!-powersq (n ./ d)
   end;

% 8 Jul 04, 22 May 05.

symbolic procedure tidy!-powersq x;
   begin scalar expts,!*precise,!*keepsqrts;
    !*keepsqrts := t;
    x := subs2q x;
    expts := find!-expts(numr x,find!-expts(denr x,nil));
    if null expts then return x;
    x := subsq(x,for each v in expts collect
		    (car v . list('expt,cadr v,cddr v)));
    x := subsq(x,for each v in expts collect
		    (cadr v
			. list('expt,car v,list('quotient,1,cddr v))));
    return x
  end;

symbolic procedure find!-expts(ff,l);
   begin scalar w;
      if domainp ff then return l;
      l := find!-expts(lc ff,find!-expts(red ff, l));
      ff := mvar ff;
      if eqcar(ff,'sqrt)
	then ff := list('expt, cadr ff,'(quotient 1 2))
       else if eqcar(ff,'expt) and eqcar(caddr ff,'quotient)
	  and numberp caddr caddr ff
	then <<w := assoc(cadr ff,l);
	       if null w
		 then <<w := cadr ff . gensym() . 1; l := w . l >>;
	       rplacd(cdr w,lcm(cddr w,caddr caddr ff))>>;
      return l
   end;

% 28 Sep 04.

symbolic procedure look_for_quad(integrand, var, zz);
 begin
  if (car zz = 'sqrt and listp cadr zz and caadr zz = 'plus) or
     (car zz = 'expt and listp cadr zz and caadr zz = 'plus and
      listp caddr zz and car caddr zz = 'quotient
         and fixp caddr caddr zz)
   then <<
    zz := simp cadr zz;
    if (cdr zz = 1) then <<
        zz := cdr coeff1(prepsq zz, var, nil);
	if length zz = 2 then return begin
          scalar a, b;
          scalar nvar, res, ss;
          a := car zz; b := cadr zz;
          if (depends(a,var) or depends(b,var)) then return nil;
          nvar := gensym();
          if !*trint then <<
                prin2 "Linear shift suggested ";
                prin2 a; prin2 " "; prin2 b; terpri();
          >>;
	  integrand := subsq(integrand,
                             list(var . list('quotient,
                                             list('difference,
                                                  list('expt,nvar,2),a),
                                                  b)));
	  integrand := multsq(integrand,
                              simp list('quotient,list('times,nvar,2),
                                        b));
          if !*trint then <<
                prin2 "Integrand is transformed by substitution to ";
                printsq integrand;
                prin2 "using substitution "; prin2 var; prin2 " -> ";
                printsq simp list('quotient,
                                 list('difference,list('expt,nvar,2),a),
                                 b);
           >>;
           res := integratesq(integrand, nvar, nil, nil);
           ss := list(nvar . list('sqrt,list('plus,list('times,var,b),
                                  a)));
           res := subsq(car res, ss) .
                  subsq(multsq(cdr res, simp list('quotient,b,
                                                  list('times,nvar,2))), ss);
           return res;
        end
	else if length zz = 3 then return begin
          scalar a, b, c;
          a := car zz; b := cadr zz; c:= caddr zz;
          if (depends(a,var) or depends(b,var) or depends(c,var)) then
                return nil;
	  a := simp list('difference, a,
                         list('times,b,b,
                         list('quotient,1,list('times,4,c))));
	  if null numr a then return nil;
          b := simp list('quotient, b, list('times, 2, c));
          c := simp c;
          return
           if minusf numr c then <<
           if minusf numr a then  begin
                            scalar !*hyperbolic;
                            !*hyperbolic := t;
                            return
                                look_for_invhyp(integrand,nil,var,a,b,c)
                        end
            else                look_for_asin(integrand,var,a,b,c)>>
          else <<
            if minusf numr a then look_for_invhyp(integrand,t,var,a,b,c)
          else                  look_for_invhyp(integrand,nil,var,a,b,c)
          >>
        end
	else if length zz = 5 then return begin
          scalar a, b, c, d, e, nn, dd, mm;
          a := car zz; b := cadr zz; c:= caddr zz;
          d := cadddr zz; e := car cddddr zz;
          if not(b = 0) or not(d = 0) then return nil;
          if (depends(a,var) or depends(c,var)) or depends(e,var) then
                return nil;
          nn := numr integrand;  dd := denr integrand;
          if denr(mm :=quotsq(nn ./ 1, !*kk2q var)) = 1 and
             even_power(numr mm, var) and even_power(dd, var) then <<
              return sqrt_substitute(numr mm, dd, var);
          >>;
          if denr(mm :=quotsq(dd ./ 1, !*kk2q var)) = 1 and
             even_power(nn, var) and even_power(numr mm, var) then <<
              return sqrt_substitute(nn, multf(dd,!*kk2f var), var);
          >>;
          return nil;
        end;
  >>>>;
  return nil
end;

% 21 Feb 06.

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 printc {"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 (argof 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,badbit,denexp,erfg,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 := zvars(getvariables x,zv,variable,t);
	   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>>;
	  erfg := erfg!*;
	  ans1 := errorset!*(list('integratesq,mkquote u,
			     mkquote variable,mkquote loglist,
			     mkquote caar nexp),
			     !*backtrace);
	  erfg!* := erfg;
	  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 << printc "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;

% 18 Oct 07.

symbolic procedure mknewsqrt u;
   begin scalar v,w;
     if not !*galois then go to new;
     v:=addf(!*p2f mksp(!*pvar,2),negf !*q2f tidysqrt simp u);
     w:=errorset!*(list('afactor,mkquote v,mkquote !*pvar),t);
     if atom w then go to new
       else w:=car w;
     if cdr w then go to notnew;
   new:
     w := mksqrt reval u;
     listofnewsqrts:=w . listofnewsqrts;
     return !*kk2f w;
   notnew:
     w:=car w;
     v:=stt(w,!*pvar);
     if car v neq 1 then errach list("Error in mknewsqrt: ",v);
     w:=addf(w,multf(cdr v,(mksp(!*pvar,car v) .* -1) .+nil));
     v:=sqrt2top(w ./ cdr v);
     w:=quotf(numr v,denr v);
     if null w
       then begin scalar oldprop;
	       oldprop := get('sqrt,'simpfn);
	       put('sqrt,'simpfn,'simpsqrt);
	       v := simp!* prepsq v;
	       put('sqrt,'simpfn,oldprop);
	       if denr v = 1 then w := numr v
	    end;
     if null w then errach list("Division failure in mknewsqrt",u);
     return w
  end;

endpatch;

patch mathpr;

% 10 Dec 04, 29 Sep 06.

symbolic procedure dflayout u;
  (begin
      scalar op, args, w;
      w := car (u := cdr u);
      u := cdr u;
      if smemq('int,w) then !*noarg := nil;
      if !*noarg and (atom w or not get(car w, 'op)) then <<
         if atom w then <<
            op := w;
	    args := assoc(op, depl!*);
            if args then args := cdr args >>
         else <<
            op := car w;
	    args := cdr w >>;
         remember!-args(op, args);
         w := op >>;
      maprin w;
      if u then <<
	 u := layout!-formula('!!dfsub!! . u, 0, nil);
         if null u then return 'failed;
         w := 1 + cddr u;
         putpline((update!-pline(0, -w, caar u) . cdar u) .
                  ((cadr u - w) . (cddr u - w))) >>
   end) where !*noarg = !*noarg;

endpatch;

patch matrix;

% 26 Jun 04.

symbolic procedure sparse_backsub(exlis,varlis);
   begin scalar d,z,c;
   if null exlis then return nil;
   d := lc car exlis;
   foreach x in exlis do
      begin scalar s,p,v,r;
      p := lc x;
      v := mvar x;
      x := red x;
      while not domainp x and mvar x member varlis do
	  <<if (c := atsoc(mvar x,z)) then
               s := addf(multf(lc x,cdr c),s)
	    else r := addf(!*t2f lt x,r);
            x := red x>>;
      s := negf quotff(addf(multf(addf(r,x),d),s),p);
      z := (v . s) . z;
      end;
   for each p in z do cdr p := cancel(cdr p ./ d);
   return z
   end;

symbolic procedure quotff(u,v);
   if null u then nil
    else (if x then x
	   else (if denr y = 1 then numr y
		  else rederr "Invalid division in backsub")
		 where y=rationalizesq(u ./ v))
	  where x=quotf(u,v);

% 23 May 06.

symbolic procedure matsm1 u;
   begin scalar x,y,z; integer n;
    a:  if null u then return z
         else if eqcar(car u,'!*div) then go to d
         else if atom car u then go to er
         else if caar u eq 'mat then go to c1
	 else if flagp(caar u,'matmapfn) and cdar u
	    and getrtype cadar u eq 'matrix
	  then x := matsm matrixmap(car u,nil)
	 else <<x := lispapply(caar u,cdar u);
	  if eqcar(x,'mat) then x := matsm x>>;
    b:  z := if null z then x
              else if null cdr z and null cdar z then multsm(caar z,x)
              else multm(x,z);
    c:  u := cdr u;
        go to a;
    c1: if not lchk cdar u then rerror(matrix,3,"Matrix mismatch");
        x := for each j in cdar u collect
                for each k in j collect xsimp k;
        go to b;
    d:  y := matsm cadar u;
        if (n := length car y) neq length y
	  then rerror(matrix,4,"Non square matrix")
	 else if (z and n neq length z)
	  then rerror(matrix,5,"Matrix mismatch")
         else if cddar u then go to h
         else if null cdr y and null cdar y then go to e;
        x := subfg!*;
        subfg!* := nil;
        if null z then z := apply1(get('mat,'inversefn),y)
         else if null(x := get('mat,'lnrsolvefn))
          then z := multm(apply1(get('mat,'inversefn),y),z)
         else z := apply2(get('mat,'lnrsolvefn),y,z);
        subfg!* := x;
        z := for each j in z collect for each k in j collect
                 <<!*sub2 := t; subs2 k>>;
        go to c;
    e:  if null caaar y then rerror(matrix,6,"Zero divisor");
        y := revpr caar y;
        z := if null z then list list y else multsm(y,z);
        go to c;
     h: if null z then z := generateident n;
        go  to c;
    er: rerror(matrix,7,list("Matrix",car u,"not set"))
   end;

symbolic procedure matrixmap(u,v);
   if flagp(car u,'matmapfn)
    then matsm!*1 for each j in matsm cadr u collect
	    for each k in j collect simp!*(car u . mk!*sq k . cddr u)
   else if flagp(car u,'matfn) then reval2(u,v)
    else typerr(car u,"matrix operator");

put('matrix,'aggregatefn,'matrixmap);

flag('(int df taylor),'matmapfn);

flag('(det trace),'matfn);

endpatch;

patch poly;

% 6 Sep 04.

symbolic procedure fctrf u;
   (begin scalar !*ezgcd,!*gcd,denom,x,y;
      if domainp u then return list u
       else if ncmp!* and not noncomfp u then ncmp!* := nil;
      !*gcd := t;
      if null !*limitedfactors and null dmode!* then !*ezgcd := t;
      if null !*mcd
        then rerror(poly,15,"Factorization invalid with MCD off")
       else if null !*exp
        then <<!*exp := t; u := !*q2f resimp !*f2q u>>;
      if dmode!* eq '!:rn!:
	then <<dmode!* := nil; alglist!* := nil . nil;
	       x := simp prepf u;
	       if atom denr x then <<denom := denr x; u := numr x>>
		else denom := 1>>;
      if null ncmp!*
	then <<x := sf2ss u;
	       if homogp x
		 then <<if !*trfac
			then prin2t
		    "This polynomial is homogeneous - variables scaled";
			y := caaar x . listsum caaadr x;
			x := fctrf1 ss2sf(car(x)
				. (reverse subs0 cadr x . 1));
			x := rconst(y,x);
			return car x . sort!-factors cdr x>>>>;
      u := fctrf1 u;
      if denom
	then <<alglist!* := nil . nil;
	       dmode!* := '!:rn!:; car u := quotf!*(car u,denom)>>;
      return  car u . sort!-factors cdr u
   end) where !*exp = !*exp, ncmp!* = ncmp!*;

% 4 Oct 05.

symbolic procedure deg(u,kern);
   <<u := simp!* u; tstpolyarg2(u,kern); numrdeg(numr u,kern)>>
     where dmode!* = gdmode!*;

symbolic procedure tstpolyarg2(u,kern);
 <<for each j in kernels numr u do
      if j=kern then nil
       else if depends(j,kern) then typerr(prepsq u,"polynomial");
   for each j in kernels denr u do
      if depends(j,kern) then typerr(prepsq u,"polynomial")>>;

% 11 Jan 07.

symbolic procedure rnfactor!: u;
   begin scalar x,y,dmode!*; integer m,n;
     x := subf(u,nil);
     if not domainp denr x then return {1,(u . 1)};
     y := factorf numr x;
     n := car y;
     dmode!* := '!:rn!:;
     y := for each j in cdr y collect
           <<n := n*(m := (lnc ckrn car j)**cdr j);
             quotfd(car j,m) . cdr j>>;
     return int!-equiv!-chk mkrn(n,denr x) . y
   end;

endpatch;

patch rlfi;

put('tex,'simpfn,'simpcar);

endpatch;

% Solve declarations.

fluid '(!*cramer bareiss!-step!-size!*);

global '(assumptions);

patch solve;

% 26 Jun 04.

symbolic procedure solvelnrsys(exlis,varlis);
   begin scalar w,x;
   if w := solvesparsecheck(exlis,varlis) then exlis := w
     else exlis := exlis . varlis;
   if null !*cramer
      and null errorp(x :=
	   errorset2{'solvebareiss,mkquote car exlis, mkquote cdr exlis}
		      where bareiss!-step!-size!* = if w then 4 else 2)
     then exlis := car x
    else exlis := solvecramer(car exlis,cdr exlis);
   return solvesyspost(exlis,varlis)
   end;

% 12 Feb 05, 5 Oct 05.

symbolic procedure solvesq (ex,var,mul);
   begin scalar r,x;
      r:= for each w in solvesq1(ex,var,mul) join
	if null cadr w
	   or eqcar(x := prepsq caar w,'root_of)
	   or numr subfx(denr ex,{caadr w . x}) then {w};
      if r and not domainp denr ex then
	   assumptions:=append(assumptions,{prepf denr ex});
      return r
  end;

% 5 Oct 05.

symbolic procedure subfx(u,v);
   begin scalar x;
      x := errorset2 {'subf,mkquote u,mkquote v};
      return if errorp x then 1 ./ 1 else car x
   end;

% 12 Feb 05

symbolic procedure polypeval u;
   begin scalar bool,v;
      v := cadr u;
      u := simpcar u;
      if cdr u neq 1 then return nil else u := kernels car u;
      while u and null bool do
	<<if v neq car u and smember(v,car u) then bool := t;
	  u := cdr u>>;
      return null bool
   end;

put('polyp,'psopfn,'polypeval);

(algebraic <<

depend(!~p,!~x);

clearrules
{root_of(~p,~x,~tg)^~n =>
	sub(x=root_of(p,x,tg),
	    -reduct(p,x)/coeffn(p,x,deg(p,x)))^(n-deg(p,x)+1)
		when fixp n and deg(p,x)>=1 and n>=deg(p,x)};

let root_of(~p,~x,~tg)^~n =>
       sub(x=root_of(p,x,tg),
           -reduct(p,x)/coeffn(p,x,deg(p,x))) ^ (n-deg(p,x)+1)
	       when polyp(p,x) and fixp n and deg(p,x)>=1 and n>=deg(p,x);

nodepend(!~p,!~x);

>>) where dmode!*=nil,!*modular=nil,!*rounded=nil,!*complex=nil;

% 30 May 05.

symbolic procedure solvenonlnrtansolve(u,x,w);
  begin scalar v,s,z,r,y;
    integer ar;
   ar:=!!arbint;
   v:=caar u;u:=prepf numr simp cdr u;
   s:=solveeval{u,'tg!-};
   !!arbint:=ar;
   for each q in cdr s do
   <<z:=reval caddr q;
     z:=reval sublis(solvenonlnrtansolve1 z,z);
     !!arbint:=ar;
     y:=solve0({'equal,{'tan,{'quotient,V,2}},z},x); 
     r:=union(y,r)>>;
   y := errorset2 {'subf,mkquote w,mkquote{x . 'pi}};
   if null errorp y and null numr y
     then <<!!arbint:=ar; r:=union(solve0({'equal,{'cos,x},-1},x),r)>>;
   return t.r end;

% 5 Oct 05.

symbolic procedure check!-solns(z,ex,var);
   begin scalar x,y;
      if not errorp (x :=
	    errorset2 {'check!-solns1,mkquote z,mkquote ex,mkquote var})
	then return car x
       else if ex = (y := (numr simp!* prepf ex where !*reduced=t))
	  or errorp (x :=
	    errorset2 {'check!-solns1,mkquote z,mkquote y,mkquote var})
	then return 'unsolved
       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
	   or
	    errorp(y := errorset2 {'check!-solns2,mkquote ex,mkquote z})
	   then <<z := nil; x := 'unsolved>>
	  else if null(y := car 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;

symbolic procedure check!-solns2(ex,z);
   if smemq('root_of,z) then rederr 'check!-solns
    else numr subf(ex,{caadar z . mk!*sq caaar z});

% 27 Jul 08.

symbolic procedure solvequartic(a4,a3,a2,a1,a0);
   begin scalar x,y,yy,cx,z,s,l,zz1,zz2,dmode,neg,!*msg,!*numval,
                a1cr,a2cr,a3cr,xcr,ycr,zcr;
      dmode := dmode!*;
      a3 := quotsq(a3,a4);
      a2 := quotsq(a2,a4);
      a1 := quotsq(a1,a4);
      a0 := quotsq(a0,a4);
      yy := subtrsq(exptsq(a3,2),multfq(4,a2));
      x := solvecubic(!*f2q 1,
                      negsq a2,
                      subs2!* subtrsq(multsq(a1,a3),multfq(4,a0)),
                      subs2!* negsq addsq(exptsq(a1,2),
                                          multsq(a0,yy)));
      cx := car x;
      for each rr in x do if one_real rr then s := append(s,list rr);
      x := if (l := length s)=1 then car s else cx;
      a3 := quotsqf(a3,2); yy := quotsqf(yy,4);
      y := addsq(yy,x);
      if l<2 then go to zz;
loop: if not pos_num negsq y then go to zz else if l=1 then
        <<x := cx; y := addsq(yy,x); go to zz>>;
      l := l-1; s := cdr s; x := car s;
      y := addsq(yy,x); go to loop;
  zz: y := sqrtq y;
      x := quotsqf(x,2);
      z := sqrtq subtrsq(exptsq(x,2),a0);
      !*numval := t;
      on rounded,complex;
      a1cr := resimp a1; a2cr := resimp a2; a3cr := resimp a3;
      xcr := resimp x; ycr := resimp y; zcr := resimp z;
      if null numr
         (zz1 :=
           resimp subtrsq(a1cr,addsq(multsq(subtrsq(a3cr,ycr),addsq(xcr,zcr)),
                 multsq(addsq(a3cr,ycr),subtrsq(xcr,zcr))))) then go to rst;
      if null numr
         (zz2 :=
           resimp subtrsq(a1cr,addsq(multsq(subtrsq(a3cr,ycr),subtrsq(xcr,zcr)),
                 multsq(addsq(a3cr,ycr),addsq(xcr,zcr)))))
         then <<neg := t; go to rst>>;
      if domainp numr zz1 and domainp numr zz2
             and numberp denr zz1 and numberp denr zz2 and
         numr simp list('sign,list('difference,list('norm,mk!*sq zz1),
           list('norm,mk!*sq zz2)))=1 then neg := t;
 rst: off rounded,complex;
      if dmode then onoff(get(dmode,'dname),t);
      if neg then z := negsq z;
      return append(solvequadratic(!*f2q 1,subtrsq(a3,y),subtrsq(x,z)),
                    solvequadratic(!*f2q 1,addsq(a3,y),addsq(x,z)))
   end;

endpatch;

patch specfn;

% 20 Feb 06.

algebraic (for all u,m let df(Jacobidn(u,m),u)
	       = -m^2 *Jacobisn(u,m)*Jacobicn(u,m));

endpatch;

patch rlisp;

!#if (member 'psl lispsystem!*)
   symbolic procedure global idlist;
      fluid idlist;
   symbolic procedure global1 id1;
      if not get(id1,'vartype) then fluid1 id1;
!#endif

endpatch;


endmodule;

end;