%----------------------------------------------------------------------------
% |
% 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);