Artifact a6377def83165b6c11712cff1108eb5a30ff4bd7ee16f26b282f5ea2695da25c:
- Executable file
r37/packages/mrvlimit/mrvlimit.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 20603) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/mrvlimit/mrvlimit.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 20603) [annotate] [blame] [check-ins using]
%---------------------------------------------------------------------------- % | % 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);