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

module TaySimp;

%*****************************************************************
%
%     The special Taylor simplification functions
%
%*****************************************************************


exports taysimpp, taysimpsq, taysimpsq!*, expttayrat, expttayrat1;

imports

% from the REDUCE kernel:
        !*f2q, !*k2q, !*p2f, !*p2q, !*t2q, addsq, apply1, denr,
        domainp, evenp, exptsq, invsq, kernp, mk!*sq, mkrn, multf,
        multpq, multsq, mvar, nth, numr, over, pdeg, prepsq, quotsq,
        reversip, sfp, simp, simp!*, tc, to, tpow,

% from the header module:
        !*q2TayExp, !*tay2f, !*tay2q, !*tayexp2q, comp!.tp!.!-p,
        cst!-taylor!*, has!-Taylor!*, find!-non!-zero,
        get!-degreelist, has!-TayVars, invert!-powerlist,
        make!-cst!-coefflis, make!-cst!-powerlist, make!-Taylor!*,
        prune!-coefflist, resimptaylor, TayCfPl, TayCfSq,
        TayCoeffList, TayFlags, TayGetCoeff, Taylor!-kernel!-sq!-p,
        Taylor!*p, Taylor!:, TayMakeCoeff, taymultcoeffs, TayOrig,
        TayTemplate, TayTpElNext, TayTpElPoint, TayTpElVars,
        TpNextList,

% from module Tayintro:
        confusion, Taylor!-error, Taylor!-error!*,

% from module Tayutils:
        addto!-all!-TayTpElOrders, get!-cst!-coeff, smallest!-increment,
        Taylor!*!-nzconstantp, Taylor!*!-zerop,

% from module Tayinterf:
        taylorexpand, taylorexpand!-sf,

% from module Taybasic:
        addtaylor, addtaylor!-as!-sq, invtaylor, makecoeffpairs,
        makecoeffs0, multtaylor, multtaylor!-as!-sq, multtaylorsq,
        quottaylor!-as!-sq;


fluid '(!*taylorautoexpand !*taylorkeeporiginal)$

comment The procedures in this module provide the higher taylor
        manipulation machinery.  Given any s.q. (s.f.,...) they
        return the equivalent Taylor kernel (disguised as a s.q.)
        if the argument contains a Taylor kernel and everything
        else may be Taylor expanded.
        Otherwise the Taylor kernels in the argument are only
        partially combined (but as far as possible);


symbolic procedure taysimpsq u;
  %
  % The argument u is any standard quotient.
  % numerator and denominator are first simplified independently,
  % then the quotient is built.
  % We have four possible cases here, as both expressions
  %  may or may not be Taylor kernels.
  %
  begin scalar nm,dd;
    dd := taysimpf denr u;
    if null numr dd then Taylor!-error('zero!-denom,'taysimpsq)
     else if Taylor!-kernel!-sq!-p dd
      then return taysimpf multf(numr u,!*tay2f invtaylor mvar numr dd);
    nm := taysimpf numr u;
    return
    if Taylor!-kernel!-sq!-p nm
      then if not has!-TayVars(mvar numr nm,dd)
             then !*tay2q resimptaylor
                            multtaylorsq(mvar numr nm,invsq dd)
            else if Taylor!*!-nzconstantp mvar numr nm
             then quotsq(get!-cst!-coeff mvar numr nm,dd)
            else if null !*taylorautoexpand or has!-Taylor!* dd
             then quotsq(nm,dd)
            else taysimpsq!*
                   quottaylor!-as!-sq(
                     nm,
                     taylorexpand(dd,TayTemplate mvar numr nm))
     else quotsq(nm,dd)
  end;

symbolic procedure taysimpsq!* u;
   %
   % Variant of taysimpsq that does not automatically expand
   %  non-Taylor expressions
   %
   taysimpsq u where !*taylorautoexpand := nil;


symbolic procedure taysimpf u;
  %
  % u is a standard form which may contain Taylor subexpressions;
  % value is a standard form
  %
  begin scalar tay,notay,x,flg;
    %
    % Remember the definition of a s.f.:
    %  it is either a domain element,
    %  or it's car is a standard term and it's cdr is a s.f.
    %
    notay := nil ./ 1;
    while u do
      %
      % Split the constituents of the s.f. into non-Taylor and
      %  Taylor parts.  Taylor s.t.'s are simplified accordingly.
      % A domain element can never be a Taylor kernel.
      %
      <<if domainp u then notay := addsq(!*f2q u,notay)
         else if not has!-Taylor!* car u
          then notay := addsq(notay,!*t2q car u)
         else <<x := taysimpt car u;
                if Taylor!-kernel!-sq!-p x
                  then if null tay then tay := mvar numr x
                    else if comp!.tp!.!-p(tay,mvar numr x)
                           then tay := addtaylor(tay,mvar numr x)
                          else <<flg := t;
                                 %
                                 % flg indicates that there are
                                 %  Taylor kernels whose templates
                                 %  are not compatible
                                 %
                                 notay := addsq(notay,x)>>
                 else notay := addsq(notay,x)>>;
         u := if domainp u then nil else cdr u>>;
    %
    % tay is now a Taylor kernel or nil.
    %
    % We first make sure that it is not actually a constant.
    %
    if not null tay and not null TayOrig tay and null numr TayOrig tay
      then return notay
    %
    % If tay is nil, return the non-taylor parts.
    %
     else if null numr notay and not null tay then return !*tay2q tay
     else if null tay or Taylor!*!-zerop tay then return notay;
    %
    % Otherwise the non-taylor parts (if the are non-nil)
    % must be expanded if !*taylorautoexpand is non-nil.
    % The only exception are terms that do not contain
    % any of the Taylor variables: these are always expanded.
    %
    if Taylor!*!-nzconstantp tay and not has!-Taylor!* notay
      then return addsq(get!-cst!-coeff tay,notay)
     else if null !*taylorautoexpand and has!-TayVars(tay,notay)
      then return addsq(!*tay2q tay,notay);
    if flg then return addsq(!*tay2q tay,notay)
     else <<
       notay := taylorexpand(notay,TayTemplate tay);
       return taysimpsq!* addtaylor!-as!-sq(notay,!*tay2q tay)>>
  end;

symbolic procedure taysimpt u;
  %
  % u is a standard term containing one or more Taylor kernels,
  % value is the simplified Taylor expression (also as a s.f.).
  %
  begin scalar rest,pow;
    %
    % Since the coefficient of a term is a s.f.
    % we call taysimpf on it.
    %
    rest := taysimpf tc u;
    if null numr rest then return rest;
    pow := tpow u;
    %
    % Then we have to distinguish three cases:
    %   the case where no Taylor kernel appears was already caught
    %   by taysimpf before taysimpt was called.
    %
    % If combination of different Taylor kernels is impossible
    %   return them separately
    %
    % Remark: the call to SMEMQLP checks if rest contains one of
    %         the Taylor variables if it is not a Taylor kernel.
    %
    return if not has!-Taylor!* pow
      then if Taylor!-kernel!-sq!-p rest
             then multpowerintotaylor(pow,mvar numr rest)
            else multpq(pow,rest)
     else <<pow := taysimpp pow;
            if not has!-Taylor!* rest and Taylor!-kernel!-sq!-p pow
              then if has!-TayVars(mvar numr pow,rest)
                 then if !*taylorautoexpand
                   then taysimpsq!* multtaylor!-as!-sq(pow,
%                          taylorexpand(rest,TayTemplate mvar numr pow))
%
% the above is not entirely correct. the expansion should be done
% in a way so that the result of the multiplication has same order 
% as pow
%
                          taylorexpand!-sf(tc u,
                                           TayTemplate mvar numr pow,
                                           nil))
                  else multsq(pow,rest)
                else !*tay2q multtaylorsq(mvar numr pow,rest)
              else multtaylor!-as!-sq(pow,rest)>>
  end;

symbolic procedure multpowerintotaylor (p, tk);
  %
  % p is a standard power, tk a Taylor kernel
  % value is p expanded to a Taylor kernel multiplied by tk
  % this requires Taylor expansion of p if it contains
  % at least one of the expansion variables
  %
  % Remark: the call to SMEMQLP checks if p contains one of
  %         the Taylor variables.
  %
  if not has!-TayVars(tk,p)
    then !*tay2q multtaylorsq(tk,!*p2q p)
   else if !*taylorautoexpand
%    then taysimpsq!*
%           multtaylor!-as!-sq(!*tay2q tk,
%                     taylorexpand(!*p2q p,TayTemplate tk))
%
% here the same comment as above applies
%
    then taysimpsq!*
           multtaylor!-as!-sq(!*tay2q tk,
                     taylorexpand!-sf(!*p2f p,TayTemplate tk,nil))
   else if Taylor!*!-nzconstantp tk
    then multpq(p,get!-cst!-coeff tk)
   else multpq(p,!*tay2q tk);


symbolic procedure taysimpp u;
  %
  % u is a standard power containing a Taylor expression,
  % value is the simplified Taylor expression, as a s.f.
  %
  % We begin by isolating base and exponent.
  % First we simplify them separately.
  % Remember that the exponent is always an integer,
  % base is a kernel.
  %
  % If the main variable of the power is a kernel made of one
  %  of the functions known to the Taylor simplifier, call
  %  the appropriate simplification function.
  %  (There is a user hook here!)
  %
  if null car u or null pdeg u then confusion 'taysimpp
   else if sfp car u then !*p2q u
%%%% taysimpsq exptsq(taysimpf car u,cdr u)
   else if not taylor!*p car u
    then ((if kernp x
             then if (x := mvar numr x) = car u then !*p2q u
                   else if cdr u=1 then !*k2q x
                   else taysimpp(x .** cdr u)
            else if cdr u=1 then x
            else taysimpsq exptsq(x,cdr u))
          where x := (taysimpmainvar car u))
   %
   % We know how to raise a Taylor series to a rational power:
   %  positive integer --> multiply
   %  negative integer --> multiply and invert
   % Zero exponent should not appear: should be already simplified
   %  to 1 by the standard simplifier
   %
   else if not fixp pdeg u or pdeg u = 0 then confusion 'taysimpp
   else if not null TayOrig car u and null numr TayOrig car u
    then (nil ./ 1)
   else !*tay2q
     if pdeg u = 1 then car u else expttayi(car u,cdr u)$


symbolic procedure taysimpmainvar u;
  if not sfp u then taysimpkernel u
   else !*f2q taysimpf u;


symbolic procedure taysimpkernel u;
  begin scalar fn, x;
    u := simp!* u;
    if not kernp u then return u
     else << x := mvar numr u;
             if atom x or Taylor!*p x then return u;
             fn := get (car x, 'taylorsimpfn);
             return if null fn then u
                     else apply1 (fn, x)>>
  end;


symbolic procedure expttayi(u,i);
  %
  % raise Taylor kernel u to integer power i
  % algorithm is a scheme that computes powers of two.
  %
  begin scalar v,flg;
    if i<0 then <<i := -i; flg := t>>;
    v := if evenp i then cst!-Taylor!*(1 ./ 1,TayTemplate u)
          else <<i := i - 1; u>>;
    while (i:=i/2)>0 do <<u := multtaylor(u,u);
                          if not evenp i then v := multtaylor(v,u)>>;
    return if flg then invtaylor v else v
  end;


comment non-integer powers of Taylor kernels;


comment The implementation of expttayrat follows the algorithm
        quoted by Knuth in the second volume of `The Art of
        Computer Programming', extended to the case of series in
        more than one variable.

        Assume you want to compute the series W(x) where

            W(x) = V(x)**alpha

        Differentiation of this equation gives

            W'(x) = alpha * V(x)**alpha * V'(x) .

        You make now the ansatz

                    -----
                    \           n
            W(x) =   >      W  x  ,
                    /        n
                    -----

        substitute this into the above equation and compare
        powers of x.  This yields the recursion formula

                       n-1
                      -----
                  1   \                  m      m
           W  = -----  >    (alpha (1 - ---) - --- ) W  V     .
            n    V    /                  n      n     m  n-m
                  0   -----
                       m=0

        The first coefficient must be calculated directly, it is

           W   = V  ** alpha .
            0     0

        To use this for series in more than one variable you have to
        calculate all partial derivatives: n and m refer then to the
        corresponding component of the multi index.  Looking closely
        one finds that there is an ambiguity: the same coefficient
        can be calculated using any of the partial derivatives.  The
        only restriction is that the corresponding component of the
        multi index must not be zero, since we have to divide by it.

        We resolve this ambiguity by simply taking the first nonzero
        component.

        We use it here only for the case that alpha is a rational
        number;


symbolic procedure expttayrat(tay,rat);
  %
  % tay is a Taylor kernel, rat is a s.q. of two integers
  % value is tay ** rat
  % algorithm as quoted by Knuth
  %
  Taylor!:
  begin scalar clist,tc,tp;
    %
    % First of all we have to find out if we can raise the leading
    %  term to the power rat.
    % If so we calculate the reciprocal of this leading coefficient
    %  and multiply all other terms with it.
    % This guarantees that the resulting Taylor kernel starts with
    %  coefficient 1.
    %
    if not Taylor!*p tay then return simp!* {'expt,tay,mk!*sq rat};
    tc := prune!-coefflist TayCoeffList tay;
    tp := TayTemplate tay;
    %
    % Find first non-zero coefficient.
    %
    if null tc
      then if minusp numr rat
             then Taylor!-error!*('not!-a!-unit,'expttayrat)
            else <<tp := for each tpel in tp collect begin scalar w;
                           w := TayTpElNext tpel * !*q2TayExp rat;
                           return {TayTpElVars tpel,
                                   TayTpElPoint tpel,
                                   w - mkrn(1,denr rat),
                                   w};
                           end;
                   clist := make!-cst!-coefflis(nil ./ 1,tp)>>
     else begin scalar c0,l,l1;
       c0 := car tc;
       l1 := for each ll in TayCfPl c0 collect
               for each p in ll collect (p * !*q2TayExp rat);
       l := invert!-powerlist TayCfPl c0;
       tp := addto!-all!-TayTpElOrders(tp,get!-degreelist l);
       l := TayMakeCoeff(l,invsq TayCfSq c0);
       %
       % We divide the rest of the kernel (without the leading term)
       %  by the leading term.
       %
       l := for each el in cdr tc collect taymultcoeffs(el,l);
       clist := expttayrat1(tp,l,rat);
       %
       % Next we multiply the resulting Taylor kernel by the leading
       %  coefficient raised to the power rat.
       %
       c0 := TayMakeCoeff(l1,simp!* {'expt,mk!*sq TayCfSq c0,
                                     {'quotient,numr rat,denr rat}});
       clist := for each el in clist collect taymultcoeffs(el,c0);
       tp := addto!-all!-TayTpElOrders(tp,get!-degreelist l1);
    end;
    %
    % Finally we fill in the original expression
    %
    return make!-Taylor!*(
             clist,
             tp,
             if !*taylorkeeporiginal and TayOrig tay
               then simp {'expt,prepsq TayOrig tay,
                          {'quotient,car rat,cdr rat}}
              else nil,
             TayFlags tay)
  end;

symbolic procedure expttayrat1(tp,tcl,rat);
  Taylor!:
   begin scalar clist,coefflis,il,l0,rat1;
     rat1 := addsq(rat,1 ./ 1);
     %
     % Now we compute the coefficients
     %
     l0 := make!-cst!-powerlist tp;
     clist := {TayMakeCoeff(l0,1 ./ 1)};
     tcl := TayMakeCoeff(l0,1 ./ 1) . tcl;
     il := smallest!-increment tcl;
     coefflis := makecoeffs0(tp,TpNextList tp,il);
     if null coefflis then return clist;
     for each cc in cdr coefflis do
       begin scalar s,pos,pp,q,n,n1;
         s := nil ./ 1;
         pos := find!-non!-zero cc;
         n := nth(nth(cc,car pos),cdr pos);
         pp := makecoeffpairs(l0,cc,l0,il);
         for each p in pp do begin scalar v,w;
           v := TayGetCoeff(cdr p,tcl);
           w := TayGetCoeff(car p,clist);
           %
           % The following line is a short cut for efficiency.
           %
           if null numr v or null numr w then return;
           w := multsq(w,v);
           n1 := nth(nth(car p,car pos),cdr pos);
           q := quotsq(!*TayExp2q(-n1),!*TayExp2q n);
           s := addsq(s,multsq(addsq(rat,multsq(q,rat1)),w))
          end;
         if not null numr s then clist := TayMakeCoeff(cc,s) . clist
       end;
     return reversip clist
   end;

endmodule;

end;