123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602 |
- module tps; % Truncated Power Series.
- % Authors: Julian Padget & Alan Barnes <barnesa@kirk.aston.ac.uk>.
- % Version 1.3 September 1990.
- %
- % Revisions:
- %
- % 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);
- % 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) . mapcar(ps!:operands ps,
- 'simp!:ps1));
- 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),
- mapcar(ps!:operands ps,
- 'ps!:find!-order));
- 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:=15);
- % 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!:,'!:bf!:,mkdmoderr('!:ps!:,'!:bf!:));
- % put('!:ps!:,'!:rn!:,mkdmoderr('!:ps!:,'!:rn!:));
- put('!:rn!:,'!:ps!:,'!*d2ps);
- put('!:ft!:,'!:ps!:,'cdr);
- put('!:bf!:,'!:ps!:,'!*d2ps);
- 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 exp1,exp2,b,psvalue;
- exp1:=rand2 a;
- if (ps!:p exp1 and constantpsp exp1) then exp1:=ps!:value exp1;
- begin scalar dmode!*, 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!:unknown!-crule(list('expt,'e,
- ps!:compile(if rand1 a='e then exp1
- else list('times, exp1,
- list('log,rand1 a)),
- d,n)),
- d,n);
- b:=ps!:compile(rand1 a,d,n);
- 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;
|