File r38/packages/int/sqrtf.red artifact dd5c2b20bc part of check-in c70d02b470

module sqrtf;   % Square root of standard forms.

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

fluid '(!*noextend zlist);

exports nrootn,domainp,minusf;  % minusdfp,sqrtdf

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

% symbolic procedure minusdfp a;
%   % Test sign of leading coedd of d.f.
%     if null a then interr "Minusdfp 0 illegal"
%     else minusf numr lc a;

% symbolic procedure sqrtdf l;
%    % Takes square root of distributive form. "Failed" usually means
%    % that the square root is not among already existing objects.
%     if null l then nil
%     else begin scalar c;
%         if lpow l=vp2 zlist then go to ok;
%         c:=sqrtsq df2q l;
%         if numr c eq 'failed
%           then return 'failed;
%         if denr c eq 'failed
%           then return 'failed;
%         return for each u in f2df numr c
%                  collect (car u).!*multsq(cdr u,1 ./ denr c);
%     ok: c:=sqrtsq lc l;
%         if  numr c eq 'failed or
%             denr c eq 'failed
%           then return 'failed
%           else return (lpow l .* c) .+nil
%     end;

% symbolic procedure sqrtsq a;
%     sqrtf numr a ./ sqrtf denr a;

symbolic procedure sqrtf p;
    begin       scalar ip,qp;
        if null p then return nil;
        ip:=sqrtf1 p;
        qp:=cdr ip;
        ip:=car ip; %respectable and nasty parts of the sqrt.
        if qp=1 then return ip; %exact root found.
         if !*noextend then return 'failed;
            % We cannot add new square roots in this case, since it is
            % then impossible to determine if one square root depends
            % on another if new ones are being added all the time.
         if zlistp qp then return 'failed;
            % Liouville's theorem tells you that you never need to add
            % new algebraics depending on the variable of integration.
        qp:=simpsqrt2 qp;
        return !*multf(ip,qp)
    end;

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

symbolic procedure sqrtf1 p;
  % Returns a . b with p=a**2*b.
    if domainp p
     then if fixp p then nrootn(p,2)
           else !*q2f simpsqrt list prepf p . 1
    else begin scalar co,pp,g,pg;
        co:=contentsmv(p,mvar p,nil); %contents of p.
        pp:=quotf(p,co); %primitive part.
        co:=sqrtf1(co); %process contents via recursion.
        g:=gcdf(pp,partialdiff(pp,mvar pp));
        pg:=quotf(pp,g);
        g:=gcdf(g,pg); %a repeated factor of pp.
        if g=1 then pg:=1 . pp
        else <<
            pg:= quotf(pp,!*multf(g,g)); %what is still left.
            pg:=sqrtf1(pg); %split that up.
            rplaca(pg,!*multf(car pg,g))>>;
                 %put in the thing found here.
        rplaca(pg, !*multf(car pg,car co));
        rplacd(pg, !*multf(cdr pg,cdr co));
        return pg
    end;

% NROOTN removed as in REDUCE base.

endmodule;

end;