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

module TayRevrt;

%*****************************************************************
%
%       Functions for reversion of Taylor kernels
%
%*****************************************************************


exports taylorrevert;

imports

% from the REDUCE kernel:
        !*a2k, !*q2a, invsq, lastpair, mvar, neq, nth, numr, over,
        reval, simp!*, typerr,

% from the header module:
        !*TayExp2q, cst!-Taylor!*, make!-cst!-coefficient,
        make!-taylor!*, multtaylorsq, prepTayExp, prune!-coefflist,
        set!-TayTemplate, TayCfPl, TayCfSq, TayCoeffList,
        TayExp!-Quotient, Taylor!:, TayMakeCoeff,
        Taylor!-kernel!-sq!-p, TayTemplate, TayTpElNext, TayTpElOrder,
        TayTpElPoint, TayTpElVars,

% from module Tayintro:
        delete!-nth, Taylor!-error,

% from module Taybasic:
        addtaylor, invtaylor, multtaylor, negtaylor,

% from module TaySimp:
        expttayrat,

% from module Tayutils:
        enter!-sorted, smallest!-increment;


symbolic procedure tayrevert (tay, okrnl, krnl);
  %
  % Reverts Taylor kernel tay that has okrnl as variable
  %  into a Taylor kernel that has krnl as variable.
  %
  % This is the driver procedure; it check whether okrnl
  %  is valid for this operation and calls tayrevert1 to do the work.
  %
  begin scalar tp, cfl, x; integer i;
    cfl := TayCoeffList tay;
    tp := TayTemplate tay;
    x := tp;
    i := 1;
    %
    % Loop through the template and make sure that the kernel
    %  okrnl is present and not part of a homogeneous template.
    %
   loop:
    if okrnl member TayTpElVars car x then <<
      if not null cdr TayTpElVars car x then <<
        Taylor!-error ('tayrevert,
                       {"Kernel", okrnl,
                        "appears in homogenous template", car x});
        return
        >>
       else goto found;
      >>
     else <<
       x := cdr x;
       i := i + 1;
       if not null x then goto loop
       >>;
    Taylor!-error
      ('tayrevert, {"Kernel", okrnl, "not found in template"});
    return;
   found:
    return tayrevert1 (tp, cfl, car x, i, okrnl, krnl)
  end;


symbolic procedure tayrevertreorder (cf, i);
  %
  % reorders coefflist cf so that
  %  a) part i of template is put first
  %  b) the resulting coefflist is again ordered properly
  %
  begin scalar cf1, pl, sq;
    for each pp in cf do <<
      pl := TayCfPl pp;
      sq := TayCfSq pp;
      pl := nth (pl, i) . delete!-nth (pl, i);
      cf1 := enter!-sorted (TayMakeCoeff (pl, sq), cf1)
      >>;
    return cf1
  end;

symbolic procedure tayrevertfirstdegreecoeffs (cf, n);
  %
  % Returns a coefflist that corresponds to the coefficient
  %  of (the first kernel in the template) ** n.
  %
  for each el in cf join
    if car car TayCfPl el = n and not null numr TayCfSq el
     then {TayMakeCoeff (cdr TayCfPl el, TayCfSq el)} else nil;

symbolic procedure tayrevert1(tp,cf,el,i,okrnl,krnl);
  %
  % This is the procedure that does the real work.
  %  tp is the old template,
  %  cf the old coefflist,
  %  el the element of the template that contains the "old" kernel,
  %  i its position in the template,
  %  okrnl the old kernel,
  %  krnl the new kernel.
  %
  Taylor!:
  begin scalar first,incr,newtp,newcf,newpoint,newel,u,u!-k,v,w,x,x1,n,
               expo,upper;
    %
    % First step: reorder the coefflist as if the okrnl appears
    %              at the beginning of the template and construct a
    %              new template by first deleting this element from it.
    %
    newcf := prune!-coefflist tayrevertreorder (cf, i);
    newtp := delete!-nth (tp, i);
    %
    % Check that the lowest degree of okrnl is -1, 0, or +1.
    %  For -1, we have a first order pole.
    %  For 0, reversion is possible only if the coefficient
    %   of okrnl is a constant w.r.t. the other Taylor variables.
    %  For +1, we use the algorithm quoted by Knuth,
    %   in: The Art of Computer Programming, vol2. p. 508.
    %
    n := car car TayCfPl car newcf;
    if n < 0
      then tayrevert1pole(tp,cf,el,i,okrnl,krnl,newcf,newtp);
    if n = 0
      then if not null newtp
             then begin scalar xx;
               xx := tayrevertfirstdegreecoeffs(newcf,0);
               if length xx > 1
                 then Taylor!-error
                       ('tayrevert,
                        "Term with power 0 is a Taylor series");
               xx := car xx;
               for each el in TayCfPl xx do
                 for each el2 in el do
                   if el2 neq 0
                     then Taylor!-error
                           ('tayrevert,
                            "Term with power 0 is a Taylor series");
               newpoint := !*q2a TayCfSq xx;
              end
            else <<newpoint := !*q2a TayCfSq car newcf;
                   newcf := prune!-coefflist cdr newcf;
                   n := car car TayCfPl car newcf>>
     else newpoint := 0;
    tp := {{krnl},newpoint,TayTpElOrder el,TayTpElNext el} . newtp;
    first := TayExp!-quotient(1,n);
    incr := car smallest!-increment newcf;
    expo := first * incr;
    if not(expo=1)
      then (<<newcf := TayCoeffList newtay;
              tp := TayTemplate newtay;
              newtp := cdr tp;
              tp := {TayTpElVars car tp,
                     reval {'expt,TayTpElPoint car tp,prepTayExp expo},
                     TayTpElOrder car tp * expo,
                     TayTpElNext car tp * expo}
                    . newtp>>
            where newtay := expttayrat(Make!-Taylor!*(newcf,tp,nil,nil),
                                       !*TayExp2q expo));
    upper := TayExp!-quotient(TayTpElNext car tp,incr) - 1;
    x := tayrevertfirstdegreecoeffs(newcf,incr);
    x1 := x := invtaylor make!-Taylor!*(x,newtp,nil,nil);
    w := for each pp in TayCoeffList x1 collect
           TayMakeCoeff({expo} . TayCfPl pp,TayCfSq pp);

    v := mkvect upper;

    for j := 2 : upper do
      putv(v,j,
           multtaylor(
               x,
               make!-Taylor!*(tayrevertfirstdegreecoeffs(newcf,j*incr),
                              newtp,nil,nil)));

    u := mkvect upper;

    putv(u,0,cst!-Taylor!*(1 ./ 1,newtp));
    for j := 2 : upper do <<
      for k := 1 : j - 2 do begin
        u!-k := getv(u,k);
        for l := k - 1 step -1 until 0 do
          u!-k := addtaylor
                   (u!-k,
                    negtaylor
                      multtaylor(getv(u,l),
                                 getv(v,k - l + 1)));
        putv(u,k,u!-k);
        end;
      u!-k := multtaylorsq(getv(v,j),!*TayExp2q j);

      for k := 1 : j - 2 do
        u!-k := addtaylor
                 (multtaylorsq(multtaylor(getv(u,k),
                                          getv(v,j - k)),
                               !*TayExp2q (j - k)),
                  u!-k);

      u!-k := negtaylor u!-k;
      putv(u,j - 1,u!-k);
%
      x1 := multtaylor(x1,x);            % x1 is now x ** j
      for each pp in TayCoeffList
             multtaylor(multtaylorsq
                           (getv(u,j - 1),
                            invsq !*TayExp2q j),x1) do
        w := enter!-sorted (TayMakeCoeff({j * expo}
                         . TayCfPl pp,TayCfSq pp),
                            w);
      >>;
%
      newtp := (car tp) . newtp;
      w := enter!-sorted(
             make!-cst!-coefficient(simp!* TayTpElPoint el,newtp),
             w);

    w := Make!-taylor!*(w,newtp,nil,nil);
    return if incr = 1 then w
            else expttayrat(w,invsq !*TayExp2q incr)
  end;

comment The mechanism for a first order pole is very simple:
        This corresponds to a first order zero at infinity,
        so we invert the original kernel and revert the result;

symbolic procedure tayrevert1pole (tp, cf, el, i, okrnl, krnl,
                                   newcf, newtp);
  begin scalar x, y, z;
    cf := TayCoeffList invtaylor make!-Taylor!*(cf,tp,nil,nil);
    x := tayrevert1 (tp, cf, el, i, okrnl, krnl);
    y := TayTemplate x;
    if TayTpElPoint car y neq 0
      then Taylor!-error ('not!-implemented,
                          "(Taylor series reversion)")
     else <<
       set!-TayTemplate (x, {{krnl}, 'infinity, TayTpElOrder car y}
                              . cdr y);
       return x >>
  end;   

comment The driver routine;

symbolic procedure TaylorRevert (u, okrnl, nkrnl);
  (if not Taylor!-kernel!-sq!-p sq
     then typerr (u, "Taylor kernel")
    else tayrevert (mvar numr sq, !*a2k okrnl, !*a2k nkrnl))
     where sq := simp!* u$
  
flag ('(TaylorRevert), 'opfn);

endmodule;

end;