Artifact 0d95a183129a1862a23f0504b156a5ed87965a24472cf7b6dfd9123a2cae6eb6:
- Executable file
r36/src/tps.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: 55472) [annotate] [blame] [check-ins using] [more...]
module tps; % Truncated Power Series. % Authors: Julian Padget & Alan Barnes <barnesa@aston.ac.uk>. % Version 1.31 March 1993. % % Revisions: % % 20 Mar 93. Bug in PS!:EXPT!-CRULE corrected. % ps!:order!-limit increased to 100. % % 16 Sep 90. Bugs in PS!:EXPT!-CRULE, PS!:EXPT!-ERULE corrected. % Treatment of expt generally improved. % PSREVERSE now only takes one argument: the series to be % inverted. PSCHANGEVAR added for those who feel the % need to change the expansion variable. % Code for pscompose and psreverse generalised to handle % the point at infinity correctly and to deal with series % with a negative order correctly. % Code for simppssum improved to detect non- % strictly increasing exponents. % % 6 Sep 90 PSSUM added to allow definition of series whose % general term is known. % % 8 Aug 90. FUNCTION changed to QUOTE to avoid FUNARG problem when % interpreting code on SLISP based systems. Need to be % careful with quoted lambdas! % % 26 Jul 90. JHD's smacro's ps!:numberp and ps!:atom added to improve % behaviour of system after ON BIGFLOAT. % Still dicky with NUMVAL ON. Need to REMPROP properties % !:bf!:, !:ft!: etc. from !:ps!: ? (AB) % 25 Jul 90. Global declaration of s added in ps:unknown-crule. Needed % in UOLISP, avoids warning messages in some other Lisps. % simp:ps1 tidied up (now uses selectors to access terms) % 24 Jan 90. Ps:evaluate corrected (missing assignment for next). % 7 Jul 89. Explicit compilation and evaluation rules for SQRT and % EXPT have been added. This avoids the infinite loops that % could sometimes occur when ps!: unknown!-crule was used to % generate recurrence relations for rational powers of a % power series. % A few bugs have been fixed notably one in difff and the % efficiency of code has been improved in several places. % Experimental code has been added to allow the domain mode % to be set to TPS by typing ON TPS. With this switch on, % quotients of power series are expanded automatically as % are (when NUMVAL is on) expressions such as sin a, where a % is a power series. % 3 Mar 89. Addition of code for reversion of series. Modification % of internal form to avoid circular lists as arguments. % Minor bug fixes. % A power series object is a tagged tuple of the following form: % % (:ps: . [<order>, % <last-term>, % <variable>, % <expansion-point>, % <value>, % <terms>, % <ps-expression]) % % <order> is the exponent of the first term of the series and is also % used to modify the index when accessing components of the % series which are addressed by power % % <last-term> the power of the last term generated in the series % largely for the benefit of the printing routine % % <variable> is the dependent variable of this expansion, needed, in % particular, for printing and when combining two series % % <expansion!-point> is self-explanatory except that % ps!:inf denotes expansion about infinity % % <value> is the originating prefix form which is needed to allow for % power series variables appearing inside other power series % expressions % % <terms> is an alist containing the terms of the series computed so % far, access is controlled using <order> as an index base. % % <ps-expression> is a power series object corresponding to the prefix % form of which the expansion was requested, the first element % of which is the ps!:operator and the rest of which are the % ps!:operands which may themselves be pso's % % In addition we have the following degenerate forms of power series % object: % (!:ps!: . <identifier>) the value of <identifier> is a vector % as above(used in automatically generated recurrence relations) % <number> % <identifier> 2nd argument of DF, INT create!-package('(tps tpscomp tpseval tpsdom tpsrev tpsfns tpselfns tpssum), '(contrib tps)); fluid '(ps!:exp!-lim knownps ps!:max!-order); flag('(!:ps!:),'noconvert); % Some structure selectors and referencers. symbolic smacro procedure rands e; cdr e; symbolic smacro procedure rand1 e; cadr e; symbolic smacro procedure rand2 e; caddr e; symbolic smacro procedure rator e; car e; symbolic smacro procedure ps!:domainp u; atom u or (car u neq '!:ps!:) and not listp u; symbolic smacro procedure constantpsp u; null ps!:depvar u; symbolic smacro procedure ps!:p u; pairp u and (car u = '!:ps!:); symbolic smacro procedure ps!:atom u; atom u or (car u neq '!:ps!: and get(car u,'dname)); symbolic smacro procedure ps!:numberp u; numberp u or (car u neq '!:ps!: and get(car u,'dname)); symbolic procedure ps!:getv(ps,i); if eqcar(ps,'!:ps!:) then if idp cdr ps then getv(eval cdr ps,i) else getv(cdr ps,i) else rerror(tps,1,list("PS:GETV: not a ps",ps)); symbolic procedure ps!:putv(ps,i,v); if eqcar(ps,'!:ps!:) then if idp cdr ps then putv(eval cdr ps,i,v) else putv(cdr ps,i,v) else rerror(tps,2,list("PS:PUTV: not a ps",ps)); symbolic procedure ps!:order ps; if ps!:atom ps then 0 else ps!:getv(ps,0); symbolic smacro procedure ps!:set!-order(ps,n); ps!:putv(ps,0,n); symbolic procedure ps!:last!-term ps; if ps!:atom ps then ps!:max!-order else ps!:getv(ps,1); symbolic (ps!:max!-order:= 2147483647); % symbolic here seems to be essential in Cambridge Lisp systems symbolic smacro procedure ps!:set!-last!-term (ps,n); ps!:putv(ps,1,n); symbolic procedure ps!:depvar ps; if ps!:atom ps then nil else ps!:getv(ps,2); symbolic smacro procedure ps!:set!-depvar(ps,x); ps!:putv(ps,2,x); symbolic procedure ps!:expansion!-point ps; if ps!:atom ps then nil else ps!:getv(ps,3); symbolic smacro procedure ps!:set!-expansion!-point(ps,x); ps!:putv(ps,3,x); symbolic procedure ps!:value ps; if ps!:atom ps then if ps then ps else 0 else ps!:getv(ps,4); symbolic smacro procedure ps!:set!-value(ps,x); ps!:putv(ps,4,x); symbolic smacro procedure ps!:terms ps; if ps!:atom ps then list (0 . ( ps . 1)) else ps!:getv(ps,5); symbolic smacro procedure ps!:set!-terms(ps,x); ps!:putv(ps,5,x); symbolic procedure ps!:expression ps; if ps!:atom ps then ps else ps!:getv(ps,6); symbolic smacro procedure ps!:set!-expression(ps,x); ps!:putv(ps,6,x); symbolic smacro procedure ps!:operator ps; car ps!:getv(ps,6); symbolic smacro procedure ps!:operands ps; cdr ps!:getv(ps,6); symbolic procedure ps!:get!-term(ps,i); (lambda(psorder,pslast); if null psorder or null pslast then nil else if i<psorder then nil ./ 1 else if i>pslast then nil else begin scalar term; term:=assoc(i-psorder, ps!:terms ps); return if term then cdr term else nil ./ 1; end) (ps!:order ps,ps!:last!-term ps); symbolic procedure ps!:term(ps,i); begin scalar term; term:=ps!:get!-term (ps,i); if term then return term; for j:=ps!:last!-term(ps)+1:i do term:= ps!:evaluate(ps,j); return term end; symbolic procedure ps!:set!-term(ps,n,x); (lambda (psorder,terms); <<if null psorder then psorder:=ps!:find!-order ps; if n < psorder then rerror(tps,3,list (n, "less than the order of ", ps)) else (lambda term; <<if (n=psorder and x = '(nil . 1)) then ps!:set!-order(ps,n+1); if term then rplacd(term,x) else if atom x or (numr x neq nil) then if terms then nconc(terms,list((n-psorder).x)) else ps!:set!-terms(ps,list((n-psorder).x)) >>) assoc(n-psorder,terms)>>) (ps!:order ps,ps!:terms ps); put('psexplim, 'simpfn, 'simppsexplim); symbolic (ps!:exp!-lim := 6); % default depth of expansion % symbolic here seems to be essential in Cambridge Lisp systems symbolic procedure simppsexplim u; begin integer n; n:=ps!:exp!-lim; if u then ps!:exp!-lim := ieval carx(u,'psexplim); return (if n=0 then nil ./ 1 else n ./ 1); end; symbolic procedure simpps a; if length a =3 then apply('simpps1,a) else rerror(tps,4, "Args should be <FORM>,<depvar>, and <point>: simpps"); put('ps,'simpfn,'simpps); symbolic procedure simpps1(form,depvar,about); if form=nil then rerror(tps,5,"Args should be <FORM>,<depvar>, and <point>: simpps") else if not kernp simp!* depvar then typerr(depvar, "kernel: simpps") else if smember(depvar,(about:=prepsq simp!* about)) then rerror(tps,6,"Expansion point depends on depvar: simpps") else begin scalar knownps; return ps!:compile(ps!:presimp form, depvar, if about='infinity then 'ps!:inf else about) ./ 1 end; put('psterm,'simpfn,'simppsterm); symbolic procedure simppsterm a; if length a=2 then apply('simppsterm1,a) else rerror(tps,7, "Args should be of form <power series>,<term>: simppsterm"); symbolic procedure simppsterm1(p,n); << n := ieval n; p:=prepsq simp!* p; if ps!:numberp p then if n neq 0 or p=0 then nil ./ 1 else p ./ 1 else if ps!:p p then << ps!:find!-order p; ps!:term(p,n)>> else typerr(p,"power series: simppsterm1") >>; put('psorder,'simpfn,'simppsorder); put('pssetorder,'simpfn,'simppssetorder); symbolic procedure simppsorder u; << u:=prepsq simp!* carx(u,'psorder); if ps!:numberp u then if u=0 then 'undefined else nil else if ps!:p u then ps!:find!-order u else typerr(u,"power series: simppsorder") >> ./ 1; symbolic procedure simppssetorder u; (lambda (psord,ps); if not ps!:p ps then typerr(ps,"power series: simppssetorder") else if not fixp psord then typerr(psord, "integer: simppssetorder") else <<u:=ps!:order ps; ps!:set!-order(ps,psord); (if u=0 then nil else u) ./ 1>>) (prepsq simp!* carx(cdr u,'pssetorder),prepsq simp!* car u); put('psexpansionpt,'simpfn,'simppsexpansionpt); symbolic procedure simppsexpansionpt u; << u:=prepsq simp!* carx(u,'psexpansionpt); if ps!:numberp u then 'undefined ./ 1 else if ps!:p u then (lambda about; if about neq 'ps!:inf then simp!* about else 'infinity ./ 1) (ps!:expansion!-point u) else typerr(u,"power series: simppsexpansionpt") >>; put('psdepvar,'simpfn,'simppsdepvar); symbolic procedure simppsdepvar u; << u:=prepsq simp!* carx(u,'psdepvar); if ps!:numberp u then 'undefined else if ps!:p u then ps!:depvar u else typerr(u,"power series: simppsdepvar") >> ./ 1; put('psfunction,'simpfn,'simppsfunction); symbolic procedure simppsfunction u; << u:=prepsq simp!* carx(u,'psfunction); if ps!:numberp u then u ./ 1 else if ps!:p u then simp!* ps!:value u else typerr(u,"power series: simppsfunction") >>; symbolic procedure ps!:presimp form; if (pairp form) and ((rator form = 'expt) or (rator form = 'int)) then list(rator form, prepsort rand1 form, prepsort rand2 form) else prepsort form; symbolic procedure prepsort u; % Improves log handling if logsort is defined. S.L. Kameny. if getd 'logsort then logsort u else prepsq simp!* u; % symbolic procedure tag!-if!-float n; % if floatp n then int!-equiv!-chk mkfloat n else n; symbolic procedure !*pre2dp u; begin scalar x; u:=simp!* u; return if fixp denr u % then if denr u = 1 and domainp(x:= tag!-if!-float numr u) % then x then if denr u = 1 and domainp(x:= numr u) then x else if fixp numr u then mkrn(numr u, denr u) end; flag('(!:ps!:),'full); put('!:ps!:,'simpfn,'simp!:ps!:); symbolic procedure simp!:ps!: ps; simp!:ps1 ps ./1; symbolic procedure simp!:ps1 ps; if ps!:atom ps or idp cdr ps then ps else <<if not atom ps!:expression ps then begin scalar i, term, simpfn; if (ps!:operator ps ='psgen) then simpfn:= 'simp!:ps1 else simpfn:= 'resimp; i:=ps!:order ps; while term:=ps!:get!-term(ps,i) do << ps!:set!-term(ps,i,apply1(simpfn,term)); i:=i+1>> end; if atom ps!:expression ps or null ps!:depvar ps or (ps!:operator ps ='ps!:summation) then nil else<<ps!:set!-expression(ps, (ps!:operator ps) . for each j in ps!:operands ps collect simp!:ps1 j); ps!:set!-value(ps,prepsq simp!* ps!:value ps)>>; ps>>; % symbolic procedure simp!:ps1 ps; % if ps!:atom ps or idp cdr ps then ps % else % <<if not atom ps!:expression ps then % begin scalar i, term, simpfn; % if (ps!:operator ps ='psgen) then % simpfn:= 'simp!:ps1 % else simpfn:= 'resimp; % i:=ps!:order ps; % while term:=ps!:get!-term(ps,i) do % << ps!:set!-term(ps,i,apply1(simpfn,term)); i:=i+1>> % end; % if atom ps!:expression ps or null ps!:depvar ps then nil % else<<ps!:set!-expression(ps, % (ps!:operator ps) . mapcar(ps!:operands ps, % 'simp!:ps1)); % ps!:set!-value(ps,prepsq simp!* ps!:value ps)>>; % ps>>; put('!:ps!:,'domain!-diff!-fn,'ps!:diff!:); put('pschangevar,'simpfn,'simppschangevar); symbolic procedure simppschangevar u; (lambda (newvar,ps, oldvar); if not ps!:p ps then typerr(ps,"power series: simppschangevar") else if not kernp newvar then typerr(prepsq newvar, "kernel: simppschangevar") else <<oldvar:=ps!:depvar ps; newvar:=prepsq newvar; if smember(newvar,ps!:value ps) and newvar neq oldvar then rerror(tps,8,"Series depends on new variable") else if oldvar then << ps!:set!-depvar(ps,newvar); ps!:set!-value(ps,subst(newvar,oldvar,ps!:value ps)); ps ./ 1>> else rerror(tps,9,"Can't change variable of constant series")>>) (simp!* carx(cdr u,'pschangevar),prepsq simp!* car u,nil); endmodule; module tpscomp; % Compile prefix expression into network of % communicating power series. % Authors: Julian Padget & Alan Barnes % The compiler is rule driven by looking for a compilation rule (crule) % property on the property list of the operator. If a rule does not % exist the expression is differentiated to get an expression which is % amenable to compilation but the process takes care to check for the % existence of cycles in the derivatives e.g. sine and cosine. % % The result is an power series object which can be evaluated by the % power series evaluator. fluid '(unknowns !*exp psintconst knownps ps!:max!-order); symbolic procedure ps!:compile(form,depvar,about); if idp form then make!-ps!-id(form,depvar,about) else if ps!:numberp form then form else if ps!:p form then if((ps!:expansion!-point form=about)and(ps!:depvar form=depvar)) then form else ps!:compile(ps!:value form,depvar,about) else if get(car form,'ps!:crule) then apply(get(car form,'ps!:crule),list(form,depvar,about)) else (if tmp then '!:ps!: . cdr tmp % ******was eval cdr tmp % get(cdr tmp,'!:ps!:) else ps!:unknown!-crule((car form) . foreach arg in cdr form collect ps!:compile(arg,depvar,about), depvar,about)) where tmp = assoc(form,knownps); symbolic procedure make!-ps!-id(id,depvar,about); begin scalar ps; ps:=make!-ps(id,id,depvar,about); if about neq 'ps!:inf then <<ps!:set!-term(ps,0,ps!:identifier!-erule(id,depvar,0,about)); ps!:set!-term(ps,1,ps!:identifier!-erule(id,depvar,1,about))>> else % expansion about infinity <<ps!:set!-order(ps,-1); ps!:set!-term(ps,-1, ps!:identifier!-erule(id,depvar,-1,about))>>; ps!:set!-last!-term(ps,ps!:max!-order); return ps end; symbolic procedure make!-ps(form,exp,depvar,about); % if f is a trivial expression (it can only ever be a number, % an identifier or a ps) then it is a 'constant' and stands for % itself, otherwise a larger ps is composed from it begin scalar ps; ps:=get('tps,'tag) . mkvect 6; ps!:set!-order(ps,0); ps!:set!-expression(ps,form); ps!:set!-value(ps,exp); ps!:set!-depvar(ps,depvar); ps!:set!-expansion!-point(ps,about); ps!:set!-last!-term(ps,-1); return ps end; %symbolic procedure ps!:plus!-crule(a,d,n); % if atom cdddr a then % make!-ps('plus . list(ps!:compile(rand1 a,d,n), % ps!:compile(rand2 a,d,n)), % ps!:arg!-values a,d,n) % else % make!-ps('plus . list(ps!:compile(rand1 a,d,n), % ps!:plus!-crule('plus . cddr a,d,n)), % ps!:arg!-values a,d,n); % put('plus,'ps!:crule,'ps!:plus!-crule); put('plus,'ps!:crule,'ps!:nary!-crule); % symbolic procedure ps!:minus!-crule(a,d,n); % make!-ps(list('minus,ps!:compile(rand1 a,d,n)), % 'minus . ps!:arg!-values a,d,n); % put('minus,'ps!:crule, 'ps!:minus!-crule); symbolic procedure ps!:unary!-crule(a,d,n); make!-ps(list(rator a,ps!:compile(rand1 a,d,n)), ps!:arg!-values a,d,n); put('minus,'ps!:crule, 'ps!:unary!-crule); symbolic procedure ps!:binary!-crule(a,d,n); make!-ps((rator a) . list(ps!:compile(rand1 a,d,n), ps!:compile(rand2 a,d,n)), ps!:arg!-values a,d,n); put('difference,'ps!:crule,'ps!:binary!-crule); symbolic procedure ps!:nary!-crule(a,d,n); if null cdddr a then make!-ps((rator a) . list(ps!:compile(rand1 a,d,n), ps!:compile(rand2 a,d,n)), ps!:arg!-values a,d,n) else make!-ps((rator a) . list(ps!:compile(rand1 a,d,n), ps!:nary!-crule((rator a) . cddr a,d,n)), ps!:arg!-values a,d,n); put('times,'ps!:crule,'ps!:nary!-crule); % put('times,'ps!:crule,'ps!:times!-crule); symbolic procedure ps!:quotient!-crule(a,d,n); % forms such as (quotient (expt <x> <y>) (expt <x> <z>)) are % detected here and transformed into (expt <x>(difference <y> <z>)) to % help avoid certain essential singularities if eqcar(rand1 a,'expt) and eqcar(rand2 a,'expt) and ((rand1 rand1 a)=(rand1 rand2 a)) then ps!:compile(list('expt, rand1 rand1 a, list('difference, rand2 rand1 a, rand2 rand2 a)),d,n) else ps!:binary!-crule(a,d,n); % make!-ps('quotient . list(ps!:compile(rand1 a,d,n), % ps!:compile(rand2 a,d,n)), % ps!:arg!-values a,d,n); put('quotient,'ps!:crule,'ps!:quotient!-crule); flag('(psintconst), 'share); algebraic (psintconst:= 0 ); symbolic procedure ps!:int!-crule(a,d,n); begin scalar r,arg1, psord; if not idp rand2 a then typerr(rand2 a, "kernel: ps!:int!-crule"); if rand2 a=d and smember(d,prepsq simp!* psintconst) then rerror(tps,10,list("psintconst depends on depvar: ",d)); arg1:=ps!:compile(prepsq simp!* rand1 a,d,n); r:= make!-ps(list('int,arg1,rand2 a),ps!:arg!-values a,d,n); psord:= ps!:find!-order arg1; if d=rand2 a then if ps!:expansion!-point(arg1) neq 'ps!:inf then if (psord < 0) and (ps!:term(arg1,-1) neq (nil ./ 1)) then rerror(tps,11,"Logarithmic Singularity") else <<ps!:set!-term(r,0,simp!* psintconst); ps!:find!-order r>> else % expansion about infinity if (psord < 2) and (ps!:term(arg1,1) neq (nil ./ 1)) then rerror(tps,12,"Logarithmic Singularity") else <<ps!:set!-term(r,0,simp!* psintconst); ps!:find!-order r>> else <<ps!:set!-term(r,0, simpint list(prepsq ps!:term(arg1,0), rand2 a)); ps!:find!-order r>>; return r; end; put('int,'ps!:crule,'ps!:int!-crule); symbolic procedure ps!:arg!-values funct; (rator funct) . (foreach arg in rands funct collect if ps!:atom arg then arg else if ps!:p arg then ps!:value arg else ps!:arg!-values arg); symbolic procedure ps!:unknown!-crule(a,d,n); % unknowns is an alist structure, the CAR of which is the % form which was differentiated and the CDR is a dotted pair whose % CDR is a gensym'ed identifier which is used to build % the cyclic structures used to represent a recurrence relation. % Minor improvements by Stan Kameny, July 1990. (lambda (dfdx,unknowns,tmp); if (tmp:=assoc(caar unknowns,cdr unknowns)) then cdr tmp else (lambda(r,s); << tmp:=ps!:arg!-values a; ps!:set!-value(r,tmp); % intern s; % apparently not needed, useful for debugging. global list s; % This is definitely needed in UOLISP. set(s,cdr r); knownps:=(tmp . s) . knownps; ps!:set!-order(r,0); % screws up if order is negative if n= 'ps!:inf then n:=0; % expansion about infinity a:=(rator a).(foreach arg in rands a collect if ps!:p arg then if ps!:find!-order arg < 0 then rerror(tps,13, "Essential Singularity") else prepsq ps!:term(arg,0) else subst(n,d,arg)); ps!:set!-term(r,0,simp!* a); ps!:set!-last!-term(r,0); r >> ) (make!-ps(list('int, ps!:compile(dfdx,d,n),d),a,d,n), cddar unknowns) ) (ps!:differentiate(a,d),(a . ('!:ps!: . gensym())) . unknowns,nil); symbolic procedure ps!:differentiate(a,v); % note the binding of !*factor to t; this ensures the result of the % differentiation will be factored, which prevents us looping % infinitely because we can't see the cycle in the derivative. % *********not necessary so removed March 1989 AB (lambda x; if eqcar(x,'df) then rerror(tps,14, list("ps:differentiate: cannot differentiate ",a)) else x) % the default method catches forms which are defined by patterns % (eg Bessel functions) % ((lambda (!*factor,!*exp); ((lambda (!*exp); if get(car a,'dfn) then prepsq diffsq(simp!* a,v) else prepsq simp!* list ('df,a,v)) % (t,nil)); (nil)); endmodule; module tpseval; % Evaluator for truncated power series. % Authors: Julian Padget & Alan Barnes % The evaluator interprets the results of the compilation phase and % is also rule driven until I get round to getting the compilation % phase to produce directly executable code % The evaluation functions live on the erule property of the name. fluid '(ps psintconst ps!:order!-limit); %symbolic (orig!*:= 0); % % symbolic here seems to be essential in Cambridge Lisp systems symbolic procedure ps!:prin!: p; (lambda (first,u,delta,symbolic!-exp!-pt,about,atinf); << if !*nat and posn!*<20 then orig!*:=posn!*; atinf:=(about='ps!:inf); ps!:find!-order p; delta:=prepf((ps!:depvar p) .** 1 .*1 .+ (negf if atinf then nil % expansion about infinity else if idp about then !*k2f about else if ps!:numberp about then !*n2f about else if (u:=!*pre2dp about) then !*n2f u else !*k2f(symbolic!-exp!-pt:= compress append(explode ps!:depvar p, explode '0)))); if symbolic!-exp!-pt then prin2!* "["; prin2!* "{"; for i:=(ps!:order p): ps!:exp!-lim do << u:=ps!:term(p,i); if not null numr u then <<if minusf numr u then <<u:=negsq u; prin2!* " - ">> else if not first then prin2!* " + "; first := nil; if posn!*>55 then <<terpri!* nil;prin2!* " ">>; if denr u neq 1 then prin2!* "("; if u neq '(1 . 1) then maprint(prepsq u,get('times,'infix)) else if i=0 then prin2!* 1; if denr u neq 1 then prin2!* ")"; if i neq 0 and u neq '(1 . 1) then prin2!* "*"; if i neq 0 then xprinf(!*p2f mksp(delta, if atinf then -i else i),nil,nil) >> >>; if first then prin2!* "0"; if posn!*>55 then terpri!* nil; u:=ps!:exp!-lim +1; if (u=1) and not atinf and (about neq 0) then prin2!* " + O" else prin2!* " + O("; xprinf(!*p2f mksp(delta,if atinf then -u else u),nil,nil); if (u=1) and not atinf and (about neq 0) then prin2!* "}" else prin2!* ")}"; if symbolic!-exp!-pt then << if posn!*>45 then terpri!* nil; prin2!* " where "; prin2!* symbolic!-exp!-pt; prin2!* " = "; maprin about; prin2!* "]" >>; terpri!* nil; >>) (t,nil,nil,nil,ps!:expansion!-point p,nil); symbolic procedure ps!:id!-order ps; (lambda u; if numberp u then u else rerror(tps,15, list("Can't find the order of ",ps!:value ps))) ps!:order ps; symbolic procedure ps!:find!-order ps; if null ps then 0 else if idp ps then ps % second arg of DF etc are identifiers else if numberp ps then 0 else if eqcar(ps,'!:ps!:) then << if idp cdr ps then ps!:id!-order ps else if atom ps!:expression ps or null ps!:depvar ps then ps!:order ps else ps!:find!-order1(ps) >> else if get(car ps,'dname) then 0 else rerror(tps,16,"Unexpected form in ps!:find!-order"); symbolic procedure ps!:find!-order1(ps); begin scalar psoperator,psord,pslast; psoperator:=ps!:operator ps; psord:=ps!:order ps; pslast:=ps!:last!-term ps; if psord leq pslast then if psoperator neq 'int then return psord else if (psord neq 0) or (pslast neq 0) then return psord; psord:=apply(get(psoperator,'ps!:order!-fn), for each j in ps!:operands ps collect ps!:find!-order j); ps!:set!-order(ps,psord); if psoperator= 'int and psord=0 then nil else ps!:set!-last!-term(ps,psord-1); if ps!:value ps =0 then <<psord:=0;ps!:set!-last!-term(ps,1000000)>> % prevents infinite loop if we have exact cancellation else while ps!:evaluate(ps,psord)=(nil ./ 1 ) do % in case we have finite # of cancellations in a sum or difference <<psord:=psord+1; if psord > ps!:order!-limit then rerror(tps,17,list("Expression ",ps!:value ps, " has zero expansion to order ", psord)) % We may not always be able to recognise zero. % Anyway give up after 15 iterations. >>; return psord end; symbolic (ps!:order!-limit:=100); % symbolic here seems to be essential in Cambridge Lisp systems put('psordlim, 'simpfn, 'simppsordlim); symbolic procedure simppsordlim u; begin integer n; n:=ps!:order!-limit; if u then ps!:order!-limit := ieval carx(u,'psordlim); return (if n=0 then nil ./ 1 else n ./ 1); end; put('times,'ps!:order!-fn, 'plus2); put('quotient,'ps!:order!-fn, 'difference); put('plus,'ps!:order!-fn, 'min2); put('difference,'ps!:order!-fn, 'min2); put('minus,'ps!:order!-fn, '(lambda (u) u)); put('int,'ps!:order!-fn,'ps!:int!-orderfn); put('df,'ps!:order!-fn,'ps!:df!-orderfn); symbolic procedure ps!:int!-orderfn(u,v); if (ps!:order ps=0) and (u geq 0) then 0 else if v=ps!:depvar ps then if ps!:expansion!-point ps neq 'ps!:inf then if u=-1 then rerror(tps,18,"Logarithmic Singularity") else u+1 else % expansion about infinity if u=1 then rerror(tps,19,"Logarithmic Singularity") else u-1 else u; symbolic procedure ps!:df!-orderfn (u,v); if v=ps!:depvar ps then if ps!:expansion!-point ps neq 'ps!:inf then if u=0 then 0 else u-1 else if u=0 then 2 else u+1 % expansion about infinity else u; symbolic procedure ps!:number!-erule(n,i); % << n:=(if numberp n then tag!-if!-float n else cdr n); <<n := if numberp n then n else cdr n; if i =0 and n neq 0 then n ./ 1 else nil ./ 1>>; symbolic procedure ps!:identifier!-erule(v,d,n,about); if v=d then if about='ps!:inf then if n=-1 then 1 ./ 1 else nil ./ 1 % expansion about infinity else if n=0 then if idp about then !*k2q about % else if ps!:numberp about then !*n2f tag!-if!-float about ./ 1 else if ps!:numberp about then !*n2f about ./ 1 else simp!* about else if n=1 then 1 ./ 1 else nil ./ 1 else if n=0 then !*k2q v else nil ./ 1; symbolic procedure ps!:evaluate(ps,n); % ps can come in two flavours (one of which is degenerate): % (i) a number % (ii) a power series object % in the last case the appropriate evaluation rule for the operator % in the ps is selected and invoked if n leq ps!:last!-term ps then ps!:get!-term(ps,n) else (lambda next; << if n < ps!:order ps then nil else << ps!:set!-term(ps,n,next:=simp!* prepsq next); ps!:set!-last!-term(ps,n) >>; next>>) apply(get(ps!:operator ps,'ps!:erule), list(ps!:expression ps,n)); symbolic procedure ps!:plus!-erule(a,n); addsq(ps!:evaluate(rand1 a,n), ps!:evaluate(rand2 a,n)); put('plus,'ps!:erule,'ps!:plus!-erule); symbolic procedure ps!:minus!-erule(a,n); negsq ps!:evaluate(rand1 a,n); put('minus,'ps!:erule,'ps!:minus!-erule); symbolic procedure ps!:difference!-erule(a,n); addsq(ps!:evaluate(rand1 a,n), negsq ps!:evaluate(rand2 a,n)); put('difference,'ps!:erule,'ps!:difference!-erule); symbolic procedure ps!:times!-erule(a,n); begin scalar aa,b,x,y,y1,z; aa:=rand1 a; b:= rand2 a; z:= nil ./ 1; x:=ps!:order(aa); y:=ps!:order(ps); % order of product ps y1 := ps!:order b; % Next "if" suggested by A.C. Norman to avoid different tan df rule % The problem with tan was in fact due to ps!:find!-order! - AB for i := 0:n-y do if n-x-i>=y1 then z:= addsq(z,multsq(ps!:evaluate(aa,i+x), ps!:evaluate(b,n-x-i))); return z end; put('times,'ps!:erule,'ps!:times!-erule); symbolic procedure ps!:quotient!-erule(a,n); begin scalar aa,b,x,y,z; aa:=rand1 a; b:=rand2 a; z:= nil ./ 1; y:=ps!:order(b); x:=ps!:order(ps); %order of quotient ps for i:=1:n-x do z:=addsq(z,multsq(ps!:evaluate(b,i+y), ps!:evaluate(ps,n-i))); return quotsq(addsq(ps!:evaluate(aa,n+y),negsq z), ps!:evaluate(b,y)) end; put('quotient,'ps!:erule,'ps!:quotient!-erule); symbolic procedure ps!:differentiate!-erule(a,n); if rand2 a =ps!:depvar rand1 a then if ps!:expansion!-point rand1 a neq 'ps!:inf then multsq((n+1) ./ 1,ps!:evaluate(rand1 a,n+1)) else multsq((1-n) ./ 1,ps!:evaluate(rand1 a,n-1)) else diffsq(ps!:evaluate(rand1 a,n),rand2 a); put('df,'ps!:erule,'ps!:differentiate!-erule); symbolic procedure ps!:int!-erule(a,n); if rand2 a=ps!:depvar ps then if n=0 then simp!* psintconst else if ps!:expansion!-point ps neq 'ps!:inf then quotsq(ps!:evaluate(rand1 a,n-1),n ./ 1) else quotsq(ps!:evaluate(rand1 a,n+1),-n ./ 1) else simpint list(prepsq ps!:evaluate(rand1 a,n),rand2 a); put('int,'ps!:erule,'ps!:int!-erule); endmodule; module tpsdom; % Domain definitions for truncated power series. % Authors: Julian Padget & Alan Barnes. fluid '(ps!:exp!-lim ps!:max!-order); global '(domainlist!*); symbolic (domainlist!*:=union('(!:ps!:),domainlist!*)); % symbolic here seems to be essential in Cambridge Lisp systems put('tps,'tag,'!:ps!:); put('!:ps!:,'dname,'tps); flag('(!:ps!:),'field); put('!:ps!:,'i2d,'i2ps); put('!:ps!:,'minusp,'ps!:minusp!:); put('!:ps!:,'plus,'ps!:plus!:); put('!:ps!:,'times,'ps!:times!:); put('!:ps!:,'difference,'ps!:difference!:); put('!:ps!:,'quotient,'ps!:quotient!:); put('!:ps!:,'zerop,'ps!:zerop!:); put('!:ps!:,'onep,'ps!:onep!:); put('!:ps!:,'prepfn,'ps!:prepfn!:); put('!:ps!:,'specprn,'ps!:prin!:); put('!:ps!:,'prifn,'ps!:prin!:); put('!:ps!:,'intequivfn,'psintequiv!:); % conversion functions put('!:ps!:,'!:mod!:,mkdmoderr('!:ps!:,'!:mod!:)); % put('!:ps!:,'!:gi!:,mkdmoderr('!:ps!:,'!:gi!:)); % put('!:ps!:,'!:rn!:,mkdmoderr('!:ps!:,'!:rn!:)); put('!:rn!:,'!:ps!:,'!*d2ps); put('!:ft!:,'!:ps!:,'cdr); put('!:gi!:,'!:ps!:,'!*d2ps); put('!:gf!:,'!:ps!:,'!*d2ps); symbolic procedure psintequiv!: u; if idp cdr u or ps!:depvar u or length (u:=ps!:expression u)>2 then nil else int!-equiv!-chk rand1 u; symbolic procedure i2ps u; if dmode!*='!:ps!: then !*d2ps u else u; symbolic procedure !*d2ps u; begin scalar ps; ps:=get('tps,'tag) . mkvect 6; ps!:set!-order(ps,0); ps!:set!-expression(ps,list ('psconstant,u)); ps!:set!-value(ps,prepsq( u ./ 1)); ps!:set!-last!-term(ps,ps!:max!-order); ps!:set!-terms(ps,list ( 0 . ( u ./ 1))); return ps end; %symbolic procedure ps!:compatiblep(u,v); % if (constantpsp u or constantpsp v ) then t % else if (ps!:depvar u) neq (ps!:depvar v) then nil % then nil % else if (ps!:expansion!-point u) neq (ps!:expansion!-point v) % else t; symbolic procedure ps!:minusp!: u; nil; symbolic procedure ps!:plus!:(u,v); ps!:operator!:('plus,u,v); symbolic procedure ps!:difference!:(u,v); ps!:operator!:('difference,u,v); symbolic procedure ps!:times!:(u,v); ps!:operator!:('times,u,v); symbolic procedure ps!:quotient!:(u,v); ps!:operator!:('quotient,u,v); symbolic procedure ps!:diff!:(u,v); (( if idp deriv then ps!:compile (deriv,ps!:depvar u,ps!:expansion!-point u) else if numberp deriv then if zerop deriv then nil else deriv else make!-ps(list('df,u,v), deriv, ps!:depvar u,ps!:expansion!-point u)) ./ 1) where (deriv = prepsq simp!* list('df, ps!:value u,v)); symbolic procedure ps!:operator!:(op,u,v); begin scalar value,x,x0,y; x:=ps!:depvar u; y:=ps!:depvar v; if null x then if null y then return << if ps!:p u then u:=rand1 ps!:expression u; if ps!:p v then v:=rand1 ps!:expression v; if op='quotient and atom u and atom v then if null u then nil else <<op:=gcdn(u,v);u:=u/op;v:=v/op; if v=1 then u else '!:rn!: . (u . v)>> else dcombine!*(u,v,op)>> else << x:= y; x0:=ps!:expansion!-point v>> else if null y then x0:=ps!:expansion!-point u else if ((x0:=ps!:expansion!-point u) neq ps!:expansion!-point v) or (x neq y) then rerror(tps,20, list("ps!:operator: incompatible power series in ", op)); value:= simp!* list(op,ps!:value u,ps!:value v); if denr value=1 and domainp numr value then return numr value; return make!-ps(list(op,u,v), prepsq value,x,x0) end; symbolic procedure ps!:zerop!: u; (numberp v and zerop v) where v=ps!:value u; symbolic procedure ps!:onep!: u; onep ps!:value u; symbolic procedure ps!:prepfn!: u; u; initdmode 'tps; endmodule; 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 PSGEN 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,21,"Wrong number of arguments to PSREVERSE"); symbolic procedure simppsrev1(series); begin scalar rev,psord, depvar,about; % if not kernp simp!* revvar then % typerr(revvar,"kernel: simppsrev"); series:=prepsq simp!* series; if not ps!:p series then rerror(tps,22, "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!:term(series,1) neq (nil ./ 1)) then about := prepsq ps!:get!-term(series,0) else if psord =-1 then about:='ps!:inf else rerror(tps,23,"Series cannot be inverted: simppsrev"); rev:=ps!:compile(list('psrev,series),depvar,about); if ps!:expansion!-point series = 'ps!:inf then return make!-ps(list('quotient,1,rev), ps!:value rev,depvar,about) ./ 1 else return rev ./ 1 end; symbolic procedure ps!:generating!-series(a,psord,inverted); begin scalar ps, depvar,point; depvar:=ps!:depvar a; point:= ps!:expansion!-point a; ps:=make!-ps(list('psgen, a,inverted),ps!:value a, depvar, point); ps!:set!-order(ps,psord); ps!:set!-last!-term(ps,psord); a:=ps!:compile(list('expt,a,if inverted then -psord else psord), depvar,point); ps!:find!-order a; ps!:set!-term(ps,psord,a); 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; rthpow:=ps!:compile(list('expt,series,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); ps!:set!-last!-term(genseries,r)>>; term:=ps!:get!-term(series,n); if null term then for j:=ps!:last!-term(series)+1:n do term:= ps!:evaluate(series,j); 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; return 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) end; symbolic procedure ps!:remove!-constant(ps); << ps:=ps!:compile(list('difference,ps,prepsq ps!:term(ps,0)), ps!:depvar ps, ps!:expansion!-point ps); ps!:find!-order ps; ps >>; % symbolic procedure ps!:reciprocal(ps); % << ps:=ps!:compile(list('quotient,1,ps), % ps!:depvar ps, % ps!:expansion!-point ps); % ps!:find!-order ps; % ps >>; put('psrev,'ps!:erule,'ps!:rev!-erule); put('psrev,'ps!:order!-fn,'ps!:rev!-orderfn); symbolic procedure ps!:rev!-orderfn u; if (u:=ps!:expansion!-point ps!:get!-rthpow(rand1 ps!:expression ps,1))=0 then 1 else if u = 'ps!:inf then 1 else 0; 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(u,v); if u=0 then 0 else ps!:find!-order(ps!:get!-rthpow(rand2 ps!:expression ps,u)); 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; return make!-ps(list('pscomp,series1, ps!:generating!-series(series2, ps!:order series1, if n1='ps!:inf then t else nil)), ps!:arg!-values a,d,n) 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; % if n=0 then return ps!:term(aa,0); 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,24, "Args should be <POWER SERIES>,<POWER SERIES>: simppscomp"); symbolic procedure simppscomp1(ps1,ps2); begin scalar x,d,n1,n; ps1:=prepsq 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:=prepsq simp!* ps2) then rerror(tps,25, "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=prepsq ps!:term(ps2,0)) then return ps!:compile(list('pscomp,ps1,ps2),d,n) ./ 1 else rerror(tps,26,"Series cannot be composed: simppscomp"); end; algebraic operator psrev,pscomp; endmodule; module tpsfns; % Expansion of elementary functions as power series using DOMAINVALCHK % TPS & NUMVAL must be on for these functions to be invoked % Example sin a where a is a power series will now be expanded % Author: Alan Barnes, March 1989 fluid '(!*numval); put('exp, '!:ps!:, 'ps!:exp!:); put('log, '!:ps!:, 'ps!:log!:); put('sin, '!:ps!:, 'ps!:sin!:); put('cos, '!:ps!:, 'ps!:cos!:); put('tan, '!:ps!:, 'ps!:tan!:); put('asin, '!:ps!:, 'ps!:asin!:); put('acos, '!:ps!:, 'ps!:acos!:); put('atan, '!:ps!:, 'ps!:atan!:); put('sinh, '!:ps!:, 'ps!:sinh!:); put('cosh, '!:ps!:, 'ps!:cosh!:); put('tanh, '!:ps!:, 'ps!:tanh!:); put('asinh, '!:ps!:, 'ps!:asinh!:); put('acosh, '!:ps!:, 'ps!:acosh!:); put('atanh, '!:ps!:, 'ps!:atanh!:); put('expt, '!:ps!:, 'ps!:expt!:); % the above is grotty but necessary as unfortunately DOMAINVALCHK % passes arglist of sin (rather than sin . arglist) to ps!:sin!: etc symbolic procedure ps!:expt!:(base,exp); begin scalar !*numval,depvar,about; % NB binding of !*numval avoids infinite loop depvar:= ps!:depvar base; if null depvar then << depvar:=ps!:depvar exp; about:= ps!:expansion!-point exp>> else about:= ps!:expansion!-point base; if null depvar then error(999,"constantps**constantps formed"); return ps!:compile(list('expt, base,exp),depvar,about) end; symbolic procedure ps!:unary!:fn(fn, arg); begin scalar !*numval; % NB binding of !*numval avoids infinite loop return ps!:compile(list(fn, arg), ps!:depvar arg, ps!:expansion!-point arg) end; symbolic procedure ps!:cos!: arg; ps!:unary!:fn('cos,arg); symbolic procedure ps!:sin!: arg; ps!:unary!:fn('sin,arg); symbolic procedure ps!:tan!: arg; ps!:unary!:fn('tan,arg); symbolic procedure ps!:log!: arg; ps!:unary!:fn('log,arg); symbolic procedure ps!:exp!: arg; ps!:unary!:fn('exp,arg); symbolic procedure ps!:cosh!: arg; ps!:unary!:fn('cosh,arg); symbolic procedure ps!:sinh!: arg; ps!:unary!:fn('sinh,arg); symbolic procedure ps!:tanh!: arg; ps!:unary!:fn('tanh,arg); symbolic procedure ps!:asin!: arg; ps!:unary!:fn('asin,arg); symbolic procedure ps!:acos!: arg; ps!:unary!:fn('acos,arg); symbolic procedure ps!:atan!: arg; ps!:unary!:fn('atan,arg); symbolic procedure ps!:asinh!: arg; ps!:unary!:fn('asinh,arg); symbolic procedure ps!:acosh!: arg; ps!:unary!:fn('acosh,arg); symbolic procedure ps!:atanh!: arg; ps!:unary!:fn('atanh,arg); endmodule; module tpselfns; % Specific compilation and evaluation rules for elementary functions % Automatically generated recurrence relations can sometimes lead to % infinite loops. % Also helps in the detection of branch point singularities % Author: Alan Barnes. March 1989 fluid '(ps!:max!-order ps); put('sqrt,'ps!:crule,'ps!:unary!-crule); put('sqrt,'ps!:order!-fn,'ps!:sqrt!-orderfn); put('sqrt,'ps!:erule,'ps!:sqrt!-erule); symbolic procedure ps!:sqrt!-orderfn u; (if v*2=u then v else rerror(tps,27,"Branch Point in Sqrt")) where v=u/2; symbolic procedure ps!:sqrt!-erule(a,n); begin scalar aa,x,y,z; aa:=rand1 a; z:= nil ./ 1; y:=ps!:order aa; x:=ps!:order(ps); %order of sqrt ps if n=x then return simpexpt(list(prepsq ps!:evaluate(aa,y), '(quotient 1 2))); for k:=1:n-x do z:=addsq(z, multsq(((lambda y; if y=0 then nil else y) (k*3-2*n+y)) ./ 1, multsq(ps!:evaluate(aa,k+y), ps!:evaluate(ps,n-k)))); return quotsq(z,multsq(2*(n-x) ./ 1,ps!:evaluate(aa,y))) end; % alternative algorithm (for order 0 only) % for i:=1:n-1 do % z:=addsq(z,multsq(multsq( i ./ 1,ps!:evaluate(ps,i)), % ps!:evaluate(ps,n-i))); % z:=multsq(z, 1 ./ (n+1)); % return quotsq(addsq(ps!:evaluate(aa,n),negsq z), % multsq(2 ./ 1,ps!:evaluate(b,x))) symbolic procedure ps!:expt!-erule(a,n); begin scalar base,exp,x,y,z,p,q; base:= rand1 a; if length a=3 then <<exp:=rand2 a;p:=exp;q:=1>> else << p:=rand2 a; q:=cadddr a; exp:=list('quotient,p,q)>>; % exp:=ps!:value rand2 a; % exponent is always p or (QUOTIENT p q) with p,q integers x:= nil ./ 1; y:=ps!:order(base); z:= ps!:order ps; % order of exponential if n=z then return simpexpt(list(prepsq ps!:evaluate(base,y),exp)) else <<for k:=1:n-z do x:=addsq(x, multsq(((lambda num; if num=0 then nil else num) (k*p+q*(k-n+z))) ./ q, multsq(ps!:evaluate(base,k+y), ps!:evaluate(ps,n-k)))); return quotsq(x,multsq((n-z) ./ 1,ps!:evaluate(base,y))) >>; end; put('expt,'ps!:erule, 'ps!:expt!-erule); put('expt,'ps!:crule,'ps!:expt!-crule); put('expt,'ps!:order!-fn,'ps!:expt!-orderfn); symbolic procedure ps!:expt!-crule(a,d,n); % we will assume that forms like (expt (expt <x> <y>) <z>) will % continue to be transformed by SIMP!* into (expt <x> (times <y> <z>)) % this is very important for the avoidance of essential singularities begin scalar exp0,eflg,exp1,exp2,b,ps,psvalue,dmode,!*msg; dmode := dmode!*; exp0:=rand1 a; eflg := evalequal(exp0,prepsq simp!* aeval 'e); if dmode then onoff(get(dmode,'dname),nil); exp1:=rand2 a; if (ps!:p exp1 and constantpsp exp1) then exp1:=ps!:value exp1; begin scalar alglist!*; exp2:=simp!* exp1; exp1:=prepsq exp2 end; if (atom numr exp2 and atom denr exp2) then <<exp1:=numr exp2; exp2:=denr exp2>> else return <<ps:=ps!:unknown!-crule(list('expt, 'e, ps!:compile(if eflg then exp1 else list('times, exp1, list('log,exp0)), d,n)), d,n); if dmode then onoff(get(dmode,'dname),t); ps>>; b:=ps!:compile(exp0,d,n); if dmode then onoff(get(dmode,'dname),t); psvalue:=ps!:arg!-values a; if exp2=1 then if exp1=nil then return if ps!:zerop!: b then rerror(tps,28,"0**0 formed: ps:expt-crule") else 1 else if exp1=1 then return b else if exp1=2 then return make!-ps(list('times,b,b),psvalue,d,n) else if exp1 = -1 then return make!-ps(list('quotient,1,b),psvalue,d,n) else return make!-ps(list('expt,b,exp1),psvalue,d,n) else return make!-ps(list('exp3,b,exp1,exp2),psvalue,d,n); end; symbolic procedure ps!:expt!-orderfn(u,v); u*rand2 ps!:expression ps; symbolic procedure ps!:exp3!-orderfn(u,v,w); begin scalar expres; expres:=ps!:expression ps; v:= rand2 expres; w:=cadddr expres; if cdr(v:=divide(u * v,w))=0 then return car v else rerror(tps,29,"Branch Point in EXPT") end; put('exp3,'ps!:order!-fn,'ps!:exp3!-orderfn); put('exp3,'ps!:erule,'ps!:expt!-erule); endmodule; module tpssum; % Written by Alan Barnes. September 1990 % Allows power series whose general term is given to be manipulated. % % pssum(<sumvar>=<lowlim>, <coeff>, <depvar>, <about>, <power>); % % <sumvar> summation variable (a kernel) % <lowlim> lower limit of summation (an integer) % <coeff> general coefficient of power series (algebraic) % <depvar> expansion variable of series (a kernel) % <about> expansion point of series (algebraic) % <power> general exponent of power series (algebraic) % <power> must be a strictly increasing function of <sumvar> % this is now partially checked by the system symbolic procedure ps!:summation!-erule(a,n); begin scalar power, coeff,sumvar,current!-index,last!-exp,current!-exp; current!-index:= rand2 a; sumvar:= rand1 a; coeff := cdddr a; power:= cadr coeff; coeff:=car coeff; last!-exp:= ieval reval subst(current!-index,sumvar,power); repeat << current!-index:=current!-index+1; current!-exp:= ieval reval subst(current!-index,sumvar,power); if current!-exp leq last!-exp then rerror(tps,30,"Exponent not strictly increasing: ps:summation"); if current!-exp < n then << ps!:set!-term(ps,current!-exp, simp!* subst(current!-index,sumvar,coeff)); ps!:set!-last!-term(ps,current!-exp); rplaca(cddr a,current!-index)>>; last!-exp:=current!-exp>> until current!-exp geq n; return if current!-exp = n then << rplaca(cddr a,current!-index); simp!* subst(current!-index,sumvar,coeff) >> else (nil ./ 1) end; put('ps!:summation, 'ps!:erule, 'ps!:summation!-erule); put('ps!:summation, 'simpfn, 'simpiden); put('pssum, 'simpfn, 'simppssum); symbolic procedure simppssum a; begin scalar !*nosubs,from,sumvar,lowlim,coeff, power,depvar,about,psord,term; if length a neq 5 then rerror(tps,31, "Args should be <FROM>,<coeff>,<depvar>,<point>,<power>: simppssum"); !*nosubs := t; % We don't want left side of eqns to change. from := reval car a; !*nosubs := nil; if not eqexpr from then errpri2(car a,t) else <<sumvar:= cadr from; if not kernp simp!* sumvar then typerr(sumvar, "kernel: simppssum"); lowlim:= ieval caddr from >>; coeff:= prepsq simp!* cadr a; a:= cddr a; depvar := car a; about:=prepsq simp!* cadr a; if about = 'infinity then about := 'ps!:inf; power:= prepsq simp!* caddr a; if not kernp simp!* depvar then typerr(depvar, "kernel: simppssum") else if depvar=sumvar then rerror(tps,32, "Summation and expansion variables are the same: simppssum") else if smember(depvar,about) then rerror(tps,33,"Expansion point depends on depvar: simppssum") else if smember(sumvar,about) then rerror(tps,34, "Expansion point depends on summation var: simppssum") else if not smember(sumvar,power) then rerror(tps,35, "Exponent does not depend on summation variable: simppssum"); lowlim:=lowlim-1; repeat << lowlim:=lowlim+1; psord:= ieval reval subst(lowlim,sumvar,power)>> until (term:=simp!* subst(lowlim,sumvar,coeff)) neq '(nil . 1); ps:=make!-ps(list('ps!:summation,sumvar,lowlim,coeff,power), list('ps!:summation,from,coeff,depvar,about,power), depvar, about); ps!:set!-order(ps,psord); ps!:set!-term(ps,psord, term); ps!:set!-last!-term(ps,psord); return (ps ./ 1) end; endmodule; end;