reduce-historical/r34/src/tps.red

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;