Artifact 48575f780adc3474ce08f69815e2d2cc64006b67b6eb798669d7c1be57640023:
- Executable file
r38/packages/misc/limits.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: 23420) [annotate] [blame] [check-ins using] [more...]
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;