reduce-historical/r38/packages/misc/limits.red

module limits;

%% A fast limit package for REDUCE for functions which are continuous
%% except for computable poles and singularities.

%% Author: Stanley L. Kameny.

%% Revised 23 Mar 1993.  Version 1.4.

%% Added capability for using either the Taylor series package or the
%% Truncated Power Series Package.

%% Added provisions for transformation of certain irrational functions
%% into rational functions before limit calculation in order to be able
%% to compute series.

%% Changed the algebraic interface so that if limit package fails, an
%% equivalent of the original expression is returned.

%% Allowed for limited recursion through limsimp.

%% Corrected several bugs.

%% Date: 10 Oct 1990. Original version.

%% The Truncated Power Series package is used for non-critical points.
%% L'Hopital's rule is used in critical cases, with preprocessing of
%% <infinity - infinity> forms and reformatting of product forms in
%% order to be able to apply l'Hopital's rule.  A limited amount of
%% bounded arithmetic is also employed where applicable.

%% This limits package makes use of the ideas embodied in the
%% limit.red package, by Ian Cohen and John Fitch, 11 July 1990
%% that is in reduce-netlib;  in fact, some code is lifted bodily.
%% The idea of using the Truncated Power Series package to compute
%% limits at non-critical points, and the substitutions used in limit!+
%% and limit!- come from there.

load!-package 'tps; %load!-package 'taylor;

lisp(ps!:order!-limit := 100);

switch usetaylor; off usetaylor;

fluid '(!*precise lhop!# lplus!# !*protfg !*msg !*rounded !*complex
	!#nnn lim00!# !*crlimtest !*lim00rec);

!*lim00rec := t;  % Default value.

global '(erfg!* exptconv!#);

global '(abslims!#);
symbolic(abslims!# := {0,1,-1,'infinity,'(minus infinity)});
 % others may be added.

fluid '(lsimpdpth); global '(ld0!#); symbolic(ld0!# := 3);

flag('(limit limit!+ limit!- limit2),'full);

symbolic
for each c in '(limit limit!+ limit!- limit2) do
  <<remflag({c},'opfn); put(c,'simpfn,'simplimit)>>;

symbolic procedure limit2(top,bot,xxx,a);
   lhopital(top,bot,xxx,a) where lhop!#=0;

symbolic procedure limit!+(ex,x,a);
<<ex := simp!* limlogsort ex;
  if a = 'infinity then rederr "Cannot approach infinity from above"
  else if a = '(minus infinity) then
     limit(prepsq subsq(ex,list(x .
        list('quotient,-1,list('expt,'!*eps!*,2)))),'!*eps!*,0)
  else limit(prepsq subsq(ex,list(x .
        list('plus,a,list('expt,'!*eps!*,2)))),'!*eps!*,0)>>;

symbolic procedure limit!-(ex,x,a);
<<ex := simp!* limlogsort ex;
  if a = 'infinity then
     limit(prepsq subsq(ex,list(x .
        list('quotient,1,list('expt,'!*eps!*,2)))),'!*eps!*,0)
  else if a = '(minus infinity) then
     rederr "Cannot approach -infinity from below"
  else limit(prepsq subsq(ex,list(x .
        list('difference,a,list('expt,'!*eps!*,2)))),'!*eps!*,0)>>;

symbolic procedure limit(ex,xxx,a); limit0(limlogsort ex,xxx,a)
  where !*combinelogs=nil,lhop!#=0,lplus!#=0,lim00!#=nil,lsimpdpth=0;

symbolic procedure limlogsort x;
   begin scalar !*precise;
      x := prepsq simp!* x;
      return if countof('log,x)>1 then logsort x else x
   end;

symbolic procedure countof(u,v);
   if u = v then 1 else if atom v then 0
      else countof(u,car v)+countof(u,cdr v);

symbolic procedure simplimit u;
   % The kludgey handling of cot needs to be fixed some day.
   begin scalar fn,exprn,var,val,old,v,!*precise,!*protfg;
     if length u neq 4
       then rerror(limit,1,
		   "Improper number of arguments to limit operator");
     fn:= car u; exprn := cadr u; var := !*a2k caddr u; val := cadddr u;
     !*protfg := t;   % ACH: I'm not sure why this is needed.
     old := get('cot,'opmtch);
     put('cot,'opmtch,
	 '(((!~x) (nil . t) (quotient (cos !~x) (sin !~x)) nil)));
     v := errorset!*({'apply,mkquote fn,mkquote {exprn,var,val}},nil);
     put('cot,'opmtch,old);
     !*protfg := nil;
     return if errorp v or (v := car v) = aeval 'failed then mksq(u,1)
	     else simp!* v
   end;

symbolic procedure limit0(exp,x,a);
   begin scalar exp1;
     exp1 := simp!* exp;
     if a = 'infinity then
        return limit00(subsq(exp1,{x . {'quotient,1,{'expt,x,2}}}),x);
     if a = '(minus infinity) then
        return limit00(subsq(exp1,{x . {'quotient,-1,{'expt,x,2}}}),x);
     return
        (<<!*protfg := t;
           y := errorset!*
             ({'subsq,mkquote(exp := simp!* exp),mkquote{(x . a)}},nil)
              where !*expandlogs=t;
           !*protfg := nil;
           if not (errorp y) and not ((y := car y) = aeval 'failed)
              then mk!*sq y
           else if neq(a,0) then limit00(subsq(exp1,{x .
              {'plus,a,x}}),x)
           else limit00(exp1,x)>> where y=nil) end;

symbolic procedure limit00(ex,x);
   begin scalar p,p1,z,xpwrlcm,lim,ls;
     if (lim := crlimitset(p := prepsq ex,x)) then go to ret;
     if not lim00!# then
       <<lim00!# := not !*lim00rec;
         p1 := factrprep prepsq ex;
         if (xpwrlcm := xpwrlcmp(p1,x)) neq 1 then
           <<ex := subsq(ex,{x . {'expt,x,xpwrlcm}});
             p1 := factrprep prepsq ex>>;
         if (z := pwrdenp(p1,x)) neq 1 then
            ex := simp!*{'expt,p1,z};
         if (lim := crlimitset(p := prepsq ex,x)) then go to ret>>;
     % tps has failed because ex has a branch point at a or is undefined
     % at a or tps itself has failed or Reduce has not recognized the
     % numeric value of an expression.
     if %xpwrlcm and xpwrlcm>1 or
     lsimpdpth>ld0!#
       then lim := aeval 'failed else
       <<lsimpdpth := lsimpdpth + 1; ls := t;
         lim := limsimp(p,x);
         if prepsq simp!* lim = 'failed and lsimpdpth=1 then
           <<exptconv!# := nil; p := expt2exp(p,x);
             if exptconv!# then lim := limsimp(p,x)>> >>;
ret: return
       <<if ls then lsimpdpth := lsimpdpth - 1;
         if not z or z = 1 or lim=0 then lim
         else if (ls := prepsq simp!* lim) = '(minus infinity)
           then if (-1)^z = 1 then aeval 'infinity else lim
         else if ls member '(infinity failed) then lim
         else mk!*sq simp!* {'expt,prepsq simp!* lim,{'quotient,1,z}}>>
    end;

symbolic procedure factrprep p;
   begin scalar !*factor;
     !*factor := t;
     return prepsq simp!* p end;

symbolic procedure expt2exp(p,x);
   if atom p then p
   else if eqcar(p,'expt)
      and not freeof(cadr p,x) and not freeof(caddr p,x) then
	<<exptconv!# := t; {'expt,'e,{'times,{'log,cadr p},caddr p}}>>
   else expt2exp(car p,x) . expt2exp(cdr p,x);

symbolic procedure xpwrlcmp(p,x);
   if atom p then 1
   else if eqcar(p,'expt) and cadr p = x then getdenom caddr p
   else if eqcar(p,'sqrt) then getdenomx(cadr p,x)
   else lcm(xpwrlcmp(car p,x),xpwrlcmp(cdr p,x));

symbolic procedure getdenomx(p,x);
   if freeof(p,x) then 1
   else if eqcar(p,'minus) then getdenomx(cadr p,x)
   else if p = x or eqcar(p,'times) and x member cdr p then 2
   else xpwrlcmp(p,x);

symbolic procedure getdenom p;
   if eqcar(p,'minus) then getdenom cadr p
   else if eqcar(p,'quotient) and numberp caddr p then caddr p
   else 1;

symbolic procedure pwrdenp(p,x);
   if atom p then 1
   else if eqcar(p,'expt) and not freeof(cadr p,x)
      then getdenom caddr p
   else if eqcar(p,'sqrt) and not freeof(cadr p,x) then 2
   else if eqcar(p,'minus) then pwrdenp(cadr p,x)
   else if car p member '(times quotient) then
      (<<for each c in cdr p do m := lcm(m,pwrdenp(c,x)); m>>
       where m=1)
   else if atom car p then 1
   else lcm(pwrdenp(car p,x),pwrdenp(cdr p,x));

symbolic procedure limitset(ex,x,a);
 if !*usetaylor then
  <<!*protfg := t;
    ex := errorset!*({'limit1t,mkquote ex,mkquote x,mkquote a},nil);
    !*protfg := nil;
    if errorp ex then nil else car ex>>
 else % use tps.
  begin scalar oldpslim;
      !*protfg := t; oldpslim := simppsexplim '(1);
      ex := errorset!*({'limit1p,mkquote ex,mkquote x,mkquote a},nil);
      !*protfg := nil; simppsexplim list car oldpslim;
      return if errorp ex then nil else car ex
  end;

symbolic procedure limit1t(ex,x,a);
   begin scalar nnn, vvv,oldklist;
     oldklist := get('taylor!*,'klist);
     ex := {ex,x,a,0};
     vvv := errorset!*({'simptaylor,mkquote ex},!*backtrace);
     put('taylor!*,'klist,oldklist);
     if errorp vvv then <<if !*backtrace then break();return nil>>
      else ex := car vvv;
     if kernp ex then ex := mvar numr ex
      else return nil;
     if not eqcar(ex,'taylor!*) then return nil
       else ex := cadr ex;
    % ex is now the list of coefs and values, but we need the lowest
    % order non-zero value, which may not be the first of these.
    % if this list is empty the result is zero
    while ex and null numr cdr car ex do ex := cdr ex;
    if null ex then return (!#nnn := 0) else
       !#nnn := nnn := caaaar ex;
     vvv := cdar ex;
     return
       if tayexp!-greaterp(nnn,0) then 0
       else if nnn=0 then mk!*sq vvv
       else if !*complex then 'infinity
       else if domainp(nnn := numr vvv) then
         (if !:minusp nnn
         then aeval '(minus infinity) else 'infinity)
       else aeval{'times,{'sign,prepsq vvv},'infinity}
   end;

symbolic procedure limit1p(ex,x,a);
   begin scalar aaa, nnn, vvv;
     aaa := mk!*sq simpps1(ex,x,a);
     !#nnn := nnn := mk!*sq simppsorder list aaa;
     vvv := simppsterm1(aaa,min(nnn,0));
     return
       if nnn>0 then 0
       else if nnn=0 then mk!*sq vvv
       else if !*complex then 'infinity
       else if domainp(nnn := car vvv) then
         (if !:minusp nnn then aeval '(minus infinity)
                          else 'infinity)
       else aeval{'times,{'sign,prepsq vvv},'infinity}
   end;

symbolic procedure crlimitset(ex,x);
 (begin scalar lim1,lim2,n1,fg,limcr,!#nnn;
    lim1 := limitset(ex,x,0);
    if null lim1 then if r and c then return nil else go to a;
    if (n1 := !#nnn) < 0 or lim1 member abslims!#
       or r and c then return lim1;
 a: if not !*crlimtest then return lim1;
    if not r then on rounded; if not c then on complex;
    if not (lim2 := limitset(ex,x,0))
      or !#nnn > n1 then <<fg := t; go to ret>>;
    if !#nnn < n1 or lim2 member abslims!# then go to ret;
   % at this point, both lim1 and lim2 have values.  If they are
   % equivalent, we want lim1; otherwise lim2.
    if (limcr := topevalsetsq lim1) and
       evalequal(prepsq simp!* lim2,prepsq limcr)
         then fg := t;
ret:if not r then off rounded; if not c then off complex;
    return if fg then lim1 else lim2 end)
  where r=!*rounded,c=!*complex,!*msg=nil;

symbolic procedure topevalsetsq u;
  <<!*protfg := t;
    if not r then on rounded; if not c then on complex;
    u := errorset!*({'simp!*,{'aeval,{'prepsq,{'simp!*,mkquote u}}}},
      nil);
    !*protfg := nil;
    if not r then off rounded;if not c then off complex;
    if errorp u then nil else car u>>
  where r=!*rounded,c=!*complex,!*msg=nil;

put('times,'limsfn,'ltimesfn);
put('quotient,'limsfn,'lquotfn);
put('plus,'limsfn,'lplusfn);
put('expt,'limsfn,'lexptfn);

symbolic procedure limsimp(ex,x);
  % called when limit1 has failed, to apply more sophisticated methods.
  % output must be aeval form.
   begin scalar y,c,z,m,ex0;
      if eqcar(ex,'minus) then <<m := t; ex := cadr ex>>;
      ex0 := ex;
      if not atom ex then  % check for plus, times, or quotient.
         <<if(z := get(y := car ex,'limsfn))
              then ex := apply(z,list(ex,x))>>
         else <<if ex eq x then ex := 0; go to ret>>;
      if y eq 'plus then go to ret;
      if y eq 'expt then if ex then return ex else ex := ex0 . 1;
      if z then<<z := car ex; c := cdr ex>>
         else <<z := prepsq !*f2q numr(ex := simp!* ex);
                c := prepsq !*f2q denr ex>>;
      ex := lhopital(z,c,x,0);
 ret: if m and prepsq simp!* ex neq 'failed then
         ex := aeval lminus2 ex;
      return ex end;

symbolic procedure lminus2 ex;
   if numberp ex then -ex
   else if eqcar(ex,'minus) then cadr ex
   else list('minus,ex);

symbolic procedure ltimesfn(ex,x); specchk(ex,1,x);

symbolic procedure lquotfn(ex,x);
 %  (if eqcar(n,'expt) and (nlim :=lexptfn(n,x))
specchk(cadr ex,caddr ex,x);

symbolic procedure lexptfn(ex,x);
   if not evalequal(cadr ex,0) and limit00(simp!* caddr ex,x)=0
     then 1;

symbolic procedure specchk(top,bot,x);
   begin scalar tlist,blist,tinfs,binfs,tlogs,blogs,tzros,bzros,
         tnrms,bnrms,m;
      if eqcar(top,'minus) then <<m := t; top := cadr top>>;
      if eqcar(bot,'minus) then <<m := not m; bot := cadr bot>>;
      tlist := limsort(timsift(top,x),x);
      blist := limsort(timsift(bot,x),x);
      tinfs := cdr(tlogs := logcomb(cadr tlist,x)); tlogs := car tlogs;
      binfs := cdr(blogs := logcomb(cadr blist,x)); blogs := car blogs;
      tzros := car tlist; tnrms := caddr tlist;
      bzros := car blist; bnrms := caddr blist;
      if tlogs and not blogs then
         <<top := triml append(tlogs,tnrms);
           bot := triml append(bzros,append(binfs,
              append(bnrms,trimq append(tinfs,tzros))))>>
      else if blogs and not tlogs then
         <<bot := triml append(blogs,bnrms);
           top := triml append(tzros,append(tinfs,
              append(tnrms,trimq append(binfs,bzros))))>>
      else
         <<top := triml append(cadr tlist,trimq bzros);
           bot := triml append(cadr blist,
              append(bnrms,trimq append(tzros,tnrms)))>>;
      if m then top := list('minus,top);
      return top . bot end;

symbolic procedure trimq l;
   if l then list list('quotient,1,
      if length l>1 then 'times . l else car l);

symbolic procedure triml l;
   if null l then 1 else if length l>1 then 'times . l else car l;

symbolic procedure limsort(ex,x);
   begin scalar zros,infs,nrms,q,s;
      for each c in ex do
         if (q := numr(s := simp!* limit00(simp!* c,x)))
              and numberp q and not zerop q then nrms := q . nrms
         else if null q or zerop q then zros := c . zros
         else if caaar q memq '(failed infinity) then infs := c.infs
         else nrms := (prepsq s) . nrms;
      return list(zros,infs,nrms) end;

symbolic procedure logcomb(tinf,x);
  % separate product list into log terms and others.
   begin scalar tlog,c,z;
      while tinf do
         <<c := car tinf; tinf := cdr tinf;
           if eqcar(c,'log)
              or eqcar(c,'expt) and eqcar(cadr c,'log)
              or eqcar(c,'plus) and
                (eqcar(cadr(c := logjoin(c,x)),'log)
                  or eqcar(cadr c,'minus) and eqcar(cadadr c,'log))
                    and freeof(cddr c,x)
           then tlog := c . tlog else z := c . z>>;
      return tlog . reversip z end;

symbolic procedure logjoin(p,x);
  % combine log terms in sum list into a single log.
   begin scalar ll,z;
      for each c in cdr p do
         if freeof(c,x) then z := c . z
         else if eqcar(c,'log) then ll := (cadr c) . ll
         else if eqcar(c,'minus) and eqcar(cadr c,'log) then
            ll := list('quotient,1,cadadr c) . ll
         else z := c . z;
      if ll then ll := list list('log,'times . ll);
      return (car p) . append(ll,reversip z) end;

symbolic procedure timsift(ex,x);
   if eqcar(ex,'times) then cdr ex
   else if eqcar(ex,'plus) then list logjoin(ex,x)
    % for plus, combine log terms, change infinity - infinity to
    % inner quotient.
   else list ex;

symbolic procedure lplusfn(ex,x);
  % combine logs and evaluate each limit term.  if infinity - infinity
  % is found, attempt conversion to quotient form for lhopital.
   begin scalar z,infs,nrms,vals,vp,vm,cz,vnix;
      lplus!# := lplus!# + 1;
     % write "lplus#=",lplus!#; terpri();
      if lplus!#>4 then return aeval 'failed;
      z := limsort(cdr ex,x); % ignore car z, a list of 0's.
      nrms := caddr z; infs := cadr z;
      if length infs>1 then
         <<infs := logjoin('plus . infs,x);
           infs := if eqcar(infs,'plus) then cdr infs else list infs>>;
     % at this point, only infs needs to be evaluated.
      vals := for each c in infs collect
         minfix prepsq simp!* limit00(simp!* c,x);
      z := infs;
      for each c in vals do
         <<cz := car z; z := cdr z;
           if c eq 'infinity then vp := cz . vp
	   else if c = '(minus infinity) then vm := cz . vm
           else if c eq 'failed then vnix := cz . vnix
           else nrms := cz . nrms>>;
      if vm and not vp or vp and not vm or length vnix = 1
         or length vm > 1 or length vp > 1 then return aeval 'failed;
      if vm then vm := qform(car vp,vm);
      if vnix then vnix := qform(car vnix,cdr vnix);
      vm := append(nrms,append(vm,vnix));
      return if null vm then 0 else
         limit00(simp!* if length vm>1 then 'plus . vm else car vm,x) 
         end;

symbolic procedure minfix v;
   if eqcar(v,'minus) and numberp cadr v then -cadr v else v;

symbolic procedure qform(a,b);
   list list('quotient,list('plus,1,
      list('quotient,if length b = 1 then car b else 'plus . b,a)),
         list ('quotient,1,a));

symbolic procedure lhopital(top,bot,xxx,a);
  begin scalar limt, limb, nvt, nvb;
     nvt := notval(limt := limfix(top,xxx,a));
     nvb := notval(limb := limfix(bot,xxx,a));
 % possibilities for lims are {failed, infinity, -infinity, bounded,
 % nonzero, zero} and each combination of cases has to be handled.
     if limt=0 and limb=0 or nvt and nvb then go to lhop;
     if specval limt or specval limb then return speccomb(limt,limb);
     if limb=0 then return aeval 'infinity;  % maybe impossible.
     return aeval list('quotient,limt,limb);
lhop: lhop!# := lhop!#+1;
    %  write "lhop#=",lhop!#; terpri();
      if lhop!#>6 then return aeval 'failed;
      return limit0(prepsq quotsq(diffsq(simp!* top,xxx),
          diffsq(simp!* bot,xxx)),xxx,a) end;

symbolic procedure notval lim;
   not lim or infinp prepsq simp!* lim;

symbolic procedure infinp x; member(x,'(infinity (minus infinity)));

symbolic procedure specval lim;
   notval lim or lim eq 'bounded;


symbolic procedure speccomb(a,b);
 aeval
  (if not a or not b or b eq 'bounded then 'failed
   else if notval b then 0
   else if notval a then
      if numberp b then
	 if b>=0 then a
         else if a eq 'infinity then '(minus infinity) else 'infinity
      else ((if c then
        <<c := prepsq c;
          if evalgreaterp(c,0) then cc := 1 else if evallessp(c,0)
            then cc := -1;
          if cc then c := if a eq 'infinity then 1 else -1;
	  if cc then
             if c*cc = 1 then 'infinity else '(minus infinity)
	    else {'times,{'sgn,b},a}>> else {'quotient,a,b})
          where c=topevalsetsq prepsq simp!* b,cc=nil)
   else 'failed);

symbolic procedure limfix(ex,x,a);
   (if val then val
       else limitest(ex,x,a))
    where val=limitset(ex,x,a);

symbolic procedure limitest(ex,x,a);
  if ex then if atom ex then if ex eq x then a else ex else
  begin scalar y,arg,val;
    if eqcar(ex,'expt) then
       if cadr ex eq 'e then ex := list('exp,caddr ex)
       else return exptest(cadr ex,caddr ex,x,a);
    if (y := get(car ex,'fixfn)) then
       <<arg := cadr ex; val := limitset(arg,x,a);
         return apply1(y,
          if val then val else limitest(arg,x,a))>>
    else if (y := get(car ex,'limcomb)) then
       return apply3(y,cdr ex,x,a) end;

symbolic procedure exptest(b,n,x,a);
   if numberp n then
      if n<0 then limquot1(1,exptest(b,-n,x,a))
      else if n=0 then 1 else
        ((if 2*y=n then limlabs limitest(b,x,a) else limitest(b,x,a))
         where y=n/2)
   else if numberp b and b>1 then limitest(list('exp,n),x,a);

symbolic procedure limlabs a;
   if null a then nil
   else if infinp a then 'infinity
   else if a eq 'bounded then 'bounded else
     begin scalar n,d; d := denr(n := simp!* a); n := numr n;
     return if null n then a else if not numberp n then nil
         else mk!*sq abs a ./ d end;

symbolic procedure limplus(exl,x,a);
   if null exl then 0
   else limplus1(mkalg limfix(car exl,x,a),limplus(cdr exl,x,a));

symbolic procedure limplus1(a,b);
   if null a or null b then nil
   else if infinp a
      then if infinp b
         then if a eq b then a else nil else a
   else if infinp b then b
   else if a eq 'bounded or b eq 'bounded then 'bounded
   else mk!*sq addsq(simp!* a,simp!* b);

symbolic procedure limtimes(exl,x,a);
   if null exl then 1
  else ltimes1(mkalg limfix(car exl,x,a),limtimes(cdr exl,x,a));

symbolic procedure mkalg x;
   minfix if eqcar(x,'!*sq) then prepsq simp!* x else x;

symbolic procedure ltimes1(a,b);
 begin scalar c;
  return if null a or null b then nil
   else if infinp a then
         if infinp b  then
           if a = b then 'infinity else '(minus infinity)
           else if b eq 'bounded or b=0 then nil
           else if (c := limposp b) eq 'failed then nil
           else if c then a else lminus1 a
   else if infinp b then
      if a eq 'bounded or a=0 then nil
      else if (c := limposp a) eq 'failed then nil
      else if c then b else lminus1 b
   else if a eq 'bounded or b eq 'bounded then 'bounded
   else mk!*sq multsq(simp!* a,simp!* b) end;

symbolic procedure limposp a;
  (if n and not numberp n then 'failed else n and n>0)
  where n=numr simp!* a;

symbolic procedure lminus(exl,x,a);
   lminus1 mkalg limfix(car exl,x,a);

symbolic procedure lminus1 a; if a then
   if a eq 'infinity then '(minus infinity)
   else if a = '(minus infinity) then 'infinity
   else if a eq 'bounded then a
   else mk!*sq negsq simp!* a;

symbolic procedure limquot(exl,x,a);
   limquot1(mkalg limfix(car exl,x,a),mkalg limfix(cadr exl,x,a));

symbolic procedure limquot1(a,b);
 begin scalar c;
  return if null a or null b then nil
   else if infinp a then
         if infinp b then nil
           else if b eq 'bounded then nil
           else if b=0 then a
           else if (c := limposp b) eq 'failed then nil
           else if c then a else lminus1 a
   else if infinp b then 0
   else if a eq 'bounded then if b=0 then nil else 'bounded
   else if b=0 or b eq 'bounded then nil
   else mk!*sq quotsq(simp!* a,simp!* b) end;

put('log,'fixfn,'fixlog);
put('sin,'fixfn,'fixsin);
put('cos,'fixfn,'fixsin);
put('sqrt,'fixfn,'fixsqrt);
put('cosh,'fixfn,'fixcosh);
put('sinh,'fixfn,'fixsinh);
put('exp,'fixfn,'fixexp);
put('plus,'limcomb,'limplus);
put('minus,'limcomb,'lminus);
put('times,'limcomb,'limtimes);
put('quotient,'limcomb,'limquot);

symbolic procedure fixlog x;
   if zerop x then '(minus infinity) else if infinp x then 'infinity;

symbolic procedure fixsqrt x;
   if zerop x then 0 else if infinp x then 'infinity;

symbolic procedure fixsin x;
   if infinp x then 'bounded;

symbolic procedure fixcosh x;
   if infinp x then 'infinity;

symbolic procedure fixsinh x;
   if infinp x then x;

symbolic procedure fixexp x;
   if x eq 'infinity then x else if x = '(minus infinity) then 0;

endmodule;

end;