reduce-historical/r38/packages/taylor/tayutils.red

module TayUtils;

%*****************************************************************
%
%       Utility functions that operate on Taylor kernels
%
%*****************************************************************


exports add!-degrees, add!.comp!.tp!., check!-for!-cst!-Taylor,
        comp!.tp!.!-p, delete!-superfluous!-coeffs, enter!-sorted,
        exceeds!-order, exceeds!-order!-variant, find!-non!-zero,
        get!-cst!-coeff, inv!.tp!., is!-neg!-pl, min2!-order,
        mult!.comp!.tp!., rat!-kern!-pow, replace!-next,
        subtr!-degrees, subtr!-tp!-order, taydegree!<,
        taydegree!-strict!<!=, taymincoeff, tayminpowerlist,
        Taylor!*!-constantp, Taylor!*!-nzconstantp, Taylor!*!-onep,
        Taylor!*!-zerop, TayTpmin2, tp!-greaterp, truncate!-coefflist,
        truncate!-Taylor!*;

imports

% from the REDUCE kernel:
        ./, gcdn, geq, lastpair, mkrn, neq, nth, numr, reversip,

% from the header module:
        !*tay2q, get!-degree, get!-degreelist, make!-cst!-powerlist,
        Make!-Taylor!*, nzerolist, prune!-coefflist, TayCfPl, TayCfSq,
        TayCoeffList, TayGetCoeff, TayFlags, Taylor!:, TayOrig,
        TayTemplate, TayTpElNext, TayTpElOrder, TayTpElPoint,
        TayTpElVars, TpDegreeList, TpNextList,

% from module Tayintro:
        confusion,

% from module TayUtils:
        TayCoeffList!-zerop;



symbolic procedure add!-degrees (dl1, dl2);
  %
  % calculates the element-wise sum of two degree lists dl1, dl2.
  %
  Taylor!:
    begin scalar dl,u,v,w;
      while dl1 do <<
        u := car dl1;
        v := car dl2;
        w := nil;
        while u do <<
          w := (car u + car v) . w;
          u := cdr u;
          v := cdr v>>;
        dl := reversip w . dl;
        dl1 := cdr dl1;
        dl2 := cdr dl2>>;
      return reversip dl
    end;

symbolic procedure subtr!-degrees(dl1,dl2);
  %
  % calculates the element-wise difference of two degree lists dl1, dl2.
  %
  Taylor!:
    begin scalar dl,u,v,w;
      while dl1 do <<
        u := car dl1;
        v := car dl2;
        w := nil;
        while u do <<
          w := (car u - car v) . w;
          u := cdr u;
          v := cdr v>>;
      dl := reversip w . dl;
      dl1 := cdr dl1;
      dl2 := cdr dl2>>;
    return reversip dl
  end;

symbolic procedure find!-non!-zero pl;
  %
  % pl is a power list.  Returns pair (n . m) of position of first
  % non zero degree.
  %
  begin scalar u; integer n, m;
    n := 1;
   loop:
    m := 1;
    u := car pl;
   loop2:
    if not (car u = 0) then return (n . m);
    u := cdr u;
    m := m + 1;
    if not null u then goto loop2;
    pl := cdr pl;
    n := n + 1;
    if null pl then confusion 'find!-non!-zero;
    goto loop
  end$


symbolic procedure lcmn(n,m);
   %
   % returns the least common multiple of two integers m,n.
   %
   n*(m/gcdn(n,m));

symbolic smacro procedure get!-denom expo;
   if atom expo then 1 else cddr expo;

symbolic procedure get!-denom!-l expol;
   <<for each el in cdr expol do
       result := lcmn(result,get!-denom el);
     result>>
      where result := get!-denom car expol;

symbolic procedure get!-denom!-ll(dl,pl);
   if null dl then nil
    else lcmn(car dl,get!-denom!-l car pl)
           . get!-denom!-ll(cdr dl, cdr pl);

symbolic procedure smallest!-increment cfl;
   if null cfl then confusion 'smallest!-increment
    else begin scalar result;
     result := for each l in TayCfPl car cfl collect get!-denom!-l l;
     for each el in cdr cfl do
       result := get!-denom!-ll(result,TayCfPl el);
     return for each n in result collect if n=1 then n else mkrn(1,n);
   end;


symbolic procedure common!-increment(dl1,dl2);
   begin scalar result,l;
     loop:
       l := lcmn(get!-denom car dl1,get!-denom car dl2);
       result := (if l=1 then l else mkrn(1,l)) . result;
       dl1 := cdr dl1;
       dl2 := cdr dl2;
       if not null dl1 then goto loop
        else if not null dl2 then confusion 'common!-increment
        else return reversip result
   end;
     
symbolic procedure min2!-order(nextlis,ordlis,dl);
  %
  % (List of Integers, List of Integers, TayPowerList) -> Boolean
  %
  % nextlis is the list of TayTpElNext numbers,
  % ordlis the list if TayTpElOrder numbers,
  % dl the degreelist of a coefficient.
  % Dcecrease the TayTpElNext number if the degree is greater than
  % the order, but smaller than the next.
  % Returns the corrected nextlis.
  %
  if null nextlis then nil
   else (Taylor!: (if dg > car ordlis then min2(dg,car nextlis)
                   else car nextlis) where dg := get!-degree car dl)
          . min2!-order(cdr nextlis,cdr ordlis,cdr dl);

symbolic procedure replace!-next(tp,nl);
   %
   % Given a template and a list of exponents, returns a template
   %  with the next part replaced.
   %
   if null tp then nil
    else {TayTpElVars car tp,
          TayTpElPoint car tp,
          TayTpElOrder car tp,
          car nl}
           . replace!-next(cdr tp,cdr nl);

symbolic procedure comp!.tp!.!-p (u, v);
  %
  % Checks templates of Taylor kernels u and v for compatibility,
  % i.e. whether variables and expansion points match.
  % Returns t if possible.
  %
  begin;
    u := TayTemplate u; v := TayTemplate v;
    if length u neq length v then return nil;
   loop:
    if not (TayTpElVars car u = TayTpElVars car v
            and TayTpElPoint car u = TayTpElPoint car v)
      then return nil;
    u := cdr u; v := cdr v;
    if null u then return t;
    goto loop
  end$

symbolic procedure add!.comp!.tp!.(u,v);
  %
  % Checks templates of Taylor kernels u and v for compatibility
  % when adding them, i.e. whether variables and expansion points
  % match.
  % Returns either a list containing a new Taylor template whose
  % degrees are the minimum of the corresponding degrees of u and v,
  % or nil if variables or expansion point(s) do not match.
  %
  Taylor!:
  begin scalar w;
    u := TayTemplate u; v := TayTemplate v;
    if length u neq length v then return nil;
    if null u then return {nil};
   loop:
    if not (TayTpElVars car u = TayTpElVars car v
            and TayTpElPoint car u = TayTpElPoint car v)
      then return nil
     else w := {TayTpElVars car u,
                TayTpElPoint car u,
                min2(TayTpElOrder car u,TayTpElOrder car v),
                min2(TayTpElNext car u,TayTpElNext car v)}
                . w;
    u := cdr u; v := cdr v;
    if null u then return {reversip w};
    goto loop
  end;

symbolic procedure taymindegreel(pl,dl);
   Taylor!:
     if null pl then nil
      else min2(get!-degree car pl,car dl)
             . taymindegreel(cdr pl,cdr dl);

symbolic procedure get!-min!-degreelist cfl;
  Taylor!:
    if null cfl then confusion 'get!-min!-degreelist
     else if null cdr cfl then get!-degreelist TayCfPl car cfl
     else taymindegreel(TayCfPl car cfl,
                        get!-min!-degreelist cdr cfl);

symbolic procedure mult!.comp!.tp!.(u,v,div!?);
  %
  % Checks templates of Taylor kernels u and v for compatibility
  % when multiplying or dividing them, i.e., whether variables and
  % expansion points match.  The difference to addition is that
  % in addition to the new template it returns two degreelists
  % and two nextlists to be used by truncate!-coefflist which
  % are made up so that the kernels have the same number of terms.
  %
  Taylor!:
  begin scalar cf1,cf2,next1,next2,ord1,ord2,w,
        !#terms!-1,!#terms!-next,dl1,dl2,mindg;
    cf1 := prune!-coefflist TayCoeffList u;
    if null cf1 then dl1 := nzerolist length TayTemplate u
     else dl1 := get!-min!-degreelist cf1;
    cf2 := prune!-coefflist TayCoeffList v;
    if null cf2 then dl2 := nzerolist length TayTemplate v
     else dl2 := get!-min!-degreelist cf2;
    u := TayTemplate u; v := TayTemplate v;
    if length u neq length v then return nil;
    if null u then return {nil,nil,nil,nil,nil};
   loop:
    if not (TayTpElVars car u = TayTpElVars car v
            and TayTpElPoint car u = TayTpElPoint car v)
      then return nil;
    mindg := if div!? then car dl1 - car dl2 else car dl1 + car dl2;
    !#terms!-1 := min2(TayTpElOrder car u - car dl1,
                       TayTpElOrder car v - car dl2);
    !#terms!-next := min2(TayTpElNext car u - car dl1,
                          TayTpElNext car v - car dl2);
    ord1 := (!#terms!-1 + car dl1) . ord1;
    ord2 := (!#terms!-1 + car dl2) . ord2;
    next1 := (!#terms!-next + car dl1) . next1;
    next2 := (!#terms!-next + car dl2) . next2;
    w := {TayTpElVars car u,TayTpElPoint car u,
          mindg + !#terms!-1,mindg + !#terms!-next}
          . w;
    u := cdr u; v := cdr v; dl1 := cdr dl1; dl2 := cdr dl2;
    if null u then return {reversip w,
                           reversip ord1,
                           reversip ord2,
                           reversip next1,
                           reversip next2};
    goto loop
  end;

symbolic procedure inv!.tp!. u;
  %
  % Checks template of Taylor kernel u for inversion. It returns a
  % template (to be used by truncate!-coefflist)
  % which is made up so that the resulting kernel has the correct
  % number of terms.
  %
  Taylor!:
  begin scalar w,cf,!#terms!-1,!#terms!-next,dl,mindg;
    cf := prune!-coefflist TayCoeffList u;
    if null cf then dl := nzerolist length TayTemplate u
     else dl := get!-degreelist TayCfPl car cf;
    u := TayTemplate u;
    if null u then return {nil,nil};
   loop:
    mindg := - car dl;
    !#terms!-1 := TayTpElOrder car u - car dl;
    !#terms!-next := TayTpElNext car u - car dl;
    w := {TayTpElVars car u,TayTpElPoint car u,mindg + !#terms!-1,
          mindg + !#terms!-next}
          . w;
    u := cdr u; dl := cdr dl;
    if null u then return {reversip w};
    goto loop
  end;

symbolic smacro procedure taycoeff!-before(cc1,cc2);
  %
  % (TayCoeff, TayCoeff) -> Boolean
  %
  % returns t if coeff cc1 is ordered before cc2
  % both are of the form (degreelist . sq)
  %
  taydegree!<(TayCfPl cc1,TayCfPl cc2);

symbolic procedure taydegree!<(u,v);
  %
  % (TayPowerList, TayPowerList) -> Boolean
  %
  % returns t if coefflist u is ordered before v
  %
  Taylor!:
  begin scalar u1,v1;
   loop:
    u1 := car u;
    v1 := car v;
   loop2:
     if car u1 > car v1 then return nil
      else if car u1 < car v1 then return t;
     u1 := cdr u1;
     v1 := cdr v1;
     if not null u1 then go to loop2;
    u := cdr u;
    v := cdr v;
    if not null u then go to loop
  end;

symbolic procedure taydegree!-strict!<!=(u,v);
  %
  % (TayPowerList, TayPowerList) -> Boolean
  %
  % returns t if every component coefflist u is less or equal than
  %  respective component of v
  %
  Taylor!:
  begin scalar u1,v1;
   loop:
    u1 := car u;
    v1 := car v;
   loop2:
     if car u1 > car v1 then return nil;
     u1 := cdr u1;
     v1 := cdr v1;
     if not null u1 then go to loop2;
    u := cdr u;
    v := cdr v;
    if not null u then go to loop;
    return t
  end;

symbolic procedure exceeds!-order(ordlis,cf);
  %
  % (List of Integers, TayPowerlist) -> Boolean
  %
  % Returns t if the degrees in coefficient cf are greater or
  % equal than those in the degreelist ordlis
  %
  if null ordlis then nil
   else Taylor!:(get!-degree car cf >= car ordlis)
          or exceeds!-order(cdr ordlis,cdr cf);

symbolic procedure exceeds!-order!-variant(ordlis,cf);
  %
  % (List of Integers, TayPowerlist) -> Boolean
  %
  % Returns t if the degrees in coefficient cf are greater or
  % equal than those in the degreelist ordlis
  %
  if null ordlis then nil
   else Taylor!:(get!-degree car cf > car ordlis)
          or exceeds!-order!-variant(cdr ordlis,cdr cf);

symbolic procedure enter!-sorted (u, alist);
  %
  % (TayCoeff, TayCoeffList) -> TayCoeffList
  %
  % enters u into the alist alist according to the standard
  % ordering for the car part
  %
  if null alist then {u}
   else if taycoeff!-before (u, car alist) then u . alist
   else car alist . enter!-sorted (u, cdr alist)$

symbolic procedure delete!-superfluous!-coeffs(cflis,pos,n);
  %
  % (TayCoeffList, Integer, Integer) -> TayCoeffList
  %
  % This procedure deletes all coefficients of a TayCoeffList cflis
  % whose degree in position pos exceeds n.
  %
  Taylor!:
  for each cc in cflis join
     (if get!-degree nth(TayCfPl cc,pos) > n then nil else {cc});

symbolic procedure truncate!-coefflist (cflis, dl);
  %
  % (TayCoeffList, List of Integers) -> TayCoeffList
  %
  % Deletes all coeffs from coefflist cflis that are equal or greater
  %  in degree than the corresponding degree in the degreelist dl.
  %
  begin scalar l;
    for each cf in cflis do
      if not exceeds!-order (dl, TayCfPl cf) then l := cf . l;
    return reversip l
  end;

symbolic procedure TayTp!-min2(tp1,tp2);
  %
  % finds minimum (w.r.t. Order and Next parts) of compatible
  %  templates tp1 and tp2
  %
  Taylor!:
    if null tp1 then nil
     else if not (TayTpElVars car tp1 = TayTpElVars car tp2 and
                  TayTpElPoint car tp1 = TayTpElPoint car tp2)
      then confusion 'TayTpmin2
     else {TayTpElVars car tp1,TayTpElPoint car tp2,
           min2(TayTpElOrder car tp1,TayTpElOrder car tp2),
           min2(TayTpElNext car tp1,TayTpElNext car tp2)}
          . TayTp!-min2(cdr tp1,cdr tp2);


symbolic procedure truncate!-Taylor!*(tay,ntp);
  %
  % tcl is a coefflist for template otp
  % truncate it to coefflist for template ntp
  %
  Taylor!:
   begin scalar nl,ol,l,tp,tcl,otp;
     tcl := TayCoeffList tay;
     otp := TayTemplate tay;
     tp := for each pp in pair(ntp,otp) collect
           {TayTpElVars car pp,
            TayTpElPoint car pp,
            min2(TayTpElOrder car pp,TayTpElOrder cdr pp),
            min2(TayTpElNext car pp,TayTpElNext cdr pp)};
     nl := TpNextList tp;
     ol := TpDegreeList tp;
     for each cf in tcl do
       if not null numr TayCfSq cf
%         then ((if not exceeds!-order(nl,pl) then l := cf . l
%                 else nl := min2!-order(nl,ol,pl))
         then ((if not exceeds!-order!-variant(ol,pl) then l := cf . l
                 else nl := min2!-order(nl,ol,pl))
                where pl := TayCfPl cf);
     return make!-Taylor!*(reversip l,replace!-next(tp,nl),
                           TayOrig tay,TayFlags tay)
   end;

symbolic procedure tp!-greaterp(tp1,tp2);
   %
   % Given two templates tp1 and tp2 with matching variables and
   %  expansion points this function returns t if the expansion
   %  order wrt at least one variable is greater in tp1 than in tp2.
   %
   if null tp1 then nil
    else Taylor!: (TayTpElOrder car tp1 > TayTpElOrder car tp2)
           or tp!-greaterp(cdr tp1,cdr tp2);

symbolic procedure subtr!-tp!-order(tp1,tp2);
   %
   % Given two templates tp1 and tp2 with matching variables and
   %  expansion points this function returns the difference in their
   %  orders.
   %
   Taylor!:
    if null tp1 then nil
     else (TayTpElOrder car tp1 - TayTpElOrder car tp2)
            . subtr!-tp!-order(cdr tp1,cdr tp2);


comment Procedures to non-destructively modify Taylor templates;

symbolic procedure addto!-all!-TayTpElOrders(tp,nl);
  Taylor!:
   if null tp then nil
    else {TayTpElVars car tp,
          TayTpElPoint car tp,
          TayTpElOrder car tp + car nl,
          TayTpElNext car tp + car nl} .
         addto!-all!-TayTpElOrders(cdr tp,cdr nl);

symbolic procedure taymincoeff cflis;
  %
  % Returns degree of first non-zero coefficient
  % or 0 if there isn't any.
  %
  if null cflis then 0
   else if not null numr TayCfSq car cflis
    then get!-degree car TayCfPl car cflis
   else taymincoeff cdr cflis;

symbolic procedure tayminpowerlist cflis;
  %
  % Returns degreelist of first non-zero coefficient of TayCoeffList
  % cflis or a list of zeroes if there isn't any.
  %
  if null cflis then confusion 'tayminpowerlist
   else tayminpowerlist1(cflis,length TayCfPl car cflis);

symbolic procedure tayminpowerlist1(cflis,l);
   if null cflis then nzerolist l
    else if null numr TayCfSq car cflis
     then tayminpowerlist1(cdr cflis,l)
    else get!-degreelist TayCfPl car cflis;

symbolic procedure get!-cst!-coeff tay;
   TayGetCoeff(make!-cst!-powerlist TayTemplate tay,TayCoeffList tay);

symbolic procedure Taylor!*!-constantp tay;
   Taylor!*!-constantp1(make!-cst!-powerlist TayTemplate tay,
                        TayCoeffList tay);

symbolic procedure Taylor!*!-constantp1(pl,tcf);
   if null tcf then t
    else if TayCfPl car tcf = pl
     then TayCoeffList!-zerop cdr tcf
    else if not null numr TayCfSq car tcf then nil
    else Taylor!*!-constantp1(pl,cdr tcf);

symbolic procedure check!-for!-cst!-Taylor tay;
   begin scalar pl,tc;
      pl := make!-cst!-powerlist TayTemplate tay;
      tc := TayCoeffList tay;
      return if Taylor!*!-constantp1(pl,tc) then TayGetCoeff(pl,tc)
              else !*tay2q tay
   end;

symbolic procedure Taylor!*!-nzconstantp tay;
   Taylor!*!-nzconstantp1(make!-cst!-powerlist TayTemplate tay,
                          TayCoeffList tay);

symbolic procedure Taylor!*!-nzconstantp1(pl,tcf);
   if null tcf then nil
    else if TayCfPl car tcf = pl
     then if null numr TayCfSq car tcf then nil
           else TayCoeffList!-zerop cdr tcf
    else if TayCfPl car tcf neq pl and
            not null numr TayCfSq car tcf
     then nil
    else Taylor!*!-nzconstantp1(pl,cdr tcf);

symbolic procedure Taylor!*!-onep tay;
   Taylor!-onep1(make!-cst!-powerlist TayTemplate tay,TayCoeffList tay);

symbolic procedure Taylor!-onep1(pl,tcf);
   if null tcf then nil
    else if TayCfPl car tcf = pl
     then if TayCfSq car tcf = (1 ./ 1)
            then TayCoeffList!-zerop cdr tcf
           else nil
    else if null numr TayCfSq car tcf
     then Taylor!*!-nzconstantp1(pl,cdr tcf)
    else nil;

symbolic procedure is!-neg!-pl pl;
  %
  % Returns t if any of the exponents in pl is negative.
  %
  Taylor!:
    if null pl then nil
     else if get!-degree car pl < 0 then t
     else is!-neg!-pl cdr pl;

symbolic procedure rat!-kern!-pow(x,pos);
   %
   % check that s.f. x is a kernel to a rational power.
   %  if pos is t allow only positive exponents.
   %  returns pair (kernel . power)
   %
   begin scalar y; integer n;
     if domainp x or not null red x or not (lc x=1) then return nil;
     n := ldeg x;
     x := mvar x;
     Taylor!:
       if eqcar(x,'sqrt) then return (cadr x . mkrn(1,2)*n)
        else if eqcar(x,'expt) and (y := simprn{caddr x})
               then if null pos or (y := car y)>0
                      then return (cadr x . (y*n))
                     else return nil
        else return (x . n)
   end;

endmodule;

end;