reduce-historical/r38/packages/int/simpsqrt.red

module simpsqrt;   % Simplify square roots.

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

fluid '(!*galois !*pvar !*tra !*trint basic!-listofallsqrts
	gaussiani basic!-listofnewsqrts intvar knowntobeindep
        listofallsqrts listofnewsqrts sqrtflag sqrtlist
        sqrt!-places!-alist varlist zlist);

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

% remd 'simpsqrt;

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

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


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


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


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


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


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


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

symbolic procedure mknewsqrt u;
   % U is prefix form.
   begin scalar v,w;
     if not !*galois then go to new;
       % no checking required.
     v:=addf(!*p2f mksp(!*pvar,2),negf !*q2f tidysqrt simp u);
     w:=errorset!*(list('afactor,mkquote v,mkquote !*pvar),t);
     if atom w then go to new
       else w:=car w; % The actual result of afactor.
     if cdr w then go to notnew;
   new:
     w := mksqrt reval u;  % Note that u need not be a canonical
			   % structure here.
     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
       % We now test to see if the quotient failure is spurious, e.g.,
       % as in int(-2x/(sqrt(2x^2+1)-2x^2+1),x); It's not clear this is
       % the right place to check though.  More information is
       % available from the earlier int-sqrt step.
       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;

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


symbolic procedure actualsqrtinner(sf,l);
   if sf =1 then 1
    else if null l
       or domainp sf or ldeg sf=1
   % Patch by A.C. Norman to prevent recursion errors.
     then actualsimpsqrt2 sf
    else begin scalar z;
       if numberp sf and (z := list('sqrt,sf)) member l
         then return !*kk2f z;
       z := argof car l;
       if z member l then z := !*kk2f car l else z := !*q2f simp z;
       if z = -1 then return actualsqrtinner(sf,cdr l);
       z:=quotf(sf,z);
       if null z then return actualsqrtinner(sf,cdr l)
        else return !*multf(!*kk2f car l,actualsimpsqrt z)
     end;


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


symbolic procedure int!-sqrt n;
   % Return sqrt of n if same is exact, or something non-numeric
   % otherwise.
    if not numberp n then 'nonnumeric
    else if n<0 then 'negative
    else if floatp n then sqrt n
    else if n<2 then n
    else int!-nr(n,(n+1)/2);


symbolic procedure int!-nr(n,root);
% root is an overestimate here. nr moves downwards to root;
 begin
    scalar w;
    w:=root*root;
    if n=w then return root;
    w:=(root+n/root)/2;
    if w>=root then return !*q2f simpsqrt list n;
    return int!-nr(n,w)
 end;


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


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


endmodule;

end;