reduce-historical/r38/packages/mrvlimit/mrvlimit.red

%----------------------------------------------------------------------------
%						                            |
% Source code for a new Exp-Log limits package in REDUCE	            |
%									    |
% Author: Neil Langmead						            |
% Place:  ZIB, Berlin							    |
% Date:   April 1997					                    |
% e/mail: langmead@zib.de						    |
%									    |
% all bugs and comments or suggestions please report to Winfried Neun,      |
% ZIB, Takustrasse 7, D-14195, Berlin Dahlem, Germany: e/mail neun@zib.de  |
%----------------------------------------------------------------------------

module mrvlimit;
create!-package ('(mrvlimit expon),nil);

global '(tracelimit!*); % for the tracing facility)

switch tracelimit;

off tracelimit; % off by default

flag('(sqchk),'opfn);

%symbolic procedure max_power(f,x);
%  if domainp f then 0 else 
%    if mvar f eq x then ldeg f else 
%          max(max_power(lc f,x),max_power(red f,x))

%put('max2,'psopfn,'max_power)
load_package limits;
load_package sets;

algebraic;
off mcd;
%in "/silo/neil/test.red";
symbolic procedure maxi(f,g);
begin scalar c;
if(freeof(f,'x)) then return g;
if(freeof(g,'x)) then return f;
if(f=g) then return  f else
<<
if(null f) then return  g else
<<     
  if(null g) then return   f
  else <<
    if (intersection(f,g) neq '(list)) then return union(f,g)
      else <<
        if(evalb('x member f)='true) then return  g 
         else <<
            if(evalb('x member g)='true) then return  f
             else  <<
              if(car f='list and cadr f='list) then % double list
                            << % only want caddr f to be given to compare 
              			c:=compare(caddr f,cadr g); 
 				%write "c is ", c; write length(c)
                                return c;
                            >>
                  else <<
                          if(car g='list and cadr g='list) then
   			  <<
                             c:=compare(cadr f,caddr g);
                             %write "c is ", c;
    			     return c;
                          >>
                          else <<
                                 c:=compare(cadr f, cadr g);
 				 %write "c is ", c;
 				 return c;
 				>>;
		        >>;

              %if(c=cadr f) then return cadr f else <<
              %         if(c=cadr g) then return (cadr g)
	      % else return union(cdr f,cdr g);
					   %   >>;
                    >>;
                >>;
            >>;
          >>;
 >>;
>>; 
end; % of max

 %max
%-------------------------------------------------------------------------

algebraic;
procedure maxi1(f,g);  lisp cadr (lisp (list('list,maxi(f,g))));

algebraic;
expr procedure compare(f,g);
begin scalar logg, logf;
   logf:=log(f);
   logg:=log(g);

if(mrv_limit(logf/logg,x,infinity)=0) then return {g} else <<
 if(mrv_limit(logg/logf,x,infinity)=0) then return {f} else
return {f,g}; >>;
end;

procedure comp(f,g); lisp('list.compare(f,g));

%----------------------------------------------------------------------------
load_package assist;

symbolic procedure mrv(li);
begin
off mcd;  on factor;
% The next line doesn't do anything in symbolic mode. Presumably li
% should be simplified in some way.  However, li is not always an
% algebraic expression. Sometimes it is a list of one, or a list of a
% number.
li:=li;
if(numberp li) then return nil else <<
if(li='(list)) then return nil else <<
if(atom li) then return lisp ('list.{li}) else <<
  if(car li='times) then <<
               if(atom cadr li and atom caddr li) then
                     << if(length(cddr li)=1) then
                        return  lisp ('list.maxi1({cadr li}, {caddr li}))
                        else return maxi1({cadr li},mrv(cddr li))
                      >>
               else return  maxi1(mrv(cadr li), mrv(cddr li))
                           >>
    else <<
     if(car li='minus) then <<
          if(atom cadr li) then return 'list.{cadr li} else return
                mrv(cadr li)  
                             >>       
                            %return mrv(append({'plus},cdr li))
      else <<
     if(car li='plus) then <<
   %if(null caddr li) then return mrv(cadr li) 
   if(length cdr li=1) then %only one argument to plus
      return mrv(cadr li) else <<
             if(atom cadr li and atom caddr li) then
               << if(length(cddr li)=1) then return 
                       lisp ('list.maxi1({cadr li},{caddr li}))
                  else return 
                  lisp ('list.maxi1({cadr li},mrv(append({'plus},cddr li))))
                >> 
                else <<
                   if(atom cadr li and pairp caddr li) then
                  return 
           maxi1('list.{cadr li}, mrv(cddr li)) % here as well
                else  <<
         if(pairp cadr li and null caddr li) then return mrv(cadr li) else
       <<
         if(pairp cadr li and atom caddr li) then 
             <<
               if(length(cdr li)>2) then % we have plus with > two args
return lisp cdr ('list.maxi1(mrv(cadr li),mrv(append({'plus},cddr li)))) %her
 else       
    return lisp cdr ('list.maxi1(mrv(cadr li), mrv(cddr li))) >>
    else
<< if(null caddr li) then return mrv(cadr li)
      else return maxi1(mrv(cadr li), mrv(append({'plus},cddr li))) >>
             >> 
                      >> 
                            >> 
                                 >>
                                      >>
           else <<
           if(car li='expt) then <<
               if(cadr li neq 'e) then 
               return  mrv(cadr li)
               else <<  %we have e to the power of something
		       if sqchk mrv_limit(caddr li,'x,'infinity)
			     eq 'infinity
			   then
                         return 
                        maxi1('list.{li},mrv(caddr li)) else
                         <<
			   if sqchk mrv_limit(caddr li,'x,'infinity)
			      = '(minus infinity)
			       then
                     return  maxi1('list.{li},'list.mrv(cddr li)) else return
                            mrv(caddr li)
                          >>
                     >>
				  >>
                else << if(car li='log) then
                        << if(atom cadr li) then return mrv(cadr li) else
                           return mrv(cdr li)
                         >> else <<
                       if(car li='sqrt) then return mrv(cdr li) else
                        return mrv(car li)
                                  >>
                      >>
                  >>
            >> 
     >> % for minus
  >> %for null
 >> % for numberp
              >>;
off mcd;
end; % of mrv

algebraic;

procedure mrv1(li);  lisp (mrv(li));
%----------------------------------------------------------------------------
% procedure to return a list of subexpressions of exp
% this will then be used for the mrv function

symbolic procedure flatten(li);
  % This procedure turns a list with possibly nested sub_lists into a single
  % List with no nested sub-lists. Easier to search this list.
  makeflat(li,nil);

symbolic procedure makeflat(li,answer);
  if li=nil then nil 
  else
    if atom li then li.answer
    else if (cdr li)=nil then makeflat(car li,answer)
       else append(makeflat(car li,answer),makeflat(cdr li,answer));

algebraic;
procedure flat(li); lisp(flatten li);
procedure mkflat(li); lisp(makeflat(li,nil));
%in "max";
%trst maxi;

symbolic procedure lim(exp,var,val);
begin scalar mrv_list, rule;
mrv_list:=mrv1(exp);
if(mrv_list='(list)) then rederr "unable to compute mrv set"
  else
    <<
       rule:=list(list ('replaceby, cdr mrv_list,'w));
       let rule;
    >>;

end;
% nedd to consider if x belongs to mrv(exp), then follow rest of alg.
algebraic;                        
expr procedure move_up(exp,x);
sub({log(x)=x,x=e^x},exp);

expr procedure move_down(exp,x);
sub({e^x=x,x=log(x)},exp);
  
%off mcd;
algebraic;
expr procedure rewrite(m);
begin scalar ans_list,summ,k,g,c,A;

ans_list:={};
g:=part(m,1); write length g;
for k:=1:arglength(m) do
<<
  summ:=length(den(part(m,k)))+length(num(part(m,k))); write summ;
  if(summ<(length(den(g))+length(num(g)))) then g:=part(m,k);
>>;
write "g is ", g;
%for each f in m do <<
  %                    if(length(f)<length(g)) then g:=f;
 %                  >>;
%			% gets smallest element in the list
%write "g is ", g;
%if(limit(g,x,infinity)=infinity) then g:=1/g; write "g is ", g;

for each f in m do
  <<
    c:=limit(log(f)/log(g),x,infinity);
    A:=e^(log(f)-c*log(g));
    f:=A*w^c;
    ans_list:=append({f}, ans_list);
  >>;
return ans_list;
end;

%trst rewrite
%rewrite({e^(-x),e^x})
%h:=e^(-x/(1+e^-x))
%rewrite({e^(x-h),e^-(x/(1+h)),h,e^x,e^(-x)})

expr procedure smallest(li);
begin scalar current,k;
current:=part(li,1);
for k:=1:arglength(li) do <<
             if(length(current)>length(part(li,k))) then
             current:=part(li,k);
                           >>;
return current;
end;

expr procedure smallest(li);
begin scalar l1,l2;
if(length li=1) then return part(li,1) else 
            <<
l1:=lngth2(part(li,1));
l2:=lngth2(part(li,2));
if(l1>l2) then return part(li,2) else
 << if(l1<l2) then return part(li,1) else return part(li,1); >>;
             >>;

end;

symbolic procedure lngth u;
 begin
     if(u='list) then return nil else
    <<
      if(atom u) then return 1 else
      <<
        if(atom car u) then return (1+lngth cdr u)
        else return lngth car u + lngth cdr u;
      >>;
    >>;
 end;

%put('lngth2,'psopfn,'lngth);
algebraic;
procedure lngth2 u; lisp lngth u;

%-------------------------------------------------------------------------

% main routine to compute limits of exp-log functions as the variable tends
% to infinity.

operator x;
operator series;
algebraic;
expr procedure mrv_limit(f,var,val);

begin scalar mrv_f,mrv1_f,w, mrv_f2,tt, lead_term, series_exp,f1, small, rule1,
const, e0,sig,rule2,k, nu,de,h,recurse;
off mcd; off factor; off rational; off exp; off precise;
%if(val neq infinity) then return sub(var=val,f);

% trigonometric expressions aren't dealt with by the algorithm
 
if(not(freeof(f,sin))) then rederr "input not an exp log function";
if(not(freeof(f,cos))) then rederr "input not an exp log function";  
if(not(freeof(f,tan))) then rederr "input not an exp log function";
if(freeof(f,var)) then % possible cases: f can be a number, an expression
		       % independent of var, or possibly not an exp log
		       % function. In all cases, return f as the answer
return f;
if(val neq infinity) then return sub(var=val,f); 

%on rational;
%on rounded;
%this checks for numbers. red doesn't recognise some objects as numbers unless
% the rounded switch is on

f1:=f; if(numberp(f1)) then return f; off rational;
off rounded;
%off mcd;
%on rational; 
%on rounded; % maybe a risk, but we'll try it
if(var neq x) then f:=sub(var=x,f) else nil;
if(numberp(f)) then return f;
%%%% special case where f=e, or pi. Don't want to use the on rounded switch
%%%% in these cases
if(f=e) then return e;
if(f=pi) then return pi;
%if(f=x) then return plus_infinity;
on factor; off mcd; mrv_f:=mrv1(f); 
%write "*********************************************************************";
if(lisp !*tracelimit) then
write "mrv_f is ", mrv_f;
%write "*********************************************************************";
off factor; 
%on mcd;

if(member(x,mrv_f)) then
  <<
    off exp; off mcd; 
    while(member(x,mrv_f)) do
      <<
        f:=move_up(f,x); %write "f is ", f;
        %mrv_f:=mrv1(f);
	mrv_f:=for k:=1:arglength(mrv_f) collect move_up(part(mrv_f,k),x);
        %write "mrv is ", mrv_f;
      >>;
    % we know that x was a member of mrv(f), so now, the mrv set will contain
    % e^x at least. Hence, write directly in terms of w=e^(-x)
    %mrv_f2:=rewrite(mrv_f); % write mrv_f2;
 
    small:=e^(-x);
    % now have the smallest element, but don't know its limiting behaviour
    % if lim is infinity, need to set w to 1/small
    rule1:={small => ww };
    let rule1;
    f:=f;
    on mcd; nu:=num(f); de:=den(f); f:=nu/de; off mcd;
    f:=f; % write f;
    clearrules rule1;
       % f now rewritten
   >>
  else <<
	 %mrv_f2:=rewrite(mrv_f); % write "f2 is ", mrv_f2;
		 % now need to rewrite f itself
	small:=smallest(mrv_f);
        h:=log(small); %write "h is ", h;
        if(mrv_limit(h,x,infinity)=infinity) then
         <<
            small:=small^-1;
            % write "small has been changed", small; 
         >>;
 
  	rule1:=
		{
		  small => ww,
                  1/small => ww^-1
 		 };
         let rule1; off mcd;
 	 %write "smallest f is ", small;
         f:=f;
         on mcd;
        % de:=den(f); nu:=num(f);
	% f:=nu/de;
         off mcd; %f:=f;
	 clearrules rule1;
	 %write "f is ", f;
         % now rewritten in terms of w, and mrv(f)=w hopefully
        >>;

% at this point, f has been rewritten. Now check lcof of series expansion


% lisp !*mcd:=nil; lisp !*factor:=nil; lisp !*exp:=t; lisp !*rational:=nil; 

off mcd; on factor; off exp; off rational;
series_exp:=taylor(f,ww,0,1);  
if(lisp !*tracelimit) then
   write "performing taylor on: ", f;

%off mcd; on exp; on factor; off rational;
if(not taylorseriesp series_exp and part(series_exp,0)=taylor)
  then rederr "could not compute the Taylor series expansion";

tt:=log(small); %write "tt is ", tt;

series_exp:=sub(log(ww)=tt,series_exp); %off mcd; off factor; off exp;
series_exp:=taylortostandard series_exp; 
if(lisp !*tracelimit) then write "series expansion is ", series_exp;
% should now have the lead term of the series expansion in terms of w
if(numberp(series_exp)) then return series_exp
else <<
       if(lisp !*tracelimit) then
       write "series is ", series_exp;
       off rational;  off mcd;  off exp; off factor;
       series_exp:=series_exp; off factor;
       const:=coeffn(series_exp,ww,0); %write "const is ", const;
       if(const neq 0) then  
       <<
	 if(numberp(const)) then return const else
         << 
	     if(lisp !*tracelimit) then
             write "performing recursion on ", const;
             off mcd; on factor; off exp;  on rational; const:=const;
         return mrv_limit(const,x,infinity); >>;
       >> 
       else
       <<
       % need to look at exponent of ww. If e0>0 then return 0, if
       % e0<0 return infinity, if e0=0 return mrv_limit(c)
       %write "series_exp is ", series_exp; 
       series exp:=series_exp; off mcd; % try it here!
       %if(lisp !*tracelimit) then
       %write "series exp is  ", series_exp;
       %  series_exp:=lisp reval series_exp;
       %off mcd;
       e0:=find_least_expt(series_exp); 
	if(lisp !*tracelimit) then write "leading exponent e0 is ", e0;
       off mcd; on factor; off exp; off rational;
       %if(part(e0,1)=expt) then 
       %<< e0:=part(e0,2);         
       %e0:=part(e0,1); if(e0=e) then return e;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% possible cases: e0:={expt_list,num_list,x_list}
% if both num_list and x_list are empty, then e0 takes the value of the
% smallest exponent in expt_list.
% if numbers in expt_list are all positive, then e0 is the value in either
% num_list, or expt_list: if num_list, then this number is returned, if
% expt_list, then we apply algorithm recusively to the value of x_exp
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       if(part(e0,1)=expt) then <<
       					e0:=part(e0,2);
       %if(not numberp e0) then return part(e0,2);
       if(e0<0 and ((lisp car series_exp) neq minus)) 
	then return infinity else 
                << %*sign(lcof(series_exp,ww)) 
          	   if(e0<0) then return -infinity
		                 
       			else 
       			<<
       			if(e0>0) then return 0 else <<
			  off mcd; off factor; off rational; on exp;
                           recurse:=lcof(series_exp,ww); return 
          		mrv_limit(recurse,x,infinity); >>;
       			>>;
                 >>;
				 >>
       else <<
       if(part(e0,1)=number) then return part(e0,2) else 
	<<if sqchk part(e0,1) = 'minus then return -infinity else nil>>
            >>;
       e0:=lpower(series_exp,ww); off mcd;
       % expr free of neg powers ? sort of
       if(numberp e0) then <<
        on mcd; on exp; e0:=lpower(num series_exp,ww)/lpower(den series_exp,ww);
        off mcd; e0:=e0;  >>;
      % if(numberp e0) then << % we have to go back to series_exp
        %  on mcd; on rational;  series_exp:=series_exp; 
         %    temp:=lpower(num(series_exp),ww)/lpower(den(series_exp),ww);
         %    off mcd; temp:=temp; e0:=temp; lisp (e0:=lisp prepsq cadr e0);
         %                  >>
      % if(numberp(e0)) then <<
       %     if(e0>0) then return infinity else <<
       %if(e0<0) then return -infinity else return 0 >>;
     %                       >>;
                               % >>;
       lisp (e0:=lisp prepsq cadr e0);
            
      % off rational; off mcd; series_exp:=series_exp;
       %write "series_exp is ", series_exp;
       if(e0=ww) then return 0 else
       <<
         if(part(e0,0)=expt) then
          <<
            if(part(e0,2)<0) then % have plus /minus infinity
               <<
                 %write "sign is ", sign(lcof(series_exp,ww));
                 %write "expt is ", part(e0,2);
                   if(sign(lcof(series_exp,ww))=-1) then
                   return -infinity else  << if
                  (sign(lcof(series_exp,ww))=1) then return infinity
                 else << 
                         rule2:= {
                                   sign(log(~x)) => 1
                                 };
                         let rule2;
                         sig:=sign(lcof(series_exp,ww));
                         clearrules rule2;
                         return infinity; 
                       >>;
                                                 >>;
               >> 
	    else
            <<
              if(part(e0,2)>0) then return 0 else 
                  % the limiting behaviour of f depends on that of c 
                << on rational;
		off mcd; return mrv_limit(lcof(series_exp,ww),x,infinity); >>; 
                >>;
             >>;
          >>;
        >>;
     >>;
off mcd;
end;
%---------------------------------------------------------------------------
%a:=log(w^-1+x)-x
%a:=a*w^-1
%a:=a/(log(w^-1+log(x)))
%taylor(a,w,0,3)
%off mcd
%my_limit(e^-x,x,infinity)
%clear ww
%off mcd
%my_limit(e^x,x,infinity)
%trst my_limit
%ex:=log(log(x)+log(log(x)))-log(log(x))
%ex:=ex/(log(log(x)+log(log(log(x)))))
%ex:=ex*log(x)
% now need to see about the recursion part. Create test file, and sort out
% mrv function, and sort out extensions.
%trst mrv
%trst my_limit

%-----------------------------------------------------------------------------
% for examples, please see the test flie included with the documentation. The
% examples labelled from ex1 to ex12 are all taken from Dominik Gruntz' Thesis
% paper. Most are complicated examples, as he himself admits. This package 
% does not use the standard limits operator in REDUCE: it has been written to
% be independent, and can be presented as a separted facility in REDUCE.

%-----------------------------------------------------------------------------
endmodule;
end;


expr procedure test(exp);
  if(freeof(exp,ww^-1) and freeof(exp,ww^-2)) then t else nil;

% we need a procedure to determine if an expression is free of negative
% powers of w

%$symbolic procedure free_of_neg_power(exp,var);
%if(freeof(exp,ww^a)) then t else nil where numberp(a) and a<0;

%procedure free(exp,var); lisp free_of_neg_power(exp,var);
end;


symbolic procedure least_exp(exp);
begin scalar expon,base,temp,expon2;
      expon:=0;
      while(exp) do <<
      if(car exp='expt) then 
      << 
      base:=cadr exp;
      temp:=caddr exp;
      if(temp<=expon) then expon:=temp;
       >> else << while exp do << exp:=cdr exp; return least_exp(exp); >>;
               >>;
 return expon;
end;

procedure least(exp); lisp least_exp(exp);