Artifact bf66c251f0ee4ac3d6906fbd1e3386bd79192fd549b8785f94821a4558ff8a6a:
- Executable file
r37/packages/tps/tpsrev.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 8352) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/tps/tpsrev.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 8352) [annotate] [blame] [check-ins using]
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;