Artifact e7d9790bfa5b65c7c8e24abacb4c25dcdb2ad25c7ae2d2b07baf9682560e75de:
- Executable file
r36/src/misc.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: 67599) [annotate] [blame] [check-ins using] [more...]
module misc; % Miscellaneous algebraic code. create!-package('(misc contfr limits pf sum),'(contrib misc)); endmodule; module contfr; % Simultaneous approximation of a real number by a % continued fraction and a rational number with optional % user controlled precision (upper bound for numerator). % Author: Herbert Melenk. symbolic procedure contfract2 (u,b1); % compute continued fraction until either numerator exceeds b1 % or approximation has reached system precision. begin scalar b0,l,a,b,g,gg,ggg,h,hh,hhh; b:= u; g:=0; gg:=1; h:=1; hh:=0; if null b1 then b0:= absf !:times(b,!:expt(10,minus precision 0)); loop: a:=rd!-fix b; ggg:=a*gg + g; hhh:=a*hh + h; if b1 and abs hhh > b1 then goto ret; g := gg; gg:=ggg; h:=hh; hh:=hhh; l:=a.l; b:=!:difference(b,a); if null b or !:lessp(absf !:difference(!:quotient(i2rd!* gg,i2rd!* hh),u) ,b0) then goto ret; b:=!:quotient(1,b); goto loop; ret: return (gg.hh) . reversip l; end; symbolic procedure !:lessp(u,v); !:minusp !:difference(u,v); symbolic procedure rd!-fix u; if atom cdr u then fix cdr u else ashift(cadr u,cddr u); symbolic procedure contfract1(u,b); begin scalar oldmode,v; if eqcar(v:=u,'!:rd!:) then goto c; oldmode := get(dmode!*,'dname).!*rounded; if car oldmode then setdmode(car oldmode,nil); setdmode('rounded,t); !*rounded := t; v:=reval u; setdmode('rounded,nil); if car oldmode then setdmode(car oldmode,t); !*rounded:=cdr oldmode; if eqcar(v,'minus) and (numberp cadr v or eqcar(cadr v,'!:rd!:)) then v:=!:minus cadr v; if fixp v then return (v . 1).{v}; if not eqcar(v,'!:rd!:) then typerr(u,"continued fraction argument"); c: return contfract2(v,b); end; symbolic procedure cont!-fract u; <<u:=contfract1(car u,if cdr u then ieval cadr u); {'list,{'quotient,caar u,cdar u}, 'list . cdr u}>>; put('continued_fraction,'psopfn,'cont!-fract); endmodule; 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 '(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; <<x := prepsq simp!* x; if countof('log,x)>1 then logsort x else x>>; 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,args,old,v; !*protfg := t; fn := car u; args := cdr u; old := get('cot,'opmtch); put('cot,'opmtch, '(((!~x) (nil . t) (quotient (cos !~x) (sin !~x)) nil))); v := errorset!*({'apply,mkquote fn,mkquote args},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 eq '(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 eq '(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; module pf; % Compute partial fractions for an expression. % Author: Anthony C. Hearn. Comment PF is the top level operator for finding the partial fractions of an expression. It returns the partial fractions as a list. The algorithms used here are relatively unsophisticated, and use the extended Euclidean algorithm to break up expressions into factors. Much more sophisticated algorithms exist in the literature; fluid '(!*exp !*limitedfactors !*gcd kord!*); symbolic operator pf; flag('(pf),'noval); % Since PF will do its own simplification. symbolic procedure pf(u,var); % Convert an algebraic expression into partial fractions. begin scalar !*exp,!*gcd,kord!*,!*limitedfactors,polypart,rfactor, u1,u2,u3,u4,var,x,xx,y; !*exp := !*gcd := t; xx := updkorder var; % Make var the main variable. x := subs2 resimp simp!* u; % To allow for OFF EXP forms. u1 := denr x; if degr(u1,var) = 0 then <<setkorder xx; return list('list,u)>>; u2 := qremsq(!*f2q numr x,!*f2q u1,var); %Extract polynomial part. if caar u2 then polypart := car u2; rfactor := 1 ./ 1; % Factor for rational part. u2 := cdr u2; u3 := fctrf u1; % Factorize denominator. x := cdr u3; u3 := car u3; % Process monomial part. while not domainp u3 do <<if mvar u3 eq var then x := (!*k2f mvar u3 . ldeg u3) . x else <<u4 := !*p2f lpow u3; rfactor := numr rfactor ./ multf(u4,denr rfactor); u1 := quotf(u1,u4)>>; u3 := lc u3>>; if u3 neq 1 then <<rfactor := numr rfactor ./ multf(u3,denr rfactor); u1 := quotf(u1,u3)>>; % Separate power factors in denominator. while length x>1 do <<u3 := exptf(caar x,cdar x); u1 := quotf(u1,u3); if degr(u3,var)=0 then rfactor := numr rfactor ./ multf(u3,denr rfactor) % then <<rfactor := numr rfactor ./ multf(u3,denr rfactor); % u2 := nil>> else <<u4 := xeucl(u1,u3,var); % Remove spurious polynomial in numerator. y := (remsq(multsq(car u4,u2),!*f2q u3,var) . car x) . y; u2 := multsq(cdr u4,u2)>>; x := cdr x>>; u3 := exptf(caar x,cdar x); if u2 = (nil ./ 1) then nil else if degr(u3,var)=0 then rfactor := numr rfactor ./ multf(u3,denr rfactor) % Remove spurious polynomial in numerator. else y := (remsq(u2,!*f2q u3,var) . car x) . y; x := nil; % Finally break down non-linear terms in denominator. for each j in y do if cddr j =1 then x := j . x else x := append(pfpower(car j,cadr j,cddr j,var),x); % Fold in rfactor. if rfactor neq '(1 . 1) then x := for each j in x collect multsq(rfactor,car j) . cdr j; x := for each j in x collect list('quotient,prepsq!* car j, if cddr j=1 then prepf cadr j else list('expt,prepf cadr j,cddr j)); if polypart then x := prepsq!* polypart . x; setkorder xx; return 'list . x end; symbolic procedure xeucl(u,v,var); % Extended Euclidean algorithm with rational coefficients. % I.e., find polynomials Q, R in var with rational coefficients (as % standard quotients) such that Q*u + R*v = 1, where u and v are % relatively prime standard forms in variable var. Returns Q . R. begin scalar q,r,s,w; q := list(1 ./ 1,nil ./ 1); r := list(nil ./ 1,1 ./ 1); if degr(u,var) < degr(v,var) then <<s := u; u := v; v := s; s := q; q := r; r := s>>; u := !*f2q u; v := !*f2q v; while numr v do <<if degr(numr v,var)=0 then w := quotsq(u,v) . (nil ./ 1) else w := qremsq(u,v,var); s := list(addsq(car q,negsq multsq(car w,car r)), addsq(cadr q,negsq multsq(car w,cadr r))); u := v; v := cdr w; q := r; r := s>>; v := lnc numr u ./ denr u; % Is it possible for this not to be % in lowest terms, and, if so, does % it matter? r := quotsq(v,u); return multsq(r,quotsq(car q,v)) . multsq(r,quotsq(cadr q,v)) end; symbolic procedure qremsq(u,v,var); % Find rational quotient and remainder (as standard quotients) % dividing standard quotients u by v wrt var. % This should really be done more directly without using quotsq. (quotsq(addsq(u,negsq x),v) . x) where x=remsq(u,v,var); symbolic procedure remsq(u,v,var); % Find rational and remainder (as a standard quotient) on % dividing standard quotients u by v wrt var. begin integer m,n; scalar x; n := degr(numr v,var); if n=0 then rederr list "Remsq given zero degree polynomial"; while (m := degr(numr u,var))>= n do <<if m=n then x := v else x := multsq(!*p2q(var.**(m-n)),v); u := addsq(u, negsq multsq(multf(lc numr u,denr v) ./ multf(lc numr v,denr u), x))>>; return u end; symbolic procedure pfpower(u,v,n,var); % Convert u/v^n into partial fractions. begin scalar x,z; while degr(numr u,var)>0 do <<x := qremsq(u,!*f2q v,var); z := (cdr x . v . n) . z; n := n-1; u := car x>>; if numr u then z := (u . v . n) . z; return z end; endmodule; module sum; % Package for summation in finite terms. % Authors: K.Yamamoto, K.Kishimoto & K.Onaga Hiroshima Univ. % Modified by: F.Kako Hiroshima Univ. % Fri Sep. 19, 1986 % Mon Sep. 7, 1987 added PROD operator (by F. Kako) % e-mail: kako@kako.math.sci.hiroshima-u.ac.jp % or D52789%JPNKUDPC.BITNET % Usage: % sum(expression,variable[,lower[,upper]]); % lower and upper are optionals. % prod(expression,variable[,lower[,upper]]); % returns product of expression. fluid '(!*trsum); % trace switch; fluid '(prod_last_attempt_rules!* sum_last_attempt_rules!*); switch trsum; symbolic procedure simp!-sum u; %ARGUMENT CAR U: expression of prefix form; % CADR U: kernel; % CADDR U: lower bound; % CADDDR U: upper bound; %value : expression of sq form; begin scalar v,y,upper,lower,lower1,dif; y := cdr u; u := simp!* car u; if null numr u then return (nil ./ 1) else if atom y then return !*p2f(car fkern(list('sum,prepsq u)) .* 1) ./ 1; if not atom cdr y then << lower := cadr y; lower1 := if numberp lower then lower - 1 else list('plus,lower,-1); upper := if not atom cddr y then caddr y else car y; dif := addsq(simp!* upper, negsq simp!* lower); if denr dif = 1 then if null numr dif then return subsq(u,list(!*a2k car y . upper)) else if fixp numr dif then dif := numr dif else dif := nil else dif := nil; if dif and dif <= 0 then return nil ./ 1; if atom cddr y then upper := nil>>; v := !*a2k car y; return simp!-sum1(u,v,y,upper,lower,lower1,dif) end; symbolic procedure simp!-sum1(u,v,y,upper,lower,lower1,dif); begin scalar w,lst,x,z,flg; lst := sum!-split!-log(u,v); w := car lst; lst := cdr lst; u := nil ./ 1; a: if null w then go to b; x := multsq(caar w, simp!-prod1(cdar w,v,y,upper,lower,lower1,dif)); u := addsq(u,simp!* list('log, prepsq x)); w := cdr w; go to a; b: if null lst then return u; flg := nil; z := car lst; if !*trsum then << prin2!* "Summation ";sqprint z;prin2!* " w.r.t "; xprinf(!*k2f v,nil,nil);terpri!* t >>; % z := reorder numr z ./ reorder denr z; w := sum!-sq(z,v); if w = 'failed then << if !*trsum then << prin2!* "UMM-SQ failed. Trying SUM-TRIG"; terpri!* t>>; w := sum!-trig(z,v); if w = 'failed then << if !*trsum then << prin2!* "SUM-TRIG failed."; terpri!* t>>; w := sum!-unknown(z,v,y,lower,dif); flg := car w; w := cdr w>>>>; if !*trsum then << prin2!* "Result = "; sqprint w; terpri!* t >>; if flg then goto c; if upper then w := addsq(subsq(w,list(v . upper)), negsq subsq(w,list(v . lower1))) else if lower then w := addsq(w , negsq subsq(w, list(v . lower1))); c: u := addsq(u,w); lst := cdr lst; goto b end; put('sum,'simpfn,'simp!-sum); %********************************************************************* % Case of trigonometric or other functions % Trigonometric functions are expressed in terms of exponetials. % Pattern matching to get the summation in closed form. %********************************************************************; global '(!*trig!-to!-exp); % variable to indicate % that the expression contains % some trig. functions. symbolic procedure sum!-trig(u,v); begin scalar lst,w; % z; !*trig!-to!-exp := nil; % trig. to exponential. u := trig!-to!-expsq(u,v); if not !*trig!-to!-exp then return 'failed; lst := sum!-term!-split(u,v); u := nil ./ 1; a: if null lst then return exp!-to!-trigsq u; % z := reorder numr car lst ./ reorder denr car lst; % w := sum!-sq(z,v); w := sum!-sq(car lst,v); if w = 'failed then return 'failed; % w := exp!-to!-trigsq w; % exponential to trig. function. u := addsq(u,w); lst := cdr lst; goto a end; sum_last_attempt_rules!* := algebraic << { sum(~f + ~g,~n,~anf,~ende) => sum(f,n,anf,ende) + sum(g,n,anf,ende) when or (part(sum(f,n,anf,ende),0) neq sum , part(sum(g,n,anf,ende),0) neq sum ), sum((~f+~g)/~dd,~n,~anf,~ende) => sum(f/dd,n,anf,ende) + sum(g/dd,n,anf,ende) when or (part(sum(f/dd,n,anf,ende),0) neq sum , part(sum(g/dd,n,anf,ende),0) neq sum ), sum(~c*~f,~n,~anf,~ende) => c* sum(f,n,anf,ende) when freeof(c,n) and c neq 1, sum(~c/~f,~n,~anf,~ende) => c* sum(1/f,n,anf,ende) when freeof(c,n) and c neq 1, sum(~c*~f/~g,~n,~anf,~ende) => c* sum(f/g,n,anf,ende) when freeof(c,n) and c neq 1} >>$ symbolic procedure sum!-unknown(u,v,y,lower,dif); begin scalar z,w; if null dif then << z := 'sum . (prepsq u . list car y); if w := opmtch z then return (nil . simp w) else if null cdr y then return (t . mksq(z,1)); z := 'sum . (prepsq u . y); let sum_last_attempt_rules!*; w:= opmtch z; rule!-list (list sum_last_attempt_rules!*,nil); return (t . if w then simp w else mksq(z,1))>>; % return (t . if w := opmtch z then simp w else mksq(z,1))>>; z := nil ./ 1; a: if dif < 0 then return (t . z); z := addsq(z,subsq(u,list(v . list('plus,lower,dif)))); dif := dif - 1; goto a end; %********************************************************************* % Summation by Gosper's algorithm. %********************************************************************; symbolic procedure sum!-sq(u,v); %Argument U : expression of s-q; % V : kernel; %value : expression of sq (result of summation.); begin scalar gn,fn,pn,rn,qn,z,k,x; if null numr u then return nil ./ 1; x := setkorder list v; z := reorder numr u; u := z ./ reorder denr u; if !*trsum then << prin2t " *** Summation by Gosper's algorithm ***"; prin2!* " A(n) = "; sqprint u;terpri!* t; terpri!* t>>; if domainp z or not (mvar z eq v) or red z then <<pn := 1; z := u>> else <<pn := !*p2f lpow z; z := lc z ./ denr u>>; z := quotsq(z,nsubsq(z,v, - 1)); gn := gcdf!*(numr z,denr z); if !*trsum then << prin2!* "A(n)/A(n-1) = ";sqprint z;terpri!* t; prin2!* "GN = ";xprinf(gn,nil,nil);terpri!* t>>; qn := quotf!*(numr z, gn); rn := quotf!*(denr z, gn); if nonpolyp(qn,v) or nonpolyp(rn,v) then go to fail; if !*trsum then << prin2!* "Initial qn, rn and pn are "; terpri!* t; prin2!* "QN = ";xprinf(qn,nil,nil);terpri!* t; prin2!* "RN = ";xprinf(rn,nil,nil);terpri!* t; prin2!* "PN = ";xprinf(pn,nil,nil);terpri!* t>>; k := compress explode '!+j; z := integer!-root(resultant(qn,nsubsf(rn,v,k),v),k); if !*trsum then << prin2 "Root of resultant(q(n),r(n+j)) are "; prin2t z >>; while z do << k := car z; gn := gcdf!*(qn,nsubsf(rn,v,k)); qn := quotf!*(qn,gn); rn := quotf!*(rn,nsubsf(gn,v, -k)); while (k := k - 1)>=0 do pn := multf(pn,nsubsf(gn,v, -k)); z := cdr z>>; if !*trsum then << prin2!* "Shift free qn, rn and pn are";terpri!* t; prin2!* "QN = ";xprinf(qn,nil,nil);terpri!* t; prin2!* "RN = ";xprinf(rn,nil,nil);terpri!* t; prin2!* "PN = ";xprinf(pn,nil,nil);terpri!* t>>; qn := nsubsf(qn,v,1); if (k := degree!-bound(pn,addf(qn,rn),addf(qn,negf rn),v)) < 0 then go to fail; if !*trsum then << prin2 "DEGREE BOUND is "; prin2t k >>; if not(fn := solve!-fn(k,pn,qn,rn,v)) then go to fail; if !*trsum then << prin2!* "FN = ";sqprint fn;terpri!* t >>; u := multsq(multsq(qn ./ 1,fn), multsq(u, 1 ./ pn)); z := gcdf!*(numr u, denr u); u := quotf!*(numr u, z) ./ quotf!*(denr u,z); if !*trsum then << prin2t " *** Gosper's algorithm completed ***"; prin2!* " S(n) = "; sqprint u;terpri!* t; terpri!* t>>; setkorder x; return (reorder numr u ./ reorder denr u); fail: if !*trsum then << prin2t " *** Gosper's algorithm failed ***"; terpri!* t>>; setkorder x; return 'failed end; %********************************************************************* % integer root isolation %********************************************************************; symbolic procedure integer!-root(u,v); % Produce a list of all positive integer root of U; begin scalar x,root,n,w; x := setkorder list v; u := reorder u; if domainp u or not(mvar u eq v) then go to a; u := numr cancel(u ./ lc u); w := u; % get trailing term; while not domainp w and mvar w eq v and cdr w do w := cdr w; if (n := degr(w,v)) > 0 then << w := lc w; while n > 0 do << root := 0 . root; n := n - 1>>>>; n := dfactors lowcoef w; % factor tail coeff. w := (v . 1) . 1; while n do << if not testdivide(u,v,car n) then << root := car n . root; u := quotf!*(u, (w . - car n))>> else n := cdr n>>; a: setkorder x; return root end; symbolic procedure lowcoef u; begin scalar lst,m; lst := dcoefl u; m := 0; a: if null lst then return m; m := gcdn(m,car lst); if m = 1 then return 1; lst := cdr lst; go to a end; symbolic procedure dcoefl u; if domainp u then if fixp u then list abs u else nil else nconc(dcoefl lc u , dcoefl red u); symbolic procedure testdivide(u,v,n); % Evaluate U at integer point (V = N); begin scalar x; a: if domainp u or not(mvar u eq v) then return addf(u,x); x := addf(multd(expt(n,ldeg u),lc u),x); if (u := red u) then go to a; return x end; %********************************************************************* %********************************************************************; symbolic procedure degree!-bound(pn,u,v,kern); % degree bound for fn; % u: q(n+1) + r(n); % v: q(n+1) - r(n); begin scalar lp,l!+, l!-, x,m,k; x := setkorder list kern; u := reorder u; v := reorder v; pn := reorder pn; l!+ := if u then degr(u,kern) else -1; l!- := if v then degr(v,kern) else -1; lp := if pn then degr(pn,kern) else -1; if l!+ <= l!- then <<k := lp - l!-;go to a>>; k := lp - l!+ + 1; if l!+ > 0 then u := lc u; if l!- > 0 then v := lc v; if l!+ = l!- + 1 and fixp(m := quotf1(multd(-2,v),u)) then k := max(m,k) else if lp = l!- then k := max(k,0); a: setkorder x; return k end; %********************************************************************* % calculate polynomial f(n) such that % p(n) - q(n+1)*f(n) + r(n)*f(n-1) = 0; %********************************************************************; symbolic procedure solve!-fn(k,pn,qn,rn,v); begin scalar i,fn,x,y,z,u,w,c,clst,flst; c := makevar('c,0); clst := list c; fn := !*k2f c; i := 0; while (i := i + 1) <= k do << c := makevar('c,i); clst := c . clst; fn := ((v . i) . !*k2f c) . fn>>; z := addf(pn, addf(negf multf(qn,fn), multf(rn,nsubsf(fn,v, - 1)))); x := setkorder (v . clst); z := reorder z; c := clst; if !*trsum then << prin2!* "C Equation is";terpri!* t; xprinf(z,nil,nil);terpri!* t >>; a: if domainp z or domainp (y := if mvar z eq v then lc z else z) then go to fail; w := mvar y; if not(w memq clst) then go to fail; if !*trsum then << prin2!* "C Equation to solve is ";xprinf(y,nil,nil);terpri!* t; prin2!* " w.r.t ";xprinf(!*k2f w,nil,nil);terpri!* t >>; u := gcdf!*(red y , lc y); u := quotf!*(negf red y, u) ./ quotf!*(lc y, u); flst := (w . u) . flst; z := subst!-cn(z,w,u); if !*trsum then << xprinf(!*k2f w,nil,nil);prin2!* " := ";sqprint u;terpri!* t >>; clst := deleteq(clst,w); if z then go to a; setkorder c; fn := reorder fn; u := 1; while not domainp fn and mvar fn memq c do << w := mvar fn; z := atsoc(w,flst); fn := subst!-cn(fn,w,if z then cdr z); if z then u := multf(u,denr cdr z); z := gcdf!*(fn,u); fn := quotf!*(fn,z); u := quotf!*(u,z)>>; setkorder x; return cancel(reorder fn ./ reorder u); fail: if !*trsum then << prin2t "Fail to solve C equation."; prin2!* "Z := ";xprinf(z,nil,nil);terpri!* t >>; setkorder x; return nil end; symbolic procedure subst!-cn(u,v,x); begin scalar z; z := setkorder list v; u := reorder u; if not domainp u and mvar u eq v then if x then u := addf(multf(lc u,reorder numr x), multf(red u,reorder denr x)) else u := red u; setkorder z; return reorder u end; symbolic procedure makevar(id,n); compress nconc(explode id, explode n); symbolic procedure deleteq(u,x); if null u then nil else if car u eq x then cdr u else car u . deleteq(cdr u, x); symbolic procedure nsubsf(u,kern,i); % ARGUMENT U : expression of sf; % KERN : kernel; % I : integer or name of integer variable; % value : expression of sf; begin scalar x,y,z,n; if null i or i = 0 then return u; x := setkorder list kern; u := reorder u; y := addf(!*k2f kern, if fixp i then i else !*k2f i); z := nil; a: if domainp u or not(mvar u eq kern) then goto b; z := addf(z,lc u); n := degr(u,kern) - degr(red u,kern); u := red u; a1: if n <= 0 then goto a; z := multf(z,y); n := n - 1; go to a1; b: z := addf(z,u); setkorder x; return reorder z end; symbolic procedure nsubsq(u,kern,i); % ARGUMENT U : expression of sq; % KERN : kernel; % I : integer or name of integer variable; % value : expression of sq; subsq(u,list(kern . list('plus, kern, i))); %********************************************************************* % dependency check %********************************************************************; symbolic procedure nonpolyp(u,v); % check U is not a polynomial of V; if domainp u then nil else (not(mvar u eq v) and depend!-p(mvar u,v)) or nonpolyp(lc u,v) or nonpolyp(red u,v); symbolic procedure depend!-sq(u,v); depend!-f(numr u,v) or depend!-f(denr u,v); symbolic procedure depend!-f(u,v); if domainp u then nil else depend!-p(mvar u,v) or depend!-f(lc u,v) or depend!-f(red u,v); symbolic procedure depend!-p(u,v); if u eq v then t else if atom u then nil else if not atom car u then depend!-f(u,v) else if car u eq '!*sq then depend!-sq(cadr u, v) else depend!-l(cdr u, v); symbolic procedure depend!-l(u,v); if null u then nil else if depend!-sq(simp car u, v) then t else depend!-l(cdr u,v); %********************************************************************* % term splitting %********************************************************************; symbolic procedure sum!-term!-split(u,v); begin scalar y,z,klst,lst,x; x := setkorder list v; z := qremf(reorder numr u, y := reorder denr u); klst := kern!-list(car z,v); lst := termlst(car z, 1 ./ 1, klst); klst := kern!-list(cdr z,v); if depend!-f(y,v) then klst := deleteq(klst,v); lst := append(lst, termlst(cdr z, 1 ./ y, klst)); setkorder x; return lst end; symbolic procedure kern!-list(u,v); % Returns list of kernels that depend on V; begin scalar x; for each j in kernels u do if depend!-p(j,v) then x := j . x; return x end; symbolic procedure termlst(u,v,klst); begin scalar x,kern,lst; if null u then return nil else if null klst or domainp u then return list multsq(!*f2q u,v); kern := car klst; klst := cdr klst; x := setkorder list kern; u := reorder u; v := reorder(numr v) ./ reorder(denr v); while not domainp u and mvar u eq kern do << lst := nconc(termlst(lc u, multsq(!*p2q lpow u, v),klst),lst); u := red u>>; if u then lst := nconc(termlst(u,v,klst),lst); setkorder x; return lst end; %********************************************************************* % Express trigonometric functions (such as sin, cos ..) % by exponentials. %********************************************************************; symbolic procedure trig!-to!-expsq(u,v); multsq(trig!-to!-expf(numr u,v), invsq trig!-to!-expf(denr u,v)); symbolic procedure trig!-to!-expf(u,v); if domainp u then u ./ 1 else addsq(multsq(trig!-to!-expp(lpow u,v), trig!-to!-expf(lc u,v)), trig!-to!-expf(red u,v)); symbolic procedure trig!-to!-expp(u,v); begin scalar !*combineexpt,w,x,z,n,wi; % We don't want to combine expt terms here, since the code % depends on the terms being separate. n := cdr u; % integer power; z := car u; % main variable; if atom z or not atom (x := car z) or not depend!-p(z,v) then return !*p2q u; if x memq '(sin cos tan sec cosec cot) then << !*trig!-to!-exp := t; w := multsq(!*k2q 'i, simp!* cadr z); w := simp!* list('expt,'e, mk!*sq w); % W := SIMP LIST('EXPT,'E, 'TIMES . ( 'I . CDR Z)); wi := invsq w; if x eq 'sin then w := multsq(addsq(w ,negsq wi), 1 ./ list(('i .** 1) .* 2)) else if x eq 'cos then w := multsq(addsq(w, wi), 1 ./ 2) else if x eq 'tan then w := multsq(addsq(w,negsq wi), invsq addsq(w,wi)) else if x eq 'sec then w := multsq(2 ./ 1, invsq addsq(w, wi)) else if x eq 'cosec then w := multsq(list(('i .** 1) .* 2), invsq addsq(w, negsq wi)) else w := multsq(addsq(w, wi), invsq addsq(w, negsq wi)) >> else if x memq '(sinh cosh tanh sech cosech coth) then << !*trig!-to!-exp := t; w := simp!* list('expt,'e,cadr z); wi := invsq w; if x eq 'sinh then w := multsq(addsq(w,negsq wi), 1 ./ 2) else if x eq 'cosh then w := multsq(addsq(w,wi), 1 ./ 2) else if x eq 'tanh then w := multsq(addsq(w,negsq wi), invsq addsq(w,wi)) else if x eq 'sech then w := multsq(2 ./ 1, invsq addsq(w, wi)) else if x eq 'cosech then w := multsq(2 ./ 1, invsq addsq(w, negsq wi)) else w := multsq(addsq(w,wi), invsq addsq(w, negsq wi)) >> else return !*p2q u; return exptsq(w,n) end; %********************************************************************* % Inverse of trig!-to!-exp. % Express exponentials in terms of trigonometric functions % (sin, cos, sinh and cosh) % Wed Dec. 17, 1986 by F. Kako; %********************************************************************; symbolic procedure exp!-to!-trigsq u; multsq(exp!-to!-trigf numr u, invsq exp!-to!-trigf denr u); symbolic procedure exp!-to!-trigf u; begin scalar v,v1,x,y,n; u := termlst1(u,1,nil ./1); v := nil; a: if null u then go to b; x := caar u; y := cdar u; u := cdr u; a1: if u and y = cdar u then << x := addf(x,caar u); u := cdr u; go to a1>>; v := (x . y) . v; go to a; b: v1 := reverse v; n := length v; u := nil ./ 1; c: if n = 0 then return u else if n = 1 then return addsq(u, multsq(!*f2q caar v, simp!* list('expt,'e,mk!*sq cdar v))); u := addsq(u,exp!-to!-trigl(caar v1,caar v,cdar v1,cdar v)); v := cdr v; v1 := cdr v1; n := n - 2; go to c end; symbolic procedure exp!-to!-trigl(a,b,c,d); % A*E**C + B*E**D % --> % ((A+B)*COSH((C-D)/2)+(A-B)*SINH((C-D)/2))*E**((C+D)/2); % A, B: sf; % C, D: sq; begin scalar x,y,z; x := !*f2q addf(a,b); y := !*f2q addf(a, negf b); z := multsq(addsq(c,negsq d), 1 ./ 2); z := real!-imag!-sincos z; return multsq(simp!* list('expt,'e, mk!*sq multsq(addsq(c,d),1 ./ 2)), addsq(multsq(x, cdr z), multsq(y, car z))) end; symbolic procedure termlst1(u,v,w); %ARGUMENT U : sf; % V : sf; % W : sq; %value : list of (sf . sq); begin scalar x,y; if null u then return nil else if domainp u then return list (multf(u,v) . w); x := mvar u; y := if atom x or not(car x eq 'expt) or not(cadr x eq 'e) then termlst1(lc u,multf(!*p2f lpow u,v),w) else termlst1(lc u,v,addsq(w,multsq(simp!* caddr x, ldeg u ./ 1))); return nconc(y,termlst1(red u,v,w)) end; %********************************************************************* % Module --- COMPLEX; % Wed Dec. 17, 1986 by F. Kako; %********************************************************************; %****************************************************************** %******* SPLIT REAL AND IMAGINARY PART ****************** %****************************************************************** symbolic procedure real!-imag!-sq u; %U is a standard quotient, %Value is the standard quotient real part and imaginary part of U. begin scalar x,y; x := real!-imag!-f numr u; y := real!-imag!-f denr u; u := addf(multf(car y, car y), multf(cdr y, cdr y)); % Re Y **2 + Im Y **2; return (cancel(addf(multf(car x, car y), multf(cdr x, cdr y)) ./ u) . cancel(addf(multf(car y, cdr x), negf multf(car x, cdr y)) ./ u)) end; symbolic procedure real!-imag!-f u; %U is a standard form. %Value is the standard form real and imag part of U. begin scalar x; if domainp u then return u . nil; x := setkorder list 'i; u := reorder u; u := if mvar u eq 'i and ldeg u = 1 then red u . lc u else u . nil; setkorder x; return (reorder car u . reorder cdr u) end; %***************************************************************** % hyperbolic functions %*****************************************************************; symbolic procedure real!-imag!-sincos u; begin scalar v,w,z; v := real!-imag!-sq u; if null cadr v then << u := prepsq u; return simp!* list('sinh,u) . simp!* list('cosh,u)>> else if null caar v then << u := prepsq cdr v; return (multsq(!*k2q 'i, simp!* list('sin,u)) . simp!* list('cos,u))>>; u := prepsq cdr v; v := prepsq car v; w := simp!* list('cos,u); u := simp!* list('sin,u); u := multsq(!*k2q 'i,u); z := simp!* list('cosh,v); v := simp!* list('sinh,v); return (addsq (multsq(w, v), multsq(u,z))) . (addsq (multsq(w,z),multsq(u,v))) end; %********************************************************************* % module --- product % package for production in finite terms % % Author: F.Kako Hiroshima Univ. % Mon Sep. 7, 1987 % % usage: % prod(expression,variable[,lower[,upper]]); % lower and upper are optionals; % % %********************************************************************; %********************************************************************* %********************************************************************; symbolic procedure simp!-prod u; %ARGUMENT CAR U: expression of prefix form; % CADR U: kernel; % CADDR U: lower bound; % CADDDR U: upper bound; %value : expression of sq form; begin scalar v,y,upper,lower,lower1,dif; y := cdr u; u := simp!* car u; if null numr u then return (1 ./ 1) else if atom y then return u; if not atom cdr y then << lower := cadr y; lower1 := if numberp lower then lower - 1 else list('plus,lower,-1); upper := if not atom cddr y then caddr y else car y; dif := addsq(simp!* upper, negsq simp!* lower); if denr dif = 1 then if null numr dif then return u else if fixp numr dif then dif := numr dif else dif := nil else dif := nil; if dif and dif <= 0 then return 1 ./ 1; if atom cddr y then upper := nil>>; v := !*a2k car y; return simp!-prod1(u,v,y,upper,lower,lower1,dif) end; symbolic procedure simp!-prod1(u,v,y,upper,lower,lower1,dif); begin scalar w,lst,x,z,flg; lst := prod!-split!-exp(u,v); w := car lst; lst := cdr lst; u := 1 ./ 1; a: if null w then go to b; x := simp!-sum1(cdar w,v,y,upper,lower,lower1,dif); u := multsq(u,simpexpt list(caar w, prepsq x)); w := cdr w; go to a; b: if null lst then return u; flg := nil; z := car lst; if !*trsum then << prin2!* "Product ";sqprint z;prin2!* " w.r.t "; xprinf(!*k2f v,nil,nil);terpri!* t >>; % z := reorder numr z ./ reorder denr z; w := prod!-sq(z,v); if w = 'failed then << if !*trsum then << prin2!* "PROD-SQ failed."; terpri!* t>>; w := prod!-unknown(z,v,y,lower,dif); flg := car w; w := cdr w>>; if !*trsum then << prin2!* "Result = "; sqprint w;terpri!* t >>; if flg then goto c; if upper then w := multsq(subsq(w,list(v . upper)), invsq subsq(w,list(v . lower1))) else if lower then w := multsq(w , invsq subsq(w, list(v . lower1))); c: u := multsq(u,w); lst := cdr lst; goto b end; put('prod,'simpfn,'simp!-prod); %********************************************************************* % Case of other functions %********************************************************************; symbolic procedure prod!-unknown(u,v,y,lower,dif); begin scalar z,w,uu; if null dif then << z := 'prod . (prepsq u . list car y); if w := opmtch z then return (nil . simp w) else if null cdr y then return (t . mksq(z,1)); load_package 'factor; % try to find factors uu := factorize prepf numr u; if length uu > 2 then << z := 'times . foreach uuu in cdr uu collect ('prod . ( prepsq multsq(if pairp uuu and eq(car uuu,'!*sq) then cadr uuu else simp uuu,1 ./ denr u)) . y); return (t . simp z) >>; z := 'prod . (prepsq u . y); % try to apply rules let prod_last_attempt_rules!*; w:= opmtch z; rule!-list (list prod_last_attempt_rules!*,nil); return (t . if w then simp w else mksq(z,1))>>; % return (t . if w := opmtch z then simp w else mksq(z,1))>>; z := 1 ./ 1; a: if dif < 0 then return (t . z); z := multsq(z,subsq(u,list(v . list('plus,lower,dif)))); dif := dif - 1; goto a end; prod_last_attempt_rules!* := algebraic << { prod(~f * ~g,~n,~anf,~ende) => prod(f,n,anf,ende) * prod(g,n,anf,ende) when g neq 1 and or(numberp prod(f,n,anf,ende), part(prod(f,n,anf,ende),0) neq prod, part(prod(g,n,anf,ende),0) neq prod), prod(~f / ~g,~n,~anf,~ende) => prod(f,n,anf,ende) / prod(g,n,anf,ende) when g neq 1 and or(numberp prod(f,n,anf,ende), % 1? part(prod(f,n,anf,ende),0) neq prod, part(prod(g,n,anf,ende),0) neq prod), prod(expt(~f,~k),~n,~anf,~ende) => (for ii:=1:k product prod(f,n,anf,ende)) when neq(part(prod(f,n,anf,ende),0),prod)} >>; %********************************************************************* % Product of rational function %********************************************************************; symbolic procedure prod!-sq(u,v); %ARGUMENT U : expression of s-q; % V : kernel; %value : expression of sq (result of product.); begin scalar gn,p1n,p2n,rn,qn,z,k,x,y; if null numr u then return 1 ./ 1; x := setkorder list v; qn := reorder numr u; rn := reorder denr u; if !*trsum then << prin2t " *** Product of A(n) = qn/rn with ***"; prin2!* "QN = ";xprinf(qn,nil,nil);terpri!* t; prin2!* "RN = ";xprinf(rn,nil,nil);terpri!* t>>; if nonpolyp(qn,v) or nonpolyp(rn,v) then go to fail; k := compress explode '!+j; z := integer!-root2(resultant(qn,nsubsf(rn,v,k),v),k); if !*trsum then << prin2 "Root of resultant(q(n),r(n+j)) are "; prin2t z >>; p2n := p1n := 1; while z do << k := car z; gn := gcdf!*(qn,nsubsf(rn,v,k)); qn := quotf!*(qn,gn); rn := quotf!*(rn,nsubsf(gn,v, -k)); if k > 0 then while (k := k - 1)>=0 do << p1n := multf(p1n,nsubsf(gn,v, -k)); if y := prod!-nsubsf(gn,v,-k) then p2n := multf(p2n,y)>> else if k < 0 then while k < 0 do << p2n := multf(p2n,nsubsf(gn,v, -k)); if y := prod!-nsubsf(gn,v,-k) then p1n := multf(p1n,y); k := k + 1>>; z := cdr z>>; if depend!-f(qn,v) or depend!-f(rn,v) then go to fail; u := multsq(p1n ./ p2n, simpexpt list(prepsq (qn ./ rn), v)); if !*trsum then << prin2t " *** Product of rational function calculated ***"; prin2!* " P(n) = "; sqprint u;terpri!* t; terpri!* t>>; setkorder x; return (reorder numr u ./ reorder denr u); return u; fail: if !*trsum then << prin2t " *** Product of rational function failed ***"; terpri!* t>>; setkorder x; return 'failed end; symbolic procedure prod!-nsubsf(u,kern,i); % ARGUMENT U : expression of sf; % KERN : kernel; % I : integer; % value : expression of sf; begin scalar x,z,n; x := setkorder list kern; u := reorder u; z := nil; a: if domainp u or not(mvar u eq kern) then goto b; z := addf(z,lc u); n := degr(u,kern) - degr(red u,kern); u := red u; if i = 0 then if n = 0 then nil else z := nil else z := multf(z,expt(i,n)); go to a; b: z := addf(z,u); setkorder x; return reorder z end; %********************************************************************* % integer (positive and negative) root isolation %********************************************************************; symbolic procedure integer!-root2(u,v); % Produce a list of all integer root of U; begin scalar x,root,n,w; x := setkorder list v; u := reorder u; if domainp u or not(mvar u eq v) then go to a; u := numr cancel(u ./ lc u); w := u; % get trailing term; while not domainp w and mvar w eq v and cdr w do w := cdr w; if (n := degr(w,v)) > 0 then << w := lc w; while n > 0 do << root := 0 . root; n := n - 1>>>>; n := dfactors lowcoef w; % factor tail coeff.; w := (v . 1) . 1; while n do << if not testdivide(u,v,car n) then << root := car n . root; u := quotf!*(u, (w . - car n))>> else if not testdivide(u,v,- car n) then << root := (- car n) . root; u := quotf!*(u, (w . car n))>> else n := cdr n>>; a: setkorder x; return root end; %********************************************************************* % log and exponential term splitting for summation and product %********************************************************************; symbolic procedure sum!-split!-log(u,v); begin scalar x,y,z,lst,llst,mlst; lst := sum!-term!-split(u,v); a: if null lst then return (llst. mlst); u := car lst; lst := cdr lst; z := numr u; if domainp z or red z or not (tdeg (z := lt z) = 1) or atom tvar z or not ((car tvar z) eq 'log) or depend!-f(tc z,v) or depend!-f(denr u,v) then <<mlst := u . mlst;go to a>>; y := reorder tc z ./ reorder denr u; z := simp!* cadr tvar z; if x := assoc(y,llst) then rplacd(x,multsq(cdr x,z)) else if x := assoc(negsq y,llst) then rplacd(x,multsq(cdr x,invsq z)) else llst := (y . z) . llst; go to a end; symbolic procedure prod!-split!-exp(u,v); begin scalar x,y,z,w,klst,lst; % lst := kernels(numr u,nil); lst := kernels numr u; % lst := kernels1denr u,lst); lst := kernels1(denr u,lst); a: if null lst then go to b; z := car lst; if not atom z and car z eq 'expt and not depend!-p(cadr z,v) and depend!-p(caddr z,v) then klst := z . klst; lst := cdr lst; go to a; b: if null klst then return (nil . list u); x := setkorder klst; z := reorder numr u; y := reorder denr u; c: if domainp z or red z or not memq(w := mvar z,klst) then go to d; v := multsq(tdeg lt z ./ 1,simp!* caddr w); w := cadr w; if u := assoc(w,lst) then rplacd(u,addsq(cdr u,v)) else lst := (w . v) . lst; z := tc lt z; go to c; d: if domainp y or red y or not memq(w := mvar y,klst) then go to e; v := multsq(tdeg lt y ./ 1,negsq simp!* caddr w); w := cadr w; if u := assoc(w,lst) then rplacd(u,addsq(cdr u,v)) else lst := (w . v) . lst; y := tc lt y; go to d; e: setkorder x; u := reorder z ./ reorder y; return (lst . list u) end; % These can be found in Abramowitz-Stegun (27.8.6 Summable Series), and % were suggested by Winfried Neun. algebraic; let {sum(sin(~n*~tt)/n,~n,1,infinity) => (pi - tt)/2, sum(sin(~n*~tt)/(~n)^3,~n,1,infinity) => pi^2*tt/6 - pi*tt^2/4 + tt^3/12, sum(sin(~n*~tt)/(~n)^5,~n,1,infinity) => pi^4*tt/90 - pi^2*tt^3/36 + pi*tt^4/48-tt^5/240}$ let {sum(cos(~n*~tt)/(~n),~n,1,infinity) => -log(2*sin(tt/2)), sum(cos(~n*~tt)/(~n)^2,~n,1,infinity) => pi^2/6 - pi*tt/2 + tt^2/4, sum(cos(~n*~tt)/(~n)^4,~n,1,infinity) => pi^4/90 - pi^2*tt^2/12 + pi*tt^3/12-tt^4/48}$ endmodule; end;