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

module jpatches;   % Routines for manipulating sf's with power folding.

% Author: James H. Davenport.

fluid '(!*noncomp sqrtflag);

exports !*addsq,!*multsq,!*invsq,!*multf,!*exptsq,!*exptf; %squashsqrtsq

% !*MULTF(A,B) multiplies the polynomials (standard forms) U and V
% in much the same way as MULTF(U,V) would, EXCEPT...
%     (1) !*MULTF inhibits the action of OFF EXP and of non-commutative
%         multiplications
%     (2) Within !*MULTF powers of square roots, and powers of
%         exponential kernels are reduced as if substitution rules
%         such as FOR ALL X LET SQRT(X)**2=X were being applied;

% Note that !*MULTF comes between MULTF and !*Q2F SUBS2F MULTF in its
% behaviour, and that it is the responsibility of the user to call it
% in sensible places where its services are needed.

% Similarly for the other functions defined here.


%symbolic procedure !*addsq(u,v);
   %U and V are standard quotients.
%  %Value is canonical sum of U and V;
%   if null numr u then v
%    else if null numr v then u
%    else if denr u=1 and denr v=1 then addf(numr u,numr v) ./ 1
%    else begin scalar nu,du,nv,dv,x;
%        x := gcdf(du:=denr u,dv:=denr v);
%        du:=quotf(du,x); dv:=quotf(dv,x);
%        nu:=numr u; nv:=numr v;
%        u:=addf(!*multf(nu,dv),!*multf(nv,du));
%        if u=nil then return nil ./ 1;
%        v:=!*multf(du,denr v);
%        return !*ff2sq(u,v)
%    end;

%symbolic procedure !*multsq(a,b);
%  begin
%    scalar n,d;
%    n:=!*multf(numr a,numr b);
%    d:=!*multf(denr a,denr b);
%    return !*ff2sq(n,d)
%  end;

%symbolic procedure !*ff2sq(a,b);
%  begin
%    scalar gg;
%    if null a then return nil ./ 1;
%    gg:=gcdf(a,b);
%    if not (gg=1) then <<
%        a:=quotf(a,gg);
%        b:=quotf(b,gg) >>;
%    if minusf b then <<
%        a:=negf a;
%        b:=negf b >>;
%    return a ./ b
%  end;

symbolic procedure !*addsq(u,v);
   %U and V are standard quotients.
   %Value is canonical sum of U and V;
   if null numr u then v
    else if null numr v then u
    else if denr u=1 and denr v=1 then addf(numr u,numr v) ./ 1
    else begin scalar du,dv,x,y,z;
        x := gcdf(du:=denr u,dv:=denr v);
        du:=quotf(du,x); dv:=quotf(dv,x);
        y:=addf(!*multf(dv,numr u),!*multf(du,numr v));
        if null y then return nil ./ 1;
        z:=!*multf(denr u,dv);
        if minusf z then <<y := negf y; z := negf z>>;
        % In this case (as opposed to ADDSQ), Y and Z may have
        % developed common factors from SQRT expansion, so a
        % gcd of Y and Z is needed.
        x := gcdf(y,z);
        return if x=1 then y ./ z else quotf(y,x) ./ quotf(z,x)
    end;

symbolic procedure !*multsq(u,v);
   %U and V are standard quotients. Result is the canonical product of
   %U and V with surd powers suitably reduced.
   if null numr u or null numr v then nil ./ 1
    else if denr u=1 and denr v=1 then !*multf(numr u,numr v) ./ 1
    else begin scalar w,x,y;
     x := gcdf(numr u,denr v);
     y := gcdf(numr v,denr u);
     w := !*multf(quotf(numr u,x),quotf(numr v,y));
     x := !*multf(quotf(denr u,y),quotf(denr v,x));
     if minusf x then <<w := negf w; x := negf x>>;
     y := gcdf(w,x);  % another factor may have been generated.
     return if y=1 then w ./ x else quotf(w,y) ./ quotf(x,y)
    end;

symbolic procedure !*invsq a;
   % Note that several examples (e.g., int(1/(x**8+1),x)) give a more
   % compact result when SQRTFLAG is true if SQRT2TOP is not called.
   if sqrtflag then sqrt2top invsq a else invsq a;

symbolic procedure !*multf(u,v);
  % U and V are standard forms
  % Value is SF for U*V;
  begin scalar x,y;
      if null u or null v then return nil
      else if u=1 then return squashsqrt v
        else if v=1 then return squashsqrt u
        else if domainp u then return multd(u,squashsqrt v)
        else if domainp v then return multd(v,squashsqrt u)
        else if !*noncomp then return multf(u,v);
      x:=mvar u;
      y:=mvar v;
      if x eq y then go to c else if ordop(x,y) then go to b;
      x:=!*multf(u,lc v);
      y:=!*multf(u,red v);
      return if null x then y
              else if not domainp lc v
                and mvar u eq mvar lc v
                and not atom mvar u
                and car mvar u memq '(expt sqrt)
               then addf(!*multf(x,!*p2f lpow v),y)
              else makeupsf(lpow v,x,y);
  b:  x:=!*multf(lc u,v);
      y:=!*multf(red u,v);
      return if null x then y
              else if not domainp lc u
                and mvar lc u eq mvar v
                and not atom mvar v
                and car mvar v memq '(expt sqrt)
               then addf(!*multf(!*p2f lpow u,x),y)
              else makeupsf(lpow u,x,y);
  c:  y:=addf(!*multf(list lt u,red v),!*multf(red u,v));
      if eqcar(x,'sqrt)
        then return addf(squashsqrt y,!*multfsqrt(x,
                        !*multf(lc u,lc v),ldeg u + ldeg v))
        else if eqcar(x,'expt) and prefix!-rational!-numberp caddr x
          then return addf(squashsqrt y,!*multfexpt(x,
                        !*multf(lc u,lc v),ldeg u + ldeg v));
      x:=mkspm(x,ldeg u + ldeg v);
      return if null x or null (u:=!*multf(lc u,lc v))
               then y
               else addf(multpf(x,u),y)
    end;

symbolic procedure makeupsf(u,x,y);
% Makes u .* x .+ y  except when u is not a valid lpow (because of
% sqrts).
   if atom car u or cdr u = 1 then addf(multpf(u,x),y)
     else if caar u eq 'sqrt then addf(!*multfsqrt(car u,x,cdr u),y)
       else if <<begin scalar v;
                    v:=car u;
                    if car v neq 'expt then return nil;
                    v:=caddr v;
                    if atom v then return nil;
                    return (car v eq 'quotient 
                            and numberp cadr v
                            and numberp caddr v)
                 end >>
         then addf(!*multfexpt(car u,x,cdr u),y)
         else addf(multpf(u,x),y);


symbolic procedure !*multfsqrt(x,u,w);
   % This code (Due to Norman a& Davenport) squashes SQRT(...)**2.
   begin scalar v;
     w:=divide(w,2);
     v:=!*q2f simp cadr x;
     u:=!*multf(u,exptf(v,car w));
     if cdr w neq 0 then u:=!*multf(u,!*p2f mksp(x,1));
     return u
   end;


symbolic procedure !*multfexpt(x,u,w);
  begin scalar expon,v;
    expon:=caddr x;
    x:=cadr x;
    w:=w * cadr expon;
    expon:=caddr expon;
    v:=gcdn(w,expon);
    w:=w/v;
    v:=expon/v;
    if not (w > 0) then rerror(int,8,"Invalid exponent")
     else if v = 1
      then return !*multf(u,exptf(if numberp x then x
                                    else if atom x then !*k2f x
                                    else !*q2f if car x eq '!*sq
                                                 then argof x
                                                else simp x,
                          w));
    expon:=0;
    while not (w < v) do <<expon:=expon + 1; w:=w-v>>;
    if expon>0 then u:=!*multf(u,exptf(!*q2f simp x,expon));
    if w = 0 then return u;
    x:=list('expt,x,list('quotient,1,v));
    return multf(squashsqrt u,!*p2f mksp(x,w))   % Cannot be *MULTF.
  end;

symbolic procedure prefix!-rational!-numberp u;
  % Tests for m/n in prefix representation.
    eqcar(u,'quotient) and numberp cadr u and numberp caddr u;

% symbolic procedure squashsqrtsq sq;
%    !*multsq(squashsqrt numr sq ./ 1,1 ./ squashsqrt denr sq);

symbolic procedure squashsqrt sf;
if (not sqrtflag) or atom sf or atom mvar sf
  then sf
  else if car mvar sf eq 'sqrt and ldeg sf > 1
    then addf(squashsqrt red sf,!*multfsqrt(mvar sf,lc sf,ldeg sf))
    else if car mvar sf eq 'expt
       and prefix!-rational!-numberp caddr mvar sf
       and ldeg sf >= caddr caddr mvar sf
      then addf(squashsqrt red sf,!*multfexpt(mvar sf,lc sf,ldeg sf))
      else (if null x then squashsqrt red sf
            else lpow sf .* x .+ squashsqrt red sf)
           where x = squashsqrt lc sf;

%remd 'simpx1;

% The following definition requires frlis!* declared global.

%symbolic procedure simpx1(u,m,n);
%   %u,m and n are prefix expressions;
%   %value is the standard quotient expression for u**(m/n);
%   begin scalar flg,z;
%        if null frlis!* or null intersection(frlis!*,flatten (m . n))
%          then go to a;
%        exptp!* := t;
%        return !*k2q list('expt,u,if n=1 then m
%                                   else list('quotient,m,n));
%    a:  if numberp m and fixp m then go to e
%         else if atom m then go to b
%         else if car m eq 'minus then go to mns
%         else if car m eq 'plus then go to pls
%         else if car m eq 'times and numberp cadr m and fixp cadr m
%                and numberp n
%          then go to tms;
%    b:  z := 1;
%    c:  if atom u and not numberp u then flag(list u,'used!*);
%        u := list('expt,u,if n=1 then m else list('quotient,m,n));
%        if not(u member exptl!*) then exptl!* := u . exptl!*;
%    d:  return mksq(u,if flg then -z else z); %u is already in lowest
%%       %terms;
%    e:  if numberp n and fixp n then go to int;
%        z := m;
%        m := 1;
%        go to c;
%    mns: m := cadr m;
%        if !*mcd then return invsq simpx1(u,m,n);
%        flg := not flg;
%        go to a;
%    pls: z := 1 ./ 1;
%    pl1: m := cdr m;
%        if null m then return z;
%        z := multsq(simpexpt list(u,
%                        list('quotient,if flg then list('minus,car m)
%                                        else car m,n)),
%                    z);
%        go to pl1;
%    tms: z := gcdn(n,cadr m);
%        n := n/z;
%        z := cadr m/z;
%        m := retimes cddr m;
%        go to c;
%    int:z := divide(m,n);
%        if cdr z<0 then z:= (car z - 1) . (cdr z+n);
%        if 0 = cdr z
%          then return simpexpt list(u,car z);
%        if n = 2
%          then return multsq(simpexpt list(u,car z),
%                             simpsqrti u);
%        return multsq(simpexpt list(u,car z),
%                        mksq(list('expt,u,list('quotient,1,n)),cdr z))
%   end;

symbolic procedure !*exptsq(a,n);
% Raises A to the power N using !*MULTSQ.
    if n=0 then 1 ./ 1
    else if n=1 then a
    else if n<0 then !*exptsq(invsq a,-n)
    else begin
      scalar q,r;
      q:=divide(n,2);
      r:=cdr q; q:=car q;
      q:=!*exptsq(!*multsq(a,a),q);
      if r=0 then return q
      else return !*multsq(a,q)
    end;


symbolic procedure !*exptf(a,n);
% Raises A to the power N using !*MULTF.
    if n=0 then 1
    else if n=1 then a
    else begin
      scalar q,r;
      q:=divide(n,2);
      r:=cdr q; q:=car q;
      q:=!*exptf(!*multf(a,a),q);
      if r=0 then return q
      else return !*multf(a,q)
    end;

endmodule;

end;