module tpsrev; % Power Series Reversion & Composition

% Author: Alan Barnes   November 1988
%
% If y is a power series in x then psreverse expresses x as a power
% series in y-y0 where y0 is zero order term of y.
% This is known as power series reversion (functional inverse)
% pscompose functionally composes two power series
%
%Two new prefix operators are introduced PSREV and PSCOMP.
%These appear in the expression part of the power series objects
%generated by calls to psreverse and pscompose respectively.
%The  argument of PSREV is the 'generating  series' of the
%series (PS1 say) to be inverted. This is a generalised power series
%object  which looks like a standard power series object except that
%each of its terms is itself a power series (rather than a standard
%quotient), the  nth term being the power series of the nth power of
%PS1. The expression part of the generating series is (PSGEN <PS1>). 
%
%When power series PS1 and PS2 are composed (i.e. PS2 is substituted
%for the expansion variable of PS1 and the result expressed as a power
%series in the expansion variable of PS2), the expression part of
%the power series object generated is
%              (PSCOMP <PS1> <generating-series of PS2>)
%The generating series should only appear inside the operators PSREV
%and PSCOMP and not at 'top level'. It cannot sensibly be printed with
%the power series print function. Special functions are needed to
%access and modify terms of the generating series, although these
%are simply defined in terms of the functions for manipulating
%standard power series objects.
%% The algorithms used are based on those described in
%Feldmar E & Kolbig K S, Computer Physics Commun.  39, 267-284 (1986).

fluid '(ps);

put('psreverse, 'simpfn, 'simppsrev);

symbolic procedure simppsrev a;
 if length a=1 then apply('simppsrev1,a)
 else rerror(tps,33,"Wrong number of arguments to PSREVERSE");

symbolic procedure simppsrev1(series);
 begin scalar rev,psord, depvar,about, knownps, ps!:level;
   ps!:level:=0;
   series:=prepsqxx simp!* series;
   if not ps!:p series then
      rerror(tps,34,
           "Argument should be a <POWER SERIES>: simppsrev");
   ps!:find!-order series;
   depvar:=ps!:depvar series;
   if (psord:=ps!:order series)=1 then
        about:=0
    else if (psord=0) and (ps!:evaluate(series,1) neq (nil ./ 1)) then
        about := prepsqxx ps!:get!-term(series,0)
    else if psord =-1 then about:='ps!:inf
    else rerror(tps,35,"Series cannot be inverted:  simppsrev");
   rev:=ps!:compile(list('psrev,series),depvar,about);
   if ps!:expansion!-point series = 'ps!:inf then <<
      rev := make!-ps(list('quotient,1,rev),
                      ps!:value rev,depvar,about);
      ps!:find!-order rev>>;
   return rev ./ 1;
 end;

symbolic procedure ps!:generating!-series(a,psord,inverted);
 begin scalar ps;
   ps:=make!-ps(list('psgen, a,inverted),ps!:value a,
                      ps!:depvar a, ps!:expansion!-point a);
   ps!:set!-order(ps,psord);
   ps!:set!-rthpow(ps,psord);
   return ps
 end;

symbolic smacro procedure ps!:get!-rthpow(genseries,r);
  ps!:get!-term(genseries,r);

symbolic procedure ps!:set!-rthpow(genseries,r);
 begin scalar rthpow, series, power;
    series:=ps!:expression genseries;
    power:= if rand2 series then -r else r;
    series:=rand1 series;
    if power = 0 then 
       rthpow := 1
    else if power=1 then 
       rthpow := series
    else << 
      if power = -1 then 
          rthpow := list('quotient, 1, series)
      else if power = 2 then
          rthpow := list('times, series, series)
      else
          rthpow := list('expt, series, power, 1);
      power := if rator rthpow = 'expt then 
                    list('expt, series, power)
               else rthpow;
      rthpow := make!-ps(rthpow, ps!:arg!-values power,
                         ps!:depvar series,ps!:expansion!-point series);
      ps!:find!-order rthpow >>;
    ps!:set!-term(genseries,r,rthpow);
    return rthpow
 end;

symbolic procedure ps!:term!-rthpow(genseries,r,n);
 begin scalar term,series;
  series:= ps!:get!-rthpow(genseries,r);
  if null series then 
     for i:=ps!:last!-term genseries +1:r do
	series:=ps!:set!-rthpow(genseries,i);
  term:=  ps!:evaluate(series,n);
  return term
 end;

put('psrev,'ps!:crule,'ps!:rev!-crule);

symbolic procedure ps!:rev!-crule(a,d,n);
  begin scalar series;
    series :=rand1 a;
    if (n neq 'ps!:inf) and (n neq 0) then
       series := ps!:remove!-constant series;
    series := make!-ps(list('psrev,
                       ps!:generating!-series(series,1,
                                              if n='ps!:inf then t
                                              else nil)),
                       list('psrev,ps!:value rand1 a, d, n),d,n);
    ps!:find!-order series;
    return series;
  end;

symbolic procedure ps!:remove!-constant(ps);
 ps!:compile(list('difference, ps,prepsqxx ps!:evaluate(ps,0)),
             ps!:depvar ps,
             ps!:expansion!-point ps);

put('psrev,'ps!:erule,'ps!:rev!-erule);
put('psrev,'ps!:order!-fn,'ps!:rev!-orderfn);


symbolic procedure ps!:rev!-orderfn ps;
   begin scalar u;
     u:=ps!:expansion!-point ps!:get!-rthpow(rand1 ps!:expression ps,1);
     return if (u=0) or (u = 'ps!:inf) then 1
            else 0
   end;


symbolic procedure ps!:rev!-erule(a,n);
 begin scalar genseries,x,z;
  z:=nil ./ 1; genseries:=rand1 a;
  if n=0 then
    if (x:=ps!:expansion!-point ps!:get!-rthpow(genseries,1))='ps!:inf
     then return (nil ./ 1)
     else return simp!* x;
  if n=1 then
      return invsq ps!:term!-rthpow(genseries,1,1);
  for i:=1:n-1 do
        z:=addsq(z,multsq(ps!:evaluate(ps,i),
                          ps!:term!-rthpow(genseries,i,n)));
  return quotsq(negsq z,ps!:term!-rthpow(genseries,n,n))
 end;

put('pscomp,'ps!:crule,'ps!:comp!-crule);
put('pscomp,'ps!:erule,'ps!:comp!-erule);
put('pscomp,'ps!:order!-fn,'ps!:comp!-orderfn);

symbolic procedure ps!:comp!-orderfn ps;
 begin scalar u;
   u:=ps!:find!-order rand1 ps!:expression ps;
   return 
     if u=0 then 0
     else 
       ps!:find!-order(ps!:get!-rthpow(rand2 ps!:expression ps,u));
 end;


symbolic procedure ps!:comp!-crule(a,d,n);
  begin scalar series1,series2,n1;
     series1:=rand1 a; series2:=rand2 a;
     n1 := ps!:expansion!-point series1;
     if (n1 neq 0) and (n1 neq 'ps!:inf) then
        series2:=ps!:remove!-constant series2;
     a:= make!-ps(list('pscomp,series1,
                       ps!:generating!-series(series2,
                                              ps!:order series1,
                                              if n1='ps!:inf then t
                                              else nil)),
                  append(ps!:arg!-values a, list(d,n)), d, n);
    ps!:find!-order a;
    return a;
  end;

symbolic procedure ps!:comp!-erule(a,n);
 begin scalar aa,genseries,z,psord1;
  z:=nil ./ 1; aa:=rand1 a; genseries:=rand2 a;
  psord1:=ps!:order aa;
  for i:=psord1:n do
        z:=addsq(z,multsq(ps!:evaluate(aa,i),
                          ps!:term!-rthpow(genseries,i,n)));
   return z
 end;


put('pscompose, 'simpfn, 'simppscomp);

symbolic procedure simppscomp a;
  if length a=2 then  apply('simppscomp1,a)
   else rerror(tps,36,
     "Args should be <POWER SERIES>,<POWER SERIES>:  simppscomp");

symbolic procedure simppscomp1(ps1,ps2);
  begin scalar x,d,n1,n, knownps, ps!:level;
    ps!:level:=0;
    ps1:=prepsqxx simp!* ps1;
    if ps!:numberp ps1 then
      return ((if zerop ps1 then nil else ps1) ./ 1);
    if not ps!:p ps1 or not ps!:p(ps2:=prepsqxx simp!* ps2) then
      rerror(tps,37,
      "Args should be <POWER SERIES>,<POWER SERIES>:  simppscomp");
    ps!:find!-order ps1;
    x:=ps!:find!-order ps2;
    d:= ps!:depvar ps2;
    n1:= ps!:expansion!-point ps1;
    n:= ps!:expansion!-point ps2;
    if (x >0 and n1 = 0) or
       (x <0 and n1 = 'ps!:inf) or
       (x=0 and n1=prepsqxx ps!:evaluate(ps2,0))
     then  return ps!:compile(list('pscomp,ps1,ps2),d,n) ./ 1
     else rerror(tps,38,"Series cannot be composed:  simppscomp");
 end;

algebraic operator psrev,pscomp;

endmodule;

end;
