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

module halfangl;  % Routines for conversion to half angle tangents.

% Author: Steve Harrington.
% Modifications by: John P. Fitch.

fluid '(!*gcd);

exports halfangle,untan;

symbolic procedure transform(u,x);
   % Transform the SQ U to remove the 'bad' functions sin, cos, cot etc.
   % in favor of half angles.
   % Do this with regard to cases like sin(x)*tan(x/2), so attempt to
   % limit times we use half angles.
   begin scalar zl,tnarg,substlist;
      zl := zlist;
      while car(tnarg := tan!-function!-in zl)
         and halfangle!-confusion(zlist,cadar tnarg)
       do <<zl := cdr tnarg;
            tnarg := car tnarg;
            if eqcar(tnarg,'tan)
              then substlist := (gensym() . tnarg) . substlist
             else substlist := (gensym() .
                                  list('quotient,1,('tan . cdr tnarg)))
                                  . substlist;
            u := subst(caar substlist,tnarg,u)>>;
      return if substlist
               then simp sublis(substlist,prepsq halfangle(u,x))
	      % simp prepsq was added so that 1/(e**x*cos(1/e**x)**2)
	      % for example returns a closed-form result.
	      else simp prepsq halfangle(u,x)
   end;

symbolic procedure tan!-function!-in zz;
   % Look at zlist for tangents or cotangents.
   <<while zz and not eqcar(car zz,'tan) and not eqcar(car zz,'cot)
	do zz := cdr zz;
     if null zz then nil . nil else zz>>;

symbolic procedure halfangle!-confusion(zz,tnarg);
   % Is there a function in the zlist with twice the tangent argument?
   <<while zz and (atom car zz
	      or not(tnarg = prepsq simp list('quotient,cadar zz,2)))
	do zz := cdr zz;
     zz>>;

symbolic procedure quotqq(u1,v1);
   multsq(u1,invsq(v1));

symbolic procedure !*subtrq(u1,v1);
   addsq(u1, negsq(v1));

symbolic procedure !*int2qm(u1);
   if u1=0 then nil . 1 else u1 . 1;

symbolic procedure halfangle(r,x);
   % Top level procedure for converting;
   % R is a rational expression to be converted,
   % X the integration variable.
   % A rational expression is returned.
   quotqq(hfaglf(numr(r),x), hfaglf(denr(r),x));

symbolic procedure hfaglf(p,x);
   % Converting polynomials,  a rational expression is returned.
   if domainp(p) then !*f2q(p)
    else subs2q addsq(multsq(exptsq(hfaglk(mvar(p),x), ldeg(p)),
                             hfaglf(lc(p),x)),
                      hfaglf(red(p),x));

symbolic procedure hfaglk(k,x);
   % Converting kernels,  a rational expression is returned.
   begin
      scalar kt;
      if atom k or not member(x,flatten(cdr(k))) then return !*k2q k;
      k := car(k) . hfaglargs(cdr(k), x);
      if cadr k eq 'pi then return !*k2q k;  % Don't consider tan(pi/2).
      kt := simp list('tan, list('quotient, cadr(k), 2));
      return if car(k) = 'sin
       then quotqq(multsq(!*int2qm(2),kt), addsq(!*int2qm(1),
                            exptsq(kt,2)))
      else if car(k) = 'cos
       then quotqq(!*subtrq(!*int2qm(1),exptsq(kt,2)),addsq(!*int2qm(1),
         exptsq(kt,2)))
      else if car(k) = 'tan
       then quotqq(multsq(!*int2qm(2),kt), !*subtrq(!*int2qm(1),
                            exptsq(kt,2)))
      else if car(k) = 'cot
       then quotqq(!*subtrq(!*int2qm(1),
                   exptsq(kt,2)),multsq(!*int2qm(2),kt))
      else if car(k) = 'sec
       then quotqq(addsq(!*int2qm(1), exptsq(kt,2)),
                   !*subtrq(!*int2qm(1),exptsq(kt,2)))
      else if car(k) = 'csc
       then quotqq(addsq(!*int2qm(1),exptsq(kt,2)),
%%%                   !*subtrq(!*int2qm(1),exptsq(kt,2)))
		    % FJW - was identical to sec!!!
		   multsq(!*int2qm(2),kt))
      else if car(k) = 'sinh then
        quotqq(!*subtrq(!*p2q mksp('expt.('e. cdr k),2),
        !*int2qm(1)), multsq(!*int2qm(2),
                             !*p2q mksp('expt . ('e . cdr(k)),1)))
      else if car(k) = 'cosh then
        quotqq(addsq(!*p2q mksp('expt.('e. cdr k),2),
        !*int2qm(1)), multsq(!*int2qm(2),
                             !*p2q mksp('expt . ('e . cdr(k)),1)))
      else if car(k) = 'tanh then
        quotqq(!*subtrq(!*p2q mksp('expt.('e. cdr k),2),
        !*int2qm(1)), addsq(!*p2q mksp ('expt.('e.cdr(k)),2),
        !*int2qm(1)))
      else if car(k) = 'coth then
        quotqq(addsq(!*p2q mksp('expt.('e.cdr(k)),2), !*int2qm(1)),
              !*subtrq(!*p2q mksp('expt.('e . cdr k),2),!*int2qm(1)))
      else if car(k) = 'acot then
        !*p2q mksp(list('atan, list('quotient, 1, cadr k)),1)
      else !*k2q(k);  % additional transformation might be added here.
   end;

symbolic procedure hfaglargs(l,x);
   % Conversion of argument list.
   if null l then nil
    else prepsq(hfaglk(car(l),x)) . hfaglargs(cdr(l),x);

symbolic procedure untanf x; 
   % This should be done by a table.
   % We turn off gcd to avoid unnecessary gcd calculations, as suggested
   % by Rainer Schoepf.
   begin scalar !*gcd,y,z,w;
      if domainp x then return x . 1; 
      y := mvar x; 
      if eqcar(y,'int) then error1();  % assume all is hopeless.
      z := ldeg x; 
      w := 1 . 1; 
      y := 
       if atom y then !*k2q y
        else if car y eq 'tan
         then begin scalar yy;
%%              printc "Recursive tan"; printc cadr y;
		yy := prepsq untan simp cadr y . nil;
%%              princ "==> "; printc yy;
                if evenp z
                then <<z := z/2; 
                       return simp list('quotient,
                                        list('plus,
                                             list('minus,
                                                  list('cos,
                                                       'times
                                                         . (2 . yy))),
                                             1),list('plus,
                                                     list('cos,
                                                          'times
                                                            . (2 . yy)),
                                              1))>>
               else if z=1
                then return simp list('quotient,
                                      list('plus,
                                           list('minus,
                                                list('cos,
                                                    'times . (2 . yy))),
                                           1),list('sin,
                                                   'times . (2 . yy)))
               else <<z := (z - 1)/2; 
                      w := 
                       simp list('quotient,
                                 list('plus,
                                      list('minus,
                                           list('cos,
                                                'times
                                                  . (2 . yy))),
                                      1),list('sin,
                                              'times
                                                . (2 . yy))); 
                      return simp list('quotient,
                                       list('plus,
                                            list('minus,
                                                 list('cos,
                                                      'times
                                                        . (2 . yy))),
                                            1),list('plus,
                                                    list('cos,
                                                         'times
                                                           . (2 . yy)),
                                                    1)) >>
        end
        else simp y;
      return addsq(multsq(multsq(exptsq(y,z),untanf lc x),w),
                   untanf red x)
   end;

% symbolic procedure untanlist(y);
%    if null y then nil
%      else (prepsq (untan(simp car y)) . untanlist(cdr y));

symbolic procedure untan(x);
   % Expects x to be canonical quotient.
   begin scalar y;
      y:=cossqchk sinsqrdchk multsq(untanf(numr x),
                                    invsq untanf(denr x));
      return if length flatten y>length flatten x then x else y
   end;

symbolic procedure sinsqrdchk(x);
   multsq(sinsqchkf(numr x), invsq sinsqchkf(denr x));

symbolic procedure sinsqchkf(x);
   begin
      scalar y,z,w;
      if domainp x then return x . 1;
      y := mvar x;
      z := ldeg x;
      w := 1 . 1;
      y := if eqcar(y,'sin) then if evenp z
       then <<z := quotient(z,2);
              simp list('plus,1,list('minus,
                                     list('expt,('cos . cdr(y)),2)))>>
      else if z = 1 then !*k2q y
      else  << z := quotient(difference(z,1),2); w := !*k2q y;
             simp list('plus,1,list('minus,
                                    list('expt,('cos . cdr(y)),2)))>>
       else !*k2q y;
      return addsq(multsq(multsq(exptsq(y,z),sinsqchkf(lc x)),w),
                   sinsqchkf(red x));
   end;

symbolic procedure cossqchkf(x);
   begin
      scalar y,z,w,x1,x2;
      if domainp x then return x . 1;
      y := mvar x;
      z := ldeg x;
      w := 1 . 1;
      x1 := cossqchkf(lc x);
      x2 := cossqchkf(red x);
      x := addsq(multsq(!*p2q lpow x,x1),x2);
      y := if eqcar(y,'cos) then if evenp z
       then <<z := quotient(z,2);
              simp list('plus,1,list('minus,
                                     list('expt,('sin . cdr(y)),2)))>>
      else if z = 1 then !*k2q y
      else  << z := quotient(difference(z,1),2); w := !*k2q y;
             simp list('plus,1,list('minus,
                                    list('expt,('sin . cdr(y)),2)))>>
       else !*k2q y;
      y := addsq(multsq(multsq(exptsq(y,z),w),x1),x2);
      return if length(y) > length(x) then x else y;
   end;

symbolic procedure cossqchk(x);
begin scalar !*gcd;
   !*gcd := t;
   return multsq(cossqchkf(numr x), invsq cossqchkf(denr x))
end;

symbolic procedure lrootchk(l,x);
   % Checks each member of list l for a root.
   if null l then nil else krootchk(car l, x) or lrootchk(cdr l, x);

symbolic procedure krootchk(f,x);
   % Checks a kernel to see if it is a root.
   if atom f then nil
    else if car(f) = 'sqrt and member(x, flatten cdr f) then t
   else if car(f) = 'expt
        and not atom caddr(f)
        and caaddr(f) = 'quotient
        and member(x, flatten cadr f)  then t
   else lrootchk(cdr f, x);

symbolic procedure rootchk1p(f,x);
   % Checks polynomial for a root.
   if domainp f then nil
    else krootchk(mvar f,x) or rootchk1p(lc f,x) or rootchk1p(red f,x);

symbolic procedure rootcheckp(f,x);
   % Checks rational (standard quotient) for a root.
   rootchk1p(numr f,x) or rootchk1p(denr f,x);

endmodule;

end;