reduce-historical/r38/packages/alg/simp.red

module simp; % Functions to convert prefix forms into canonical forms.

% Author: Anthony C. Hearn.

% Modifications by: J.H. Davenport, F. Kako, S. Kameny, E. Schruefer and
%                   Francis J. Wright.

% Copyright (c) 1998, Anthony C. Hearn.  All rights reserved.

fluid '(!*allfac !*div);

fluid '(!*asymp!* !*complex !*exp !*gcd !*ifactor !*keepsqrts !*mcd
	!*mode !*modular !*notseparate !*numval !*precise !*rationalize
	!*reduced !*resimp !*sub2 !*uncached alglist!* dmd!* dmode!*
	varstack!* !*combinelogs !*expandexpt !*msg frlis!* subfg!*
	!*norationalgi factorbound!* ncmp!* powlis1!* !*nospurp
	!*ncmp);

global '(!*match
         den!*
%        exptl!*   No-one else refers to this variable - just slows us
         initl!*
         mul!*
         simpcount!*
         simplimit!*
         tstack!*
         ws);

switch expandexpt; % notseparate;

!*expandexpt := t;

% The NOTSEPARATE switch inhibits an expression such as x^(4/3) to
% become x*x^(1/3).  At the present time, no one is using this.

factorbound!* := 10000;   % Limit for factoring with IFACTOR off.

% !*KEEPSQRTS uses SQRT rather than EXPT for square roots.
% Normally set TRUE in the integrator, false elsewhere.

put('ifactor,'simpfg,'((t (rmsubs))));

put('alglist!*,'initl,'(cons nil nil));

put('simpcount!*,'initl,0);

initl!* := union('(alglist!* simpcount!*),initl!*);

simplimit!* := 1000;

symbolic procedure noncom u;
   % Declare vars u to be noncom.
   <<rmsubs(); for each j in u do noncom1 j>>;

symbolic procedure noncom1 u;
   <<!*ncmp := t; flag(list u,'noncom)>>;

put('noncom,'stat,'rlis);

symbolic procedure simp!* u;
   begin scalar !*asymp!*,x;
        if eqcar(u,'!*sq) and caddr u and null !*resimp
          then return cadr u;
        x := mul!* . !*sub2;    % Save current environment.
        mul!* := nil;
        u:= simp u;
	if !*nospurp then mul!* := union(mul!*,'(isimpq));
        for each j in mul!* do u:= apply1(j,u);
        mul!* := car x;
        u := subs2 u;
        if !*combinelogs then u := clogsq!* u;
	% Must be here, since clogsq!* can upset girationalizesq!:.
	% For defint, it is necessary to turn off girationalizesq - SLK.
        if dmode!* eq '!:gi!: and not !*norationalgi
          then u := girationalize!: u
	  else if !*rationalize then u := rationalizesq u
	 else u := rationalizei u;
        !*sub2 := cdr x;
        % If any leading terms have cancelled, a gcd check is required.
        if !*asymp!* and !*rationalize then u := gcdchk u;
        return u
   end;

symbolic procedure rationalizei u;
   % Remove overall factor of i in denominator.
   begin scalar v,w;
      if domainp (v := denr u) or not smemq('i,v) then return u;
      v := reordsq u where kord!* = 'i . kord!*;
      return if lpow (w := denr v) = '(i . 1) and null red w
	       then negf multf(!*k2f 'i,reorder numr v) ./ reorder lc w
	      else u
   end;

symbolic procedure subs2 u;
   begin scalar xexp,v,w,x;
    if null subfg!* then return u
     else if !*sub2 or powlis1!* then u := subs2q u;
    u := exptchksq u;
    x := get('slash,'opmtch);
    if null (!*match or x) or null numr u then return u
     else if null !*exp
      then <<xexp:= t; !*exp := t; v := u; w := u := resimp u>>;
    u := subs3q u;
    if xexp then <<!*exp := nil; if u=w then u := v>>;
    if x then u := subs4q u;
    return u
   end;

symbolic procedure simp u;
   (begin scalar x,y;
    % This case is sufficiently common it is done first.
    if fixp u
      then if u=0 then return nil ./ 1
	    else if not dmode!* then return u ./ 1
	    else nil
     else if u member varstack!* then recursiveerror u;
    varstack!* := u . varstack!*;
    if simpcount!*>simplimit!*
      then <<simpcount!* := 0;
	     rerror(alg,12,"Simplification recursion too deep")>>
     else if eqcar(u,'!*sq) and caddr u and null !*resimp
      then return cadr u
     else if null !*uncached and (x := assoc(u,car alglist!*))
      then return <<if cadr x then !*sub2 := t; cddr x>>;
    simpcount!* := simpcount!*+1; % undone by returning through !*SSAVE.
    if atom u then return !*ssave(simpatom u,u)
     else if not idp car u or null car u
      then if atom car u then typerr(car u,"operator")
	    else if idp caar u and (x := get(caar u,'name))
	     then return !*ssave(u,u)     %%% not yet correct
	    else if eqcar(car u,'mat)
		and numlis(x := revlis cdr u) and length x=2
	     then return !*ssave(simp nth(nth(cdar u,car x),cadr x),u)
	    else errpri2(u,t)
     else if flagp(car u,'opfn)
      then if null(y := getrtype(x := opfneval u))
	     then return !*ssave(simp_without_resimp x,u)
	    else if y eq 'yetunknowntype and null getrtype(x := reval x)
	     then return simp x
	    else typerr(u,"scalar")
     else if x := get(car u,'psopfn)
      then if getrtype(x := apply1(x,cdr argnochk u))
	     then typerr(u,"scalar")
        else if x=u then return !*ssave(!*k2q x,u)
	else return !*ssave(simp_without_resimp x,u)
     % Note in above that the psopfn MUST return a *sq form,
     % otherwise an infinite recursion occurs.
     else if x := get(car u,'polyfn)
      then return
	<<argnochk u;
	  !*ssave(!*f2q lispapply(x,
			    for each j in cdr u collect !*q2f simp!* j),
		  u)>>
     else if get(car u,'opmtch)
        and not(get(car u,'simpfn) eq 'simpiden)
	and (x := opmtchrevop u)
      then return !*ssave(simp x,u)
     else if x := get(car u,'simpfn)
      then return !*ssave(apply1(x,
				 if x eq 'simpiden or flagp(car u,'full)
				   then argnochk u
				 else cdr argnochk u),
			  u)
     else if (x := get(car u,'rtype)) and (x := get(x,'getelemfn))
      then return !*ssave(simp apply1(x,u),u)
     else if flagp(car u,'boolean) or get(car u,'infix)
      then typerr(if x := get(car u,'prtch) then x else car u,
              "algebraic operator")
     else if flagp(car u,'nochange)
      then return !*ssave(simp lispeval u,u)
     else if get(car u,'psopfn) or get(car u,'rtypefn)
      then typerr(u,"scalar")
     else <<redmsg(car u,"operator");
        mkop car u;
	varstack!* := delete(u,varstack!*);
        return !*ssave(simp u,u)>>;
   end) where varstack!* = varstack!*;

symbolic procedure opmtchrevop u;
   % The following structure is designed to make index mu; p1.mu^2;
   % work.  It also introduces a redundant revlis in most cases.
   if null !*val or smemq('cons,u) then opmtch u
    else opmtch(car u . revlis cdr u);

symbolic procedure simp_without_resimp u;
   simp u where !*resimp := nil;

put('array,'getelemfn,'getelv);

put('array,'setelemfn,'setelv);

symbolic procedure getinfix u;
   %finds infix symbol for U if it exists;
   begin scalar x; return if x := get(u,'prtch) then x else u end;

symbolic procedure !*ssave(u,v);
   % We keep !*sub2 as well, since there may be an unsubstituted
   % power in U.
   begin
      if not !*uncached
    then rplaca(alglist!*,(v . (!*sub2 . u)) . car alglist!*);
      simpcount!* := simpcount!*-1;
      return u
   end;

symbolic procedure numlis u;
   null u or (numberp car u and numlis cdr u);

symbolic procedure simpatom u;
%  if null u then typerr("NIL","algebraic identifier")
   if null u then nil ./ 1   % Allow NIL as default 0.
    else if numberp u
     then if u=0 then nil ./ 1
	   else if not fixp u then ('!:rd!: . cdr fl2bf u) ./ 1
             % we assume that a non-fixp number is a float.
           else if flagp(dmode!*,'convert) and u neq 1 % Don't convert 1
            then !*d2q apply1(get(dmode!*,'i2d),u)
           else u ./ 1
    else if stringp u then typerr(list("String",u),"identifier")
    else if flagp(u,'share) then
      <<(if x eq u then mksq(u,1) else simp x) where x=lispeval u>>
    else begin scalar z;
      if z := get(u,'idvalfn) then return apply1(z,u)
       else if !*numval and dmode!* and flagp(u,'constant)
          and (z := get(u,dmode!*))
	  and not errorp(z := errorset!*(list('lispapply,mkquote z,nil),
			 nil))
        then return !*d2q car z
       else if getrtype u then typerr(u,'scalar)
       else return mksq(u,1) end;

flag('(e pi),'constant);

symbolic procedure mkop u;
   begin scalar x;
    if null u then typerr("Local variable","operator")
     else if (x := gettype u) eq 'operator
      then lprim list(u,"already defined as operator")
     else if x and not(x memq '(fluid global procedure))
      then typerr(u,'operator)
%    else if u memq frlis!* then typerr(u,"free variable")
     else put(u,'simpfn,'simpiden)
   end;

symbolic procedure operatorp u;
    gettype u eq 'operator;

symbolic procedure simpcar u;
   simp car u;

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

symbolic procedure share u;
   begin scalar y;
      for each v in u do
         if not idp v then typerr(v,"id")
         else if flagp(v,'share) then nil
         else if flagp(v,'reserved) then rsverr v
         else if (y := getrtype v) and y neq 'list
	    then rerror(alg,13,list(y,v,"cannot be shared"))
         else
          % if algebraic value exists, transfer to symbolic.
	 <<if y then remprop(v,'rtype);
	   if y := get(v,'avalue)
	     then <<setifngfl(v,cadr y); remprop(v,'avalue)>>
          % if no algebraic value but symbolic value, leave unchanged.
	    else if not boundp v then setifngfl(v,v);
          % if previously unset, set symbolic self pointer.
	   flag(list v,'share)>>
   end;

symbolic procedure boundp u;
   % Determines if the id u has a value.
   % NB:  this function must be redefined in many systems (e.g., CL).
   null errorp errorset!*(u,nil);

symbolic procedure setifngfl(v,y);
   <<if not globalp v then fluid list v; set(v,y)>>;

rlistat '(share);

flag('(ws !*mode),'share);

flag('(share),'eval);


% ***** SIMPLIFICATION FUNCTIONS FOR EXPLICIT OPERATORS - EXP *****

symbolic procedure simpexpon u;
   % Exponents must not use non-integer arithmetic unless NUMVAL is on,
   % in which case DOMAINVALCHK must know the mode.
   simpexpon1(u,'simp!*);

symbolic procedure simpexpon1(u,v);
   if !*numval and (dmode!* eq '!:rd!: or dmode!* eq '!:cr!:) 
     then apply1(v,u)
    else begin scalar dmode!*,alglist!*; return apply1(v,u) end;

symbolic procedure simpexpt u;
   % We suppress reordering during exponent evaluation, otherwise
   % internal parts (as in e^(a*b)) can have wrong order.
   begin scalar expon;
      expon := simpexpon carx(cdr u,'expt) where kord!*=nil;
      % We still need the right order, else
      % explog := {sqrt(e)**(~x*log(~y)/~z) => y**(x/z/2)};
      % on ezgcd,gcd; let explog; fails.
      expon := simpexpon1(expon,'resimp);
      return simpexpt1(car u,expon,nil)
   end;

symbolic procedure simpexpt1(u,n,flg);
   % FLG indicates whether we have done a PREPSQ SIMP!* U or not: we
   % don't want to do it more than once.
   begin scalar !*allfac,!*div,m,x,y;
      if onep u then return 1 ./ 1;
      !*allfac := t;
      m := numr n;
      if m=1 and denr n=1 then return simp u;
     % this simplifies e^(n log x) -> x^n  for all n,x.
      if u eq 'e and domainp denr n and not domainp m and ldeg m=1
         and null red m and eqcar(mvar m,'log) then return
	    simpexpt1(prepsq!* simp!* cadr mvar m,lc m ./ denr n,nil);
      if not domainp m or not domainp denr n
	then return simpexpt11(u,n,flg);
      x := simp u;
      if null m
	then return if null numr x then rerror(alg,14,"0**0 formed")
		     else 1 ./ 1;
      % We could use simp!* here, except it messes up the handling of
      % gamma matrix expressions.
%     if denr x=1 and not domainp numr x and not(denr n=1)
%       then <<y := sqfrf numr x;
%%      then <<y := fctrf numr x;
%%              if car y=1 then y := cdr y
%%               else if minusp car y then y := {1};
%              if length y>1 then return simpexptfctr(y,n)>>;
      return if null numr x
	       then if domainp m and minusf m
		      then rerror(alg,15,"Zero divisor")
		     else nil ./ 1
              else if atom m and denr n=1 and domainp numr x
                 and denr x=1
               then if atom numr x and m>0 then !*d2q(numr x**m)
		     else <<x := !:expt(numr x,m) ./ 1;
			    %remove rationals where possible.
			    if !*mcd then resimp x else x>>
              else if y := domainvalchk('expt,list(x,n)) then y
              else if atom m and denr n=1
	       then <<if not(m<0) then exptsq(x,m)
                       else if !*mcd then invsq exptsq(x,-m)
                       else multf(expf(numr x,m),mksfpf(denr x,-m))
                               ./ 1>>     % This uses OFF EXP option.
                      % There may be a pattern matching problem though.
     % We need the subs2 in the next line to take care of power and
     % product simplification left over from the call of simp on u.
	      else simpexpt11(if flg then u else prepsq!* subs2!* x,n,t)
   end;

symbolic procedure simpexptfctr(u,n);
   begin scalar x;
     x := 1 ./ 1;
     for each j in u do
	 x:= multsq(simpexpt1(prepf car j,multsq(cdr j ./ 1,n),nil),x);
     return x
   end;

symbolic procedure simpexpt11(u,n,flg);
   % Expand exponent to put expression in canonical form.
   begin scalar x;
      return if domainp denr n
		or not(car(x := qremf(numr n,denr n)) and cdr x)
	       then simpexpt2(u,n,flg)
	      else multsq(simpexpt1(u,car x ./ 1,flg),
			  simpexpt1(u,cdr x ./ denr n,flg))
   end;

symbolic procedure simpexpt2(u,n,flg);
   % The "non-numeric exponent" case.  FLG indicates whether we have
   % done a PREPSQ SIMP!* U or not: we don't want to do it more than
   % once.
   begin scalar m,n,x,y;
    if u=1 then return 1 ./ 1;
%  The following is now handled in mkrootsq.
%    else if fixp u and u>0 and (u<factorbound!* or !*ifactor)
%      and (length(x := zfactor u)>1 or cdar x>1)
%     then <<y := 1 ./ 1;
%            for each j in x do
%               y := multsq(simpexpt list(car j,
%                                         prepsq multsq(cdr j ./ 1,n)),
%                           y);
%            return y>>;
    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
%               and numberp numr n and not evenp numr n
%              then u := cadr u % list('abs,cadr u)
	       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)
     % We need the !*precise check for, say, sqrt((1+a)^2*y*z).
     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>>
        % For a product under *precise we isolate positive factors.
     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
%    The next lines did not allow, e.g., sqrt(a/b) => sqrt(a)/sqrt(b).
%    when precise is on and there is a risk of
%    E.g., sqrt(a/b) neq sqrt(a)/sqrt(b) when a=1, b=-1.
%    We allow however the denominator to be a positive number.
	and (not !*precise
%               or alg_constant_exptp(cadr u,n)
%               or alg_constant_exptp(caddr u,n)
		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})>>
     % Special case of (-expression)^(1/2).
%    else if car u eq 'minus
%            and (n = '(1 . 2) or n = '((!:rd!: . 0.5) . 1)
%                 or n = '((!:rd!: 5 . -1) . 1)
%                 or n = '((!:rn!: 1 . 2) . 1))
%     then return simptimes list('i,list('expt,cadr u,prepsq n))>>;
%    else if car u eq 'minus and numberp m and denr n=1
%     then return multsq(simpexpt list(-1,m),
%                simpexpt list(cadr u,m))>>;
     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 <<% Don't expand say e and pi, since whole expression is not
	     % numerical.
	    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;
       % "power" is not unique here.
    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)  % To make pattern matching work.
   end;

symbolic procedure posnump u;
   % True if u is a positive number. Test is naive but correct.
   if atom u then (numberp u and u>0) or u memq '(e pi)
    else if car u memq '(expt plus quotient sqrt times)
     then posnumlistp cdr u
    else nil;

symbolic procedure posnumlistp u;
   null u or posnump car u and posnumlistp cdr u;

% symbolic procedure alg_constant_exptp(u,v);
%    % U an expression, v a standard quotient.
%    alg_constantp u and alg_constantp car v and alg_constantp cdr v;

% symbolic procedure alg_constantp u;
%    % True if u is an algebraic constant whose surd is unique.
%    if atom u then numberp u
%    else if car u memq
%        '(difference expt plus minus quotient sqrt times)
%      then alg_constant_listp cdr u
%     else nil;

% symbolic procedure alg_constant_listp u;
%    null u or alg_constantp car u and alg_constant_listp cdr u;

put('expt,'simpfn,'simpexpt);

symbolic procedure split!-sign u;
  % U is a list of factors. Split into positive, negative
  % and unknown sign part. Nil if no sign is known.
  begin scalar p,n,w,s;
    for each f in u do
      if 1=(s:=sign!-of f) then p:=f.p else if -1=s then n:=f.n
          else w:=f.w;
    if null p and null n then return nil; 
    return p.n.w; 
  end;
     
symbolic procedure conv2gid(u,d);
   if null u or numberp u or eqcar(u,'!:gi!:) then d
    else if domainp u
     then if eqcar(u,'!:crn!:) then lcm(d,lcm(cdadr u,cdddr u))
	   else if eqcar(u,'!:rn!:) then lcm(d,cddr u) else d
    else conv2gid(lc u,conv2gid(red u,d));

symbolic procedure conv2gi2 u;
   if null u then u
   else if numberp u then u * den!*
   else if eqcar(u,'!:gi!:) then '!:gi!:.((den!**cadr u).(den!**cddr u))
   else if eqcar(u,'!:crn!:)
    then <<u := cdr u;
	   u:= '!:gi!: . ((den!*/cdar u*caar u).(den!*/cddr u*cadr u))>>
   else if eqcar(u,'!:rn!:) then den!*/cddr u*cadr u
   else if domainp u then rerror(alg,16,list("strange domain",u))
   else lpow u .* conv2gi2(lc u) .+ conv2gi2(red u);

symbolic procedure simpx1(u,m,n);
   % U,M and N are prefix expressions.
   % Value is the standard quotient expression for U**(M/N).
   % FLG is true if we have seen a "-" in M.
    begin scalar flg,x,z;
      % Check for imaginary result.
      if eqcar(u,'!*minus!*)
	 then if m=1 and fixp n and remainder(n,2)=0
	     or n=1 and eqcar(m,'quotient) and cadr m=1 and fixp caddr m
	       and remainder(caddr m,2)=0
		 then return multsq(simp list('expt,'i,
					      list('quotient,1,n/2)),
			  simpexpt list(cadr u,list('quotient,m,n)))
      % and for negative result.
		else if m=1 and fixp n          % n must now be odd.
		 then return negsq
			  simpexpt list(cadr u,list('quotient,m,n));
    if numberp m and numberp n
       or null(smemqlp(frlis!*,m) or smemqlp(frlis!*,n))
      then go to a;
    % exptp!* := t;
    return mksq(list('expt,u,if n=1 then m
                   else list('quotient,m,n)),1);
    a:
    if numberp m then
        if minusp m then <<m := -m; go to mns>>
           else if fixp m then
		   if fixp n then <<
                      if flg then m := -m;
		      z := m;
		      if !*mcd and (fixp u or null !*notseparate)
			then <<z := z-n*(m := m/n);
			       if z<0 then <<m := m-1; z := z+n>>>>
		       else m := 0;
		      x := simpexpt list(u,m);
		      if z=0 then return x
                      else if n=2 and !*keepsqrts
		       then <<x := multsq(x,apply1(get('sqrt,'simpfn),
						   list u));
			      % z can be 1 or -1. I'm not sure if other
			      % values can occur.
			      if z<0 then <<x := invsq x; z := -z>>;
			      return exptsq(x,z)>>
      % Note the indirect call: the integrator rebinds this property.
      % JHD understands this interaction - don't change without
      % consulting him.  Note that, since KEEPSQRTS is true, SIMPSQRT
      % won't recurse on SIMPEXPT1.
		      else return
			      multsq(x,exptsq(simprad(simp!* u,n),z))>>
		   else <<z := m; m := 1>>
                else z:=1
     else if atom m then z:=1
     else if car m eq 'minus then <<m := cadr m; go to mns>>
     else if car m eq 'plus and !*expandexpt then <<
         z := 1 ./ 1;
         for each x in cdr m do
             z := multsq(simpexpt list(u,
                         list('quotient,if flg then list('minus,x)
                                               else x,n)),
                           z);
         return z >>
%%   else if car m eq 'times and fixp cadr m and numberp n
%%    then <<
%%      z := gcdn(n,cadr m);
%%      n := n/z;
%%      z := cadr m/z;
%%      m := retimes cddr m >>
%% BEGIN modification by Francis J. Wright:
     else if car m eq 'times and fixp cadr m
      then <<
	if numberp n
	  then <<z := gcdn(n,cadr m); n := n/z; z := cadr m/z>>
	 else z := cadr m;
        % retimes seems to me to be overkill here, so try just ...
	m := if cdddr m then 'times . cddr m else caddr m>>
%% END   modification by FJW.
     else if car m eq 'quotient and n=1 and !*expandexpt
      then <<n := caddr m; m := cadr m; go to a>>
     else z := 1;
     if idp u and not flagp(u,'used!*) then flag(list u,'used!*);
        if u = '(minus 1)
               and n=1
               and null numr simp list('difference,m,'(quotient 1 2))
         then <<u := simp 'i; return if flg then negsq u else u>>;
    u := list('expt,u,if n=1 then m else list('quotient,m,n));
    return mksq(u,if flg then -z else z); %U is already in lowest terms;
    mns: %if numberp m and numberp n and !*rationalizeflag
     %  then return multsq(simpx1(u,n-m,n),invsq simp u) else
    % return invsq simpx1(u,m,n)
    if !*mcd then return invsq simpx1(u,m,n);
    flg := not flg;
    go to a;
   end;

symbolic procedure expf(u,n);
   %U is a standard form. Value is standard form of U raised to
   %negative integer power N. MCD is assumed off;
   %what if U is invertable?;
   if null u then nil
    else if u=1 then u
    else if atom u then mkrn(1,u**(-n))
    else if domainp u then !:expt(u,n)
    else if red u then mksp!*(u,n)
    else (lambda x; if x>0 and sfp mvar u
             then multf(exptf(mvar u,x),expf(lc u,n))
            else mvar u .** x .* expf(lc u,n) .+ nil)
     (ldeg u*n);

% ******* The "radical simplifier" section ******

symbolic procedure simprad(u,n);
   % Simplifies radical expressions.
   if !*reduced then multsq(radfa(numr u,n),invsq radfa(denr u,n))
     else begin scalar iflag,x,y,z;
       if !*rationalize then << % Move all radicands into numerator.
	  y:=list(denr u,1); % A partitioned expression.
          u:=multf(numr u, exptf(denr u,n-1)) ./ 1 >>
         else y := radf(denr u,n);
       if n=2 and minusf numr u % Should this be 'evenp n'?
	 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
      % This test allows transformations like sqrt(2/3)=>sqrt(2)/sqrt(3)
      % whereas we really don't want to do this for symbolic elements
      % since we can introduce paradoxes that way.
	 then z := multsq(mkrootsq(prepf numr z,n),
			  invsq mkrootsq(prepf denr z,n))
	else <<if iflag
		 then <<iflag := nil; % Absorb the "i" in square root.
			z := negsq z>>;
	       z := mkrootsq(prepsq z,n)>>;
       z := multsq(multsq(if !*precise and evenp n 
			    then car x ./ 1   % mkabsf0 car x
                           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);
    % if !*precise and evenp n then y := mkabsf0 car x ./ 1 else
      y := exptf(car x,m) ./ 1;
      return multsq(exptsq(mkrootlsq(cdr x,n),m),y)
   end;

symbolic procedure mkrootlsq(u,n);
   % U is a list of prefix expressions, N an integer.
   % Value is standard quotient for U**(1/N);
   % NOTE we need the REVAL call so that PREPSQXX is properly called on
   % the argument for consistency with the pattern matcher.  Otherwise
   % for all x,y let sqrt(x)*sqrt(y)=sqrt(x*y); sqrt(30*(l+1))*sqrt 5;
   % goes into an infinite loop.
   if null u then !*d2q 1
    else if null !*reduced then mkrootsq(reval retimes u,n)
    else mkrootlsq1(u,n);

symbolic procedure mkrootlsq1(u,n);
   if null u then !*d2q 1
    else multsq(mkrootsq(car u,n),mkrootlsq1(cdr u,n));

symbolic procedure mkrootsq(u,n);
   % U is a prefix expression, N an integer.
   % Value is a standard quotient for U**(1/N).
   if u=1 then !*d2q 1
    else if n=2 and (u= -1 or u= '(minus 1)) then simp 'i
    else if eqcar(u,'expt) and fixp caddr u
     then exptsq(mkrootsq(cadr u,n),caddr u)
    else begin scalar x,y;
	    if fixp u and not minusp u
		      and (length(x := 
                            zfactor1(u,u<factorbound!* or !*ifactor))>1
                           or cdar x>1)
	      then return mkrootsql(x,n);
	    x := if n=2 then mksqrt u
		  else list('expt,u,list('quotient,1,n));
	    if y := opmtch x then return simp y
	     else return mksq(x,1)
         end;

symbolic procedure mkrootsql(u,n);
   if null u then !*d2q 1
    else if cdar u>1
     then multsq(exptsq(mkrootsq(caar u,n),cdar u),mkrootsql(cdr u,n))
    else multsq(mkrootsq(caar u,n),mkrootsql(cdr u,n));


comment The following four procedures return a partitioned root
    expression, which is a dotted pair of integral part (a standard
    form) and radical part (a list of prefix expressions). The whole
    structure represents U**(1/N);

symbolic procedure check!-radf!-sign(rad,result,n);
   % Changes the sign of result if result**n = -rad. rad and result are
   % s.f.'s, n is an integer.
   (if evenp n and s = -1 or
       not evenp n and numberp s and 
	((numberp s1 and s neq s1)
	   where s1 = reval {'sign,mk!*sq !*f2q rad})
      then negf result
     else result)
    where s = reval{'sign,mk!*sq !*f2q result};

symbolic procedure radf(u,n);
   % U is a standard form, N a positive integer. Value is a partitioned
   % root expression for U**(1/N).
   begin scalar ipart,rpart,x,y,z,!*gcd,!*mcd;
      if null u then return list u;
      !*gcd := !*mcd := t;  % mcd cannot be off in this code.
      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 if !*precise
%                    then !*p2f mksp(numr
%                    then exptf(numr
%                                    simp list('abs,if sfp mvar u
%                                                     then prepf mvar u
%                                                    else mvar u),
                    else check!-radf!-sign(!*p2f(mvar u .** pdeg car y),
                                           !*p2f(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);   % We need *exp on here.
       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
		   <<%if !*precise and evenp n
		      % then y := !*kk2f {'abs,prepf y};
		     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;  % was aconc(rpart,z).
      return ipart . rpart
   end;

symbolic procedure radf1(u,n);
   %U is a form_power list, N a positive integer. Value is a
   %partitioned root expression for U**(1/N);
   begin scalar ipart,rpart,x;
      ipart := 1;
      for each z in u do
     <<x := divide(cdr z,n);
       if not(car x=0)
            then ipart := multf(
                 check!-radf!-sign(!*p2f z,exptf(car z,car x),n),ipart);
          if not(cdr x=0)
            then rpart := mkexpt(prepsq!*(car z ./ 1),cdr x)
                   . rpart>>;
      return ipart . rpart
   end;

symbolic procedure radd(u,n);
   %U is a domain element, N an integer.
   %Value is a partitioned root expression for U**(1/N);
   begin scalar bool,ipart,x;
      if not atom u then return list(1,prepf u);
%      then if x := integer!-equiv u then u := x
%            else return list(1,prepf u);
      if u<0 and evenp n then <<bool := t; u := -u>>;
      x := nrootnn(u,n);
      if bool then if !*reduced and n=2
             then <<ipart := multd(car x,!*k2f 'i);
                x := cdr x>>
            else <<ipart := car x; x := -cdr x>>
       else <<ipart := car x; x := cdr x>>;
      return if x=1 then list ipart else list(ipart,x)
   end;

% symbolic procedure iroot(m,n);
%    %M and N are positive integers.
%   %If M**(1/N) is an integer, this value is returned, otherwise NIL;
%   begin scalar x,x1,bk;
%      if m=0 then return m;
%      x := 10**iroot!-ceiling(lengthc m,n);   %first guess;
%   a: x1 := x**(n-1);
%      bk := x-m/x1;
%      if bk<0 then return nil
%       else if bk=0 then return if x1*x=m then x else nil;
%      x := x - iroot!-ceiling(bk,n);
%      go to a
%   end;

 symbolic procedure iroot(n,r);
    % N, r are integers; r >= 1.  If n is an exact rth power then its
    % rth root is returned, otherwise NIL.
    begin scalar tmp;
       tmp := irootn(n,r);
       return if tmp**r = n then tmp else nil
  end;

symbolic procedure iroot!-ceiling(m,n);
   %M and N are positive integers. Value is ceiling of (M/N) (i.e.,
   %least integer greater or equal to M/N);
   (lambda x; if cdr x=0 then car x else car x+1) divide(m,n);

symbolic procedure mkexpt(u,n);
   if n=1 then u else list('expt,u,n);

% The following definition is due to Eberhard Schruefer.

symbolic procedure nrootn(n,x);
   % N is an integer, x a positive integer. Value is a pair
   % of integers r,s such that r*s**(1/x)=n**(1/x).
   begin scalar fl,r,s,m,signn;
     r := 1;
     s := 1;
     if n<0 then <<n := -n; if evenp x then signn := t else r := -1>>;
     fl := zfactor n;
     for each j in fl do
         <<m := divide(cdr j,x);
           r := car j**car m*r;
           s := car j**cdr m*s>>;
     if signn then s := -s;
     return r . s
   end;

% symbolic procedure nrootn(n,x);
%   % N is an integer, X a positive integer. Value is a pair
%   % of integers I,J such that I*J**(1/X)=N**(1/X).
%   begin scalar i,j,r,signn;
%      r := 1;
%      if n<0 then <<n := -n; if evenp x then signn := t else r := -1>>;
%      j := 2**x;
%      while remainder(n,j)=0 do <<n := n/j; r := r*2>>;
%      i := 3;
%      j := 3**x;
%      while j<=n do
%         <<while remainder(n,j)=0 do <<n := n/j; r := r*i>>;
%           if remainder(i,3)=1 then i := i+4 else i := i+2;
%           j := i**x>>;
%      if signn then n := -n;
%      return r . n
%   end;

% ***** simplification functions for other explicit operators *****

symbolic procedure simpiden u;
   % Convert the operator expression U to a standard quotient.
   % Note: we must use PREPSQXX and not PREPSQ* here, since the REVOP1
   % in SUBS3T uses PREPSQXX, and terms must be consistent to prevent a
   % loop in the pattern matcher.
   begin scalar bool,fn,x,y,z;
    fn := car u; u := cdr u;
    % Allow prefix ops with names of symbolic functions.
    if (get(fn,'!:rn!:) or get(fn,'!:rd!:)) and (x := valuechk(fn,u))
      then return x;
    % Keep list arguments in *SQ form.
    if u and eqcar(car u,'list) and null cdr u
      then return mksq(list(fn,aeval car u),1);
    x := for each j in u collect aeval j;
    u := for each j in x collect
              if eqcar(j,'!*sq) then prepsqxx cadr j
               else if numberp j then j
               else <<bool := t; j>>;
%   if u and car u=0 and (flagp(fn,'odd) or flagp(fn,'oddreal))
    if u and car u=0 and flagp(fn,'odd)
         and not flagp(fn,'nonzero)
      then return nil ./ 1;
    u := fn . u;
    if flagp(fn,'noncom) then ncmp!* := t;
    if null subfg!* then go to c
     else if flagp(fn,'linear) and (z := formlnr u) neq u
      then return simp z
     else if z := opmtch u then return simp z;
 %   else if z := get(car u,'opvalfn) then return apply1(z,u);
 %    else if null bool and (z := domainvalchk(fn,
 %                for each j in x collect simp j))
 %     then return z;
    c:  if flagp(fn,'symmetric) then u := fn . ordn cdr u
     else if flagp(fn,'antisymmetric)
      then <<if repeats cdr u then return (nil ./ 1)
          else if not permp(z:= ordn cdr u,cdr u) then y := t;
         % The following patch was contributed by E. Schruefer.
         fn := car u . z;
         if z neq cdr u and (z := opmtch fn)
           then return if y then negsq simp z else simp z;
         u := fn>>;
%    if (flagp(fn,'even) or flagp(fn,'odd))
%       and x and minusf numr(x := simp car x)
%     then <<if flagp(fn,'odd) then y := not y;
%   if (flagp(fn,'even) or flagp(fn,'odd) or flagp(fn,'oddreal)
%          and x and not_imag_num car x)
    if (flagp(fn,'even) or flagp(fn,'odd))
	 and x and minusf numr(x := simp car x)
     then <<if not flagp(fn,'even) then y := not y;
        u := fn . prepsqxx negsq x . cddr u;
        if z := opmtch u
          then return if y then negsq simp z else simp z>>;
    u := mksq(u,1);
    return if y then negsq u else u
   end;

switch rounded;

symbolic procedure not_imag_num a;
 % Tests true if a is a number that is not a pure imaginary number.
 % Rebinds sqrtfn and *keepsqrts to make integrator happy.
   begin scalar !*keepsqrts,!*msg,!*numval,dmode,sqrtfn;
      dmode := dmode!*;
      !*numval := t;
      sqrtfn := get('sqrt,'simpfn);
      put('sqrt,'simpfn,'simpsqrt);
      on rounded,complex;
      a := resimp simp a;
      a := numberp denr a and domainp numr a and numr repartsq a;
      off rounded,complex;
      if dmode then onoff(get(dmode,'dname),t);
      put('sqrt,'simpfn,sqrtfn);
      return a
   end;

flagop even,odd,nonzero;

symbolic procedure domainvalchk(fn,u);
   begin scalar x;
      if (x := get(dmode!*,'domainvalchk)) then return apply2(x,fn,u);
      % The later arguments tend to be smaller ...
      u := reverse u;
  a:  if null u then return valuechk(fn,x)
       else if denr car u neq 1 then return nil;
      x := mk!*sq car u . x;
      u := cdr u;
      go to a
   end;

symbolic procedure valuechk(fn,u);
   begin scalar n;
      if (n := get(fn,'number!-of!-args)) and length u neq n
         or not n
            and u and cdr u and (get(fn,'!:rd!:) or get(fn,'!:rn!:))
       then rerror(alg,17,list("Wrong number of arguments to",fn));
      u := opfchk!!(fn . u);
      if u then return znumrnil
          ((if eqcar(u,'list) then list((u . 1) . 1) else u) ./ 1)
   end;

symbolic procedure znumrnil u; if znumr u then nil ./ 1 else u;

symbolic procedure znumr u;
   null (u := numr u) or numberp u and zerop u
   or not atom u and domainp u and
     (y and apply1(y,u) where y=get(car u,'zerop));

symbolic procedure opfchk!! u;
   begin scalar fn,fn1,sf,sc,int,ce; fn1 := fn := car u; u := cdr u;
     % first save fn and check to see whether fn is defined.
     % Integer functions are defined in !:rn!:,
     % real functions in !:rd!:, and complex functions in !:cr!:.
      fn := if flagp(fn,'integer) then <<int := t; get(fn,'!:rn!:)>>
         else if !*numval and dmode!* memq '(!:rd!: !:cr!:)
            then get(fn,'!:rd!:);
      if not fn then return nil;
      sf := if int then 'simprn
         else if (sf := get(fn,'simparg)) then sf else 'simprd;
     % real function fn is defined.  now check for complex argument.
      if int or not !*complex then go to s; % the simple case.
     % mode is complex, so check for complex argument.
     % list argument causes a slight complication.
      if eqcar(car u,'list)
         then if (sc := simpcr revlis cdar u) and eqcar(sc,nil)
            then go to err else go to s;
      if not (u := simpcr revlis u) then return nil
     % if fn1 = 'expt, then evaluate complex function only; else
     % if argument is real, evaluate real function, but if error
     % occurs, then evaluate complex function.
      else if eqcar(u,nil) or
	  fn1 eq 'expt and rd!:minusp caar u then u := cdr u
	 else <<ce := cdr u; u := car u; go to s>>;
     % argument is complex or real function failed.
     % now check whether complex fn is defined.
 evc: if fn := get(fn1,'!:cr!:) then go to a;
 err: rerror(alg,18,list(fn1,"is not defined as complex function"));
   s: if not (u := apply1(sf, revlis u)) then return nil;
   a: u := errorset!*(list('apply,mkquote fn,mkquote u),nil);
      if errorp u then
         if ce then <<u := ce; ce := nil; go to evc>> else return nil
       else return if int then intconv car u else car u
   end;

symbolic procedure intconv x;
   if null dmode!* or dmode!* memq '(!:rd!: !:cr!:) then x
   else apply1(get(dmode!*,'i2d),x);

symbolic procedure simpcr x;
 % Returns simprd x if all args are real, else nil . "simpcr" x.
  if atom x then nil else
   <<(<<if not errorp y then z := car y;
        y := simplist x where dmode!* = '!:cr!:;
        if y then z . y else z>>)
   where z=nil,y=errorset!*(list('simprd,mkquote x),nil)>>;

symbolic procedure simprd x;
   % Converts any argument list that can be converted to list of rd's.
   if atom x then nil else <<simplist x where dmode!* = '!:rd!:>>;

symbolic procedure simplist x;
   begin scalar fl,c; c := get(dmode!*,'i2d);
     x := for each a in x collect (not fl and
	<<if null (a := mconv numr b) then a := 0;
	  if numberp a then a := apply1(c,a)
	     else if not(domainp a and eqcar(a,dmode!*)) then fl := t;
	  if not fl and
	     (numberp(b := mconv denr b) and (b := apply1(c,b))
	     or domainp b and eqcar(b,dmode!*))
		then apply2(get(dmode!*,'quotient),a,b) else fl := t>>
	where b=simp!* a);
    if not fl then return x
   end;

symbolic procedure mconv v; <<dmconv0 dmode!*; mconv1 v>>;

symbolic procedure dmconv0 dmd;
   dmd!* := if null dmd then '!:rn!:
        else if dmd eq '!:gi!: then '!:crn!: else dmd;

symbolic procedure dmconv1 v;
   if null v or eqcar(v,dmd!*) then v
   else if atom v then if flagp(dmd!*,'convert)
      then apply1(get(dmd!*,'i2d),v) else v
   else if domainp v then apply1(get(car v,dmd!*),v)
   else lpow v .* dmconv1(lc v) .+ dmconv1(red v);

symbolic procedure mconv1 v;
   if domainp v then drnconv v
   else lpow v .* mconv1(lc v) .+ mconv1(red v);

symbolic procedure drnconv v;
   if null v or numberp v or eqcar(v,dmd!*) then v else
   <<(if y and atom y then apply1(y,v) else v)
   where y=get(car v,dmd!*)>>;

% Absolute Value Function.

symbolic procedure simpabs u;
   if null u or cdr u then mksq('abs . revlis u, 1)  % error?.
    else begin scalar x;
      u := car u;
      if numberp u then return abs u ./ 1
       else if x := sign!-abs u then return x;
      u := simp!* u;
      return if null numr u then nil ./ 1
        else quotsq(simpabs1 numr u, simpabs1 denr u);
    end;

symbolic procedure simpabs1 u;
   % Currently abs(sqrt(2)) does not simplify, whereas it clearly
   % should simplify to just sqrt(2).  The facts that abs(i) -> 1 and
   % abs(sqrt(-2)) -> abs(sqrt(2)) imply that REDUCE regards abs as
   % the complex modulus function, in which case I think it is always
   % correct to commute abs and sqrt.  However, I will do this only if
   % the result is a simplification.  FJW, 18 July 1998
   begin scalar x,y,w;
      x:=prepf u; u := u ./ 1;
      if eqcar(x,'minus) then x:=cadr x;
      % FJW: abs sqrt y -> sqrt abs y if abs y simplifies.
      if eqcar(x,'sqrt) then
	 return !*kk2q if eqcar(y:=reval('abs.cdr x), 'abs)
	 then {'abs, x} else {'sqrt, y};
%%    if eqcar(x,'times) and (y:=split!-sign cdr x) then
%%    <<w:=simp!* retimes car y; u:=quotsq(u,w);
%%      if cadr y then
%%       <<y:=simp!* retimes cadr y; u:=quotsq(u,y);
%%         w:=multsq(negsq y,w)>>
%%    >>;
      if eqcar(x,'times) then
      begin scalar abslist, noabs;
	 for each fac in cdr x do
	    % FJW: abs sqrt y -> sqrt abs y if abs y simplifies.
	    if eqcar(fac,'sqrt)
	       and not eqcar(y:=reval('abs.cdr fac), 'abs)
	    then noabs := {'sqrt, y} . noabs
	    else abslist := fac . abslist;
	 abslist := reversip abslist;
	 if noabs then
	    u := quotsq(u, noabs := simp!*('times . reversip noabs));
	 if (y:=split!-sign abslist) then
	 <<w:=simp!* retimes car y; u:=quotsq(u,w);
	    if cadr y then
	    <<y:=simp!* retimes cadr y; u:=quotsq(u,y);
	       w:=multsq(negsq y,w)>>;
	    if noabs then w := multsq(noabs, w)
	 >>
	 else w := noabs
      end;
      if numr u neq 1 or denr u neq 1 then
	 u:=quotsq(mkabsf1 absf numr u,mkabsf1 denr u);
      if w then u:=multsq(w,u);
      return u
   end;

%symbolic procedure rd!-abs u;
%   % U is a prefix expression.  If it represents a constant, return the
%   % abs of u.
%   (if !*rounded or not constant_exprp u then nil
%    else begin scalar x,y,dmode!*;
%	    setdmode('rounded,t) where !*msg := nil;
%	    x := aeval u;
%	    if evalnumberp x 
%	      then if null !*complex or 0=reval {'impart,x}
%		     then y := if evalgreaterp(x,0) then u
%				else if evalequal(x,0) then 0
%				else {'minus,u};
%	    setdmode('rounded,nil) where !*msg := nil;
%	    return if y then simp y else nil
%       end) where alglist!*=alglist!*;

symbolic procedure sign!-abs u;
   % Sign based evaluation of abs - includes the above rd!-abs
   % method as sub-branch.
    <<if not numberp n then nil else
      simp if n<0 then {'minus,u} else if n=0 then 0 else u
      >> where n=sign!-of u;

symbolic procedure constant_exprp u;
   % True if u evaluates to a constant (i.e., number).
   if atom u
     then numberp u or flagp(u,'constant) or u eq 'i and idomainp()
    else (flagp(car u,'realvalued)
	     or flagp(car u,'alwaysrealvalued)
	     or car u memq '(plus minus difference times quotient)
	     or get(car u,'!:rd!:)
	     or !*complex and get(car u,'!:cr!:))
	  and not atom cdr u
	  and constant_expr_listp cdr u;

symbolic procedure constant_expr_listp u;
   % True if all members of u are constant_exprp.
   % U can be a dotted pair as well as a list.
   if atom u
     then null u or numberp u or flagp(u,'constant)
	or u eq 'i and idomainp()
    else constant_exprp car u and constant_expr_listp cdr u;

symbolic procedure mkabsf0 u; simp{'abs,mk!*sq !*f2q u};

symbolic procedure mkabsf1 u;
   if domainp u then mkabsfd u
    else begin scalar x,y,v;
	   x := comfac!-to!-poly comfac u;
	   u := quotf1(u,x);
	   y := split!-comfac!-part x;
	   x := cdr y;
	   y := car y;
	   if positive!-sfp u then <<y := multf(u,y); u := 1>>;
	   u := multf(u,x);
	   v := lnc y;
	   y := quotf1(y,v);
	   v := multsq(mkabsfd v,y ./ 1);
	   return if u = 1 then v
		   else multsq(v,simpiden list('abs,prepf absf u))
	end;

symbolic procedure mkabsfd u;
   if null get('i,'idvalfn) then absf u ./ 1
    else (simpexpt list(prepsq nrm,'(quotient 1 2))
	  where nrm = addsq(multsq(car us,car us),
			     multsq(cdr us,cdr us))
	  where us = splitcomplex u);

symbolic procedure positive!-sfp u;
   if domainp u
      then if get('i,'idvalfn)
	      then !:zerop impartf u and null !:minusp repartf u
	    else null !:minusp u
    else positive!-powp lpow u and positive!-sfp lc u
	 and positive!-sfp red u;

symbolic procedure positive!-powp u;
   not atom car u and caar u memq '(abs norm);

% symbolic procedure positive!-powp u;
%    % This definition allows for the testing of positive valued vars.
%    if atom car u then flagp(car u, 'positive)
%     else ((if x then apply2(x,car u,cdr u) else nil)
%           where x = get(caar u,'positivepfn));

symbolic procedure split!-comfac!-part u;
   split!-comfac(u,1,1);

symbolic procedure split!-comfac(u,v,w);
   if domainp u then multd(u,v) . w
    else if red u then
       if positive!-sfp u then multf(u,v) . w
	else v . multf(u,w)
    else if mvar u eq 'i then split!-comfac(lc u,v,w)
    else if positive!-powp lpow u
      then split!-comfac(lc u,multpf(lpow u,v),w)
    else split!-comfac(lc u,v,multpf(lpow u,w));

put('abs,'simpfn,'simpabs);

symbolic procedure simpdiff u;
   <<ckpreci!# u; addsq(simpcar u,simpminus cdr u)>>;

put('difference,'simpfn,'simpdiff);

symbolic procedure simpminus u;
   negsq simp carx(u,'minus);

put('minus,'simpfn,'simpminus);

symbolic procedure simpplus u;
   begin scalar z;
     if length u=2 then ckpreci!# u;
     z := nil ./ 1;
  a: if null u then return z;
     z := addsq(simpcar u,z);
     u := cdr u;
     go to a
   end;

put('plus,'simpfn,'simpplus);

symbolic procedure ckpreci!# u;
   % Screen for complex number input.
   !*complex
      and (if a and not b then ckprec2!#(cdar u,cadr u)
	    else if b and not a then ckprec2!#(cdadr u,car u))
	    where a=timesip car u,b=timesip cadr u;

symbolic procedure timesip x; eqcar(x,'times) and 'i memq cdr x;

symbolic procedure ckprec2!#(im,rl);
   % Strip im and rl to domains.
   <<im := if car im eq 'i then cadr im else car im;
     if eqcar(im,'minus) then im := cadr im;
     if eqcar(rl,'minus) then rl := cadr rl;
     if domainp im and domainp rl and not(atom im and atom rl)
        then ckprec3!#(!?a2bf im,!?a2bf rl)>>;

remflag('(!?a2bf),'lose);   % Until things stabilize.

symbolic smacro procedure make!:ibf (mt, ep);
   '!:rd!: . (mt . ep);

symbolic smacro procedure i2bf!: u; make!:ibf (u, 0);

symbolic procedure !?a2bf a;
   % Convert decimal or integer to bfloat.
   if atom a then if numberp a then i2bf!: a else nil
    else if eqcar(a,'!:dn!:) then a;

symbolic procedure ckprec3!#(x,y);
  % if inputs are valid, check for precision increase.
   if x and y then
   precmsg max(length explode abs cadr x+cddr x,
               length explode abs cadr y+cddr y);

symbolic procedure simpquot q;
   (if null numr u
      then if null numr v then rerror(alg,19,"0/0 formed")
	    else rerror(alg,20,"Zero divisor")
     else if dmode!* memq '(!:rd!: !:cr!:) and domainp numr u
        and domainp denr u and domainp denr v
            and !:onep denr u and !:onep denr v
      then (if null numr v then nil else divd(numr v,numr u)) ./ 1 
     else <<q := multsq(v,simprecip cdr q);
	    if !*modular and null denr q
	      then rerror(alg,201,"Zero divisor");
	    q>>)
   where v=simpcar q,u=simp cadr q;

put('quotient,'simpfn,'simpquot);

symbolic procedure simprecip u;
   if null !*mcd then simpexpt list(carx(u,'recip),-1)
    else invsq simp carx(u,'recip);

put('recip,'simpfn,'simprecip);

symbolic procedure simpset u;
  begin scalar x;
     x := prepsq simp!* car u;
     if null x % or not idp x
       then typerr(x,"set variable");
     let0 list(list('equal,x,mk!*sq(u := simp!* cadr u)));
     return u
  end;

put ('set, 'simpfn, 'simpset);

symbolic procedure simpsqrt u;
   if u=0 then nil ./ 1 else
   if null !*keepsqrts
     then simpexpt1(car u, simpexpon '(quotient 1 2), nil)
    else begin scalar x,y;
       x := xsimp car u;
       return if null numr x then nil ./ 1
	       else if denr x=1 and domainp numr x and !:minusp numr x
		then if numr x=-1 then simp 'i
		      else multsq(simp 'i,
				  simpsqrt list prepd !:minus numr x)
	       else if y := domainvalchk('sqrt,list x) then y
	       else simprad(x,2)
     end;

symbolic procedure xsimp u; expchk simp!* u;

symbolic procedure simptimes u;
   begin scalar x,y;
    if null u then return 1 ./ 1;
    if tstack!* neq 0 or null mul!* then go to a0;
    y := mul!*;
    mul!* := nil;
    a0: tstack!* := tstack!*+1;
    x := simpcar u;
    a:    u := cdr u;
    if null numr x then go to c
     else if null u then go to b;
    x := multsq(x,simpcar u);
    go to a;
    b:    if null mul!* or tstack!*>1 then go to c;
    x:= apply1(car mul!*,x);
    alglist!* := nil . nil;   % since we may need MUL!* set again.
    mul!*:= cdr mul!*;
    go to b;
    c:    tstack!* := tstack!*-1;
    if tstack!* = 0 then mul!* := y;
    return x;
   end;

put('times,'simpfn,'simptimes);

symbolic procedure resimp u;
   % U is a standard quotient.
   % Value is the resimplified standard quotient.
   resimp1 u where varstack!*=nil;

symbolic procedure resimp1 u;
   begin
      u := quotsq(subf1(numr u,nil),subf1(denr u,nil));
      !*sub2 := t;
      return u
   end;

symbolic procedure simp!*sq u;
   if cadr u and null !*resimp then car u else resimp1 car u;

put('!*sq,'simpfn,'simp!*sq);

endmodule;

end;