Artifact 8cd4749f47012c308951a50d412ab0c7b81acbfb6eff1c70bc67f9a7754ffc4b:
- Executable file
r37/packages/taylor/tayexpnd.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: 23561) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/taylor/tayexpnd.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: 23561) [annotate] [blame] [check-ins using]
module tayexpnd; %***************************************************************** % % The Taylor expansion machinery % %***************************************************************** exports taylorexpand; imports % from the REDUCE kernel: !*k2q, !*p2q, .*, .+, ./, aeval, addsq, apply1, denr, dependsl, dfn_prop, diffsq, domainp, eqcar, error1, errorp, errorset!*, exptsq, kernp, lastpair, lc, let, lpow, lprim, mk!*sq, mkquote, mksq, multsq, mvar, nconc!*, neq, nlist, nth, numr, operator, prepsq, quotsq, red, rederr, setcar, sfp, simp!*, subsq, subtrsq, % from the header module: !*tay2q, cst!-Taylor!*, has!-Taylor!*, make!-cst!-coefficient, make!-Taylor!*, prune!-coefflist, set!-TayOrig, TayCfPl, TayCfSq, TayCoeffList, TayFlags, Taylor!*p, Taylor!-kernel!-sq!-p, Taylor!-trace, Taylor!-trace!-mprint, Taylor!:, TayMakeCoeff, TayOrig, TayTemplate, TayTpElNext, TayTpElOrder, TayTpElPoint, TayTpElVars, TayTpVars, TayVars, % from module TayBasic: addtaylor!*, multtaylor, multtaylor!*, quottaylor!-as!-sq, % from module TayConv: prepTaylor!*, % from module TayFns: inttaylorwrttayvar, taycoefflist!-union, % from module TayInterf: taylor1, % from module TayIntro: nzerolist, replace!-nth, smemberlp, Taylor!-error, Taylor!-error!*, var!-is!-nth, % from module TaySimp: expttayrat, taysimpsq, taysimpsq!*, % from module TayUtils: add!.comp!.tp!., addto!-all!-taytpelorders, enter!-sorted, mult!.comp!.tp!., subtr!-tp!-order, tayminpowerlist, taytp!-min2, tp!-greaterp, truncate!-Taylor!*; fluid '(!*backtrace !*taylor!-assoc!-list!* !*tayexpanding!* !*taylorkeeporiginal !*taylornocache !*tayrestart!* !*trtaylor); global '(!*sqvar!* erfg!*); symbolic smacro procedure !*tay2q!* u; ((u . 1) .* 1 .+ nil) ./ 1; switch taylornocache; symbolic procedure init!-taylor!-cache; !*taylor!-assoc!-list!* := nil . !*sqvar!*; put('taylornocache,'simpfg,'((t (init!-taylor!-cache)))); symbolic init!-taylor!-cache(); symbolic procedure taylor!-add!-to!-cache(krnl,tp,result); if null !*taylornocache then car !*taylor!-assoc!-list!* := ({krnl,sublis({nil . nil},tp)} . result) . car !*taylor!-assoc!-list!*; fluid '(!*taylor!-max!-precision!-cycles!*); symbolic(!*taylor!-max!-precision!-cycles!* := 5); symbolic procedure taylorexpand(sq,tp); begin scalar result,oldklist,!*tayexpanding!*,!*tayrestart!*,ll; integer cycles; ll := tp; oldklist := get('taylor!*,'klist); !*tayexpanding!* := t; restart: !*tayrestart!* := nil; result := errorset!*({'taylorexpand1,mkquote sq,mkquote ll,'t}, !*trtaylor); put('taylor!*,'klist,oldklist); if null errorp result then <<result := car result; if cycles>0 and Taylor!-kernel!-sq!-p result then result := !*tay2q truncate!-Taylor!*( mvar numr result,tp); return result>>; if null !*tayrestart!* then error1(); erfg!* := nil; Taylor!-trace {"Failed with template",ll}; cycles := cycles + 1; if cycles > !*taylor!-max!-precision!-cycles!* then Taylor!-error('max_cycles,cycles - 1); ll := addto!-all!-TayTpElOrders(ll,nlist(2,length ll)); Taylor!-trace {"Restarting with template",ll}; goto restart end; symbolic procedure taylorexpand1(sq,ll,flg); % % sq is a s.q. that is expanded according to the list ll % which has the form % ((var_1 var0_1 deg1) (var_2 var0_2 deg_2) ...) % flg if true indicates that the expansion order should be % automatically increased if the result has insufficient order. % begin scalar oldresult,result,lll,lmin,dorestart,nl; integer count; lll := ll; if null cadr !*taylor!-assoc!-list!* then !*taylor!-assoc!-list!* := nil . !*sqvar!*; restart: count := count + 1; if count > !*taylor!-max!-precision!-cycles!* or oldresult and TayTemplate oldresult = TayTemplate result then Taylor!-error('max_cycles,count - 1); oldresult := result; if denr sq = 1 then result := taysimpsq!* taylorexpand!-sf(numr sq,lll,t) % else if not dependsl(denr sq,TayTpVars lll) then begin scalar nn; % nn := taylorexpand!-sf(numr sq,lll,t); % if Taylor!-kernel!-sq!-p nn % then result := !*tay2q multtaylorsq( % truncate!-Taylor!*(mvar numr nn,lll), % 1 ./ denr sq) % else result := taysimpsq!* quotsq(nn,1 ./ denr sq) % end % else if numr sq = 1 then begin scalar dd; % dd := taylorexpand!-sf(denr sq,lll,nil); % if null numr dd % then Taylor!-error!*('zero!-denom,'taylorexpand) % else if Taylor!-kernel!-sq!-p dd % then if Taylor!*!-zerop mvar numr dd % then <<!*tayrestart!* := t; % Taylor!-error('zero!-denom, % 'taylorexpand)>> % else result := !*tay2q invtaylor mvar numr dd % else result := taysimpsq!* invsq dd % end % else if not dependsl(numr sq,TayTpVars lll) then begin scalar dd; % dd := taylorexpand!-sf(denr sq,lll,nil); % if null numr dd % then Taylor!-error!*('zero!-denom,'taylorexpand) % else if Taylor!-kernel!-sq!-p dd % then if Taylor!*!-zerop mvar numr dd % then <<!*tayrestart!* := t; % Taylor!-error('zero!-denom, % 'taylorexpand)>> % else result := !*tay2q multtaylorsq( % truncate!-Taylor!*( % invtaylor mvar numr dd, % lll), % numr sq ./ 1) % else result := taysimpsq!* quotsq(numr sq ./ 1,dd) % end else begin scalar nn,dd; dd := taylorexpand!-sf(denr sq,lll,nil); if null numr dd then Taylor!-error!*('zero!-denom,'taylorexpand) else if not Taylor!-kernel!-sq!-p dd then return result := taysimpsq!* quotsq(taylorexpand!-sf(numr sq,lll,t),dd); lmin := prune!-coefflist TayCoeffList mvar numr dd; if null lmin then Taylor!-error!*('zero!-denom,'taylorexpand); lmin := tayminpowerlist lmin; nn := taylorexpand!-sf( numr sq, addto!-all!-TayTpElOrders(lll,lmin),t); if not (Taylor!-kernel!-sq!-p nn and Taylor!-kernel!-sq!-p dd) then result := taysimpsq!* quotsq(nn,dd) else result := quottaylor!-as!-sq(nn,dd); end; if not Taylor!-kernel!-sq!-p result then return if not smemberlp(TayTpVars ll,result) then !*tay2q cst!-Taylor!*(result,ll) else result; result := mvar numr result; dorestart := nil; Taylor!: begin scalar ll1; ll1 := TayTemplate result; for i := (length ll1 - length ll) step -1 until 1 do ll := nth(ll1,i) . ll; if flg then <<nl := subtr!-tp!-order(ll,ll1); for each o in nl do if o>0 then dorestart := t>> end; if dorestart % if tp!-greaterp(ll,TayTemplate result) then <<lll := addto!-all!-TayTpElOrders(lll,nl); Taylor!-trace {"restarting (prec loss):", "old =",ll, "result =",result, "new =",lll}; goto restart>>; result := truncate!-Taylor!*(result,ll); if !*taylorkeeporiginal then set!-TayOrig(result,sq); return !*tay2q result end; symbolic procedure taylorexpand!-sf(sf,ll,flg); begin scalar lcof,lp,next,rest,x,dorestart,lll,xcl,nl,tps; integer l,count; lll := ll; restart: count := count + 1; if count > !*taylor!-max!-precision!-cycles!* then Taylor!-error('max_cycles,count - 1); x := nil ./ 1; rest := sf; while not null rest do << if domainp rest then <<next := !*tay2q!* cst!-Taylor!*(rest ./ 1,lll); rest := nil>> else << lp := taylorexpand!-sp(lpow rest,lll,flg); if lc rest=1 then next := lp else << lcof := taylorexpand!-sf(lc rest,lll,flg); if Taylor!-kernel!-sq!-p lcof and Taylor!-kernel!-sq!-p lp and (tps := mult!.comp!.tp!.(mvar numr lcof,mvar numr lp,nil)) then next := !*tay2q!* multtaylor!*(mvar numr lcof,mvar numr lp,tps) else next := taysimpsq!* multsq(lcof,lp); >>; rest := red rest>>; if null numr x then x := next else if Taylor!-kernel!-sq!-p x and Taylor!-kernel!-sq!-p next and (tps := add!.comp!.tp!.(mvar numr x,mvar numr next)) then x := !*tay2q!* addtaylor!*(mvar numr x,mvar numr next,car tps) else x := taysimpsq!* addsq(x,next)>>; if not Taylor!-kernel!-sq!-p x then return x; if null flg then << xcl := prune!-coefflist TayCoeffList mvar numr x; if null xcl then << lll := addto!-all!-TayTpElOrders(lll,nlist(2,length lll)); Taylor!-trace {"restarting (no coeffs)...(",count,")"}; goto restart >> else Taylor!: <<l := tayminpowerlist xcl; dorestart := nil; nl := for i := 1:length lll collect if nth(l,i)>0 then <<dorestart := t; 2*nth(l,i)>> else 0; if dorestart then <<flg :=t; lll := addto!-all!-TayTpElOrders(lll,nl); Taylor!-trace {"restarting (no cst trm)...(",count,"):", "result =",mvar numr x, "new =",lll}; goto restart>>>>>>; return x end; symbolic procedure taylorexpand!-sp(sp,ll,flg); Taylor!: begin scalar fn,krnl,pow,sk,vars; % if (sk := assoc({sp,ll},car !*taylor!-assoc!-list!*)) % then return cdr sk; Taylor!-trace "expanding s.p."; Taylor!-trace!-mprint mk!*sq !*p2q sp; vars := TayTpVars ll; krnl := car sp; pow := cdr sp; if idp krnl and krnl memq vars then <<pow := 1; sk := !*tay2q!* make!-pow!-Taylor!*(krnl,cdr sp,ll); % taylor!-add!-to!-cache(sp,TayTemplate mvar numr sk,sk) >> else if sfp krnl then sk := taylorexpand!-sf(krnl,ll,flg) else if (sk := assoc({sp,ll},car !*taylor!-assoc!-list!*)) then <<pow := 1; sk := cdr sk>> else if not(pow=1) and (sk := assoc({krnl . 1,ll},car !*taylor!-assoc!-list!*)) then sk := cdr sk else <<sk := if idp krnl then if dependsl(krnl,vars) then taylorexpand!-diff(krnl,ll,flg) else !*tay2q!* cst!-Taylor!*(simp!* krnl,ll) else if Taylor!*p krnl then if smemberlp(vars,TayVars krnl) then taylorexpand!-samevar(krnl,ll,flg) else taylorexpand!-taylor(krnl,ll,flg) else if not idp car krnl then taylorexpand!-diff(krnl,ll,flg) % else if (fn := get(car krnl,'taylordef)) % then else if null(fn := get(car krnl,'taylorsimpfn)) then taylorexpand!-diff(krnl,ll,flg) else begin scalar res,args,flg,!*taylorautoexpand; args := for each el in cdr krnl collect if not dependsl(el,vars) then el else <<flg := t; prepsq taylorexpand(simp!* el,ll)>>; if has!-Taylor!* args then res := apply1(fn,car krnl . args) else if null flg then res := !*tay2q!* cst!-Taylor!*(mksq(krnl,1),ll) else res := mksq(krnl,1); return res end; if Taylor!-kernel!-sq!-p sk then taylor!-add!-to!-cache( krnl . 1,TayTemplate mvar numr sk,sk)>>; if not(pow = 1) then <<if not Taylor!-kernel!-sq!-p sk then sk := taysimpsq exptsq(sk,pow) else <<sk := mvar numr sk; sk := !*tay2q!* if pow=2 then multtaylor(sk,sk) else expttayrat(sk,pow ./ 1)>>; if Taylor!-kernel!-sq!-p sk then taylor!-add!-to!-cache( sp,TayTemplate mvar numr sk,sk)>>; Taylor!-trace "result of expanding s.p. is"; Taylor!-trace!-mprint mk!*sq sk; return sk end; symbolic procedure make!-pow!-Taylor!*(krnl,pow,ll); Taylor!: begin scalar pos,var0,nxt,ordr,x; integer pos1; pos := var!-is!-nth(ll,krnl); pos1 := cdr pos; pos := car pos; var0 := TayTpElPoint nth(ll,pos); ordr := TayTpElOrder nth(ll,pos); nxt := TayTpElNext nth(ll,pos); % if ordr < pow % then ll := replace!-nth(ll,pos, ({TayTpElVars tpel, TayTpElPoint tpel, max2(pow,ordr),max2(pow,ordr)+nxt-ordr} where tpel := nth(ll,pos))); % if ordr < 1 then return % make!-Taylor!*(nil,replace!-nth(ll,pos, % ({TayTpElVars tpel, % TayTpElPoint tpel, % TayTpElOrder tpel, % 1} % where tpel := nth(ll,pos))), % nil,nil) % else if var0 = 0 or var0 eq 'infinity then return make!-Taylor!*( {make!-var!-coefflist(ll,pos,pos1,pow, var0 eq 'infinity)}, ll, if !*taylorkeeporiginal then mksq(krnl,pow), nil); x := make!-Taylor!*( {make!-cst!-coefficient(simp!* var0,ll), make!-var!-coefflist(ll,pos,pos1,1,nil)}, ll, nil, nil); if not (pow=1) then x := expttayrat(x,pow ./ 1); if !*taylorkeeporiginal then set!-TayOrig(x,mksq(krnl,pow)); return x end; symbolic procedure make!-var!-coefflist(tp,pos,pos1,pow,infflg); TayMakeCoeff(make!-var!-powerlist(tp,pos,pos1,pow,infflg),1 ./ 1); symbolic procedure make!-var!-powerlist(tp,pos,pos1,pow,infflg); if null tp then nil else ((if pos=1 then for j := 1:l collect if j neq pos1 then 0 else if infflg then -pow else pow else nzerolist l) where l := length TayTpElVars car tp) . make!-var!-powerlist(cdr tp,pos - 1,pos1,pow,infflg); symbolic procedure taylorexpand!-taylor(tkrnl,ll,flg); begin scalar tay,notay,x; notay := nil ./ 1; for each cc in TayCoeffList tkrnl do << % x := taylorexpand1(TayCfSq cc,ll,t); x := taylorexpand1(TayCfSq cc,ll,flg); if Taylor!-kernel!-sq!-p x then tay := nconc(tay, for each cc2 in TayCoefflist mvar numr x collect TayMakeCoeff( append(TayCfPl cc,TayCfPl cc2), TayCfSq cc2)) else Taylor!-error('expansion,"(possbile singularity)")>>; %notay := aconc!*(notay,cc)>>; return if null tay then nil ./ 1 else !*tay2q!* make!-Taylor!*(tay, append(TayTemplate tkrnl,ll), nil,nil); end; comment The cache maintained in !*!*taylorexpand!-diff!-cache!*!* is the key to handle the case of a kernel whose derivative involves the kernel itself. It is guaranteed that in every recursive step the expansion order is smaller than in the previous one; fluid '(!*!*taylorexpand!-diff!-cache!*!*); symbolic procedure taylorexpand!-diff(krnl,ll,flg); begin scalar result; % % We use a very simple strategy: if we know a partial derivative % of the kernel, we pass the problem to taylorexpand!-diff1 which % will try to differentiate the kernel, expand the result and % integrate again. % % NOTE: THE FOLLOWING test can be removed, but needs more checking % removing it seems to slow down processing % if null atom krnl and get(car krnl,dfn_prop krnl) then (result := errorset!*({'taylorexpand!-diff1, mkquote krnl,mkquote ll,mkquote flg}, !*backtrace) where !*!*taylorexpand!-diff!-cache!*!* := !*!*taylorexpand!-diff!-cache!*!*); % % If this fails we fall back to simple differentiation and % substitution at the expansion point. % if result and not errorp result then result := car result else if !*tayrestart!* then error1() % propagate else <<result := !*k2q krnl; for each el in ll do % % taylor1 is the function that does the real work % result := !*tay2q!* taylor1(result, TayTpElVars el, TayTpElPoint el, TayTpElOrder el)>>; if !*taylorkeeporiginal and Taylor!-kernel!-sq!-p result then set!-TayOrig(mvar numr result,!*k2q krnl); return result end; symbolic procedure taylorexpand!-diff1(krnl,ll,flg); <<if y and TayTpVars ll = TayTpVars(y := cdr y) and not tp!-greaterp(y,ll) then ll := for each el in y collect {TayTpElVars el,TayTpElPoint el, TayTpElOrder el - 1,TayTpElNext el - 1}; !*!*taylorexpand!-diff!-cache!*!* := (krnl . ll) . !*!*taylorexpand!-diff!-cache!*!*; for each el in ll do result := taylorexpand!-diff2(result,el,nil); result>> where result := !*k2q krnl, y := assoc(krnl,!*!*taylorexpand!-diff!-cache!*!*); symbolic procedure taylorexpand!-diff2(sq,el,flg); begin scalar l,singlist,c0,tay,l0,tp,tcl,sterm; singlist := nil ./ 1; tp := {el}; % % We check whether there is a rule for taylorsingularity. % sterm := simp!* {'taylorsingularity,mvar numr sq, 'list . TayTpElVars el,TayTpElPoint el}; if kernp sterm and eqcar(mvar numr sterm,'taylorsingularity) then sterm := nil % failed else sq := subtrsq(sq,sterm); if TayTpElOrder el > 0 then << l := for each var in TayTpElVars el collect begin scalar r; r := taylorexpand1(diffsq(sq,var), {{TayTpElVars el, TayTpElPoint el, TayTpElOrder el - 1, TayTpElNext el - 1}},flg); if Taylor!-kernel!-sq!-p r then return inttaylorwrttayvar(mvar numr r,var) else Taylor!-error('expansion,"(possible singularity)"); end; tcl := TayCoeffList!-union for each pp in l collect TayCoeffList cdr pp; for each pp in l do if car pp then singlist := addsq(singlist,car pp); if not null numr singlist then Taylor!-error('expansion,"(possible singularity)")>>; % % If we have a special singularity, then set c0 to zero. % if not null sterm then c0 := nil ./ 1 else c0 := subsq(sq,for each var in TayTpElVars el collect (var . TayTpElPoint el)); l0 := {nzerolist length TayTpElVars el}; if null numr c0 then nil else if not Taylor!-kernel!-sq!-p c0 then tcl := TayMakeCoeff(l0,c0) . tcl else <<c0 := mvar numr c0; tp := nconc!*(TayTemplate c0,tp); for each pp in TayCoeffList c0 do tcl := enter!-sorted(TayMakeCoeff(append(TayCfPl pp,l0), TayCfSq pp), tcl)>>; if not null l then for each pp in l do tp := TayTp!-min2(tp,TayTemplate cdr pp); tay := !*tay2q!* make!-Taylor!*(tcl,tp, if !*taylorkeeporiginal then sq else nil,nil); if not null numr singlist then tay := addsq(singlist,tay); if null sterm then return tay else return taysimpsq!* addsq(taylorexpand(sterm,tp),tay) end; algebraic operator taylorsingularity; if null get('psi,'simpfn) then algebraic operator psi; algebraic let { taylorsingularity(dilog(~x),~y,~y0) => pi^2/6 - log(x)*log(1-x), taylorsingularity(ei(~x),~y,~y0) => log(x) - psi(1) }; symbolic procedure taylorexpand!-samevar(u,ll,flg); Taylor!: begin scalar tpl; tpl := TayTemplate u; for each tpel in ll do begin scalar tp,varlis,mdeg,n; integer pos; varlis := TayTpElVars tpel; pos := car var!-is!-nth(tpl,car varlis); tp := nth(tpl,pos); if length varlis > 1 and not (varlis = TayTpElVars tp) then Taylor!-error('not!-implemented, "(homogeneous expansion in TAYLORSAMEVAR)"); n := TayTpElOrder tpel; if TayTpElPoint tp neq TayTpElPoint tpel then u := taylor1(if not null TayOrig u then TayOrig u else simp!* prepTaylor!* u, varlis,TayTpElPoint tpel,n); mdeg := TayTpElOrder tp; if n=mdeg then nil else if n > mdeg % % further expansion required % then if null flg then nil else Taylor!-error('expansion, "Cannot expand further... truncated.") else u := make!-Taylor!*( for each cc in TayCoeffList u join if nth(nth(TayCfPl cc,pos),1) > n then nil else list cc, replace!-nth(tpl,pos,{varlis,TayTpElPoint tpel,n,n+1}), TayOrig u,TayFlags u) end; return !*tay2q!* u end; endmodule; end;