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

module tidysqrt;  % General tidying up of square roots.

% Authors: Mary Ann Moore and Arthur C. Norman.
% Modifications by J.H. Davenport.

exports sqrt2top,tidysqrt;

%symbolic procedure tidysqrtdf a;
%    if null a then nil
%    else begin    scalar tt,r;
%        tt:=tidysqrt lc a;
%        r:=tidysqrtdf red a;
%        if null numr tt then return r;
%        return ((lpow a) .* tt) .+ r
%    end;
%
symbolic procedure tidysqrt q;
    begin    scalar n1,dd;
        n1:=tidysqrtf numr q;
        if null n1 then return nil ./ 1; %answer is zero.
        dd:=tidysqrtf denr q;
        return multsq(n1,invsq dd)
    end;

symbolic procedure tidysqrtf p;
%Input - standard form.
%Output - standard quotient.
%Simplifies sqrt(a)**n with n>1.
    if domainp p then p ./ 1
    else begin    scalar v,w;
        v:=lpow p;
        if car v='i then v:=mksp('(sqrt -1),cdr v); %I->sqrt(-1);
        if eqcar(car v,'sqrt) and not onep cdr v
          then begin    scalar x;
             %here we have a reduction to apply.
            x:=divide(cdr v,2); %halve exponent.
            w:=exptsq(simp cadar v,car x); %rational part of answer.
            if cdr x neq 0
              then w := multsq(w,((mksp(car v,1) .* 1) .+ nil) ./ 1);
            %the next line allows for the horrors of nested sqrts.
            w:=tidysqrt w
            end
        else w:=((v .* 1) .+ nil) ./ 1;
        v:=multsq(w,tidysqrtf lc p);
        return addsq(v,tidysqrtf red p)
    end;


symbolic procedure multoutdenr q;
  % Move sqrts in a sq to the numerator.
    begin  scalar n,d,root,conj;
        n:=numr q;
        d:=denr q;
        while (root:=findsquareroot d) do <<
          conj:=conjugatewrt(d,root);
          n:=!*multf(n,conj);
          d:=!*multf(d,conj) >>;
        while (root:=findnthroot d) do <<
          conj:=conjugateexpt(d,root,kord!*);
          n:=!*multf(n,conj);
          d:=!*multf(d,conj) >>;
        return (n . d);
        end;

symbolic procedure conjugateexpt(d,root,kord!*);
  begin scalar ord,ans,repl,xi;
  ord:=caddr caddr root; % the denominator of the exponent;
  ans:=1;
  kord!*:= (xi:=gensym()) . kord!*;
  % XI is an ORD'th root of unity;
  for i:=1:ord-1 do <<
    ans:=!*multf(ans,numr subf(d,
                   list(root . list('times,root,list('explt,xi,i)))));
    while (mvar ans eq xi) and ldeg ans > ord do
      ans:=addf(red ans,(xi) to (ldeg ans - ord) .* lc ans .+ nil);
    if (mvar ans eq xi) and ldeg ans = ord then
      ans:=addf(red ans,lc ans) >>;
  if (mvar ans eq xi) and ldeg ans = ord-1 then <<
    repl:=-1;
    for i:=1:ord-2 do
      repl:=(xi) to i .* -1 .+ repl;
    ans:=addf(red ans,!*multf(lc ans,repl)) >>;
  if not domainp ans and mvar ans eq xi
    then interr "Conjugation failure";
  return ans;
  end;

symbolic procedure sqrt2top q;
begin
  scalar n,d;
  n:=multoutdenr q;
  d:=denr n;
  n:=numr n;
  if d eq denr q
    then return q;%no change.
  if d iequal 1
    then return (n ./ 1);
  q:=gcdcoeffsofsqrts n;
  if q iequal 1
    then if minusf d
      then return (negf n ./ negf d)
      else return (n ./ d);
  q:=gcdf(q,d);
  n:=quotf(n,q);
  d:=quotf(d,q);
  if minusf d
    then return (negf n ./ negf d)
    else return (n ./ d)
    end;

%symbolic procedure denrsqrt2top q;
%begin
%  scalar n,d;
%  n:=multoutdenr q;
%  d:=denr n;
%  n:=numr n;
%  if d eq denr q
%    then return d; % no changes;
%  if d iequal 1
%    then return 1;
%  q:=gcdcoeffsofsqrts n;
%  if q iequal 1
%    then return d;
%  q:=gcdf(q,d);
%  if q iequal 1
%    then return d
%    else return quotf(d,q)
%  end;

symbolic procedure findsquareroot p;
  % Locate a sqrt symbol in poly p.
    if domainp p then nil
    else begin scalar w;
        w:=mvar p; %check main var first.
        if atom w
          then return nil; %we have passed all sqrts.
        if eqcar(w,'sqrt) then return w;
        w:=findsquareroot lc p;
        if null w then w:=findsquareroot red p;
        return w
    end;

symbolic procedure findnthroot p;
   nil;   % Until corrected.

% symbolic procedure x!-findnthroot p;
%     % Locate an n-th root symbol in poly p.
%     if domainp p then nil
%     else begin scalar w;
%         w:=mvar p; %check main var first.
%         if atom w
%           then return nil; %we have passed all sqrts.
%         if eqcar(w,'expt) and eqcar(caddr w,'quotient) then return w;
%         w:=findnthroot lc p;
%         if null w then w:=findnthroot red p;
%         return w
%     end;

symbolic procedure conjugatewrt(p,var);
  % Var -> -var in form p.
    if domainp p then p
    else if mvar p=var then begin
        scalar x,c,r;
        x:=tdeg lt p; %degree
        c:=lc p; %coefficient
        r:=red p; %reductum
        x:=remainder(x,2); %now just 0 or 1.
        if x=1 then c:=negf c; %-coefficient.
        return (lpow p .* c) .+ conjugatewrt(r,var) end
    else if ordop(var,mvar p) then p
    else (lpow p .* conjugatewrt(lc p,var)) .+
        conjugatewrt(red p,var);

symbolic procedure gcdcoeffsofsqrts u;
if atom u
  then if numberp u and minusp u
    then -u
    else u
  else if eqcar(mvar u,'sqrt)
    then begin
      scalar v;
      v:=gcdcoeffsofsqrts lc u;
      if v iequal 1
        then return v
        else return gcdf(v,gcdcoeffsofsqrts red u)
      end
    else begin
      scalar root;
      root:=findsquareroot u;
      if null root
        then return u;
      u:=makemainvar(u,root);
      root:=gcdcoeffsofsqrts lc u;
      if root iequal 1
        then return 1
        else return gcdf(root,gcdcoeffsofsqrts red u)
      end;

endmodule;

end;