reduce-historical/r38/packages/taylor/tayfns.red

module TayFns;

%*****************************************************************
%
%       Simplification functions for special functions
%
%*****************************************************************


exports taysimpexpt, taysimpatan, taysimplog, taysimpexp,
        taysimptan, taysimpsin, taysimpsinh, taysimpasin;

imports

% from the REDUCE kernel: 
        !*f2q, !:minusp, addsq, aeval, denr, domainp, eqcar, evenp,
        freeof, invsq, kernp, lastpair, let, lprim, lnc, mk!*sq, mksq,
        multsq, mvar, negsq, neq, nlist, nth, numr, over, prepd,
        prepsq, quotsq, retimes, reval, reversip, simp, simp!*,
        simplogi, simplogsq, subs2!*, subsq, subtrsq,

% from the header module:
        !*tay2q, !*TayExp2q, constant!-sq!-p, cst!-Taylor!*,
        find!-non!-zero, get!-degree, has!-TayVars,
        make!-cst!-powerlist, make!-Taylor!*, prune!-coefflist,
        set!-TayCoeffList, set!-TayFlags, set!-TayOrig, TayCfPl,
        TayCfSq, TayCoeffList, TayFlags, TayGetCoeff, Taylor!*p,
        Taylor!-kernel!-sq!-p, Taylor!:, TayMakeCoeff, TayOrig,
        TayTemplate, TayTpElNext, TayTpElOrder, TayTpElPoint,
        TayTpElVars, TayTpVars, TayVars, TpNextList,

% from the module Tayintro:
        confusion, delete!-nth, delete!-nth!-nth, replace!-nth,
        replace!-nth!-nth, Taylor!-error, Taylor!-error!*,
        var!-is!-nth,

% from the module Tayutils:
        addto!-all!-TayTpElORders, get!-cst!-coeff, is!-neg!-pl,
        smallest!-increment, subtr!-degrees, Taylor!*!-constantp,
        Taylor!*!-zerop,

% from the module Taybasic:
        addtaylor, addtaylor!-as!-sq, invtaylor, makecoeffs0,
        makecoeffpairs, makecoeffpairs1, multtaylor,
        multtaylor!-as!-sq, multtaylorsq, negtaylor, negtaylor1,
        quottaylor,

% from the module TayExpnd:
        taylorexpand,

% from the module Taysimp:
        expttayrat, taysimpsq, taysimpsq!*,

% from the module Taydiff:
        difftaylorwrttayvar,

% from the module TayConv:
        prepTayCoeff, prepTaylor!*,

% from the module Tayfrontend:
        taylorcombine, taylortostandard;


fluid '(!*taylorkeeporiginal !*!*taylor!-epsilon!*!* frlis!*);


symbolic procedure taysimpexpt u;
  %
  % Argument is of the form ('expt base exponent)
  % where both base and exponent (but a least one of them)
  % may contain Taylor kernels given as prefix forms.
  % Value is the equivalent Taylor kernel.
  %
  if not (car u eq 'expt) or cdddr u then confusion 'taysimpexpt
   else if cadr u eq 'e then taysimpexp {'exp, caddr u}
   else begin scalar bas, expn;
     bas := taysimpsq simp!* cadr u;
     expn := taysimpsq simp!* caddr u;
     if constant!-sq!-p bas then return taysimpexp
         {'exp,mk!*sq multsq(simp!*{'log,mk!*sq bas},expn)};
     if null kernp bas
       then if not(denr bas = 1)
              then return mksq({'expt,prepsq bas,prepsq expn},1)
             else if domainp numr bas
              then return taysimpexp {'exp,
                  prepsq multsq(simp!* {'log,prepd numr bas},expn)}
             else return mksq({'expt,prepsq bas,prepsq expn},1);
     if fixp numr expn and fixp denr expn
       then return !*tay2q expttayrat(mvar numr bas,expn);
     if denr expn = 1 and eqcar(numr expn,'!:rn!:)
       then return !*tay2q expttayrat(mvar numr bas,cdr numr expn);
     bas := mvar numr bas;
     return
       if Taylor!*p bas
         then if Taylor!-kernel!-sq!-p expn
                then if TayTemplate bas = TayTemplate mvar numr expn
                       then taysimpexp {'exp,
                                        mk!*sq taysimpsq
                                          multtaylor!-as!-sq(
                                              expn,
                                              taysimplog {'log,bas})}
                      else mksq({'expt,bas,mvar numr expn},1)
               else if not has!-TayVars(bas,expn)
                then begin scalar logterm;
                  logterm := taysimplog{'log,bas};
                  return
                    if Taylor!-kernel!-sq!-p logterm
                      then taysimpexp{'exp,
                                      multtaylorsq(mvar numr logterm,
                                                   expn)}
                     else taysimpsq
                            simp!* {'exp,mk!*sq multsq(logterm,expn)}
                 end
               else mksq({'expt,bas,mk!*sq expn},1)
        else if Taylor!-kernel!-sq!-p expn
               then if not has!-TayVars(mvar numr expn,bas)
                      then taysimpexp{'exp,
                                      multtaylorsq(mvar numr expn,
                                                   simp!*{'log,bas})}
                     else if Taylor!*!-constantp mvar numr expn
                      then taylorexpand(
                             simp!* {'expt,bas,
                                     prepTaylor!* mvar numr expn},
                             TayTemplate mvar numr expn)
                     else mksq({'expt,bas,mk!*sq expn},1)
        else mksq({'expt,bas,mk!*sq expn},1);
  end;

put('expt,'taylorsimpfn,'taysimpexpt);


symbolic procedure TayCoeffList!-union u;
  if null cdr u then car u
   else TayCoeffList!-union2 (car u, TayCoeffList!-union cdr u)$

symbolic procedure TayCoeffList!-union2 (x, y);
  %
  % returns union of TayCoeffLists x and y
  %
  << for each w in y do
       if null (assoc (car w, x)) then x := w . x;
     x >>$

symbolic procedure inttaylorwrttayvar(tay,var);
  %
  % integrates Taylor kernel tay wrt variable var
  %
  inttaylorwrttayvar1(TayCoeffList tay,TayTemplate tay,var)$

symbolic procedure inttaylorwrttayvar1(tcl,tp,var);
  %
  % integrates Taylor kernel with TayCoeffList tcl and template tp
  %  wrt variable var
  %
  Taylor!:
  begin scalar tt,u,w,singlist,csing; integer n,n1,m;
    u := var!-is!-nth(tp,var);
    n := car u;
    n1 := cdr u;
    tt := nth(tp,n);
    u := for each cc in tcl join <<
           m := nth(nth(TayCfPl cc,n),n1);
           if TayTpElPoint nth(tp,n) eq 'infinity
             then <<
               if m=1 then <<singlist :=
                               TayMakeCoeff(
                                 delete!-nth!-nth(TayCfPl cc,n,n1),
                                 TayCfSq cc) . singlist;nil>>
                else {TayMakeCoeff(
                        replace!-nth!-nth(TayCfPl cc,n,n1,m-1),
                        multsq(TayCfSq cc,invsq !*TayExp2q(-m + 1)))}>>
            else <<
               if m=-1 then <<singlist :=
                                TayMakeCoeff(
                                  delete!-nth!-nth(TayCfPl cc,n,n1),
                                  TayCfSq cc) . singlist;nil>>
               else {TayMakeCoeff(
                       replace!-nth!-nth(TayCfPl cc,n,n1,m+1),
                       multsq(TayCfSq cc,invsq !*TayExp2q(m + 1)))}>>>>;
    w := {TayTpElVars tt,TayTpElPoint tt,
          if var member TayTpElVars tt
            then if TayTpElPoint tt eq 'infinity
                   then TayTpElOrder tt - 1
                  else TayTpElOrder tt + 1
           else TayTpElOrder tt,
          if var member TayTpElVars tt
            then if TayTpElPoint tt eq 'infinity
                   then TayTpElNext tt - 1
                  else TayTpElNext tt + 1
           else TayTpElOrder tt};
    if singlist then begin scalar tpel;
        tpel := nth(tp,n);
        singlist := reversip singlist;
        if TayCfPl car singlist = '(nil) % no Taylor left
          then csing := TayCfSq car singlist
         else csing := !*tay2q
                         make!-Taylor!*(
                           singlist,
                           replace!-nth(
                             tp,n,
                             {delete!-nth(TayTpElVars tpel,n1),
                              TayTpElPoint tpel,
                              TayTpElOrder tpel,
                              TayTpElNext tpel}),
                           nil,nil);
        csing := multsq(csing,simp!* {'log,nth(TayTpElVars tpel,n1)})
       end;

    return (csing . make!-Taylor!*(u,replace!-nth(tp,n,w),nil,nil))
%
% The following is not needed yet
%
%     return make!-Taylor!*(
%              u,
%              replace!-nth(TayTemplate tay,n,w),
%              if !*taylorkeeporiginal and TayOrig tay
%                then simp {'int,mk!*sq TayOrig tay,var}
%               else nil,
%              TayFlags u)
  end;


comment The inverse trigonometric and inverse hyperbolic functions
        of a Taylor kernel are calculated by first computing the
        derivative(s) with respect to the Taylor variable(s) and
        integrating the result.  The derivatives can easily be
        calculated by the manipulation functions defined above.

        The method is best illustrated with an example.  Let T(x)
        be a Taylor kernel depending on one variable x.  To compute
        the inverse tangent T1(x) = atan(T(x)) we calculate the
        derivative

                                T'(x)
                    T1'(x) = ----------- .
                                      2
                              1 + T(x)

        (If T and T1 depend on more than one variable replace
        the derivatives by gradients.)
        This is integrated again with the integration constant
        T1(x0) = atan(T(x0)) yielding the desired result.
        If there is more than variable we have to find the
        potential function T1(x1,...,xn) corresponding to
        the vector grad T1(x1,...,xn) which is always possible
        by construction.

        The prescriptions for the eight functions asin, acos, asec,
        acsc, asinh, acosh, asech, and acsch can be put together
        in one procedure since the expressions for their derivatives
        differ only in certain signs.  The same remark applies to
        the four functions atan, acot, atanh, and acoth.

        These expressions are:

         d                 1
         -- asin x = ------------- ,
         dx           sqrt(1-x^2)

         d                -1
         -- acos x = ------------- ,
         dx           sqrt(1-x^2)

         d                 1
         -- asinh x = ------------- ,
         dx            sqrt(1+x^2)

         d                 1
         -- acosh x = ------------- ,
         dx            sqrt(x^2-1)

         d               1
         -- atan x = --------- ,
         dx           1 + x^2

         d               -1
         -- acot x = --------- ,
         dx           1 + x^2

         d                1
         -- atanh x = --------- ,
         dx            1 - x^2

         d                1
         -- acoth x = --------- ,
         dx            1 - x^2

        together with the relations

                       1
         asec x = acos - ,
                       x

                       1
         acsc x = asin - ,
                       x

                         1
         asech x = acosh - ,
                         x

                         1
         acsch x = asinh -
                         x .


        This method has one drawback: it is applicable only when T(x0)
        is a regular point of the function in question. E.g., if T(x0)
        = 0, then asech(T(x)) cannot be calculated by this method, as
        asech has a logarithmic singularity at 0. This means that the
        constant term of the series cannot be determined by computing
        asech(T(x0)).  In that case, we use the following relations
        instead:


         asin z = -i log(i z + sqrt(1 - z^2)),


         acos z = -i log(z + sqrt(z^2 - 1)),

                    1        1 + i z
         atan z = ----- log ---------,
                   2 i       1 - i z

                   -1        i z + 1
         acot z = ----- log ---------,
                   2 i       i z - 1


         asinh z = log(z + sqrt(1 + z^2)),


         acosh z = log(z + sqrt(z^2 - 1)),

                    1       1 + z
         atanh z = --- log -------,
                    2       1 - z

                    1       z + 1
         acoth z = --- log -------.
                    2       z - 1


         These formulas are applied at the following points:

           infinity for all functions,

           +i/-i for atan and acot,

           +1/-1 for atanh and acoth.


         There are still some branch points, where the calculation is
         not always possible:

           +1/-1 for asin and acos, and consequently for asec and acsc,

           +i/-i for asinh, acosh, asech and acsch.

         In these cases, the above formulas are applied as well, but
         simplification of the square roots and the logarithm may lead
         to a rather complicated result;



symbolic procedure taysimpasin u;
  if not (car u memq '(asin acos acsc asec asinh acosh acsch asech))
     or cddr u
    then confusion 'taysimpasin
   else Taylor!:
     begin scalar l,l0,c0,v,tay0,tay,tay2,tp,singlist;
     tay0 := taysimpsq simp!* cadr u;
     if not Taylor!-kernel!-sq!-p tay0
       then return mksq({car u,mk!*sq tay0},1);
     tay0 := mvar numr tay0; % asin's argument
     l0 := make!-cst!-powerlist TayTemplate tay0;
     c0 := TayGetCoeff(l0,TayCoeffList tay0);
     if car u memq '(asec acsc asech acsch)
       then if null numr c0 then return taysimpasec!*(car u,tay0)
             else tay := invtaylor tay0
      else tay := tay0;
     tp := TayTemplate tay;
     l := prune!-coefflist TayCoeffList tay;
     if null l then return !*tay2q cst!-Taylor!*(simp!* {car u,0},tp); 
     if is!-neg!-pl TayCfPl car l then return taysimpasin!*(car u,tay);
     tay2 := multtaylor(tay,tay);
     if car u memq '(asin acos acsc asec)
       then tay2 := negtaylor tay2;
     tay2 := addtaylor(
               cst!-Taylor!*(
                 !*f2q(if car u memq '(acosh asech) then -1 else 1),
                 tp),
               tay2);
     if Taylor!*!-zerop tay2
       then Taylor!-error!*('branch!-point,car u)
      else if null numr TayGetCoeff(l0,TayCoeffList tay2)
       then return taysimpasin!*(car u,tay);

     tay2 := invtaylor expttayrat(tay2,1 ./ 2);
     if car u memq '(acos asec) then tay2 := negtaylor tay2;
     l := for each krnl in TayVars tay collect
            inttaylorwrttayvar(
              multtaylor(difftaylorwrttayvar(tay,krnl),tay2),
              krnl);
     v := TayCoeffList!-union
            for each pp in l collect TayCoeffList cdr pp;
     singlist := nil ./ 1;
     for each pp in l do
       if car pp then singlist := addsq(singlist,car pp);
     %
     % special treatment for zeroth coefficient
     %
     c0 := simp {car u,mk!*sq c0};
     v := TayMakeCoeff(l0,c0) . v;
     tay := make!-Taylor!*(
              v,
              tp,
              if !*taylorkeeporiginal and TayOrig tay
                then simp {car u,mk!*sq TayOrig tay}
               else nil,
              TayFlags tay);
     if null numr singlist then return !*tay2q tay;
     if !*taylorkeeporiginal and TayOrig tay
       then set!-TayOrig(tay,subtrsq(TayOrig tay,singlist));
     return addsq(singlist,!*tay2q tay)
  end;

symbolic procedure taysimpasec!*(fn,tay);
   begin scalar result,tay1,tay2,i1;
     i1 := simp 'i;
     if fn memq '(asin acsc) then tay := multtaylorsq(tay,i1);
     tay1 := cst!-Taylor!*(1 ./ 1,TayTemplate tay);
     tay2 := multtaylor(tay,tay);
     if fn memq '(asec asech) then tay2 := negtaylor tay2;
     result := taysimplog {'log,
                           addtaylor(
                             expttayrat(addtaylor(tay2,tay1),1 ./ 2),
                             tay1)};
     tay1 := taysimplog {'log,tay};
     if fn memq '(asin acos asec acsc)
       then <<result := taysimpsq negsq multsq(result,i1);
              result := addtaylor!-as!-sq(result,multsq(i1,tay1))>>
      else result := addtaylor!-as!-sq(result,
                                       negsq taysimplog {'log,tay});
     return taysimpsq!* result
   end;

symbolic procedure taysimpasin!*(fn,tay);
   begin scalar result,tay1;
     if fn memq '(asin acsc)
       then tay := multtaylorsq(tay,simp 'i);
     tay1 := cst!-Taylor!*(
               (if fn memq '(asin asinh acsc acsch)
                  then 1
                 else -1) ./ 1,
                TayTemplate tay);
     result := taysimplog {'log,
                           addtaylor(
                             expttayrat(addtaylor(multtaylor(tay,tay),
                                                  tay1),
                                        1 ./ 2),
                             tay)};
     if fn memq '(asin acos asec acsc)
       then result := quotsq(result,simp 'i);
     return taysimpsq!* result
   end;


put('asin,'taylorsimpfn,'taysimpasin);
put('acos,'taylorsimpfn,'taysimpasin);
put('asec,'taylorsimpfn,'taysimpasin);
put('acsc,'taylorsimpfn,'taysimpasin);
put('asinh,'taylorsimpfn,'taysimpasin);
put('acosh,'taylorsimpfn,'taysimpasin);
put('asech,'taylorsimpfn,'taysimpasin);
put('acsch,'taylorsimpfn,'taysimpasin);


symbolic procedure taysimpatan u;
  if not (car u memq '(atan acot atanh acoth)) or cddr u
    then confusion 'taysimpatan
   else begin scalar l,l0,c0,v,tay,tay2,tp,singlist;
     tay := taysimpsq simp!* cadr u;
     if not Taylor!-kernel!-sq!-p tay
       then return mksq({car u,mk!*sq tay},1);
     tay := mvar numr tay; % atan's argument
     tp := TayTemplate tay;
     l0 := make!-cst!-powerlist tp;
     l := prune!-coefflist TayCoeffList tay;
     if null l then return !*tay2q cst!-Taylor!*(simp!* {car u,0},tp); 
     if is!-neg!-pl TayCfPl car l then return taysimpatan!*(car u,tay);
     c0 := get!-cst!-coeff tay;
     if car u memq '(atan acot)
       then c0 := subs2!* multsq(c0,simp 'i);
     if c0 = (1 ./ 1) or c0 = (-1 ./ 1)
       then return taysimpatan!*(car u,tay);
     tay2 := multtaylor(tay,tay);
     if car u memq '(atanh acoth) then tay2 := negtaylor tay2;
     tay2 := invtaylor addtaylor(cst!-Taylor!*(1 ./ 1,tp),tay2);
     if car u eq 'acot then tay2 := negtaylor tay2;
     l := for each krnl in TayVars tay collect
            inttaylorwrttayvar(
              multtaylor(difftaylorwrttayvar(tay,krnl),tay2),
              krnl);
     v := TayCoeffList!-union
            for each pp in l collect TayCoeffList cdr pp;
     singlist := nil ./ 1;
     for each pp in l do
       if car pp then singlist := addsq(singlist,car pp);
     %
     % special treatment for zeroth coefficient
     %
     c0 := simp {car u,
                 mk!*sq TayGetCoeff(l0,TayCoeffList tay)};
     v := TayMakeCoeff (l0,c0) . v;
     tay := make!-Taylor!*(
              v,
              tp,
              if !*taylorkeeporiginal and TayOrig tay
                then simp {car u,mk!*sq TayOrig tay}
               else nil,
              TayFlags tay);
     if null numr singlist then return !*tay2q tay;
     if !*taylorkeeporiginal and TayOrig tay
       then set!-TayOrig(tay,subtrsq(TayOrig tay,singlist));
     return addsq(singlist,!*tay2q tay)
  end;

symbolic procedure taysimpatan!*(fn,tay);
   begin scalar result,tay1;
     if fn memq '(atan acot)
       then tay := multtaylorsq(tay,simp 'i);
     tay1 := cst!-Taylor!*(1 ./ 1,TayTemplate tay);
     tay := quottaylor(addtaylor(tay1,tay),
                       addtaylor(tay1,negtaylor tay));
     result := multsq(taysimplog {'log,tay},1 ./ 2);
     if fn eq 'atan then result := quotsq(result,simp 'i)
      else if fn eq 'acot then result := multsq(result,simp 'i);
     return taysimpsq!* result
   end;



put('atan,'taylorsimpfn,'taysimpatan);
put('acot,'taylorsimpfn,'taysimpatan);
put('atanh,'taylorsimpfn,'taysimpatan);
put('acoth,'taylorsimpfn,'taysimpatan);


comment For the logarithm and exponential we use the extension of
        an algorithm quoted by Knuth who shows how to do this for
        series in one expansion variable.

        We extended this to the case of several variables which is
        straightforward except for one point, see below.
        Knuth's algorithm works as follows: Assume you want to compute
        the series W(x) where

            W(x) = log V(x)

        Differentiation of this equation gives

                    V'(x)
            W'(x) = ----- ,   or V'(x) = W'(x)V(x) .
                     V(x)

        You make now the ansatz

                    -----
                    \           n
            W(x) =   >      W  x  ,
                    /        n
                    -----

        substitute this into the above equation and compare
        powers of x.  This yields the recursion formula

                               n-1
                 V            -----
                  n       1   \
           W  = ---- - ------  >    m W  V     .
            n    V      n V   /        m  n-m
                  0        0  -----
                               m=0

        The first coefficient must be calculated directly, it is

           W   = log V  .
            0         0

        To use this for series in more than one variable you have to
        calculate all partial derivatives: n and m refer then to the
        corresponding component of the multi index.  Looking closely
        one finds that there is an ambiguity: the same coefficient
        can be calculated using any of the partial derivatives.  The
        only restriction is that the corresponding component of the
        multi index must not be zero, since we have to divide by it.

        We resolve this ambiguity by simply taking the first nonzero
        component.

        The case of the exponential is nearly the same: differentiation
        gives

            W'(x) = V'(x) W(x) ,

        from which we derive the recursion formula

                   n-1
                  -----
                  \     n-m
            W  =   >    --- W  V     , W  = exp V  .
             n    /      m   m  n-m     0        0
                  -----
                   m=0

        ;


symbolic procedure taysimplog u;
  %
  % Special Taylor expansion function for logarithms
  %
  if not (car u eq 'log) or cddr u then confusion 'taysimplog
   else Taylor!:
    begin scalar a0,clist,coefflis,il,l0,l,tay,tp,csing,singterm;
    u := simplogi cadr u;
    if not kernp u then return taysimpsq u;
    u := mvar numr u;
    if not (car u eq 'log) then confusion 'taysimplog;
    u := taysimpsq simp!* cadr u;
    if not Taylor!-kernel!-sq!-p u then return mksq({'log,mk!*sq u},1);
    tay := mvar numr u;
    tp := TayTemplate tay;
    l0 := make!-cst!-powerlist tp;
    %
    % The following relies on the standard ordering of the
    % TayCoeffList.
    %
    l := prune!-coefflist TayCoeffList tay;
    if null l then Taylor!-error!*('not!-a!-unit,'taysimplog);
    %
    % The assumption at this point is that the first term is the one
    % with lowest degree, i.e. dividing by this term yields a series
    % which starts with a constant term.
    %
    if TayCfPl car l neq l0 then
      <<csing := TayCfPl car l;
        l := for each pp in l collect begin scalar newpl;
                 newpl := subtr!-degrees(TayCfPl pp,csing);
                 if is!-neg!-pl newpl
                   then Taylor!-error('not!-a!-unit,'taysimplog)
                  else return TayMakeCoeff(newpl,TayCfSQ pp);
               end;
        tp := addto!-all!-TayTpElOrders(
                tp,
                for each nl in csing collect
                  - get!-degree nl);
        singterm := simp!* retimes preptaycoeff(csing,tp);
        if !:minusp lnc numr TayCfSq car l
          then <<singterm := negsq singterm;
                 l := negtaylor1 l>>>>;
    a0 := TayGetCoeff(l0,l);
    clist := {TayMakeCoeff(l0,simplogi mk!*sq a0)};
    il := if null l then nlist(1,length tp)
           else smallest!-increment l;
    coefflis := makecoeffs0(tp,TpNextList tp,il);
    if not null coefflis
      then for each cc in cdr coefflis do
             begin scalar s,pos,pp,n,n1;
               s := nil ./ 1;
               pos := find!-non!-zero cc;
               n := nth(nth(cc,car pos),cdr pos);
               pp := makecoeffpairs(l0,cc,TayCfPl car l,il);
               for each p in pp do <<
                 n1 := nth(nth(car p,car pos),cdr pos);
                 s := addsq(s,
                            multsq(!*TayExp2q n1,
                                   multsq(TayGetCoeff(car p,clist),
                                          TayGetCoeff(cdr p,l))))>>;
%               for each p in pp addsq
%                 multsq(!*TayExp2q nth(nth(car p,car pos),cdr pos),
%                        multsq(TayGetCoeff(car p,clist),
%                               TayGetCoeff(cdr p,l)));
               s := subtrsq(TayGetCoeff(cc,l),quotsq(s,!*TayExp2q n));
               if not null numr s
                 then clist := TayMakeCoeff(cc,quotsq(s,a0)) . clist;
             end;
    tay := make!-Taylor!*(
             reversip clist,
             tp,
             if !*taylorkeeporiginal and TayOrig tay
               then simplogi mk!*sq TayOrig tay
              else nil,
             TayFlags tay);
    if null csing then return !*tay2q tay;
    singterm := simplogsq singterm;
    if Taylor!*!-zerop tay then return singterm;
    if !*taylorkeeporiginal and TayOrig tay 
      then set!-TayOrig(tay,subtrsq(TayOrig tay,singterm));
    return addsq(singterm,!*tay2q tay)
  end;

put('log,'taylorsimpfn,'taysimplog);


symbolic procedure taysimpexp u;
  %
  % Special Taylor expansion function for exponentials
  %
  if not (car u eq 'exp) or cddr u then confusion 'taysimpexp
   else Taylor!:
    begin scalar a0,clist,coefflis,il,l0,l,lm,lp,tay,tp;
    u := taysimpsq simp!* cadr u;
    if not Taylor!-kernel!-sq!-p u
      then return mksq ({'exp,mk!*sq u},1);
    tay := mvar numr u;
    tp := TayTemplate tay;
    l0 := make!-cst!-powerlist tp;
    %
    % The following relies on the standard ordering of the
    % TayCoeffList.
    %
    l := prune!-coefflist TayCoeffList tay;
    if null l then return !*tay2q cst!-Taylor!*(1 ./ 1,tp);
    for each pp in l do
      if is!-neg!-pl TayCfPl pp then lm := pp . lm
       else if not null numr TayCfSq pp then lp := pp . lp;
    lm := reversip lm;
    l := reversip lp;

    if lm
      then lm := simp!* {'exp,
                         preptaylor!* make!-Taylor!*(lm,tp,nil,nil)};

    if null l then return lm;

    a0 := TayGetCoeff(l0,l);
    clist := {TayMakeCoeff(l0,simp!* {'exp,mk!*sq a0})};
    il := smallest!-increment l;
    coefflis := makecoeffs0(tp,TpNextList tp,il);

    if not null coefflis
      then for each cc in cdr coefflis do
             begin scalar s,pos,pp,n,n1;
               s := nil ./ 1;
               pos := find!-non!-zero cc;
               n := nth(nth(cc,car pos),cdr pos);
               pp := makecoeffpairs(l0,cc,l0,il);
               for each p in pp do <<
                 n1 := nth(nth(car p,car pos),cdr pos);
                 s := addsq(s,
                            multsq(!*TayExp2q(n - n1),
                                   multsq(TayGetCoeff(car p,clist),
                                          TayGetCoeff(cdr p,l))))>>;
               s := subs2!* quotsq(s,!*TayExp2q n);
               if not null numr s
                 then clist := TayMakeCoeff(cc,s) . clist
             end;

    clist := reversip clist;

    u := !*tay2q
           make!-Taylor!*(
             clist,
             tp,
             if !*taylorkeeporiginal and TayOrig tay
               then simp {'exp,mk!*sq TayOrig tay}
              else nil,
             TayFlags tay);
    return if null lm then u else multsq(u,lm)
  end;

put('exp,'taylorsimpfn,'taysimpexp);


comment The algorithm for the trigonometric functions is also
        derived from their differential equation.
        The simplest case is that of tangent whose equation is

                            2
           tan'(x) = 1 + tan (x) .          (*)

        For the others we have

                               2
           cot'(x) = - (1 + cot (x)),

                              2
           tanh'(x) = 1 - tanh (x),

                              2
           coth'(x) = 1 - coth (x) .



        Let T(x) be a Taylor series,

                  -----
                  \         N
           T(x) =  >    a  x
                  /      N
                  -----
                   N=0

        Now, let

                              -----
                              \         N
           T1(x) = tan T(x) =  >    b  x
                              /      N
                              -----
                               N=0

        from which we immediately deduce that b  = tan a .
                                               0        0

        From (*) we get
                              2
           T1'(x) = (1 + T1(x) ) T'(x) ,

        or, written in terms of the series:

        Inserting this into (*) we get

          -----              /     /  -----       \ 2 \  -----
          \           N-1    |     |  \         N |   |  \           L
           >    N b  x    =  | 1 + |   >    b  x  |   |   >    L a  x
          /        N         |     |  /      N    |   |  /        L
          -----              \     \  -----       /   /  -----
           N=1                         N=0                L=1

        We perform the square on the r.h.s. using Cauchy's rule
        and obtain:


           -----
           \           N-1
            >    N b  x    =
           /        N
           -----
            N=1

                              N
               /     -----  -----            \  -----
               |     \      \              N |  \           L
               | 1 +  >      >    b    b  x  |   >    L a  x  .
               |     /      /      N-M  M    |  /        L
               \     -----  -----            /  -----
                      N=0    M=0                 L=1

        Expanding this once again with Cauchy's product rule we get

           -----
           \           N-1
            >    N b  x    =
           /        N
           -----
            N=1

                                L-1     N
           -----      /        -----  -----                    \
           \      L-1 |        \      \                        |
            >    x    | L a  +  >      >    b    b  (L-N) a    | .
           /          |    L   /      /      N-M  M        L-N |
           -----      \        -----  -----                    /
            L=1                 N=0    M=0

        From this we immediately deduce the recursion relation

                      L-1                 N
                     -----              -----
                     \       L-N        \
           b  = a  +  >     ----- a      >    b    b  .  (**)
            L    L   /        L    L-N  /      N-M  M
                     -----              -----
                      N=0                M=0

        This relation is easily generalized to the case of a
        series in more than one variable, where the same comments
        apply as in the case of log and exp above.

        For the hyperbolic tangent the relation is nearly the same.
        Since its differential equation has a `-' where that of
        tangent has a `+' we simply have to do the same substitution
        in the relation (**).  For the cotangent we get an additional
        overall minus sign.

        There is one additional problem to be handled: if the constant
        term of T(x), i.e. T(x0) is a pole of tangent. This can be
        solved quite easily: for a small quantity TAYEPS we calculate

             Te(x) = T(x) + TAYEPS .

        and perform the above calculation for Te(x). Then we recover
        the desired result via the relation

            tan T(x)  = tan (Te(x) - TAYEPS)

                           tan Te(x) - tan TAYEPS
                      = ---------------------------- .
                         1 + tan Te(x) * tan TAYEPS

        For the other functions we have similar relations:

            tanh T(x) = tanh (Te(x) - TAYEPS)

                           tanh Te(x) - tanh TAYEPS
                      = ------------------------------ ,
                         1 - tanh Te(x) * tanh TAYEPS

            cot T(x)  = cot (Te(x) - TAYEPS)

                         1 + cot Te(x) * cot TAYEPS
                      = ---------------------------- ,
                           cot TAYEPS - cot Te(x)

            coth T(x) = coth (Te(x) - TAYEPS)

                         1 - coth Te(x) * coth TAYEPS
                      = ------------------------------ .
                           coth Te(x) - coth TAYEPS

        We know that this result is independent of TAYEPS, and we can
        thus evaluate it for any value of TAYEPS.

        Since we further know that T(x0) is a pole of the function in
        question, we can express tan (T(x0) + TAYEPS) as

                                1
                        - ------------ ,
                           tan TAYEPS

        and similarly all the other possible expressions involving
        TAYEPS can be written in terms of tan TAYEPS and tanh TAYEPS,
        respectively. This makes it possible to just substitute any
        occurrence of tan TAYEPS or tanh TAYEPS by zero, which greatly
        simplifies the final calculation

        ;


!*!*taylor!-epsilon!*!* := compress '(t a y e p s);

symbolic procedure taysimptan u;
  %
  % Special Taylor expansion function for circular and hyperbolic
  %  tangent and cotangent
  %
  if not (car u memq '(tan cot tanh coth)) or cddr u
    then confusion 'taysimptan
   else Taylor!:
    begin scalar a,a0,clist,coefflis,il,l0,l,poleflag,tay,tp;
    tay := taysimpsq simp!* cadr u;
    if not Taylor!-kernel!-sq!-p tay
      then return mksq({car u,mk!*sq tay},1);
    tay := mvar numr tay;
    tp := TayTemplate tay;
    l0 := make!-cst!-powerlist tp;
    %
    % First we get rid of possible zero coefficients.
    %
    l := prune!-coefflist TayCoeffList tay;
%    if null l then return !*tay2q cst!-Taylor!*(simp!* {car u,0},tp);
    %
    % The following relies on the standard ordering of the
    % TayCoeffList.
    %
    if not null l and is!-neg!-pl TayCfPl car l
      then Taylor!-error('essential!-singularity,car u);
    a0 := TayGetCoeff(l0,l);
    il := if null l then nlist(1,length tp)
           else smallest!-increment l;
    %
    %%% handle poles of function
    %
    a := quotsq(a0,simp 'pi);
    if car u memq '(tanh coth) then a := subs2!* multsq(a,simp 'i);
    if car u memq '(tan tanh) and
       denr(a := multsq(a,simp '2)) = 1 and
       fixp (a := numr a) and not evenp a
       or car u memq '(cot coth) and denr a = 1 and
          (null (a := numr a) or fixp a)
      then <<
        %
        % OK, now we are at a pole, so we add a small quantity, compute
        %  the series and use the addition formulas to get the real
        %  result.
        %
        poleflag := t;
        a0 := if car u eq 'tan
                then negsq invsq simp!* {'tan,!*!*taylor!-epsilon!*!*}
               else if car u eq 'tanh
                then invsq simp!* {'tanh,!*!*taylor!-epsilon!*!*}
               else if car u eq 'cot
                then invsq simp!* {'tan,!*!*taylor!-epsilon!*!*}
               else invsq simp!* {'tanh,!*!*taylor!-epsilon!*!*};
        clist := {TayMakeCoeff(l0,a0)};
        >>
     else clist := {TayMakeCoeff(l0,simp!* {car u,mk!*sq a0})};
    %
    coefflis := makecoeffs0(tp,TpNextList tp,il);

    if not null coefflis
      then for each cc in cdr coefflis do
             begin scalar cf,s,pos,x,y,n,n1;
               s := nil ./ 1;
               pos := find!-non!-zero cc;
               n := nth(nth(cc,car pos),cdr pos);
               for each p in makecoeffpairs(l0,cc,l0,il) do <<
                 x := reversip makecoeffpairs1(l0,car p,l0,il);
                 y := nil ./ 1;
                 for each z in x do
                   y := addsq(y,multsq(TayGetCoeff(car z,clist),
                                       TayGetCoeff(cdr z,clist)));
                 n1 := nth(nth(car p,car pos),cdr pos);
                 s := addsq(s,
                            multsq(!*TayExp2q(n - n1),
                                   multsq(y,TayGetCoeff(cdr p,l))))>>;
               cf := quotsq(s,!*TayExp2q n);
               if car u memq '(tanh coth) then cf := negsq cf;
               cf := addsq(TayGetCoeff(cc,l),cf);
               if null numr cf then return;  % short cut for efficiency
               if car u eq 'cot then cf := negsq cf;
               clist := TayMakeCoeff(cc,cf) . clist
             end;
    a := make!-Taylor!*(reversip clist,tp,nil,nil);
    %
    % Construct ``real'' series in case of pole
    %
    if poleflag then begin scalar x;
       x := if car u eq 'cot
              then cst!-Taylor!*(
                     invsq simp {'tan,!*!*taylor!-epsilon!*!*},tp)
             else if car u eq 'coth
              then cst!-Taylor!*(
                     invsq simp {'tanh,!*!*taylor!-epsilon!*!*},tp)
             else cst!-Taylor!*(
                     simp {car u,!*!*taylor!-epsilon!*!*},tp);

       if car u eq 'tan then
         a := quottaylor(addtaylor(a,negtaylor x),
                         addtaylor(cst!-taylor!*(1 ./ 1,tp),
                                   multtaylor(a,x)))
        else if car u eq 'cot then
         a := quottaylor(addtaylor(multtaylor(a,x),
                                   cst!-taylor!*(1 ./ 1,tp)),
                         addtaylor(x,negtaylor a))
        else if car u eq 'tanh then
         a := quottaylor(addtaylor(a,negtaylor x),
                         addtaylor(cst!-taylor!*(1 ./ 1,tp),
                                   negtaylor multtaylor(a,x)))
        else if car u eq 'coth then
         a := quottaylor(addtaylor(multtaylor(a,x),
                                   cst!-taylor!*(-1 ./ 1,tp)),
                         addtaylor(x,negtaylor a));
      
        if not (a freeof !*!*taylor!-epsilon!*!*)
          then set!-TayCoeffList(a,
                  for each pp in TayCoeffList a collect
                    TayMakeCoeff(TayCfPl pp,
                                 subsq(TayCfSq pp,
                                       {!*!*taylor!-epsilon!*!* . 0})));
      end;
    %
    if !*taylorkeeporiginal and TayOrig tay
      then set!-TayOrig(a,simp {car u,mk!*sq TayOrig tay});
    set!-Tayflags(a,TayFlags tay);
    return !*tay2q a
  end;

put('tan,'taylorsimpfn,'taysimptan);
put('cot,'taylorsimpfn,'taysimptan);
put('tanh,'taylorsimpfn,'taysimptan);
put('coth,'taylorsimpfn,'taysimptan);



comment For the circular sine and cosine and their reciprocals
        we calculate the exponential and use it via the formulas


                     exp(+I*x) - exp(-I*x)
            sin x = ----------------------- ,
                               2

                     exp(+I*x) + exp(-I*x)
            cos x = ----------------------- ,
                               2

        etc. To avoid clumsy expressions we separate the constant term
        of the argument,

               T(x) = a0 + T1(x),

        and use the addition theorems which give

                       1
            sin T(x) = - exp(+I*T1(x)) * (sin a0 - I*cos a0) +
                       2

                       1
                       - exp(-I*T1(x)) * (sin a0 + I*cos a0) ,
                       2
       
                       1
            cos T(x) = - exp(+I*T1(x)) * (cos a0 + I*sin a0) +
                       2

                       1
                       - exp(-I*T1(x)) * (cos a0 - I*sin a0) .
                       2

        ;


symbolic procedure taysimpsin u;
  %
  % Special Taylor expansion function for circular sine, cosine,
  %  and their reciprocals
  %
  if not (car u memq '(sin cos sec csc)) or cddr u
    then confusion 'taysimpsin
   else Taylor!:
    begin scalar l,tay,result,tp,i1,l0,a0,a1,a2;
    tay := taysimpsq simp!* cadr u;
    if not Taylor!-kernel!-sq!-p tay
      then return mksq({car u,mk!*sq tay},1);
    tay := mvar numr tay;
    tp := TayTemplate tay;
    l0 := make!-cst!-powerlist tp;
    l := prune!-coefflist TayCoeffList tay;
%    if null l then return !*tay2q cst!-Taylor!*(simp!* {car u,0},tp);
%    if is!-neg!-pl TayCfPl car l
%      then Taylor!-error('essential!-singularity,car u);
    a0 := TayGetCoeff(l0,l);
    %
    % make constant term to 0
    %
    i1 := simp 'i;
    if not null numr a0
      then tay := addtaylor(tay,cst!-Taylor!*(negsq a0,tp));
    result := taysimpexp{'exp,multtaylor(tay,cst!-Taylor!*(i1,tp))};
    a1 := simp!* {'sin,mk!*sq a0} . simp!* {'cos,mk!*sq a0};
    if car u memq '(sin csc) then <<
      a2 := addsq(car a1,multsq(i1,cdr a1));
      a1 := addsq(car a1,negsq multsq(i1,cdr a1));
      >>
     else <<
      a2 := addsq(cdr a1,negsq multsq(i1,car a1));
      a1 := addsq(cdr a1,multsq(i1,car a1));
     >>;
    result := multsq(addsq(multsq(result,a1),
                           multsq(taysimpsq!* invsq result,a2)),
                     1 ./ 2);
    if car u memq '(sec csc) then result := invsq result;
    return taysimpsq!* result
  end;


put('sin,'taylorsimpfn,'taysimpsin);
put('cos,'taylorsimpfn,'taysimpsin);
put('sec,'taylorsimpfn,'taysimpsin);
put('csc,'taylorsimpfn,'taysimpsin);



comment For the hyperbolic sine and cosine and their reciprocals
        we calculate the exponential and use it via the formulas


                     exp(+x) - exp(-x)
           sinh x = ------------------- ,
                             2

                     exp(+x) + exp(-x)
           cosh x = ------------------- ,
                             2

        etc. To avoid clumsy expressions we separate the constant term
        of the argument,

               T(x) = a0 + T1(x),

        and use the addition theorems which give

                        1
            sinh T(x) = - exp(+T1(x)) * (sinh a0 + cosh a0) +
                        2

                        1
                        - exp(-T1(x)) * (sinh a0 - cosh a0) ,
                        2
       
                        1
            cosh T(x) = - exp(+T1(x)) * (cosh a0 + sinh a0) +
                        2

                        1
                        - exp(-T1(x)) * (cosh a0 - sinh a0) .
                        2
        ;


symbolic procedure taysimpsinh u;
  %
  % Special Taylor expansion function for circular sine, cosine,
  %  and their reciprocals
  %
  if not (car u memq '(sinh cosh sech csch)) or cddr u
    then confusion 'taysimpsin
   else Taylor!:
    begin scalar l,tay,result,tp,l0,a0,a1,a2;
    tay := taysimpsq simp!* cadr u;
    if not Taylor!-kernel!-sq!-p tay
      then return mksq({car u,mk!*sq tay},1);
    tay := mvar numr tay;
    tp := TayTemplate tay;
    l0 := make!-cst!-powerlist tp;
    l := prune!-coefflist TayCoeffList tay;
%    if null l then return !*tay2q cst!-Taylor!*(simp!* {car u,0},tp);
%    if is!-neg!-pl TayCfPl car l
%      then Taylor!-error('essential!-singularity,car u);
    a0 := TayGetCoeff(l0,l);
    %
    % make constant term to 0
    %
    if not null numr a0
      then tay := addtaylor(tay,cst!-Taylor!*(negsq a0,tp));
    result := taysimpexp{'exp,tay};
    a1 := simp!* {'sinh,mk!*sq a0} . simp!* {'cosh,mk!*sq a0};
    if car u memq '(sinh csch) then <<
      a2 := addsq(car a1,cdr a1);
      a1 := subtrsq(car a1,cdr a1);
      >>
     else <<
      a2 := addsq(cdr a1,car a1);
      a1 := subtrsq(cdr a1,car a1);
     >>;
    result := multsq(addsq(multsq(result,a2),
                           multsq(taysimpsq!* invsq result,a1)),
                     1 ./ 2);
    if car u memq '(sech csch) then result := invsq result;
    return taysimpsq!* result
  end;


put('sinh, 'taylorsimpfn, 'taysimpsinh);
put('cosh, 'taylorsimpfn, 'taysimpsinh);
put('sech, 'taylorsimpfn, 'taysimpsinh);
put('csch, 'taylorsimpfn, 'taysimpsinh);


comment Support for the integration of Taylor kernels.
        Unfortunately, with the current interface, only Taylor kernels
        on toplevel can be treated successfully.

        The way it is down means stretching certain interfaces beyond
        what they were designed for...but it works!

        First we add a rule that replaces a call to INT with a Taylor
        kernel as argument by a call to TAYLORINT, then we define
        REVALTAYLORINT as simplification function for that;


algebraic let {

  int(symbolic algebraic(taylor!*(~x,~y,~z,~w)),~u)
      => taylorint(~x,~y,~z,~w,~u),

  int(log(~u)^~~n * symbolic algebraic(taylor!*(~x,~y,~z,~w)),~u)
      => log(u)^n * int(symbolic('(taylor!* x y z w)),u)
          - int(log(u)^(n-1)
             * taylorcombine(int(symbolic('(taylor!* x y z w)),u)/u),u),

  int(~x,~y) => taylorint1(~x,~y)
                  when not symbolic algebraic(~x freeof 'Taylor!*)
};


put('taylorint1, 'psopfn, 'revaltaylorint1);

symbolic procedure revaltaylorint1 x;
  begin scalar u,v;
    u := car x;
    v := cadr x;
    if Taylor!*p u then return revaltaylorint append(cdr u,{v});
    u := reval taylorcombine u;
    if Taylor!*p u then return revaltaylorint append(cdr u,{v});
    if not atom u 
      then if car u memq '(plus minus difference)
             then return
                    reval (car u . for each term in cdr u collect
                                     {'int,term,v});
    lprim "Converting Taylor kernels to standard representation";
    return aeval {'int,taylortostandard car x,v}
  end;

put('taylorint, 'psopfn, 'revaltaylorint);

symbolic procedure revaltaylorint u;
  begin scalar taycfl, taytp, tayorg, tayflgs, var;
    taycfl := car u;
    taytp := cadr u;
    tayorg := caddr u;
    tayflgs := cadddr u;
    var := car cddddr u;
    if null (var member TayTpVars taytp)
      then return mk!*sq !*tay2q
            make!-Taylor!*(
              for each pp in taycfl collect 
                TayCfPl pp . simp!* {'int,mk!*sq TayCfSq pp,var},
              taytp,
              if not null tayorg
                then simp!* {'int,mk!*sq tayorg,var}
               else nil,
            nil)
     else return mk!*sq
            ((if null car result then !*tay2q cdr result
               else addsq(car result,!*tay2q cdr result))
             where result := inttaylorwrttayvar1(taycfl,taytp,var))
   end;

endmodule;

end;