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;