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

module TayIntrf;

%*****************************************************************
%
%      The interface to the REDUCE simplificator
%
%*****************************************************************


exports simptaylor, simpTaylor!*, taylorexpand$

imports

% from the REDUCE kernel:
        !*f2q, aconc!*, denr, depends, diffsq, eqcar, kernp, lastpair,
        leq, lprim, mkquote, mksq, multsq, mvar, neq, nth, numr, over,
        prepsq, revlis, reversip, simp, simp!*, subs2, subsq, typerr,

% from the header module:
        !*tay2q, get!-degree, has!-Taylor!*, has!-TayVars,
        make!-Taylor!*, multintocoefflist, resimptaylor, TayCfPl,
        TayCfSq, TayCoeffList, TayFlags, TayMakeCoeff, TayOrig,
        TayTemplate, TayTpElOrder, TayTpElPoint,
        Taylor!-kernel!-sq!-p, taymincoeff,

% from module Tayintro:
        replace!-nth, Taylor!-error, var!-is!-nth,

% from module TayExpnd:
        taylorexpand,

% from module Tayutils:
        delete!-superfluous!-coeffs,

% from module taybasic:
        invtaylor1, quottaylor1,

% from module Tayconv:
        prepTaylor!*;


fluid '(!*backtrace !*precise !*taylorkeeporiginal !*taylorautocombine
        frlis!* subfg!*);

global '(mul!*);

comment The following statement forces all expressions to be
        re-simplified if the switch `taylorautocombine' is set to on,
        unfortunately, this is not always sufficient;

put ('taylorautocombine, 'simpfg, '((t (rmsubs))));


symbolic procedure simptaylor u;
  %
  % (PrefixForm) -> s.q.
  %
  % This procedure is called directly by the simplifier.
  % Its argument list must be of the form
  %     (exp, [var, var0, deg, ...]),
  % where [...] indicate one or more occurences.
  % This means that exp is to be expanded w.r.t var about var0
  % up to degree deg.
  % var may also be a list of variables, which means that the
  % the expansion takes place in a homogeneous way.
  % If var0 is the special atom infinity var is replaced by 1/var
  % and the result expanded about 0.
  %
  % This procedure returns the input if it cannot expand the expression.
  %
  if remainder(length u,3) neq 1
    then Taylor!-error('wrong!-no!-args,'taylor)
   else if null subfg!* then mksq('taylor . u,1)
   else begin scalar !*precise,arglist,degree,f,ll,result,var,var0;
     %
     % Allow automatic combination of Taylor kernels.
     %
     if !*taylorautocombine and not ('taysimpsq memq mul!*)
       then mul!* := aconc!*(mul!*,'taysimpsq);
     %
     % First of all, call the simplifier on exp (i.e. car u),
     %
     f := simp!* car u;
     u := revlis cdr u; % reval instead of simp!* to handle lists
     arglist := u;
     %
     % then scan the rest of the argument list.
     %
     while not null arglist do
       << var := car arglist;
          var := if eqcar(var,'list) then cdr var else {var};
          % In principle one should use !*a2k
          % but typerr (maprin) does not correctly print the atom nil
          for each el in var collect begin
            el := simp!* el;
            if kernp el then return mvar numr el
             else typerr(prepsq el,'kernel)
            end;
          var0 := cadr arglist;
          degree := caddr arglist;
          if not fixp degree
            then typerr(degree,"order of Taylor expansion");
          arglist := cdddr arglist;
          ll := {var,var0,degree,degree + 1} . ll>>;
     %
     % Now ll is a Taylor template, i.e. of the form
     %  ((var_1 var0_1 deg1 next_1) (var_2 var0_2 deg_2 next_2) ...)
     %
     result := taylorexpand(f,reversip ll);
     %
     % If the result does not contain a Taylor kernel, return the input.
     %
     return if has!-Taylor!* result then result
             else mksq('taylor . prepsq f . u,1)
   end;

put('taylor,'simpfn,'simptaylor)$


%symbolic procedure taylorexpand (f, ll);
%  %
%  % f is a s.q. that is expanded according to the list ll
%  %  which has the form
%  %  ((var_1 var0_1 deg1) (var_2 var0_2 deg_2) ...)
%  %
%  begin scalar result;
%    result := f;
%    for each el in ll do
%      %
%      % taylor1 is the function that does the real work
%      %
%      result := !*tay2q taylor1 (result, car el, cadr el, caddr el);
%      if !*taylorkeeporiginal then set!-TayOrig (mvar numr result, f);
%      return result
%  end$


symbolic procedure taylor1 (f, varlis, var0, n);
  %
  % Taylor expands s.q. f w.r.t. varlis about var0 up to degree n.
  % var is a list of kernels, which means that the expansion
  % takes place in a homogeneous way if there is more than one
  % kernel.
  % If var0 is the special atom infinity all kernels in varlis are
  % replaced by 1/kernel.  The result is then expanded about 0.
  %
  taylor1sq (if var0 eq 'infinity
               then subsq (f,
                           for each krnl in varlis collect
                             (krnl . list ('quotient, 1, krnl)))
              else f,
             varlis, var0, n)$

symbolic procedure taylor1sq (f, varlis, var0, n);
  %
  % f is a standard quotient, value is a Taylor kernel.
  %
  % First check if it is a Taylor kernel
  %
  if Taylor!-kernel!-sq!-p f
    then if has!-TayVars(mvar numr f,varlis)
           %
           % special case: f has already been expanded w.r.t. var.
           %
           then taylorsamevar (mvar numr f, varlis, var0, n)
          else begin scalar y, z;
            f := mvar numr f;
            %
            % taylor2 returns a list of the form
            %  ((deg1 . coeff1) (deg2 . coeff2) ... )
            % apply this to every coefficient in f.
            % car cc is the degree list of this coefficient,
            % cdr cc the coefficient itself.
            % Finally collect the new pairs
            %  (degreelist . coefficient)
            %
            z :=
              for each cc in TayCoeffList f join
                for each cc2 in taylor2 (TayCfSq cc, varlis, var0, n)
                  collect
                    TayMakeCoeff (append (TayCfPl cc, TayCfPl cc2),
                                  TayCfSq cc2);
            %
            % Append the new list to the Taylor template and return.
            %
            y := append(TayTemplate f,list {varlis,var0,n,n+1});
            return make!-Taylor!* (z, y, TayOrig f, TayFlags f)
            end
   %
   % Last possible case: f is not a Taylor expression.
   % Expand it.
   %
   else make!-Taylor!* (
                 taylor2 (f, varlis, var0, n),
                 list {varlis,var0,n,n+1},
                 if !*taylorkeeporiginal then f else nil,
                 nil)$

symbolic procedure taylor2 (f, varlis, var0, n);
  begin scalar result,oldklist;
    oldklist := get('Taylor!*,'klist);
    result := errorset (list ('taylor2e,
                               mkquote f,
                               mkquote varlis,
                               mkquote var0,
                               mkquote n),
                        nil, !*backtrace);
    put('Taylor!*,'klist,oldklist);
    if atom result
      then Taylor!-error ('expansion, "(possible singularity!)")
     else return car result
  end$

symbolic procedure taylor2e (f, varlis, var0, n);
  %
  % taylor2e expands s.q. f w.r.t. varlis about var0 up to degree n.
  % var is a list of kernels, which means that the expansion takes
  % place in a homogeneous way if there is more than one kernel.
  % The case that var0 is the special atom infinity has to be taken
  % care of by the calling function(s).
  % Expansion is carried out separately for numerator and
  % denominator.  This approach has the advantage of not producing
  % complicated s.q.'s which usually appear if an s.q. is
  % differentiated.  The procedure is (roughly) as follows:
  %  if the denominator of f is free of var
  %   then expand the numerator and divide,
  %  else if the numerator is free of var expand the denominator,
  %   take the reciprocal of the Taylor series and multiply,
  %  else expand both and divide the two series.
  % This fails if there are nontrivial dependencies, e.g.,
  %  if a variable is declared to depend on a kernel in varlis.
  % It is of course necessary that the denominator yields a unit
  %  in the ring of Taylor series. If not, an error will be
  %  signalled by invtaylor or quottaylor, resp.
  %
  if cdr varlis then taylor2hom (f, varlis, var0, n)
   else if denr f = 1 then taylor2f (numr f, car varlis, var0, n, t)
   else if not depends (denr f, car varlis)
    then multintocoefflist (taylor2f (numr f, car varlis, var0, n, t),
                            1 ./ denr f)
   else if numr f = 1
    then delete!-superfluous!-coeffs
           (invtaylor1 ({varlis,var0,n,n+1},
                        taylor2f (denr f, car varlis, var0, n, nil)),
            1, n)
   else if not depends (numr f, car varlis)
    then multintocoefflist
           (invtaylor1 ({varlis,var0,n,n+1},
                        taylor2f (denr f, car varlis, var0, n, nil)),
            !*f2q numr f)
   else begin scalar denom; integer n1;
     denom := taylor2f (denr f, car varlis, var0, n, nil);
     n1 := n + taymincoeff denom;
     return
       delete!-superfluous!-coeffs
         (quottaylor1 ({varlis,var0,n1,n1+1},
                       taylor2f (numr f, car varlis, var0, n1, t),
                       denom),
          1, n)
  end$

symbolic procedure taylor2f (f, var, var0, n, flg);
  %
  % taylor2f is the procedure that does the actual expansion
  % of the s.f. f.
  % It returns a list of the form
  %  ((deglis1 . coeff1) (deglis2 . coeff2) ... )
  % For the use of the parameter `flg' see below.
  %
  begin scalar x, y, z; integer k;
    %
    % Calculate once what is needed several times.
    % var0 eq 'infinity is a special case that has already been taken
    % care of by the calling functions by replacing var by 1/var.
    %
    if var0 eq 'infinity then var0 := 0;
    x := list (var . var0);
    y := simp list ('difference, var, var0);
    %
    % The following is a first attempt to treat expressions
    % that have poles at the expansion point.
    % Currently nothing more than factorizable poles, i.e.
    % factors in the denominator are handled.
    % We use the following trick to calculate enough terms: If the
    % first l coefficients of the Taylor series are zero we replace n
    % by n + 2l.  This is necessary since we separately expand
    % numerator and denominator of an expression.  If the expansion of
    % both parts starts with, say, the quadratic term we have to
    % expand them up to order (n+2) to get the correct result up to
    % order n. However, if the numerator starts with a constant term
    % instead, we have to expand up to order (n+4).  It is important,
    % however, to drop the additional coefficients as soon as they are
    % no longer needed.  The parameter `flg' is used here to control
    % this behaviour.  The calling function must pass the value `t' if
    % it wants to inhibit the calculation of additional coefficients.
    % This is currently the case if the s.q. f has a denominator that
    % may contain the expansion variable.  Otherwise `flg' is used to
    % remember if we already encountered a non-zero coefficient.
    %
    f := !*f2q f;
    z := subs2 subsq (f, x);
    if null numr z and not flg then n := n + 1 else flg := t;
    y := list TayMakeCoeff ((list list 0), z);
    k := 1;
    while k <= n and not null numr f do
      << f := multsq (diffsq (f, var), 1 ./ k);
                                             % k is always > 0!
         % subs2 added to simplify expressions involving roots
         z := subs2 subsq (f, x);
         if null numr z and not flg then n := n + 2 else flg := t;
         if not null numr z then y := TayMakeCoeff(list list k, z) . y;
         k := k + 1 >>;
    return reversip y
  end;

symbolic procedure taylor2hom (f, varlis, var0, n);
  %
  % Homogeneous expansion of s.q. f wrt the variables in varlis,
  % i.e. the sum of the degrees of the kernels is varlis is <= n
  %
  if null cdr varlis then taylor2e (f, list car varlis, var0, n)
   else for each u in taylor2e (f, list car varlis, var0, n) join
     for each v in taylor2hom (cdr u,
                               cdr varlis,
                               var0,
                               n - get!-degree TayCfPl car u)
           collect list (car TayCfPl car u . TayCfPl car v) . cdr v$

symbolic procedure taylorsamevar (u, varlis, var0, n);
  begin scalar tp; integer mdeg, pos;
    if cdr varlis
      then Taylor!-error ('not!-implemented,
                          "(homogeneous expansion in TAYLORSAMEVAR)");
    tp := TayTemplate u;
    pos := car var!-is!-nth (tp, car varlis);
    tp := nth (tp, pos);
    if TayTpElPoint tp neq var0
      then return taylor1 (if not null TayOrig u then TayOrig u
                            else simp!* prepTaylor!* u,
                           varlis, var0, n);
    mdeg := TayTpElOrder tp;
    if n=mdeg then return u
     else if n > mdeg
      %
      % further expansion required
      %
      then << lprim "Cannot expand further... truncated.";
              return u >> ;
    return make!-Taylor!* (
        for each cc in TayCoeffList u join
          if nth (nth (TayCfPl cc, pos), 1) > n
            then nil
           else list cc,
        replace!-nth(TayTemplate u,pos,
                      {varlis,TayTpElPoint tp,n,n+1}),
        TayOrig u, TayFlags u)
  end$


comment The `FULL' flag causes the whole term (including the
        Taylor!* symbol) to be passed to SIMPTAYLOR!* ;

symbolic procedure simpTaylor!* u;
  if TayCoefflist u memq frlis!* or eqcar(TayCoefflist u,'!~)
    then !*tay2q u
   else <<
     %
     % Allow automatic combination of Taylor kernels.
     %
     if !*taylorautocombine and not ('taysimpsq memq mul!*)
       then mul!* := aconc!* (mul!*, 'taysimpsq);
     !*tay2q resimptaylor u >>$

flag ('(Taylor!*), 'full)$

put ('Taylor!*, 'simpfn, 'simpTaylor!*)$

comment The following is necessary to properly process Taylor kernels
        in substitutions;

flag ('(Taylor!*), 'simp0fn);

endmodule;

end;